控制系统数字仿真与CAD习题

1、控制系统计算机仿真及辅助设计课程上机作业1月5日上机作业一：2-1 思考题：（1）数学模型的微分方程，状态方程，传递函数，零极点增益和部分分式五种形式，各有什么特点？（2）数学模型各种形式之间为什么要互相转换？（3）控制系统建模的基本方法有哪些？他们的区别和特点是什么？（4）控制系统计算机仿真中的“实现问题”是什么含意？（5）数值积分法的选用应遵循哪几条原则？（1）微分方程是直接描述系统输入和输出量之间的制约关系，是连续控制系统其他数学模型表达式的基础。状态方程能够反映系统内部各状态之间的相互关系，适用于多输入多输出系统。传递函数是零极点形式和部分分式形式的基础。零极点增益形式可用于分析系统的

2、稳定性和快速性。利用部分分式形式可直接分析系统的动态过程。（2）不同的控制系统的分析和设计方法，只适用于特定的数学模型形式。（3）控制系统的建模方法大体有三种：机理模型法，统计模型法和混合模型法。机理模型法就是对已知结构，参数的物理系统运用相应的物理定律或定理，经过合理的分析简化建立起来的各物理量间的关系。该方法需要对系统的内部结构和特性完全的了解，精度高。统计模型法是采用归纳的方法，根据系统实测的数据，运用统计规律和系统辨识等理论建立的系统模型。该方法建立的数学模型受数据量不充分，数据精度不一致，数据处理方法的不完善，很难在精度上达到更高的要求。混合法是上述两种方法的结合。（4）“实现问题”

3、就是根据建立的数学模型和精度，采用某种数值计算方法，将模型方程转换为适合在计算机上运行的公式和方程，通过计算来使之正确的反映系统各变量动态性能，得到可靠的仿真结果。（5）数值积分法应该遵循的原则是在满足系统精度的前提下，提高数值运算的速度和并保证计算结果的稳定。2-2.用matlab语言求下列系统的状态方程、传递函数、零极点增益、和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：(1) G（s）= (2) =y=0 2 0 2 X解：（1）程序代码如下：clear;clc;num=1 7 24 24;den=1 10 35 50 24;Z,P,K=tf2zp(num,den);A,B

4、,C,D=zp2ss(Z,P,K);R,P1,H=residue(num,den);G1=tf(num,den)G2=zpk(Z,P,K)G3=ss(A,B,C,D)程序的输出结果如下： Transfer function: s3 + 7 s2 + 24 s + 24-s4 + 10 s3 + 35 s2 + 50 s + 24 Zero/pole/gain:(s+1.539) (s2 + 5.461s + 15.6)-(s+4) (s+3) (s+2) (s+1)a = x1 x2 x3 x4 x1 -3 -1.414 0 0 x2 1.414 0 0 0 x3 1 1.088 -7 -3.

5、464 x4 0 0 3.464 0 b = u1 x1 1 x2 0 x3 0 x4 0 c = x1 x2 x3 x4 y1 1 1.088 -1.539 1.038 d = u1 y1 0 Continuous-time model.（2）程序代码如下：clear;clc;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;B=4 2 2 0;C=0 2 0 2;D=0;num,den=ss2tf(A,B,C,D);Z,P,K=ss2zp(A,B,C,D);R,P1,H

6、=residue(num,den);G1=tf(num,den)G2=zpk(Z,P,K)G3=ss(A,B,C,D)程序输出结果如下：Transfer function: 4 s3 + 14 s2 + 22 s + 15-s4 + 4 s3 + 6.25 s2 + 5.25 s + 2.25 Zero/pole/gain:4 (s+1.5) (s2 + 2s + 2.5)- (s+1.5)2 (s2 + s + 1) a = x1 x2 x3 x4 x1 2.25 -5 -1.25 -0.5 x2 2.25 -4.25 -1.25 -0.25 x3 0.25 -0.5 -1.25 -1 x4

