《控制系统计算机仿真》上机报告.docx
兰州理工大学控制系统计算机仿真上机报告院系:申,信学院班级:基地二班姓名:文志浩学号:08220104时间:2011 年5月9日电气工程与信息工程学院2-2用MATLAB语言求下列系统的状态方程、传递函数、零极点增益和部分分式 形式的模型参数,并分别写出其相应的数学模型表达式:()-+ 75- + 245 + 24y = o 254 +1053 +3552 +505 + 24*2.25一5-1.25-0.5''4'2.25-4.25-1.25-0.25X +20.25-0.5-1.25一 121.25-1.75-0.25-0.750(2) X = u0 2x解:G(s) =_ 53+752+ 245 + 24 _54+153+3552 +505 + 241)传递函数转化为状态方程:» num=7,24,24num =72424» den= 10,35,50,24den =10355024»A,B,C,D=tf2ss(num,den)A =-3.5000-5.0000-2.40001.00000001.00000B =100c =0.70002.40002.4000D =0状态方程为:2)传递函数转换为零极点增益:» num=7,24,24num =72424» den= 10,35,50,24den =10355024» Z,P,K=tf2zp(num,den)Z =-1.7143+ 0.6999i-1.7143-0.6999iP =-1.2973+ 0.9838i-1.2973-0.9838i-0.90530.7000零极点增益方程为:3)传递函数转换为部分分时形式:» num=7,24,24num =72424» den= 10,35,50,24den =10355024» R,P,H=residue(num,den)R =-0.0071 -0.2939i-0.0071 + 0.2939i0.7141P =-1.2973+ 0.9838i-1.2973-0.9838i-0.9053H =4）部分分式形式为:"2.25-5-1.25-0.5''42.25 -4.25-1.25-0.252(2) X =0.25 -0.5-1.25一 1X +2U1.25 -1.75-0.25-0.750y = 0 2 0 2X解:1）状态方程转换为传递函数为:» A=2.25 -5 -1.25 -0.5225 -4.25-1.25 -0.25;0.25 -0.5 -1.25-1;1.25-1.75 -0.25 -0.75A =2.2500-5.0000-1.2500-0.50002.2500-4.2500-1.2500-0.25000.2500-0.5000-1.2500-1.00001.2500-1.7500-0.2500-0.7500» B=4;2;2;B =4220» C=0 2 0 2C =0202» D=0D =0» num,den=ss2tf(A,B,C,D) num =04.000014.000022.000015.0000den =1.00004.00006.25005.25002.2500传递函数为: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.75A =2.2500-5.0000-1.2500-0.50002.2500-4.2500-1.2500-0.25000.2500-0.5000-1.2500-1.00001.2500-1.7500-0.2500-0.7500» B=4;2;2;0B =422 0» C=0 2 0 2C =0202» D=0D =0»Z,RK=ss2zp(A,B,C,D)Z =-1.0000+ 1.2247i-1.0000- 1.2247i-1.5000P =-0.5000 + 0.8660i-0.5000 - 0.8660i-1.5000-1.5000K =4.0000零极点增益方程为:3)转换成部分分式形式:» num=0 4 14 22 15num =04142215» den=l 4 6.25 5.25 2.25den =1.00004.00006.25005.25002.2500» R P H=residue(num,den)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用殴拉法求下列系统的输出响应y(。在0 4, W1上,厶= 0.1时的数值解。y=-2y, y(0) = 1要求保留4位小数,并将结果以图形的方式与真解y(r) = e-”比较。t=0:0.1:l;h=0.1;y(1)=1；for i=l：10y(i+1)=y(i)+h*(-2*y (i)；endplot (t,y, *r ')hold onm=exp(-2*t)plot(t,m,bo')2-5用四阶龙格一库塔梯形法求解2-3的数值解,并通过与真值及殴拉法的比较,分析其精度。h=0.1；for i=l：10kl=-2*y(i)k2= -2*(y(i)+kl*h/2)k3=-2*(y(i)+k2*h/2)k4=-2*(y(i)+h*k3)y(i+1)=y(i)+(kl+2*k2+2*k3+k4)*h/6 endplot (t,y, o1)hold onm=exp(-2*t)plot(t,m,r*')lea=y-mplot(t,lea, 1g')hold off4-2设典型闭环结构控制系统如下图所示,当阶跃输入幅值/? = 20时,用 sp4_l. m求取输出y (t)的响应。解 K=：La= 1 54 204.2 213.8 63.5;b= 1875000 1562000;X0= 0,0,0,01 ;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：n+l)；A=rot90(rot90(eye(n-1, n)；-fliplr(A);B= zeros(1,n-1),1'；ml=length(b)；C= fliplr(b),zeros(l,n-ml)；Ab=A-B*C*V;X=X0'；y=0；t=T0;N=round(Tf-TO)/h);for i=l：NKl=Ab*X+B*R;K2=Ab*(X+h*Kl/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« ,y'plot(t,y)4-5下图中,若各环节传递函数已知为:G（s） =1 + 0.015G(s) =1 + 0.1550.210.1（G式s）=繭'G，（s）=时7，G5（s）=雨而，G6（5）=1 + 0.1750.08550.0044 /, 旦1 + 0.015G7（5）= 0.212；试列写链接矩阵W、w。和非零元素阵Wu,将程序sp4一2完善后,应用此程序求输出九的响应曲线。解: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;l,0.01,0.0044,0;W= 0 00000 ;1 0000 -l; 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=l/0/0/0/0/0/n=6 /Y0 = l/Yt0= 0 0 0 0 0 0/h=0.01/Ll=10/T0 = 0/Tf=20/nout=4/A=diag(P(:,1)/B=diag(P(:,2);C=diag(P(：,3)/D=diag(P(:,4)；Q=B-D*W/Qn=inv(Q)/R=C*W-A/V1=C*WO/Ab=Qn*R/bl=Qn*Vl/Y=YtO'/y=Y(nout)/t=T0/N=round(Tf-TO)/(h*Ll)/ for i=l：Nfor j =1:LI/Kl=Ab*Y+bl*Y0K2=Ab*(Y+h*Kl/2)+bl*Y0K3=Ab*(Y+h*K2/2)+bl*Y0K4=Ab*(Y+h*K3)+bl*Y0Y=Y+h*(K1+2*K2+2*K3+K4)/6/ endy=y,Y(nout)/t=t,t(i)+h*Ll/ end f ,y' plot (t,y)4-7用离散相似法仿真程序sp4_3.m重求上题输出H的数据与曲线,并与四阶 龙格一库塔法比较精度。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.01z0.1, 0;1,0.01,0.0044,0;W= 0 00000 ;1 0000 -1; 0 1 0 0 -1 0; 0 0 1 -0.212 0 0 ;0 0 0 1 00;00010 0;W0= l;0;0;0;0;0;n=6;Y0 = l;Yt0= 0 0 0 0 0 0 ;h= 0.01；Ll=10;T0 = 0;Tf=40;nout=5；A=diag(P(：,1)；B=diag(P(：,2)；C=diag(P(：,3)；D=diag(P(：,4)；for i=l：6if(A(i,i)= = 0);FI (i,i)=l;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)=l;FID(i,i)=0;if(D(i,i)=0)；FID(i,i)=D(i,i)/B(i,i)；elseendelseFI(i,i)=exp(-h*A(i,i)/B(iz i)； FIM(i,i)=(l-FI(iri)*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)=l;FID(i,i)=0； if(D(i,i)-=0);F1C(i,i)=C(i,i)/D(i,i)-A(i,i)/B(i,i);FlD(i/i)=D(i,i)/B(i,i)；elseendendendY=zeros(6,1)；X=Y；y=0；Uk=zeros(6,1)；Ub=Uk； t=T0:h*Ll:Tf;N=length(t);for k=l：N-1for 1=1：LIUb=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;endy=y,Y(nout)；endplot(t,y,'r*1)hold onQ=B-D*W;Qn=inv(Q);R=C*W-A;V1=C*WO；Ab=Qn*R；bl=Qn*Vl;Y=YtO'；y=Y(nout);t=T0；N=round(Tf-T0)/(h*Ll);for i=l：Nfor j =1：LI;Kl=Ab*Y+bl*Y0K2=Ab*(Y+h*Kl/2)+bl*Y0K3=Ab*(Y+h*K2/2)+bl*Y0K4=Ab*(Y+h*K3)+bl*Y0Y=Y+h*(K1+2*K2+2*K3+K4)/6;endy=y,Y(nout)；t=t,t(i)+h*Ll；endf ,y'plot (t,y)P=1,0.01,1,0;,0.085,1,0.17;。,0.051,1,0.15;1,.:15,0.21,0;!,0.01,0.1, 0;l,0.01,0.0044,0;W= 0 00000 /1 0000 -1; 0 1 0 0 -1 0/ 0 0 1 -0.212 0 0 ;0 0 0 1 00;00010 0;W0=l;0;0;0;0;0;n=6;Y0 = l;YtO= 0 0 0 0 0 0 ;h=0.01;Ll=10;T0 = 0;Tf=40;nout=5;A=diag (P (：, 1) ) ；B=diag (P (： ,2)；C=diag(P(:,3)；D=diag(P(:,4);for i=l ： 6if(A(i,i)=0);FI(i,i)=l;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)=l;FID(i,i)=0;if(D(i,i)=0)；FID(i,i)=D(i,i)/B(i,i);elseendelseFI(i,i)=exp(-h*A(i,i)/B(i,i);FIM(i,i)=(l-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)=l;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);elseendendendY=zeros(6,1)；X=Y；y=0；Uk=zeros(6,1)；Ub=Uk；t=T0：h*Ll：Tf；N=length(t)；for k=l：N-1for 1=1：LIUb=Uk；Uk=W*Y+WO*YO;Udot= (Uk-Ub) /h； Uf=2*Uk-Ub;X=FI'*X+FIM'*Uk+FIJ1*Udot;Y=FIC *X+FID'*Uf;endy=y,Y(nout)；endplot(t,y,'r*')hold onQ=B-D*W;Qn=inv(Q);R=C*W-A;V1=C*WO;Ab=Qn*R；bl=Qn*Vl；Y=YtO'；y=Y(nout);t=T0；N=round(Tf-TO)/(h*Ll)；for i=l：Nfor j =1:LI;Kl=Ab*Y+bl*YOK2=Ab*(Y+h*Kl/2)+bl*Y0K3=Ab*(Y+h*K2/2)+bl*Y0K4=Ab*(Y+h*K3)+bl*Y0Y=Y+h*(K1+2*K2+2*K3+K4)/6;endy=y,Y(nout)；t=t,t(i)+h*Ll；endIt' ,y'plot(t,y)4-8求下图非线性系统的输出响应y(t),并与无非线性环节情况进行比较。r(t)=10%filename：satu.mfunction Uc=satu(Ur# SI) if (abs(Ur)>=S1)if(Ur>0)Uc=Sl;else Uc=-Sl；end else Uc=Ur； end%filename：dead.m function Uc=dead(Ur,SI) if(abs(Ur)>=S1if(Ur>0)Uc=Ur-Sl； else Uc=Ur+Sl； endelse Uc=0；end%filename：bcaklash.mfunctionUc,Ubb=backlash(Urb,Ur,Ucb, SI) if(Ur>Urb)if(Ur-Sl)>=Ucb) Uc=Ur-Sl；else Uc=Ucb； end else if(Ur<Urb) if(Ur+Sl)<=Ucb) Uc=Ur+Sl;else Uc=Ucb； endendUbb=Ur；P= 0.1 1 0.5 l;0 1 20 0;2 1 1 0;10 110 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;Ll=30;n=4 ；T0 = 0;Tf=20;nout=4；Y0=10;sp4_4;plot(t*g.) hold onZ=4 0 0 0 ；S= 5 0 0 0;sp4_4;plot(try, 'r* ) %filename：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=l：mif(WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,l)rWIJ(k,2)=WIJ(k/3)； end； end； for i=l：nif(A(i,i)=O)；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)=l;FID(i,i)=0;if(D(i,i)=0);FID(i,i)=D(i,i)/B(i,i)；else end elseFI(i,i)=exp(-h*A(i,i)/B(i,i)； FIM(i,i)=(l-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)=l;FID(i,i)=0;if(D(i,i)-=0);FIC(iri)=C(i,i)/D(i,i) -A(i,i)/B(ifi); FID(i,i)=D(i,i)/B(i/i);else end end endY=zeros(nr1)；X=Y；y=0;Uk=zeros(n,1)；Ubb=Uk；t=T0：h*Ll：Tf；N=length(t)；for k=l：N-1for 1=1：LIUb=Uk；Uk=W*Y+W0*Y0； for i=l：nif(Z(i)=0)if(Z(i)=l)Uk(i, i)=satu(Uk(i,i)z 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)=backlash(Ubb(i,i),Uk(iri),Ub(i,i),S(i)； end end end Udot=(Uk-Ub)/h； Uf=2*Uk-Ub;X=FI'*X+FIM1*Uk+FIJ*Udot;Yb=Y;Y=FIC'*X+FID* *Uf;for i=l：nif(Z(i)=0)if(Z(i)=4)Y(i,i)=satu(Y(izi),S(i); endif (Z(i)= = 5)Y(i,i)=dead(Y(i,i),S(i)； endif(Z(i)=6)Y(i,i),Ubb(iri)=backlash(Ubb(i,i),Y(i,i),Yb(i,i),S(i)； endendendendy=y,Y(nout)；end
