neu控制系统仿真CAD作业

1、控制系统计算机辅助设计 第一部分1.M文件：function dx=rossler(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(rossler,0,100,x0); plot(t,y) figure plot3(y(:,1),y(:,2),y(:,3),grid2.m文件function c,ceq=opt_con1(x)ceq=;c=4*x(1)*x(1)+x(2)*x(2)-4;-x(1);-x(2);指令 y=(x)x(1)*x(1)-2*x(1)+2*x(2); x0=0;0;

2、 xm=0;0; xM=; A=;B=;Aeq=;Beq=; x, f_opt=fmincon(y,x0,A,B,Aeq,Beq,xm,xM,opt_con1)x = 1.0000 0f_opt =-13.（a） s=tf(s);G=(s3+4*s+2)/(s3*(s2+2)*(s2+1)3+2*s+5) Transfer 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 funct

3、ion: z2 + 0.568-z3 - 1.2 z2 + 1.19 z - 0.99Sampling time: 0.14. 手工运算整理得令x1=y ; x2=dx1/dt ; x3=dx2/dt ; dx3/dt=2*u-5*x1-4*x2-13*x3;A=0,1,0;0,0,1;-5,-4,-13；即可得出状态方程：dx1/dt ;dx2/dt ;dx3/dt =A*x1;x2;x3+0;0;2;将上述状态方程输入到MATLAB命令窗口；调用tf（G），zpk（G）指令即可得出传递函数和零极点模型 A=0 1 0;0 0 1;-5,-4,-13;B=0;0;2;C=1 0 0;D=0;

4、 G=ss(A,B,C,D)a = x1 x2 x3 x1 0 1 0 x2 0 0 1 x3 -5 -4 -13b = u1 x1 0 x2 0 x3 2 c = x1 x2 x3 y1 1 0 0d = u1 y1 0 Continuous-time model. tf(G) Transfer function: 2-s3 + 13 s2 + 4 s + 5 zpk(G)Zero/pole/gain: 2-(s+12.72) (s2 + 0.2836s + 0.3932)当然可以由微分方程直接得出传递函数；用拉普拉斯变换法可直接输入： G=tf(2,1 13 4 5)Transfer fu

5、nction: 2-s3 + 13 s2 + 4 s + 55. 采用z变换方法；设采样周期为1s z=tf(z,1);G=(z+2)/(z2+z+0.016) Transfer function: z + 2-z2 + z + 0.016 Sampling time: 16. syms G J Gc Kp Ki s G=(s+1)/(J*s2+2*s+5);Gc=(Kp*s+Ki)/s; GG=feedback(G*Gc,1) GG = (s+1)*(Kp*s+Ki)/(J*s3+2*s2+5*s+Kp*s2+s*Ki+Kp*s+Ki)7.(a)G=tf(211.87,317.64,conv

6、(conv(1,20,1,94.34),1,0.1684);Gc=tf(169.6,400,1,4,0);H=tf(1,0.01,1); GG=feedback(G*Gc,H)Transfer function: 359.3 s3 + 3.732e004 s2 + 1.399e005 s + -0.01s6 + 2.185 s5 + 142.1 s4 + 2444 s3 + 4.389e004 s2 + 1.399e005 s + zpk(GG) Zero/pole/gain: 35933.152 (s+100) (s+2.358) (s+1.499)-0.02(s2 + 3.667s + 3

7、.501) (s2 + 11.73s + 339.1) (s2 + 203.1s + 1.07e004)(b)z=tf(z,1);G=(35786.7*z-1+)/(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: - z6 + 1.844e004 z5 + 1.789e004 z4-1.144e005 z6 + 2.876e004 z5 + 274.2 z4 + 782.4 z3 + 47.52 z2 + 0.5 zSampling time: 1 z

8、pk(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: 18. g1=tf(1,1,1);gc1=tf(1,0,1,0,2);h1=tf(4,2,1,2,1);G=feedback(g1*gc1,h1) Transfer function: s3 + 2 s2 + s-s5 + 3 s4 + 5 s3 + 11 s2 + 8 s + 2 g2=tf(1,1,0,0);h2=50;Gc=feedback(

9、g2,h2) Transfer function: 1-s2 + 50 H=tf(1,0,2,1,0,0,14);GG=3*feedback(G*Gc,H)Transfer function: 3 s6 + 6 s5 + 3 s4 + 42 s3 + 84 s2 + 42 s-s10 + 3 s9 + 55 s8 + 175 s7 + 300 s6 + 1323 s5 + 2656 s4 + 3715 s3 + 7732 s2 + 5602 s + 14009.num=conv(conv(1,1,1,1),1,2,400);den=conv(conv(conv(1,5,1,5),1,3,100

10、),1,3,2500);G=tf(num,den); step(G) grid; hold on; G1=c2d(G,0.01); step(G1) hold on; G2=c2d(G,0.1);step(G2) hold on; G3=c2d(G,1);step(G3)10. （a） G=tf(1,1,2,1,2);eig(G)ans = -2.0000 0.0000 + 1.0000i 0.0000 - 1.0000i临界稳定（b） G=tf(1,6,3,2,1,1);eig(G)ans =-0.4949 + 0.4356i - 0.4949 - 0.4356i 0.2449 + 0.56

11、88i 0.2449 - 0.5688i不稳定（c） G=tf(1,4,1,-3,-1,2);eig(G)ans = -0.8057 + 0.4664i -0.8057 - 0.4664i 0.6807 + 0.3372i 0.6807 - 0.3372i不稳定11（a） z=tf(z,1);H=(-3*z+2)/(z3-0.2*z2-0.25*z+0.05);abs(eig(H)ans = 0.5000 0.5000 0.2000稳定（b） z=tf(z,1);H=(3*z-0.39*z-0.09)/(z4-1.7*z3+1.04*z2+0.268*z+0.024);abs(eig(H)an

12、s = 1.1939 1.1939 0.1298 0.1298不稳定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; rank(ctrb(A,B) ans = 5 完全可控13显然，这个自治的微分方程组的解析解为x(t) = x(0) * expm(A*t) A=-5 2 0 0;0 -4 0 0 ;-3 2 -4 -1;-3 2 0 -4;x0=1 2 0 1; syms s t;x=expm(A*t)*x0x = -3*exp(-5*t)+4*exp(

13、-4*t) 2*exp(-4*t) 18*exp(-4*t)-18*exp(-5*t)-18*t*exp(-4*t)+4*t2*exp(-4*t) 10*exp(-4*t)-9*exp(-5*t)-8*t*exp(-4*t)由于符号运算不能图解，将x分成四个表达式，用离散化的t向量为为自变量，绘图 t=0:0.02:2; x1= -3*exp(-5*t)+4*exp(-4*t); x2=2*exp(-4*t);x3=18*exp(-4*t)-18*exp(-5*t)-18*t.*exp(-4*t)+4*t.2.*exp(-4*t);x4=10*exp(-4*t)-9*exp(-5*t)-8*t

14、.*exp(-4*t); plot(t,x1,t,x2,t,x3,t,x4) ， grid数值解法： f=(t,x)-5 2 0 0; 0 -4 0 0 ; -3 2 -4 -1; -3 2 0 -4 * x;t,y=ode45(f,0,2,1 2 0 1); plot(t,y)可见，数值解与解析解很接近。14.（a） G=tf(conv(1,6,1,-6),conv(conv(conv(1,0,1,3),1,4-4*j),1,4+4*j) rlocus(G)增益Gain从0+无穷，在右半平面都存在不稳定极点，因此不论开环增益为多大，系统始终不稳定。不妨当增益为1时，求闭环系统的特征根。 GG

15、=feedback(G,1); eig(GG)ans = -3.9404 - 4.2387i -3.9404 + 4.2387i -3.4323 0.3132 （右半平面的极点，不稳定） （b） G=tf(1,2,2,1,1,14,8);rlocus(G) 从根轨迹可以看出，所有根轨迹均在s左半平面，因此无论增益为多大，闭环系统始终稳定15. s=tf(s);G=(s-1)/(s+1)5;G.ioDelay=2; G=pademod(G,1,2)；rlocus(G) 为了方便看清晰临界稳定增益，此处截取了部分根轨迹，从图中可见，当增益Gain s=tf(s);G=8*(s+1)/s2/(s+1

16、5)/(s2+6*s+10);nyquist(G) bode(G), grid; gm,pm,wg,wp=margin(G)gm = 30.4686 pm = 4.2340 wg =1.5811 Wp = 0.2336 nichols(G) 以下是局部Nichols图，上面反映了相位裕度和幅值裕度相位裕量和幅值裕量均大于零，因此单位反馈构成的闭环系统稳定 step(feedback(G,1)（b） z=tf(z,1);H=0.45*(z+1.31)*(z+0.45)*(z-0.957)/z/(z-1)/(z-0.368)/(z-0.99); nyquist(H) bode(H)，grid on

17、 gm,pm,wg,wp=margin(H)Warning: The closed-loop system is unstable.gm = 0.5434 pm = -26.6555 wg = 0.8620 wp =1.1975相位裕量小于零，因此单位反馈构成的闭环系统不稳定。 step(feedback(H,1)17G=100*tf(1/2.5,1,conv(1/0.5,1,0,1/50,1); Gc=1000*tf(conv(1,1,1,2.5),conv(1,0.5,1,50); nichols(G*Gc) gm,pm,wg,wp=margin(G)gm = Inf pm =63.295

18、3 wg = Inf wp = 18.8688 由以上图形和数据，可判断出单位负反馈系统稳定以下是单位阶跃响应： GG=feedback(G*Gc,1);step(GG) 第二部分2. sim(simu2,100),plot(tout,yout) f=(t,y)y(2);y(3);y(4);-2*y(1)-4*y(2)-63*y(3)-5*y(4)+exp(-3*t)+exp(-5*t)*sin(4*t+pi/3); t,y=ode45(f,0,100,1 0.5 0.5 0.2);plot(t,y(:,1)3Simulink仿真模型 plot(tout,yout)4换成一阶微分方程组的形式如

19、下： dx1/dt=t*sin(x2*exp(-2.3*x4);dx2/dt=x1;dx3/dt=sin(x2*exp(-2.3*x4);dx4/dt=x3;5系统仿真模型A,B,C,D=linmod(test5);G=ss(A,B,C,D);subplot(221),step(G),grid,subplot(222),bode(G),grid,subplot(223),nyquist(G),grid,subplot(224),nichols(G),grid阶跃响应与频率响应曲线如下：6 G=210*tf(1 1.5,conv(conv(1 1.75,1 16),conv(1 1.5+3*j,

20、1 1.5-3*j); GG=feedback(G*Gc,1); step(feedback(G,1) gm,gp,wm,wp=margin(G)gm =4.8921 gp =60.0634 wm =7.9490 wp =3.9199 step(GG) gm1,gp1,wm1,wp1=margin(GG)Warning: The closed-loop system is unstable. In warning at 26 In lti.margin at 67gm1 =Inf gp1 =-43.3034 wm1 =NaN wp1 = 24.60597. A=0 1 0 0;0 0 1 0;-3 1 2 3;2 1 0 0;B=1 2 3 4;0 1 2 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