7、 1.25 -1.75 -0.25 -0.75 b = u1 x1 4 x2 2 x3 2 x4 0 c = x1 x2 x3 x4 y1 0 2 0 2 d = u1 y1 0 Continuous-time model.2-3.用欧拉法求下面系统的输出响应y(t)在0t1上，h=0.1时的数值。 要求保留4位小数，并将结果与真解比较。解：程序代码如下：clear;clc;h=0.1;y=1;disp('欧拉法求得值为：');for i=0:h:1 m=y; disp(y); y=m-m*h;enddisp('真解为：');for i=0:h:1 y=exp(

8、-i); disp(y);end输出结果如下：欧拉法值10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.36792-4用二阶龙格库塔法求解2-3的数值解，并于欧拉法求得的结果比较。解：程序代码如下：clear;clc;h=0.1;y=1;disp('二阶龙格库塔法求得值为：');for i=0:h:1 disp(y); k1=-y; k2=-(y+k1*h); y=y+(k1+k2)*h/2

9、;enddisp('欧拉法求得值为：'); y=1;for i=0:h:1 m=y; disp(y); y=m-m*h;end程序输出结果如下：欧拉法值10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487龙格库塔法值10.90500.81900.74120.67080.60710.54940.49720.45000.40720.36852-5用四阶龙格-库塔法求解题2-3数值解，并与前两题结果相比较。解：程序代码如下：clear;clc;h=0.1;y=1;disp('四阶龙格库塔法求得值为：'

10、;);for i=0:h:1 disp(y); k1=-y; k2=-(y+k1*h/2); k3=-(y+k2*h/2); k4=-(y+k3*h); y=y+(k1+2*k2+2*k3+k4)*h/6;end程序结果如下：欧拉法值10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487二阶龙格库塔法值10.90500.81900.74120.67080.60710.54940.49720.45000.40720.3685四阶龙格库塔法值10.90480.81870.74080.67030.60650.54880.49660.

11、44930.40660.36792-6已知二阶系统状态方程为写出取计算步长为h时，该系统状态变量X=的四阶龙格-库塔法递推关系式。2-7单位反馈系统的开环传递函数已知如下用matlab语句 、函数求取系统闭环零极点，并求取系统闭环状态方程的可控标准型实现。解：程序代码如下：clear;clc;num=5 100;den1=1 0;den2=1 4.6;den3=1 3.4 16.35;den4=0 0 0 5 100;den=conv(den1,den2);den=conv(den,den3);den=den+den4Z,P,K=tf2zp(num,den);P1=roots(den)Z1=r

12、oots(num)ZP输出结果如下：den = 1.0000 8.0000 31.9900 80.2100 100.0000P1 = -0.9987 + 3.0091i -0.9987 - 3.0091i -3.0013 + 0.9697i -3.0013 - 0.9697iZ1 = -20Z = -20P = -0.9987 + 3.0091i -0.9987 - 3.0091i -3.0013 + 0.9697i -3.0013 - 0.9697i2-8用matlab语言编制单变量系统三阶龙格-库塔法求解程序，程序入口要求能接收状态方程各系数阵（A,B,C,D）,和输入阶跃函数r(t)=R

13、*1(t);程序出口应给出输出量y（t）的动态响应数值解序列。解：函数程序代码如下：function y = hs( A,B,C,D,R,T,h ) y=0; r=R; x=0;0;0;0; N=T/h; for t=1:1:N k1=A*x+B*R; k2=A*(x+h*k1/3)+B*R; k3=A*(x+2*h*k2/3)+B*R; x=x+h*(k1+3*k3)/4; y=C*x+D*R;end 将程序保存为名为hs的文件，用以求解2-1求解出的A,B,C,D进行验证如下：>> hs(A,B,C,D,1,1,1)ans = 2.66672-10.用式（2-34）梯形法求解试

