控制系统计算机仿真作业

兰州理工大学《控制系统计算机仿真》上机报告Ⅰ院系：电气工程与信息工程学院班级：14级自动化3班姓名：**学号： **********时间： 2017 年 6 月 15 日电气工程与信息工程学院《控制系统计算机仿真》上机实验任务书Ⅰ（2017）一、上机实验内容及要求1.matlab软件要求利用课余时间熟悉掌握matlab软件的基本数值运算、基本符号运算、基本程序设计方法及常用的图形命令操作。

2.各章节仿真实验内容及要求具体实验内容及要求请详见上机实验报告。

二、上机实验时间安排及相关事宜1.依据课程教学大纲要求，上机实验学时共16学时，学生须在每次上机之前做好相应的准备工作，以确保在有限的机时内完成仿真实验要求的内容；2.实验完成后按规定完成相关的仿真实验报告；3.仿真实验报告请按有关样本制作并A4打印，侧面装订，作为成绩评定的一部分。

自动化系《控制系统计算机仿真》课程组2017年3月一、Matlab 基础操作1-1用MATLAB 语言求下列系统的状态方程、传递函数、零极点增益和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：（1）2450351024247)(23423+++++++=s s s s s s s s G程序如下:num=[7,24,24]den=[10,35,50,24][A,B,C,D]=tf2ss(num,den)系统的状态方程：A =-3.5000 -5.0000 -2.40001.0000 0 00 1.0000 0B =1C =0.7000 2.4000 2.4000D =零极点增益形式：[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）u X X ⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------------=⋅022475.025.075.125.1125.15.025.025.025.125.425.25.025.1525.2 []X y 2020=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 =422>> C=[0 2 0 2]C =0 2 0 2>> D=0D =零极点增益形式：>> [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=[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.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 编程实现下列系统的输出响应)(t y 在10≤≤t 上，1.0=h 时的数值解。

y y -='，1)0(=y要求保留4位小数，并将结果以图形的方式与真解t e t y -=)(比较。

t=0:0.1:1t =Columns 1 through 90 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000Columns 10 through 110.9000 1.0000h=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));endplot(t,y,'r')hold onm=exp(-1*t)m =Columns 1 through 91.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493Columns 10 through 110.4066 0.3679plot(t,m,'bo')00.10.20.30.40.50.60.70.80.911-3 用四阶龙格－库塔梯形法matlab编程实现1-2题的数值解，要求以图形的方式通过与真值及殴拉法的比较，分析其精度。

h=0.1;y(1)=1;t=0:0.1:1t =Columns 1 through 90 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000Columns 10 through 110.9000 1.0000for 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/6endk1 =-1k2 =-0.9500k3 =-0.9525k4 =-0.9047y =1.0000 0.9048 k1 =-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.7408 k1 =-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.6065 k1 =-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 91.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493Column 100.4066k1 =-0.4066k2 =-0.3862k3 =-0.3873k4 =-0.3678y =Columns 1 through 91.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493Columns 10 through 110.4066 0.3679plot(t,y,'o')hold onm=exp(-1*t)plot(t,m,'r*')lea=y-mplot(t,lea,'g')hold offm =Columns 1 through 81.0000 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966Columns 9 through 110.4493 0.4066 0.3679lea =1.0e-006 *Columns 1 through 80 0.0820 0.1483 0.2013 0.2429 0.2747 0.2983 0.3149Columns 9 through 110.3256 0.3315 0.33321-4采用matlab 语言编程实现∑==6302i i k 。

程序:disp('y=') h=1; y=0;for i=0:1:63 y=y+2^i; enddisp(y); 运行结果： y=1.8447e+019 y=1.8447e+0191-5编写matlab 的M-函数，以实现∑==mi i k 02。

要求在函数中给出必要的解释和说明，同时检测输入和返回变量的个数。

function sum=zuoye5(m) sum=0;0.10.20.30.40.50.60.70.80.91format long for i=0:m t=1;for j=1:i t=t*2; end sum=sum+t; end二、控制系统分析2-1设典型闭环结构控制系统如下图所示，当阶跃输入幅值20=R 时，用sp3_1.m 求取输出y(t)的响应。

进一步考虑：当反馈通道为ss H 1.011)(+=时，如何通过对sp3_1.m 程序的修改，以实现此时所对应的y(t)。

解：程序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];y(t)r(t)_2[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];[t',y']plot(t,y,'r-');title('反馈通道为惯性环节时的阶跃响应曲线'); hold on;grid on2-2 下图中，若各环节传递函数已知为：s s G 01.011)(1+=，sss G 085.017.01)(2+=，s s G 01.011)(3+=，s s s G 051.015.01)(4+=，s s G 0067.0170)(5+=，s s G 15.0121.0)(6+=，s s G 130)(7=，s s G 01.011.0)(8+=，ss G 01.010044.0)(9+=， 212.0)(10=s G ；试列写链接矩阵W 、W 0和非零元素阵W IJ ，将程序sp4_2完善后，应用此程序求输出7y 的响应曲线。