西电电院自动控制技术上机报告.doc

《西电电院自动控制技术上机报告.doc》由会员分享，可在线阅读，更多相关《西电电院自动控制技术上机报告.doc（40页珍藏版）》请在装配图网上搜索。

自动控制技术上机实验报告 班级：021215 学号：021214 姓名： 时域分析 程序源代码： close all; clear all; ft = 30; M=1; B=5; K=20; %系统参数 t0=0; tfinal = 5; tspan = [t0 tfinal]; %设置仿真开始和结束时间 x0 = [0,0]; %系统初始值，零初始条件 options = odeset(AbsTol,[1e-6;1e-6]); %设置仿真计算精度 [t,x] = ode113(xt4odefile,tspan,x0,options); %微分方程求解，计算位移x(:,1)和速度x(:,2) a = 1/M*(ft-B*x(:,2)-K*x(:,1)); %计算加速度 i = 1; while (abs(a(i))>0.0001|(abs(x(i,2))>0.0001)) i = i+1; end disp(稳态时系统的位移、速度和加速度及对应的时间分别为： ); %显示计算结果 result = sprintf(位移 d=%6.4f\n,x(i,1)); disp(result); result = sprintf(速度 v=%8.6f\n,x(i,2)); disp(result); result = sprintf(加速度 a=%9.6f\n,a(i)); disp(result); result = sprintf(时间 t=%4.2f\n,t(i)); disp(result); d = x(:,1); subplot(1,3,1),plot(t,d); %绘制时间-位移曲线 xlabel(时间（秒）);ylabel(位移（米）); title(时间-位移曲线);grid; v = x(:,2); subplot(1,3,2),plot(t,v); %绘制时间-速度曲线 xlabel(时间（秒）);ylabel(速度（米/秒）); title(时间-速度曲线);grid; subplot(1,3,3),plot(d,v); %绘制位移-速度曲线 xlabel(位移（米）);ylabel(速度（米/秒）); title(位移-速度曲线);grid; 运行结果： >> EX11 稳态时系统的位移、速度和加速度及对应的时间分别为： 位移 d=1.5000 速度 v=-0.000086 加速度 a=-0.000084 时间 t=4.46 源程序代码： close all; clear all; num = [0,2,5,7]; den = [1,6,10,6]; [z1,p1,k1] = tf2zp(num,den) [r2,p2,k2] = residue(num,den) 运行结果： >> EX12 z1 = -1.2500 + 1.3919i -1.2500 - 1.3919i p1 = -3.7693 + 0.0000i -1.1154 + 0.5897i -1.1154 - 0.5897i k1 = 2 r2 = 2.2417 + 0.0000i -0.1208 - 1.0004i -0.1208 + 1.0004i p2 = -3.7693 + 0.0000i -1.1154 + 0.5897i -1.1154 - 0.5897i k2 = [] 源程序代码： clc a=[6.3223 18 12.811]; b=[1 6 1.322 18 12.811]; sys=tf(a,b); t=0:.005:35; step(sys) title(系统的单位阶跃响应) [y,t]=step(a,b,t); r=1; while(y(r)<1.001) r=r+1;end; rise_time=(r-1)*0.005 %求系统上升时间 [ymax,tp]=max(y); peak_time=(tp-1)*0.005 k=dcgain(sys); max_overshoot=(ymax-k)/k s=7001; while(y(s)>0.98&y(s)<1.02) s=s-1;end setting_time=(s-1)*0.005 运行结果： rise_time = 0.7050 peak_time = 34.6100 max_overshoot = 2.6819e+06 setting_time = 35 源程序代码： num=[1 1];den=[1 5 6 0]; sys1=tf(num,den) subplot(3,1,1) step(sys1) subplot(3,1,2) rlocus(num,den) title(系统的根轨迹曲线) r=rlocus(num,den,20.575) a=[1];b=[1 0.8989]; sys2=tf(a,b) sys=series(sys1,sys2) %串联 subplot(3,1,3) step(sys) 运行结果： rise_time = 0.7050 peak_time = 34.6100 max_overshoot = 2.6819e+06 setting_time = 35 >> EX14 sys1 = s + 1 ----------------- s^3 + 5 s^2 + 6 s Continuous-time transfer function. r = -2.0505 + 4.3225i -2.0505 - 4.3225i -0.8989 + 0.0000i sys2 = 1 ---------- s + 0.8989 Continuous-time transfer function. sys = s + 1 ------------------------------------- s^4 + 5.899 s^3 + 10.49 s^2 + 5.393 s Continuous-time transfer function. 频域分析： 源程序代码： clc num=[0.01,0.0001,0.01]; dun=[0.25,0.01,1,0,0]; sys=tf(num,dun) figure(1) bode(sys) figure(2) sys2=feedback(sys,1) bode(sys2) 运行结果： sys = 0.01 s^2 + 0.0001 s + 0.01 -------------------------- 0.25 s^4 + 0.01 s^3 + s^2 Continuous-time transfer function. sys2 = 0.01 s^2 + 0.0001 s + 0.01 ------------------------------------------------ 0.25 s^4 + 0.01 s^3 + 1.01 s^2 + 0.0001 s + 0.01 Continuous-time transfer function. 源程序代码： close all; clear all; num = [0 20 20 10]; %开环传递函数分子 den = conv([1 1 0],[1 10]); %开环传递函数的分母 nyquist(num,den) %v = [-2 3 -3]; axis([-2 2 -3 3]) grid x = -pi:0.01:pi; plot(x,sin(x)), grid on 运行结果： 源程序代码： close all; clear all; num = [2000,2000]; %开环传递函数的分子 den = conv([1 0.5 0],[1 14 400]); %开环传递函数的分母 nichols(num,den) %绘制nichols图 v = [-270 -90 -40 40]; axis(v) ngrid %标出nichols图线 运行结果： 源程序代码： num=[0 2000 2000]; den=conv([1 0.5 0],[1 14 400]); h=tf(num,den); [gm,pm,wg,wc]=margin(h); gm,pm,wg,wc 运行结果： >> EX24 gm = 2.7493 pm = 73.3527 wg = 19.8244 wc = 5.3477 源程序代码： clc num=[1]; den=[0.5 1.5 1 0]; sys=tf(num,den) sys2=feedback(sys,1) bode(sys2) [gm pm wg wp]=margin(sys) 运行结果： ------------------------- 0.5 s^3 + 1.5 s^2 + s + 1 Continuous-time transfer function. gm = 3.0000 pm = 32.6133 wg = 1.4142 wp = 0.7494 现代控制理论 源程序代码： close all; clear all; %3.1_A num = [1 2 3]; %传递函数分子多项式的系数 den = [1 3 3 1]; %传递函数分母多项式的系数 [A,B,C,D] = tf2ss(num,den) %3.1_B z = [-1;-3]; %传递函数的零点 p = [0;-2;-4;-6;]; %传递函数的极点 k = 4; [A,B,C,D] = zp2ss(z,p,k) %3.1_C A = [0,1;-1,-2]; B = [0;1]; C = [1,3]; D = [1]; [num,den] = ss2tf(A,B,C,D) printsys(num,den,s) [z,p,k] = ss2zp(A,B,C,D) 运行结果： >> EX31 A = -3 -3 -1 1 0 0 0 1 0 B = 1 0 0 C = 1 2 3 D = 0 A = -10.0000 -4.8990 0 0 4.8990 0 0 0 -6.0000 -4.2866 -2.0000 0 0 0 1.0000 0 B = 1 0 1 0 C = 0 0 0 4 D = 0 num = 1.0000 5.0000 2.0000 den = 1 2 1 num/den = s^2 + 5 s + 2 ------------- s^2 + 2 s + 1 z = -0.4384 -4.5616 p = -1 -1 k = 1 源程序代码： close all; clear all; A1 = [0,1;-1,-2]; B1 = [0;1]; C1 = [1,3]; D1 = [1]; A2 = [0,1;-1,-3]; B2 = [0;1]; C2 = [1,4]; D2 = [0]; [A,B,C,D] = series(A1,B1,C1,D1,A2,B2,C2,D2) [A,B,C,D] = parallel(A1,B1,C1,D1,A2,B2,C2,D2) [A,B,C,D] = feedback(A1,B1,C1,D1,A2,B2,C2,D2) [A,B,C,D] = feedback(A1,B1,C1,D1,A2,B2,C2,D2,+1) 运行结果： >> EX32 A = 0 1 0 0 -1 -3 1 3 0 0 0 1 0 0 -1 -2 B = 0 1 0 1 C = 1 4 0 0 D = 0 A = 0 1 0 0 -1 -2 0 0 0 0 0 1 0 0 -1 -3 B = 0 1 0 1 C = 1 3 1 4 D = 1 A = 0 1 0 0 -1 -2 -1 -4 0 0 0 1 1 3 -2 -7 B = 0 1 0 1 C = 1 3 -1 -4 D = 1 A = 0 1 0 0 -1 -2 1 4 0 0 0 1 1 3 0 1 B = 0 1 0 1 C = 1 3 1 4 D = 1 源程序代码： clc; close all; clear all; A = [0 -2;1 -3]; t = 0.2; Phi = expm(A*t); %求状态转移矩阵 B = [2;0]; C = [0 3]; D = [0]; x0 = [1 1]; t = [0 0.2]; u = 0*t; [y,x] = lsim(A,B,C,D,u,t,x0) %求系统响应 运行结果： y = 3.0000 2.0110 x = 1.0000 1.0000 0.6703 0.6703 源程序代码： clc A=[-3 1;1 -3]; B=[1 1;1 1]; C=[1 1;1 -1]; D=[0 0;0 0]; N=size(A); n=N(1); [num,den]=ss2tf(A,B,C,D,2); disp(可控矩阵：) S=ctrb(A,B) f=rank(S) if(f==n) disp(系统是可控的) else disp(系统是不可控的) end disp() disp(可观测矩阵：) V=obsv(A,C) m=rank(V) if(f==n) disp(系统时可观测的) else disp(系统是不可观测的) end 运行结果： 可控矩阵： S = 1 1 -2 -2 1 1 -2 -2 f = 1 系统是不可控的 可观测矩阵： V = 1 1 1 -1 -2 -2 -4 4 m = 2 系统是不可观测的 源程序代码： clc A=[0 1;-2 -3]; B=[0;1]; C=[2 0];D=0; P_S=[-1 -2]; k=acker(A,B,P_S); P_O=[-3 -3]; h=(acker(A,C,P_O)); A1=[A,-B*k;h*C A-B*k-h*C]; B1=[B;B]; C1=[C,zeros(1,2)]; D1=0; sys=ss(A1,B1,C1,D1) tf(sys) 运行结果： 校正设计 源程序代码： ts=0.001; sys=tf([400],[1,50,0]); dsys=c2d(sys,ts,z); [num,den]=tfdata(dsys,v); u_1=0.0;u_2=0.0; y_1=0.0;y_2=0.0; x=[0,0,0]; error_1=0; error_2=0; for k=1:1:1000 time(k)=k*ts; kp=8;ki=0.1;kd=10; rin(k)=0.5*sin(2*pi*k*ts); du(k)=kp*x(1)+kd*x(2)+ki*x(3); u(k)=u_1+du(k); yout(k)=-den(2)*y_1-den(3)*y_2+num(2)*u_1+num(3)*u_2; error(k)=rin(k)-yout(k); u_2=u_1;u_1=u(k); y_2=y_1;y_1=yout(k); x(1)=error(k)-error_1; x(2)=error(k)-2*error_1+error_2; x(3)=error(k); error_2=error_1; error_1=error(k); end figure(1); plot(time,rin,b,time,yout,r),grid on gtext(rin\rightarrow) gtext(\leftarrowyout) title(系统输出曲线) xlabel(time(s)),ylabel(rin,yout); figure(2); plot(time,error,r),grid on title(误差曲线) xlabel(time(s));ylabel(error); 注：输入信号为正弦，采样信号为1ms。 运行结果： 源程序代码： clc; clear all; ts=20; sys=tf([1],[60,1],inputdelay,80); dsys=c2d(sys,ts,z); [num,den]=tfdata(dsys,v); u_1=0;u_2=0;u_3=0;u_4=0;u_5=0; y_1=0;y_2=0;y_3=0; error_1=0; ei=0; for k=1:1:200 time(k)=k*ts; rin(k)=40; kp=0.80; ki=0.005; kd=3.0; yout(k)=-den(2)*y_1+num(2)*u_5; error(k)=rin(k)-yout(k); M=2; if M==1 if abs(error(k))>=30&abs(errpr(k))<=40 beta=0.3; elseif abs(error(k))>=20&abs(error(k))<=30 beta=0.6; elseif abs(error(k))>=10&abs(error(k))<=20 beta=0.9; else beta=1.0; end elseif M==2 beta=1.0; end ei=ei+error(k)*ts; u(k)=kp*error(k)+kd*(error(k)-error_1)/ts+beta*ki*ei; if u(k)>=110 u(k)=110; end if u(k)<=-110 u(k)=-110; end u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=u(k); y_3=y_2;y_2=y_1;y_1=yout(k); error_1=error(k); end figure(1); plot(time,rin,b,time,yout,r); title(不加积分作用效果) xlabel(time(s));ylabel(rin,yout); figure(2); plot(time,u,r); xlabel(time(s));ylabel(u); 运行结果： M=2时，不加积分效果： M=1时，加积分效果： 源程序代码： clc; clear all; ts=0.001; sys=tf(5.235e005,[1,87.35,1.047e004,0]); dsys=c2d(sys,ts,z); [num,den]=tfdata(dsys,v); u_1=0.0;u_2=0.0;u_3=0.0; y_1=0;y_2=0;y_3=0; x=[0,0,0]; error_1=0; um=6; kp=0.85;ki=9.0;kd=0.0; rin=30; for k=1:1:800 time(k)=k*ts; u(k)=kp*x(1)+kd*x(2)+ki*x(3); if u(k)>=um u(k)=um; end if u(k)<=-um y(k)=-um; end yout(k)=-den(2)*y_1-den(3)*y_2-den(4)*y_3+num(2)*u_1+num(3)*u_2+num(4)*u_3; error(k)=rin-yout(k); M=1; if M==1 if u(k)>=num if error(k)>0 alpha=0; else alpha=1; end elseif u(k)<=-um if error(k)>0 alpha=1; else alpha=0; end else alpha=1; end elseif M==2 alpha=1; end u_3=u_2;u_2=u_1;u_1=u(k); y_3=y_2;y_2=y_1;y_1=yout(k); error_1=error(k); x(1)=error(k); x(2)=(error(k)-error_1)/ts; x(3)=x(3)+alpha*error(k)*ts; xi(k)=x(3); end figure(1); subplot(3,1,1); plot(time,rin,b,time,yout,r); xlabel(time(s));ylabel(Controller?output); subplot(3,1,3); plot(time,xi,r); xlabel(time(s));ylabel(Integration); 运行结果： M=2时，无抗积分饱和： M=1时，有抗积分饱和： 源程序代码： clc; clear all; ts=0.001; sys=tf(5.235e005,[1,87.35,1.047e004,0]); dsys=c2d(sys,ts,z); [num,den]=tfdata(dsys,v); u_1=0;u_2=0;u_3=0;u_4=0;u_5=0; y_1=0;y_2=0;y_3=0; yy_1=0; error_1=0;error_2=0;ei=0; sys1=tf([1],[0.04,1]); dsys1=c2d(sys1,ts,tucsin); [num1,den1]=tfdata(dsys1,v); f_1=0; for k=1:1:2000 time(k)=k*ts; rin(k)=1; yout(k)=-den(2)*y_1-den(3)*y_2-den(4)*y_3+num(2)*u_1+num(3)*u_2+num(4)*u_3; D(k)=0.50*rands(1); yyout(k)=yout(k)+D(k); filty(k)=-den1(2)*f_1+num1(1)*(yyout(k)+yy_1); error(k)=rin(k)-filty(k); if abs(error(k))<=0.20 ei=ei+error(k)*ts; else ei=0; end kp=0.50;ki=0.10;kd=0.020; u(k)=kp*error(k)+ki*ei+kd*(error(k)-error_1)/ts; M=2; if M==1 u(k)=u(k); elseif M==2 if abs(error(k))<=0.10 u(k)=0; end end if u(k)>=10 u(k)=10; end if u(k)<=-10 u(k)=-10; end rin_1=rin(k); u_3=u_2;u_2=u_1;u_1=u(k); y_3=y_2;y_2=y_1;y_1=yout(k); f_1=filty(k); yy_1=yyout(k); error_2=error_1; error_1=error(k); end figure(1); subplot(2,1,1); plot(time,rin,r,time,filty,b); xlabel(time(s));ylabel(rin,yout); subplot(2,1,2); plot(time,u,r); xlabel(time(s));ylabel(u); figure(2); plot(time,D,r); xlabel(time(s));ylabel(Disturbance?signal); 运行结果： M=1时，不带死区结果： M=2时，带死区结果：