东北大学自动化控制系统计算机辅助设计实验

1、东北大学11控制系统计算机辅助设计MATLAB语言与应用控制系统计算机辅助设计第一部分1 function>> x0=0;0;0;>> t,y=ode45('lorenzeq',0,100,x0);>> plot(t,y)>> figure;plot3(y(:,1),y(:,2),y(:,3),grid,2function>> Aeq=;Beq=;A=;B=;>> xm=0;0; xM=;x0=0;0;>> f1=inline('x(1)2-2*x(1)+x(2)');>&

2、gt; 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+

3、0.99)Transfer function:-4>> 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 0x2 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.>&

4、gt; 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 fu

5、nction: 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-Sampling time: 16functionifend>> syms J Kp Ki s;>> gc=(Kp*s+Ki)/s;>> g=

6、(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

7、=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)（

8、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 Sampli

9、ng 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

10、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

11、 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 >>step(G1)>> G2=c2d(G,0.1)Transfer function:0.0003982 z5 - 0.0003919 z4 - 0.000336 z3 + 0.0007842 z2 -z6 - 2.644 z5 + 4.044 z4 - 3.94 z3 + 2.549 z2 - 1.056 z >>s

12、tep(G2)>> G3=c2d(G,1)Transfer function:8.625e-005 z5 - 4.48e-005 z4 + 6.545e-006 z3 + 1.211e-005 z2 -z6 - 0.0419 z5 - 0.07092 z4 - 0.0004549 z3 + 0.002495 z2 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);>

13、;> 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

14、*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

15、 0 0 -10b = u1 x1 0 x2 0 x3 0 x4 0 x5 30c = x1 x2 x3 x4 x5y1 0 0 0 0 0d = u1 y1 0Continuous-time model.>> pzmap(G)系统所有极点均在虚轴左侧，所以系统稳定或>> eig(G)ans =极点均在左半轴，所以系统稳定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=

16、1 2 0 0;>> t,x=ode45(f,0,t_final,x0);>> plot(t,x)解析解functionifif.elseendendif.forendend>> 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

17、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);>&

18、gt; rlocus(g)（b）>> s=tf('s');>> G=(s2+2*s+2)/(s4+s3+14*s2+8*s);>> rlocus(G)15>> 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'

19、;);>> 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 =pm =wg =wp =>> GG=feedback(G,1)Transfer function: 8 s + 8-s5 + 21 s4 + 100 s3 + 150 s2 + 8 s + 8>> pzmap(GG)通过以上的分析，可以看

20、出该系统是稳定的。采用阶跃响应进行验证如下图：>> 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)>> nyquist(G);grid>> bode(G)>> nichols(G),grid&g

21、t;> gm,pm,wg,wp=margin(G)gm =pm =wg =wp =>> 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) (s

22、2 + 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);>&