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

1、东北大学11控制系统计算机辅助设计MATLAB语言与应用控制系统计算机辅助设计第一部分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

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

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

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

5、+ 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

6、 + 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.389e

7、004 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: -10844

8、4 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/

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

10、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.01step(G1) G2=c2d(G,0.1)Transfer function:0.0003982 z5 - 0.0003919 z4 - 0.000336 z3 + 0.0007842 z2 - 0.0

11、00766 z + 0.0003214-z6 - 2.644 z5 + 4.044 z4 - 3.94 z3 + 2.549 z2 - 1.056 z + 0.2019Sampling time: 0.1step(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

12、 + 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.3

13、9*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 -

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

15、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=A

16、a 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

17、=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,xx0ans = -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

18、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 s=tf(s); G=(s-1)/(s+1)5; G.ioDelay=2Transfer function: s - 1exp(-2*s) * - s5 + 5 s4 + 10 s3 + 10

19、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)通过以上的分析

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)Sampling time: 0.1 nyquist(G);grid bode(G) nichols(G),grid gm,pm,wg,wp=margin(G)gm = 0.9578pm = -1.7660wg = 10.4

21、641wp = 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;-50Gc=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

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