东北大学 自动化 控制系统计算机辅助设计实验.docx
控制系统计算机辅助设计RTO11RTORTO控制系统计算机辅助设计第一部分1 function dx=lorenzeq(t,x)dx=-x(2)-x(3);x(1)+0.2*x(2);0.2+(x(1)-5.7)*x(3);>> x0=0;0;0;>> t,y=ode45('lorenzeq',0,100,x0);>> plot(t,y)>> figure;plot3(y(:,1),y(:,2),y(:,3),grid,2function y,yeq=f2a(x)yeq=;y=4*x(1)2+x(2)2-4;>> Aeq=;Beq=;A=;B=;>> xm=0;0; xM=;x0=0;0;>> f1=inline('x(1)2-2*x(1)+x(2)');>> x,f=fmincon(f1,x0,A,B,Aeq,Beq,xm,xM,'f2a');x',fans = 1.0000 0f =-13（a）> >s=tf('s');G=(s3+4*s+2)/s3/(s2+2)/(s2+1)3+2*s+5Transfer function: s3 + 4 s + 2-s11 + 5 s9 + 9 s7 + 2 s6 + 12 s5 + 4 s4 + 12 s3（b）>> z=tf('z',0.1);>> H=(z2+0.568)/(z-1)/(z2-0.2*z+0.99)Transfer function: z2 + 0.568-z3 - 1.2 z2 + 1.19 z - 0.99Sampling time: 0.14>> A=0 1 0;0 0 1;-5 -4 -13;>> B=0;0;2;>> C=1 0 0;0 0 0;0 0 0;>> D=0;>> G=ss(A,B,C,D);>> Ga = x1 x2 x3 x1 0 1 0 x2 0 0 1 x3 -5 -4 -13b = u1 x1 0 x2 0 x3 2c = x1 x2 x3 y1 1 0 0 y2 0 0 0 y3 0 0 0d = u1 y1 0 y2 0 y3 0Continuous-time model.>> G=tf(G)Transfer function from input to output. 2 #1: - s3 + 13 s2 + 4 s + 5 #2: 0 #3: 0>> GG=zpk(G)Zero/pole/gain from input to output. 2 #1: - (s+12.72) (s2 + 0.2836s + 0.3932) #2: 0 #3: 0根据微分方程也可以直接写出传递函数模型：>> num=2;>> den=1,13,4,5;>> G=tf(num,den);>> GTransfer function: 2-s3 + 13 s2 + 4 s + 5>> GG=zpk(G)Zero/pole/gain: 2-(s+12.72) (s2 + 0.2836s + 0.3932)5>> num=1,2;>> den=1,1,0.16;>> H=tf(num,den,'Ts',1);>> HTransfer function: z + 2-z2 + z + 0.16Sampling time: 16function H=feedback(G1,G2,key)if nargin=2;key=-1;end,H=G1/(sym(1)-key*G1*G2);H=simple(H);>> syms J Kp Ki s;>> gc=(Kp*s+Ki)/s;>> g=(s+1)/(J*s2+2*s+5);>> gg=feedback(g*gc,1)>> gg=feedback(g*gc,1)gg =(Ki + Kp*s)*(s + 1)/(J*s3 + (Kp + 2)*s2 + (Ki + Kp + 5)*s + Ki)7（a）>> s=tf('s');>> G=(211.87*s+317.64)/(s+20)/(s+94.34)/(s+0.1684);>> Gc=(169.6*s+400)/s/(s+4);>> Hs=1/(0.01*s+1);>> GG=feedback(G*Gc,Hs)Transfer function: 359.3 s3 + 3.732e004 s2 + 1.399e005 s + 127056-0.01 s6 + 2.185 s5 + 142.1 s4 + 2444 s3 + 4.389e004 s2 + 1.399e005 s + 127056>> zpk(GG)Zero/pole/gain: 35933.152 (s+100) (s+2.358) (s+1.499)-(s2 + 3.667s + 3.501) (s2 + 11.73s + 339.1) (s2 + 203.1s + 1.07e004)（b）>> z=tf('z');>> G=(35786.7*z-1+108444)/(z-1+4)*(z-1+20)*(z-1+74.04);>> Gc=1/(z-1-1);>> H=1/(0.5*z-1-1);>> GG=feedback(G*Gc,H)Transfer function: -108444 z6 + 1.844e004 z5 + 1.789e004 z4-1.144e005 z6 + 2.876e004 z5 + 274.2 z4 + 782.4 z3 + 47.52 z2 + 0.5 z Sampling time: unspecified>> zpk(GG)Zero/pole/gain: -0.94821 z4 (z-0.5) (z+0.33)-z (z+0.3035) (z+0.04438) (z+0.01355) (z2 - 0.11z + 0.02396)Sampling time: unspecified8>> s=tf('s');c1=feedback(1/(s+1)*s/(s2+2),(4*s+2)/(s+1)2);c2=feedback(1/s2,50);G=feedback(c1*c2,(s2+2)/(s3+14)Transfer function: s6 + 2 s5 + s4 + 14 s3 + 28 s2 + 14 s-s10 + 3 s9 + 55 s8 + 175 s7 + 300 s6 + 1323 s5 + 2656 s4 + 3715 s3 + 7732 s2 + 5602 s + 1400 9>> s=tf('s');>> G=(s+1)2*(s2+2*s+400)/(s+5)2/(s2+3*s+100)/(s2+3*s+2500);>> G1=c2d(G,0.01)Transfer function:4.716e-005 z5 - 0.0001396 z4 + 9.596e-005 z3 + 8.18e -005 z2 - 0.0001289 z + 4.355e-005 -z6 - 5.592 z5 + 13.26 z4 - 17.06 z3 + 12.58 z2 - 5.032 z + 0.8521Sampling time: 0.01>>step(G1)>> G2=c2d(G,0.1)Transfer function:0.0003982 z5 - 0.0003919 z4 - 0.000336 z3 + 0.0007842 z2 - 0.000766 z + 0.0003214-z6 - 2.644 z5 + 4.044 z4 - 3.94 z3 + 2.549 z2 - 1.056 z + 0.2019Sampling time: 0.1>>step(G2)>> G3=c2d(G,1)Transfer function:8.625e-005 z5 - 4.48e-005 z4 + 6.545e-006 z3 + 1.211e-005 z2 - 3.299e-006 z + 1.011e-007-z6 - 0.0419 z5 - 0.07092 z4 - 0.0004549 z3 + 0.002495 z2 - 3.347e-005 z + 1.125e-007 Sampling time: 1>> step(G3)10（a）>> G=1/(s3+2*s2+s+2);>>pzmap(G)系统极点均在虚轴左侧，系统稳定（b）>> G=1/(6*s4+3*s3+2*s2+s+1);>> pzmap(G)系统极点在虚轴右侧侧，系统不稳定（c）>> G=1/(s4+s3-3*s2-s+2);>> pzmap(G)系统极点在虚轴右侧侧，系统不稳定11（a）>> z=tf('z',0.1);>> H=(-3*z+2)/(z3-0.2*z2-0.25*z+0.05);>> pzmap(H)系统所有极点均在单位圆内，所以系统稳定（b）>> z=tf('z',0.1);>> H=(3*z2-0.39*z-0.09)/(z4-1.7*z3+1.04*z2+0.268*z+0.024);>> pzmap(H)系统所有极点不全单位圆内，所以系统不稳定12>> A=-0.2 0.5 0 0 0;0 -0.5 1.6 0 0;0 0 -14.3 85.8 0;0 0 0 -33.3 100;0 0 0 0 -10;>> B=0;0;0;0;30;>> C=0 0 0 0 0;>> G=ss(A,B,C,0)a = x1 x2 x3 x4 x5 x1 -0.2 0.5 0 0 0 x2 0 -0.5 1.6 0 0 x3 0 0 -14.3 85.8 0 x4 0 0 0 -33.3 100 x5 0 0 0 0 -10b = u1 x1 0 x2 0 x3 0 x4 0 x5 30c = x1 x2 x3 x4 x5 y1 0 0 0 0 0d = u1 y1 0Continuous-time model.>> pzmap(G)系统所有极点均在虚轴左侧，所以系统稳定或>> eig(G)ans = -0.2000 -0.5000 -14.3000 -33.3000 -10.0000极点均在左半轴，所以系统稳定13数值解：>> f=(t,x)-5*x(1)+2*x(2);-4*x(2);-3*x(1)+2*x(2)-4*x(3)-x(4);-3*x(1)+2*x(2)-4*x(4);>> t_final=10;>> x0=1 2 0 0;>> t,x=ode45(f,0,t_final,x0);>> plot(t,x)解析解function Ga,Xa=ss_augment(G,cc,dd,X)G=ss(G);Aa=G.a;Ca=G.c;Xa=X;Ba=G.b;D=G.d;if (length(dd)>0&sum(abs(dd)>1e-5), if (abs(dd(4)>1e-5), Aa=Aa dd(2)*Ba,dd(3)*Ba;. zeros(2,length(Aa),dd(1),-dd(4);dd(4),dd(1); Ca=Ca dd(2)*D dd(3)*D;Xa=Xa;1;0;Ba=Ba;0;0; else, Aa=Aa dd(2)*B;zeros(1,length(Aa) dd(1); Ca=Ca dd(2)*D;Xa=Xa;1;Ba=B;0; endendif (length(cc)>0&sum(abs(cc)>1e-5),M=length(cc); Aa=Aa Ba zeros(length(Aa),M-1);zeros(M-1,length(Aa)+1). eye(M-1);zeros(1,length(Aa)+M); Ca=Ca D zeros(1,M-1);Xa=Xa;cc(1);ii=1; for i=2:M,ii=ii*i;Xa(length(Aa)+i)=cc(i)*ii; end,endGa=ss(Aa,zeros(size(Ca'),Ca,D);>> cc=2;dd=-3,0,2,2;x0=1;2;0;1;>> A=-5,2,0,0;0,-4,0,0;-3,2,-4,-1;-3,2,0,-4;B=0;0;0;1;C=2 1 0 0;D=0;>> G=ss(A,B,C,D);>> Ga,xx0=ss_augment(G,cc,dd,x0);Ga.a,xx0'ans = -5 2 0 0 0 0 0 0 -4 0 0 0 0 0 -3 2 -4 -1 0 0 0 -3 2 0 -4 0 2 1 0 0 0 0 -3 -2 0 0 0 0 0 2 -3 0 0 0 0 0 0 0 0ans = 1 2 0 1 1 0 2>> syms t;y=Ga.c*expm(Ga.a*t)*xx0;>> latex(y);>> latex(y)ans =-6,e-5,t+10,e-4,t14（a）>> s=tf('s');>> g=(s+6)*(s-6)/s/(s+3)/(s+4+4i)/(s+4-4i);>> rlocus(g)（b）>> s=tf('s');>> G=(s2+2*s+2)/(s4+s3+14*s2+8*s);>> rlocus(G)K<5.6215>> s=tf('s');>> G=(s-1)/(s+1)5;>> G.ioDelay=2Transfer function: s - 1exp(-2*s) * - s5 + 5 s4 + 10 s3 + 10 s2 + 5 s + 1>> rlocus(pade(G,2)16（a）>> s=tf('s');>> G=8*(s+1)/s2/(s+15)/(s2+6*s+10);>> nyquist(G),>> nyquist(G),grid>> bode(G)>> figure;nichols(G),grid>> gm,pm,wg,wp=margin(G)gm = 30.4686pm = 4.2340wg = 1.5811wp =0.2336>> GG=feedback(G,1)Transfer function: 8 s + 8-s5 + 21 s4 + 100 s3 + 150 s2 + 8 s + 8>> pzmap(GG)通过以上的分析，可以看出该系统是稳定的。采用阶跃响应进行验证如下图：>> pzmap(GG)>> step(GG)（b）>> Z=-1.31;-0.054;0.957;>> P=0;1;0.368;0.99;>> G=zpk(Z,P,0.45,'Ts',0.1)Zero/pole/gain:0.45 (z+1.31) (z-0.957) (z+0.054)- z (z-1) (z-0.99) (z-0.368)Sampling time: 0.1>> nyquist(G);grid>> bode(G)>> nichols(G),grid>> gm,pm,wg,wp=margin(G)gm = 0.9578pm = -1.7660wg = 10.4641wp = 10.7340>> GG=feedback(G,1);>> pzmap(GG)>> step(GG)17>> z=-2.5;p=0;-0.5;-50;>> G=zpk(z,p,100/2.5*0.5*50);>> z=-1;-2.5;p=-0.5;-50>>Gc=zpk(z,p,1000);>> GG=feedback(G*Gc,1) Zero/pole/gain: 1000000 (s+1) (s+2.5)2-(s+1) (s2 + 4.99s + 6.239) (s2 + 95.01s + 1.002e006)>> bode(GG)通过bode图可以看出该系统是稳定的。验证如下：>> step(GG)第二部分2>> Y=dsolve('D4y+5*D3y+6*D2y+4*Dy+2*y=exp(-3*t)+exp(-5*t)*sin(4*t+pi/3)',.)'y(0)=1','Dy(0)=1/2','D2y(0)=1/2','D3y(0)=0.2');>> latex(Y)3Simulink仿真图形：系统输出曲线：系统误差曲线：4x1(t)=sin(x2(t)*exp(-2.3*x4(t)x2(t)=x1(t)x3=sin(x2*exp(-2.3*x4)x4=x35Simulink仿真图形：阶跃响应曲线：6>> s=tf('s');>>g=210*(s+1.5)/(s+1.75)*(s+16)*(s+1.5+3*j)*(s+1.5-3*j);>> gc=52.5*(s+1.5)/(s+14.86);>> step(feedback(g,1)>> step(feedback(g*gc,1)通过两个图比较可以看出，原系统超调量过大，震荡严重，加入控制器后，系统变得不稳定。7>> A=0 1 0 0;0 0 1 0;-3 1 2 3;2 1 0 0;>> B=1 0;2 1;3 2;4 3;>> Q=diag(1,2,3,4);>> R=eye(2);>> K,S=lqr(A,B,Q,R)K = -0.0978 1.2118 1.8767 0.7871 -3.8819 -0.4668 2.6713 1.0320S = 5.4400 0.6152 -2.3163 0.0452 0.6152 1.8354 -0.0138 -0.7582 -2.3163 -0.0138 1.9214 -0.3859 0.0452 -0.7582 -0.3859 0.8540>> eig(A-B*K)ans = -12.2563 -1.6786 + 0.9981i -1.6786 - 0.9981i -1.4627 >> step(ss(A-B*K,B,eye(4),zeros(4,2)8>> A=-0.2 0.5 0 0 0;0 -0.5 1.6 0 0;0 0 -14.3 85.8 0;0 0 0 -33.3 100;0 0 0 0 -10;>> B=0;0;0;0;30;>> C=1 0 0 0 0;>> D=0;>> G=ss(A,B,C,D);>> pole(G)ans = -0.2000 -0.5000 -14.3000 -33.3000 -10.0000>> zero(G)ans = Empty matrix: 0-by-1>> P=-1 -2 -3 -4 -5;>> K=place(A,B,P)K = Columns 1 through 4 0.0004 0.0004 -0.0035 0.3946 Column 5 -1.4433>> eig(A-B*K)ans = -5.0000 -4.0000 -3.0000 -2.0000 -1.00009>> open_system(simulink)PID控制参数如下：0.40725 0.19096 0.21986 7.9397 0.9011 78.237724