14、验方程，分析对计算步长h有何限制，说明h对数值稳定性的影响。1月6日上机作业二：3-1 求解下列线性方程，并进行解得验证：(1),(2) 由A*X=B得：X=AB解：（1）程序代码如下：clear;clc;A=7 2 1 -2 9 15 3 -2 -2 -2 11 5 1 3 2 13;B=4 7 -1 0;x=inv(A)*B输出结果如下：x = 0.4979 0.1445 0.0629 -0.0813（2）程序代码如下：clear;clc;A= 5 7 6 5 1 7 10 8 7 2 6 8 10 9 3 5 7 9 10 4 1 2 3 4 5;B=24 96 34 136 36 14

15、4 35 140 15 60;x=inv(A)*B输出结果如下：x = 1.0000 4.0000 1.0000 4.0000 1.0000 4.0000 1.0000 4.00001.0000 4.00003-2进行下列计算，给出不使用for和while等循环语句的计算方法。（1） (2)求出y=x*sin(x) 在0<x<100条件下的每个峰值解：（1）程序代码如下：clear;clc;a1=1;q=2;n=64;sum=a1*(1-qn)/(1-q)输出结果如下：sum = 1.8447e+019（2）先画出0<x<100条件下y的图形代码如下：x=0:0.1:1

16、00;y=x.*sin(x);plot(x,y)xlabel('x');ylabel('y');图形如下：可以看出sin(x)在峰值时y也是在峰值位置，即每个点的峰值为3-3绘制下面的图形。（1）sin(1/t),-1<t<1 (2) -1<t<1解：（1）代码如下：t=-1:0.01:1;y=sin(1./t);plot(t,y);xlabel('t');ylabel('y');输出图形如下：（2）代码如下：clear;clc;t=-1:0.01:1;y=1-(cos(7*t).3;plot(t,y);x

17、label('t');ylabel('y');输出图形如下：3-4已知元件的实验数据如下，拟合这一数据，并尝试给出其特性方程。 X 0.0100 1.0100 2.0100 3.0100 4.0100 Y 2.5437 7.8884 9.6242 11.6071 11.9727 X 5.0100 6.0100 7.0100 8.0100 9.0100 y 13.2189 14.2679 14.6134 15.4045 15.0805 解：程序代码如下：clear;clc;x= 0.0100 1.0100 2.0100 3.0100 4.0100 5.0100 6

18、.0100 7.0100 8.0100 9.0100;y= 2.5437 7.8884 9.6242 11.6071 11.9727 13.2189 14.2679 14.6134 15.4045 15.0805 ;P=polyfit(x,y,3);xi=0:0.01:9.01;yi=polyval(P,xi);plot(x,y,'o',xi,yi,'*');plot(x,y,xi,yi);xlabel('x');ylabel('y');输出结果如下：3-5分别使用解微分方程方法、控制工具箱、simulink求解具有如下闭环传递函

19、数的系统的阶跃响应。解：Simulink中建立如下系统阶跃响应图像如下：输入以下代码：num=6 26 6 20;den=1 3 4 2 2;step(num,den,10);输出波形如下3-6已知系统的闭环传递函数，试分析该系统的稳定性。解：画出系统零极点图代码如下：num=6 26 6 20;den=1 8 36 40 10;pzmap(num,den)输出结果为：可以看出系统的几点都在左半平面中，所以该系统为稳定系统。3-8某小功率随动系统动态结构如下图所示，已知： 若系统输入分别为，适用simulink分析系统的输出分别如何？（1）输入为1（t）:输出为：（2）输入为 t 时： 输出：

20、（3） 输入为l(t)-l(t-1.5)时：1月7日上机作业三：4-2设典型闭环结构控制系统如图4-47所示，当阶跃输入幅值 时，用sp4_1.m求取输出的响应。解：代码为：a=0.016 0.864 3.27 3.42 1;b=30 25;X0=0 0 0 0; %系统状态向量初值为零V=2; %反馈系数 n=4;T0=0;Tf=10;h=0.01;R=20 ; sp4_1;plot(t,y)输出结果为：4-8求图4-49非线性系统的输出响应y(t),并与无非线性环节情况进行比较。解：程序代码为：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=0Tf=20;nout=4;Y0=10;sp4_4;plot(t,y,'r') hold on