控制新版系统计算机仿真作业

兰州理工大学《控制系统计算机仿真》上机汇报Ⅰ院系：电气工程和信息工程学院班级：14级自动化3班姓名：孙悦学号：时间：年6月15日电气工程和信息工程学院《控制系统计算机仿真》上机试验任务书Ⅰ（）一、上机试验内容及要求matlab软件要求利用课余时间熟悉掌握matlab软件基础数值运算、基础符号运算、基础程序设计方法及常见图形命令操作。2.各章节仿真试验内容及要求具体试验内容及要求请详见上机试验汇报。二、上机试验时间安排及相关事宜依据课程教学纲领要求，上机试验课时共16课时，学生须在每次上机之前做好对应准备工作，以确保在有限机时内完成仿真试验要求内容；试验完成后按要求完成相关仿真试验汇报；仿真试验汇报请按相关样本制作并A4打印，侧面装订，作为成绩评定一部分。自动化系《控制系统计算机仿真》课程组3月一、Matlab基础操作1-1用MATLAB语言求下列系统状态方程、传输函数、零极点增益和部分分式形式模型参数，并分别写出其对应数学模型表示式：（1）程序以下:num=[7,24,24]den=[10,35,50,24][A,B,C,D]=tf2ss(num,den)系统状态方程：A=-3.5000-5.0000-2.40001.00000001.00000B=100C=0.70002.40002.4000D=0零极点增益形式：[Z,P,K]=tf2zp(num,den)Z=-1.7143+0.6999i-1.7143-0.6999iP=-1.2973+0.9838i-1.2973-0.9838i-0.9053K=0.7000部分分式：[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=[]（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.75]A=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;0]B=4220>>C=[0202]C=0202>>D=0D=0零极点增益形式：>>[Z,P,K]=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+0.0000i-1.5000-0.0000iK=4.0000传输函数形式：>>num=[04142215]num=04142215>>den=[146.255.252.25]den=1.00004.00006.25005.25002.2500部分分式：>>[R，P，H]=residue(num,den)R=4.0000-0.00000.0000-2.3094i0.0000+2.3094iP=-1.5000-1.5000-0.5000+0.8660i-0.5000-0.8660iH=[]1-2用殴拉法matlab编程实现下列系统输出响应在上，时数值解。，要求保留4位小数，并将结果以图形方法和真解比较。t=0:0.1:1t=Columns1through900.10000.0.30000.40000.50000.60000.70000.8000Columns10through110.90001.0000h=0.1;y(1)=1;t=0:0.1:1;h=0.1;y(1)=1;fori=1:10y(i+1)=y(i)+h*(-1*y(i));endplot(t,y,'r')holdonm=exp(-1*t)m=Columns1through91.00000.90480.81870.74080.67030.60650.54880.49660.4493Columns10through110.40660.3679plot(t,m,'bo')1-3用四阶龙格－库塔梯形法matlab编程实现1-2题数值解，要求以图形方法经过和真值及殴拉法比较，分析其精度。h=0.1;y(1)=1;t=0:0.1:1t=Columns1through900.10000.0.30000.40000.50000.60000.70000.8000Columns10through110.90001.0000fori=1:10k1=-1*y(i)k2=-1*(y(i)+k1*h/2)k3=-1*(y(i)+k2*h/2)k4=-1*(y(i)+h*k3)y(i+1)=y(i)+(k1+2*k2+2*k3+k4)*h/6endk1=-1k2=-0.9500k3=-0.9525k4=-0.9047y=1.00000.9048k1=-0.9048k2=-0.8596k3=-0.8619k4=-0.8187y=1.00000.90480.8187k1=-0.8187k2=-0.7778k3=-0.7798k4=-0.7407y=1.00000.90480.81870.7408k1=-0.7408k2=-0.7038k3=-0.7056k4=-0.6703y=1.00000.90480.81870.74080.6703k1=-0.6703k2=-0.6368k3=-0.6385k4=-0.6065y=1.00000.90480.81870.74080.67030.6065k1=-0.6065k2=-0.5762k3=-0.5777k4=-0.5488y=1.00000.90480.81870.74080.67030.60650.5488k1=-0.5488k2=-0.5214k3=-0.5227k4=-0.4965y=1.00000.90480.81870.74080.67030.60650.54880.4966k1=-0.4966k2=-0.4718k3=-0.4730k4=-0.4493y=1.00000.90480.81870.74080.67030.60650.54880.49660.4493k1=-0.4493k2=-0.4269k3=-0.4280k4=-0.4065y=Columns1through91.00000.90480.81870.74080.67030.60650.54880.49660.4493Column100.4066k1=-0.4066k2=-0.3862k3=-0.3873k4=-0.3678y=Columns1through91.00000.90480.81870.74080.67030.60650.54880.49660.4493Columns10through110.40660.3679plot(t,y,'o')holdonm=exp(-1*t)plot(t,m,'r*')lea=y-mplot(t,lea,'g')holdoffm=Columns1through81.00000.90480.81870.74080.67030.60650.54880.4966Columns9through110.44930.40660.3679lea=1.0e-006*Columns1through800.08200.14830.0.24290.27470.29830.3149Columns9through110.32560.33150.33321-4采取matlab语言编程实现。程序:disp('y=')h=1;y=0;fori=0:1:63y=y+2^i;enddisp(y);运行结果：y=1.8447e+019y=1.8447e+0191-5编写matlabM-函数，以实现。要求在函数中给出必需解释和说明，同时检测输入和返回变量个数。functionsum=zuoye5(m)sum=0;formatlongfori=0:mt=1;forj=1:it=t*2;endsum=sum+t;end二、控制系统分析2-1设经典闭环结构控制系统以下图所表示，当阶跃输入幅值时，用sp3_1.m求取输出y(t)响应。深入考虑：当反馈通道为时，怎样经过对sp3_1.m程序修改，以实现此时所对应y(t)。y(t)y(t)r(t)_2解：程序1：clearall;closeall;a=[0.0160.8643.273.421];b=[3025];V=2;R=20;X0=[0000];n=4;T0=0;Tf=10;h=0.01;b=b/a(1);a=a/a(1);A=a(2:n+1);A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];B=[zeros(1,n-1),1]';m1=length(b);C=[fliplr(b),zeros(1,n-m1)];Ab=A-B*C*V;X=X0';y=0;t=T0;N=round(Tf-T0)/h;fori=1:NK1=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+K2*2+K3*2+K4)/6;y=[y,C*X];t=[t,t(i)+h];end[t',y']plot(t,y,'b-')title('反馈系数为2时阶跃响应曲线');holdon;gridon程序2：clearall;closeall;a=[0.0160.8643.273.421];b=[3025];a1=conv(a,[0.11]);b1=b;b1=b1/a1(1);a1=a1/a1(1);m1=length(b1);b=b/a(1);a=a/a(1);m=length(b);V=1;R=20;T0=0;Tf=15;h=0.01;X10=[00000];n1=5;X0=[0000];n=4;A1=a1(2:n1+1);A=a(2:n+1);C1=[fliplr(b1),zeros(1,n1-m1)];C=[fliplr(b),zeros(1,n-m)];A1=[rot90(rot90(eye(n1-1,n1)));-fliplr(A1)];A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];B1=[zeros(1,n1-1),1]';B=[zeros(1,n-1),1]';Ab=A1-B1*C1*V;X1=X10';y1=0;t=T0;N=round(Tf-T0)/h;fori=1:N;K1=Ab*X1+B1*R;K2=Ab*(X1+h*K1/2)+B1*R;K3=Ab*(X1+h*K2/2)+B1*R;K4=Ab*(X1+h*K3)+B1*R;X1=X1+h*(K1+K2*2+K3*2+K4)/6;y1=[y1,C1*X1];t=[t,t(i)+h];endX=X0';y=0;t=T0;fori=1:NK1=A*X+B*(R-y1(i));K2=A*(X+h*K1/2)+B*(R-y1(i));K3=A*(X+h*K2/2)+B*(R-y1(i));K4=A*(X+h*K3)+B*(R-y1(i));X=X+h*(K1+K2*2+K3*2+K4)/6;y=[y,C*X];t=[t,t(i)+h];end[t',y']plot(t,y,'r-');title('反馈通道为惯性步骤时阶跃响应曲线');holdon;gridon2-2下图中，若各步骤传输函数已知为：，，，，，，，，，；试列写链接矩阵W、W0和非零元素阵WIJ，将程序sp4_2完善后，应用此程序求输出响应曲线。GG2(s)G5(s)G10(s)G4(s)y0_y7_G1(s)G3(s)G6(s)G7(s)G8(s)G9(s)解：程序clearall;closeall;a=[101011011]';b=[0.010.0850.010.0510.00670.1510.010.01]';c=[1111700.211300.10.0044]';d=[00.1700.1500000]';P=[a,b,c,d];WIJ=[101;211;29-1;321;431;48-1;541;651;67-0.212;761;861;971];n=9;Y0=1;Yt0=[000000000];h=0.001;L1=10;T0=0;Tf=3;nout=7;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);fork=1:mif(WIJ(k,2)==0)W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);endendQ=B-D*W;R=C*W-A;V1=C*W0;Ab=Q\R;b1=Q\V1;Y=Yt0';y=Y(nout);t=T0;N=round((Tf-T0)/(h*L1));fori=1:Nforj=1:L1K1=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+K2*2+K3*2+K4)/6;endy=[y,Y(nout)];t=[t,t(i)+h*L1];end[t',y']plot(t,y,'r');title('单位阶跃响应曲线');gridon2-3用离散相同法仿真程序sp4_3.m重求上题输出数据和曲线，并和四阶龙格－库塔法比较精度，同时分析两种方法对仿真步长选择区分。解：程序clearall;closeall;a=[101011011]';b=[0.010.0850.010.0510.00670.1510.010.01]';c=[1111700.211300.10.0044]';d=[00.1700.1500000]';P=[a,b,c,d];WIJ=[101;211;29-1;321;431;48-1;541;651;67-0.212;761;861;971];n=9;Y0=1;h=0.001;L1=10;T0=0;Tf=3;nout=7;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1));W0=zeros(n,1);W=zeros(n,n);fork=1:mif(WIJ(k,2)==0)W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);endendfori=1:nif(A(i)==0);FI(i)=1;FIM(i)=h*C(i)/B(i);FIJ(i)=h*h*C(i)/B(i)/2;FIC(i)=1;FID(i)=0;if(D(i)~=0)FID(i)=D(i)/B(i);elseendelseFI(i)=exp(-h*A(i)/B(i));FIM(i)=(1-FI(i))*C(i)/A(i);FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i);FIC(i)=1;FID(i)=0;if(D(i)~=0)FIM(i)=(1-FI(i))*D(i)/A(i)FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)/A(i)FIC(i)=C(i)/D(i)-A(i)/B(i);FID(i)=D(i)/B(i);elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;t=T0:h*L1:Tf;N=length(t);fork=1:N-1forl=1:L1Ub=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)];end[t',y']plot(t,y)2-4求下图非线性系统输出响应y(t)，并和无非线性步骤情况进行比较。y(t)y(t)e(t)r(t)=105-5P=[0.110.51;01200;2110;10110];WIJ=[101;14-1;211;321;431];Z=[0000];S=[0000];h=0.01;L1=25;n=4;T0=0;Tf=20;nout=4;Y0=10;sp4_4;plot(t,y,'r')holdonZ=[4000];S=[5000];sp4_4;plot(t,y,'g')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);fork=1:mif(WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);end;end;fori=1:nif(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);elseendelseFI(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);elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;t=T0:h*L1:Tf;N=length(t);fork=1:N-1fori=1:L1Ub=Uk;Uk=W*Y+W0*Y0;fori=1:nif(Z(i)~=0)if(Z(i)==1)Uk(i,i)=satu(Uk(i,i),S(i));endif(Z(i)==2)Uk(i,i)=dead(Uk(i,i),S(i));endif(Z(i)==3)[Uk(i,i),Ubb(i,i)]=backlash(Ubb(i,i),Uk(i,i),Ub(i,i),s(i));endendendUdot=(Uk-Ub)/h;Uf=2*Uk-Ub;X=FI'*X+FIM'*Uk+FIJ'*Udot;Yb=Y;Y=FIC'*X+FID'*Uf;fori=1;nif(Z(i)~=0)if(Z(i)==4)Y(i,i)=satu(Y(i,i),S(i));endif(Z(i)==5)Y(i,i)=dead(Y(i,i),S(i));endif(Z(i)==6)[Y(i,i),Ubb(i,i)]=backlash(Ubb(i,i),Y(i,i),Yb(i,i),s(i));endendendendy=[y,Y(nout)];endbacklash函数:function[Uc,Ubb]=backlash(Urb,Ur,Ucb,S1)if(Ur>Urb)if((Ur-S1)>=Ucb)Uc=Ur-S1;elseUc=Ucb;endelseif(Ur<Urb)if((Ur+S1)<=Ucb)Uc=Ur+S1;elseUc=Ucb;endelseUc=Ucb;endendUbb=Ur;Satu函数: