控制系统计算机仿真上机报告.doc

*理 工 大 学控制系统计算机仿真上机报告院系： 班级： 姓名： 学号： 时间： 年 月 日控制系统计算机仿真上机报告2-2用MATLAB语言求下列系统的状态方程、传递函数、零极点增益和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：（1）（2） 解：（1）源程序：num=1 7 24 24;den=1 10 35 50 24;G=tf(num,den)Z,P,K=tf2zp(num,den)A,B,C,D=tf2ss(num,den)R,P,H=residue(num,den)结果：Transfer function: s3 + 7 s2 + 24 s + 24-s4 + 10 s3 + 35 s2 + 50 s + 24Z = -2.7306 + 2.8531i -2.7306 - 2.8531i -1.5388 P = -4.0000 -3.0000 -2.0000 -1.0000K = 1A = -10 -35 -50 -24 1 0 0 0 0 1 0 0 0 0 1 0B = 1 0 0 0C = 1 7 24 24D = 0R = 4.0000 -6.0000 2.0000 1.0000P = -4.0000 -3.0000 -2.0000 -1.0000H = （2）源程序：A=2.25 -5 -1.25 -0.5;2.25 -4.25 -1.25 -0.25;0.25 -0.5 -1.25 -1;1.25 -1.75 -0.25 -0.75B=4;2;2;0C=0 2 0 2D=0num,den=ss2tf(A,B,C,D)Z,P,K=tf2zp(num,den)R,P,H=residue(num,den)结果：A = 2.2500 -5.0000 -1.2500 -0.5000 2.2500 -4.2500 -1.2500 -0.2500 0.2500 -0.5000 -1.2500 -1.0000 1.2500 -1.7500 -0.2500 -0.7500B = 4 2 2 0C = 0 2 0 2D = 0num = 0 4.0000 14.0000 22.0000 15.0000den = 1.0000 4.0000 6.2500 5.2500 2.2500Z = -1.0000 + 1.2247i -1.0000 - 1.2247i -1.5000 P = -1.5000 -1.5000 -0.5000 + 0.8660i -0.5000 - 0.8660iK = 4.0000R = 4.0000 0.0000 0.0000 - 2.3094i 0.0000 + 2.3094iP = -1.5000 -1.5000 -0.5000 + 0.8660i -0.5000 - 0.8660iH = 2-3 用殴拉法求下列系统的输出响应在上，时的数值解。，要求保留4位小数，并将结果以图形的方式与真解比较。解：源程序：t=0:0.1:1;h=0.1;y(1)=1;for i=1:1:10; y(i+1)=y(i)+h*(-2*y(i); Y=y,y(i)endplot(t,y,r*);hold onm=exp(-2*t);plot(t,m)结果：Y = Columns 1 through 8 1.0000 0.8000 0.6400 0.5120 0.4096 0.3277 0.2621 0.2097Columns 9 through 12 0.1678 0.1342 0.1074 0.1342分析：经由图像比较，用欧拉法求得的解与真值大体一致。2-5 用四阶龙格库塔梯形法求解2-3的数值解，并通过与真值及殴拉法的比较，分析其精度。解：源程序：t=0:0.1:1;h=0.1;y(1)=1;for i=1:1:10; k1=-2*(y(i); k2=-2*(y(i)+h/2*k1); k3=-2*(y(i)+h/2*k2); k4=-2*(y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4); Y=y,y(i)endplot(t,y,r*);hold onm=exp(-2*t);plot(t,m)结果：Y = Columns 1 through 8 1.0000 0.8187 0.6703 0.5488 0.4493 0.3679 0.3012 0.2466 Columns 9 through 12 0.2019 0.1653 0.1353 0.1653分析：经由图像比较，用四阶龙格库塔梯形法求解题2-3的数值解与真值一致，比用欧拉法求得的解更接近于真值。4-2设典型闭环结构控制系统如下图所示，当阶跃输入幅值时，用sp4_1.m求取输出y(t)的响应。y(t)r(t)_0.002解：源程序：b=1.875e+6,1.562e+6;a=1,54,204.2,213.8,63.5;X0=0,0,0,0;V=0.002;n=4;T0=0;Tf=10;h=0.01;R=20;b=b/a(1);a=a/a(1);A=a(2:5);A=rot90(rot90(eye(3,4);-fliplr(A);B=zeros(1,3),1;m1=length(b);C=fliplr(b),zeros(1,4-m1);Ab=A-B*C*V;X=X0;y=0;t=T0;N=round(Tf-T0)/h;for i=1:N; K1=Ab*X+B*R K2=Ab*(X+h*K1/2)+B*R; K3=Ab*(X+h*K2/2)+B*R; K4=Ab*(X+h*K3)+B*R; X=X+h*(K1+2*K2+2*K3+K4)/6; y=y,C*X; t=t,t(i)+h;endt,yplot(t,y)b=1.875e+6,1.562e+6;a=1,54,204.2,213.8,63.5;X0=0,0,0,0;V=0.002;n=4;T0=0;Tf=10;h=0.01;R=20;b=b/a(1);a=a/a(1);A=a(2:5);A=rot90(rot90(eye(3,4);-fliplr(A);B=zeros(1,3),1;m1=length(b);C=fliplr(b),zeros(1,4-m1);Ab=A-B*C*V;X=X0;y=0;t=T0;N=round(Tf-T0)/h;for i=1:N; K1=Ab*X+B*R K2=Ab*(X+h*K1/2)+B*R; K3=Ab*(X+h*K2/2)+B*R; K4=Ab*(X+h*K3)+B*R; X=X+h*(K1+2*K2+2*K3+K4)/6; y=y,C*X; t=t,t(i)+h;endt,yplot(t,y)结果：分析：由结论及图形反映了开环放大洗漱K与反馈系数V对结果的影响，K值增大会使系统震荡加剧，甚至可能造成不稳定。4-5 下图中，若各环节传递函数已知为：，但；试列写链接矩阵W、W0和非零元素阵WIJ，将程序sp4_2完善后，应用此程序求输出的响应曲线。G2(s)G7(s)G3(s)y0_y4_G1(s)G4(s)G5(s)G6(s)解：源程序：P=1,0.01,1,0; 0,0.085,1,0.017; 0,0.051,1,0.15; 1,0.15,0.21,0; 1,0.01,0.1,0; 1,0.01,0.0044,0; ;WIJ=1 0 1; 2 1 1; 2 6 -1; 3 2 1; 3 5 -1; 4 3 1; 4 4 -0.212; 5 4 1; 6 4 1 ;n=6;Y0=1;Yt0=0 0 0 0 0 0 ;h=0.01;L1=10;T0=0;Tf=40;nout=5;A=diag(P(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:m if(WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); end;end;Q=B-D*W;Qn=inv(Q);R=C*W-A;V1=C*W0;Ab=Qn*R;b1=Qn*V1;Y=Yt0;y=Y(nout);t=T0;N=round(Tf-T0)/(h*L1);for i=1:N for j=1:L1; K1=Ab*Y+b1*Y0 K2=Ab*(Y+h*K1/2)+b1*Y0 K3=Ab*(Y+h*K2/2)+b1*Y0 K4=Ab*(Y+h*K3)+b1*Y0 Y=Y+h*(K1+2*K2+2*K3+K4)/6; end y=y,Y(nout); t=t,t(i)+h*L1;endt,yplot(t,y)结果：4-7用离散相似法仿真程序sp4_3.m重求上题输出的数据与曲线，并与四阶龙格库塔法比较精度。解：源程序：P=1,0.01,1,0;0,0.085,1,0.17;0,0.051,1,0.15;1,0.15,0.21,0;1,0.01,0.1,0;1,0.01,0.0044,0;W=0 0 0 0 0 0 ;1 0 0 0 0 -1; 0 1 0 0 -1 0; 0 0 1 -0.212 0 0 ;0 0 0 1 0 0 ;0 0 0 1 0 0;W0=1;0;0;0;0;0;n=6;Y0=1;Yt0=0 0 0 0 0 0 ;h=0.01;L1=10;T0=0;Tf=40;nout=4;A=diag(P(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);for i=1:6 if(A(i,i)=0); FI(i,i)=1; FIM(i,i)=h*C(i,i)/B(i,i); FIJ(i,i)=h*h*C(i,i)/B(i,i)/2; FIC(i,i)=1;FID(i,i)=0; if(D(i,i)=0); FID(i,i)=D(i,i)/B(i,i); else end else FI(i,i)=exp(-h*A(i,i)/B(i,i); FIM(i,i)=(1-FI(i,i)*C(i,i)/A(i,i); FIJ(i,i)=h*C(i,i)/A(i,i)-FIM(i,i)*B(i,i)/A(i,i); FIC(i,i)=1;FID(i,i)=0; if(D(i,i)=0); FIC(i,i)=C(i,i)/D(i,i)-A(i,i)/B(i,i); FID(i,i)=D(i,i)/B(i,i); else end end end Y=zeros(6,1);X=Y;y=0;Uk=zeros(6,1);Ub=Uk; t=T0:h*L1:Tf;N=length(t);for k=1:N-1 for l=1:L1 Ub=Uk; Uk=W*Y+W0*Y0; Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI*X+FIM*Uk+FIJ*Udot; Y=FIC*X+FID*Uf; end y=y,Y(nout);endplot(t,y,r*)结果：分析：用离散相似法与四阶龙格库塔法比较来看，两者精度很接近，并且与四阶龙格-库塔法相比，离散相似法更容易推广到解决非线性系统的仿真问题。4-8求下图非线性系统的输出响应y(t)，并与无非线性环节情况进行比较。y(t)e(t)r(t)=105-5解：源程序：主程序：P=0.1 1 0.5 1;0 1 20 0;2 1 1 0;10 1 1 0; WIJ=1 0 1;1 4 -1 ;2 1 1;3 2 1;4 3 1;Z=0 0 0 0;S=0 0 0 0; h=0.01;L1=25; n=4; T0=0;Tf=20; nout=4; Y0=10;sp4_4;plot(t,y,r)hold onZ=4 0 0 0;S=5 0 0 0;sp4_4;plot(t,y,g)子程序sp4_4.m文件：A=diag(P(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:m if(WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); end;end;for i=1:n if(A(i,i)=0); FI(i,i)=1; FIM(i,i)=h*C(i,i)/B(i,i); FIJ(i,i)=h*h*C(i,i)/B(i,i)/2; FIC(i,i)=1;FID(i,i)=0; if(D(i,i)=0); FID(i,i)=D(i,i)/B(i,i); else end else FI(i,i)=exp(-h*A(i,i)/B(i,i); FIM(i,i)=(1-FI(i,i)*C(i,i)/A(i,i); FIJ(i,i)=h*C(i,i)/A(i,i)-FIM(i,i)*B(i,i)/A(i,i); FIC(i,i)=1;FID(i,i)=0; if(D(i,i)=0); FIC(i,i)=C(i,i)/D(i,i)-A(i,i)/B(i,i); FID(i,i)=D(i,i)/B(i,i); else end endendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;t=T0:h*L1:Tf;N=length(t);for k=1:N-1 for i=1:L1 Ub=Uk; Uk=W*Y+W0*Y0; for i=1:n if(Z(i)=0) if(Z(i)=1) Uk(i,i)=satu(Uk(i,i),S(i); end if(Z(i)=2) Uk(i,i)=dead(Uk(i,i),S(i); end if(Z(i)=3) Uk(i,i),Ubb(i,i)=b