控制系统计算机仿真作业(29页).doc
-兰 州 理 工 大 学控制系统计算机仿真上机报告院系： 电气工程与信息工程学院 班级： 14级自动化3班 姓名： 孙悦 学号： 1405220323 时间： 2017 年 6 月 15 日电气工程与信息工程学院控制系统计算机仿真上机实验任务书（2017）一、上机实验内容及要求1. matlab软件要求利用课余时间熟悉掌握matlab软件的基本数值运算、基本符号运算、基本程序设计方法及常用的图形命令操作。2.各章节仿真实验内容及要求 具体实验内容及要求请详见上机实验报告。二、上机实验时间安排及相关事宜1. 依据课程教学大纲要求，上机实验学时共16学时，学生须在每次上机之前做好相应的准备工作，以确保在有限的机时内完成仿真实验要求的内容；2. 实验完成后按规定完成相关的仿真实验报告；3. 仿真实验报告请按有关样本制作并A4打印，侧面装订，作为成绩评定的一部分。自动化系控制系统计算机仿真课程组2017年3月-第 27 页-一、Matlab基础操作1-1用MATLAB语言求下列系统的状态方程、传递函数、零极点增益和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：（1）程序如下: num=7,24,24den=10,35,50,24 A,B,C,D=tf2ss(num,den)系统的状态方程：A = -3.5000 -5.0000 -2.4000 1.0000 0 0 0 1.0000 0B = 1 0 0C = 0.7000 2.4000 2.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.9053 K = 0.7000部分分式： R,P,H=residue(num,den)R = -0.0071 - 0.2939i -0.0071 + 0.2939i 0.7141 P = -1.2973 + 0.9838i -1.2973 - 0.9838i -0.9053 H = （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.5000 2.2500 -4.2500 -1.2500 -0.2500 0.2500 -0.5000 -1.2500 -1.0000 1.2500 -1.7500 -0.2500 -0.7500>> B=4;2;2;0B = 4 2 2 0>> C=0 2 0 2C = 0 2 0 2>> D=0D = 0零极点增益形式：>> Z,P,K=ss2zp(A,B,C,D)Z = -1.0000 + 1.2247i -1.0000 - 1.2247i -1.5000 P = -0.5000 + 0.8660i -0.5000 - 0.8660i -1.5000 + 0.0000i -1.5000 - 0.0000iK = 4.0000传递函数形式：>> num=0 4 14 22 15num = 0 4 14 22 15>> den=1 4 6.25 5.25 2.25den = 1.0000 4.0000 6.2500 5.2500 2.2500部分分式：>> R，P，H=residue(num,den)R = 4.0000 -0.0000 0.0000 - 2.3094i 0.0000 + 2.3094iP = -1.5000 -1.5000 -0.5000 + 0.8660i -0.5000 - 0.8660iH = 1-2 用殴拉法matlab编程实现下列系统的输出响应在上，时的数值解。，要求保留4位小数，并将结果以图形的方式与真解比较。t=0:0.1:1t = Columns 1 through 9 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 10 through 11 0.9000 1.0000 h=0.1;y(1)=1; t=0:0.1:1; h=0.1; y(1)=1; for i=1:10y(i+1)=y(i)+h*(-1*y(i);end plot(t,y,'r') hold on m=exp(-1*t)m = Columns 1 through 9 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 Columns 10 through 11 0.4066 0.3679 plot(t,m,'bo')1-3 用四阶龙格库塔梯形法matlab编程实现1-2题的数值解，要求以图形的方式通过与真值及殴拉法的比较，分析其精度。 h=0.1; y(1)=1; t=0:0.1:1t = Columns 1 through 9 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 10 through 11 0.9000 1.0000 for i=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/6 endk1 = -1k2 = -0.9500k3 = -0.9525k4 = -0.9047y = 1.0000 0.9048k1 = -0.9048k2 = -0.8596k3 = -0.8619k4 = -0.8187y = 1.0000 0.9048 0.8187k1 = -0.8187k2 = -0.7778k3 = -0.7798k4 = -0.7407y = 1.0000 0.9048 0.8187 0.7408k1 = -0.7408k2 = -0.7038k3 = -0.7056k4 = -0.6703y = 1.0000 0.9048 0.8187 0.7408 0.6703k1 = -0.6703k2 = -0.6368k3 = -0.6385k4 = -0.6065y = 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065k1 = -0.6065k2 = -0.5762k3 = -0.5777k4 = -0.5488y = 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488k1 = -0.5488k2 = -0.5214k3 = -0.5227k4 = -0.4965y = 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966k1 = -0.4966k2 = -0.4718k3 = -0.4730k4 = -0.4493y = 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493k1 = -0.4493k2 = -0.4269k3 = -0.4280k4 = -0.4065y = Columns 1 through 9 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 Column 10 0.4066k1 = -0.4066k2 = -0.3862k3 = -0.3873k4 = -0.3678y = Columns 1 through 9 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 Columns 10 through 11 0.4066 0.3679 plot(t,y,'o') hold on m=exp(-1*t) plot(t,m,'r*') lea=y-m plot(t,lea,'g') hold offm = Columns 1 through 8 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 Columns 9 through 11 0.4493 0.4066 0.3679lea = 1.0e-006 * Columns 1 through 8 0 0.0820 0.1483 0.2013 0.2429 0.2747 0.2983 0.3149 Columns 9 through 11 0.3256 0.3315 0.33321-4采用matlab语言编程实现。程序:disp('y=')h=1; y=0;for i=0:1:63 y=y+2i; end disp(y); 运行结果：y= 1.8447e+019y= 1.8447e+0191-5编写matlab的M-函数，以实现。要求在函数中给出必要的解释和说明，同时检测输入和返回变量的个数。function sum=zuoye5(m)sum=0;format longfor i=0:m t=1; for j=1:i t=t*2; endsum=sum+t;end二、控制系统分析2-1设典型闭环结构控制系统如下图所示，当阶跃输入幅值时，用sp3_1.m求取输出y(t)的响应。进一步考虑：当反馈通道为时，如何通过对sp3_1.m程序的修改，以实现此时所对应的y(t)。y(t)r(t)_2解：程序1：clear all;close all;a=0.016 0.864 3.27 3.42 1;b=30 25;V=2;R=20;X0=0 0 0 0;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;for i=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时的阶跃响应曲线'); hold on; grid on程序2：clear all;close all;a=0.016 0.864 3.27 3.42 1;b=30 25;a1=conv(a,0.1 1);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=0 0 0 0 0;n1=5;X0=0 0 0 0;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;for i=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;for i=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;endt',y'plot(t,y,'r-');title('反馈通道为惯性环节时的阶跃响应曲线');hold on;grid on2-2 下图中，若各环节传递函数已知为：， ；试列写链接矩阵W、W0和非零元素阵WIJ，将程序sp4_2完善后，应用此程序求输出的响应曲线。G2(s)G5(s)G10(s)G4(s)y0_y7_G1(s)G3(s)G6(s)G7(s)G8(s)G9(s)解：程序clear all;close all;a=1 0 1 0 1 1 0 1 1'b=0.01 0.085 0.01 0.051 0.0067 0.15 1 0.01 0.01'c=1 1 1 1 70 0.21 130 0.1 0.0044'd=0 0.17 0 0.15 0 0 0 0 0'P=a,b,c,d;WIJ=1 0 1;2 1 1;2 9 -1;3 2 1;4 3 1; 4 8 -1;5 4 1;6 5 1;6 7 -0.212;7 6 1; 8 6 1;9 7 1;n=9;Y0=1;Yt0=0 0 0 0 0 0 0 0 0;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);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); endendQ=B-D*W;R=C*W-A;V1=C*W0;Ab=QR;b1=QV1;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+K2*2+K3*2+K4)/6; end y=y,Y(nout); t=t,t(i)+h*L1;endt',y'plot(t,y,'r');title('单位阶跃响应曲线');grid on2-3用离散相似法仿真程序sp4_3.m重求上题输出的数据与曲线，并与四阶龙格库塔法比较精度，同时分析两种方法对仿真步长选取的区别。解：程序clear all;close all;a=1 0 1 0 1 1 0 1 1'b=0.01 0.085 0.01 0.051 0.0067 0.15 1 0.01 0.01'c=1 1 1 1 70 0.21 130 0.1 0.0044'd=0 0.17 0 0.15 0 0 0 0 0'P=a,b,c,d;WIJ=1 0 1;2 1 1;2 9 -1;3 2 1;4 3 1; 4 8 -1;5 4 1;6 5 1;6 7 -0.212;7 6 1; 8 6 1;9 7 1;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);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); endend for i=1:n if(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); else end else FI(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); else end endendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,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);endt',y'plot(t,y)2-4求下图非线性系统的输出响应y(t)，并与无非线性环节情况进行比较。y(t)e(t)r(t)=105-5P=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 on Z=4 0 0 0; S=5 0 0 0; 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); 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 end end Y=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)=backlash(Ubb(i,i),Uk(i,i),Ub(i,i),s(i); end end end Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'*X+FIM'*Uk+FIJ'*Udot; Yb=Y; Y=FIC'*X+FID'*Uf; for i=1;n if(Z(i)=0) if(Z(i)=4) Y(i,i)=satu(Y(i,i),S(i); end if(Z(i)=5) Y(i,i)=dead(Y(i,i),S(i); end if(Z(i)=6) Y(i,i),Ubb(i,i)=backlash(Ubb(i,i),Y(i,i),Yb(i,i),s(i); end end end end y=y,Y(nout); end backlash函数: functionUc,Ubb=backlash(Urb,Ur,Ucb,S1) if(Ur>Urb) if(Ur-S1)>=Ucb) Uc=Ur-S1; else Uc=Ucb; end else if(Ur<Urb) if(Ur+S1)<=Ucb) Uc=Ur+S1; else Uc=Ucb; end else Uc=Ucb; end end Ubb=Ur;Satu函数: function Uc=satu(Ur,S1) if(abs(Ur)>=S1) if(Ur>0) Uc=S1; else Uc=-S1; end else Uc=Ur; end dead函数: function Uc=dead(Ur,S1) if(abs(Ur)>=S1) if(Ur>0) Uc=Ur-S1; else Uc=Ur+S1; end else Uc=0; end2-5 已知系统，当分别取-1、0、1时，判别系统的能控性与能观测性。a=-1时num=-1den=10,27,18A,B,C,D=tf2ss(num,den)P=ctrb(A,B)N=rank(P)Q=obsv(A,C)N=rank(Q)A = -2.7000 -1.8000 1.0000 0B = 1 0C = 0 -0.1000D = 0P = 1.0000 -2.7000 0 1.0000N = 2Q = 0 -0.1000 -0.1000 0N = 2a=0时A = -2.7000 -1.8000 1.0000 0B = 1 0C = 0 0D = 0P = 1.0000 -2.7000 0 1.0000N = 2Q = 0 0 0 0N = 0a=1时A = -2.7000 -1.8000 1.0000 0B = 1 0C = 0 0.1000D = 0P = 1.0000 -2.7000 0 1.0000N = 2Q = 0 0.1000 0.1000 0N = 22-6 已知系统的状态方程为，采用状态反馈，将系统的几点配置到-1、-2、-3，求状态反馈阵K。A=-2 -1 1;1 0 1;-1 0 1B=1 1 1'P=-1,-2,-3K=acker(A,B,P)K = -1 2 4
