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

资源描述

兰州理工大学《控制系统计算机仿真》上机汇报Ⅰ 院系：电气工程和信息工程学院班级： 14级自动化3班姓名：孙悦学号：时间：年 6 月 15 日电气工程和信息工程学院《控制系统计算机仿真》上机试验任务书Ⅰ（）一、上机试验内容及要求 1. matlab软件要求利用课余时间熟悉掌握matlab软件基础数值运算、基础符号运算、基础程序设计方法及常见图形命令操作。 2.各章节仿真试验内容及要求具体试验内容及要求请详见上机试验汇报。二、上机试验时间安排及相关事宜 1. 依据课程教学纲领要求，上机试验课时共16课时，学生须在每次上机之前做好对应准备工作，以确保在有限机时内完成仿真试验要求内容； 2. 试验完成后按要求完成相关仿真试验汇报； 3. 仿真试验汇报请按相关样本制作并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.4000 1.0000 0 0 0 1.0000 0 B = 1 0 0 C = 0.7000 2.4000 2.4000 D = 0 零极点增益形式： [Z,P,K]=tf2zp(num,den) Z = -1.7143 + 0.6999i -1.7143 - 0.6999i P = -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.75] 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.7500 >> B=[4;2;2;0] B = 4 2 2 0 >> C=[0 2 0 2] C = 0 2 0 2 >> D=0 D = 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.0000i K = 4.0000 传输函数形式： >> num=[0 4 14 22 15] num = 0 4 14 22 15 >> den=[1 4 6.25 5.25 2.25] den = 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.3094i P = -1.5000 -1.5000 -0.5000 + 0.8660i -0.5000 - 0.8660i H = [] 1-2 用殴拉法matlab编程实现下列系统输出响应在上，时数值解。，要求保留4位小数，并将结果以图形方法和真解比较。 t=0:0.1:1 t = Columns 1 through 9 0 0.1000 0. 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:10 y(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:1 t = Columns 1 through 9 0 0.1000 0. 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 10 through 11 0.9000 1.0000 for i=1:10 k1=-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 end k1 = -1 k2 = -0.9500 k3 = -0.9525 k4 = -0.9047 y = 1.0000 0.9048 k1 = -0.9048 k2 = -0.8596 k3 = -0.8619 k4 = -0.8187 y = 1.0000 0.9048 0.8187 k1 = -0.8187 k2 = -0.7778 k3 = -0.7798 k4 = -0.7407 y = 1.0000 0.9048 0.8187 0.7408 k1 = -0.7408 k2 = -0.7038 k3 = -0.7056 k4 = -0.6703 y = 1.0000 0.9048 0.8187 0.7408 0.6703 k1 = -0.6703 k2 = -0.6368 k3 = -0.6385 k4 = -0.6065 y = 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 k1 = -0.6065 k2 = -0.5762 k3 = -0.5777 k4 = -0.5488 y = 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 k1 = -0.5488 k2 = -0.5214 k3 = -0.5227 k4 = -0.4965 y = 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 k1 = -0.4966 k2 = -0.4718 k3 = -0.4730 k4 = -0.4493 y = 1.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 k1 = -0.4493 k2 = -0.4269 k3 = -0.4280 k4 = -0.4065 y = 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.4066 k1 = -0.4066 k2 = -0.3862 k3 = -0.3873 k4 = -0.3678 y = 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 off m = 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.3679 lea = 1.0e-006 * Columns 1 through 8 0 0.0820 0.1483 0. 0.2429 0.2747 0.2983 0.3149 Columns 9 through 11 0.3256 0.3315 0.3332 1-4采取matlab语言编程实现。程序: disp('y=') h=1; y=0; for i=0:1:63 y=y+2^i; end disp(y); 运行结果： y= 1.8447e+019 y= 1.8447e+019 1-5编写matlabM-函数，以实现。要求在函数中给出必需解释和说明，同时检测输入和返回变量个数。 function sum=zuoye5(m) sum=0; format long for i=0:m t=1; for j=1:i t=t*2; end sum=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: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+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]; end X=X0';y=0;t=T0; for i=1:N K1=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('反馈通道为惯性步骤时阶跃响应曲线'); hold on;grid on 2-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); end end Q=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)); 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]; end [t',y'] plot(t,y,'r'); title('单位阶跃响应曲线'); grid on 2-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); end end 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 end end Y=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)]; end [t',y'] plot(t,y) 2-4求下图非线性系统输出响应y(t)，并和无非线性步骤情况进行比较。 y(t) e(t) r(t)=10 5 -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 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函数: function[Uc,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; end 2-5 已知系统，当分别取-1、0、1时，判别系统能控性和能观察性。 a=-1时 num=[-1] den=[10,27,18] [A,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 0 B = 1 0 C = 0 -0.1000 D = 0 P = 1.0000 -2.7000 0 1.0000 N = 2 Q = 0 -0.1000 -0.1000 0 N = 2 a=0时 A = -2.7000 -1.8000 1.0000 0 B = 1 0 C = 0 0 D = 0 P = 1.0000 -2.7000 0 1.0000 N = 2 Q = 0 0 0 0 N = 0 a=1时 A = -2.7000 -1.8000 1.0000 0 B = 1 0 C = 0 0.1000 D = 0 P = 1.0000 -2.7000 0 1.0000 N = 2 Q = 0 0.1000 0.1000 0 N = 2 2-6 已知系统状态方程为，采取状态反馈，将系统几点配置到-1、-2、-3，求状态反馈阵K。 A=[-2 -1 1;1 0 1;-1 0 1] B=[1 1 1]' P=[-1,-2,-3] K=acker(A,B,P) K = -1 2 4

展开阅读全文