控制系统仿真及CAD试题.doc

控制系统仿真及CAD试题（研2010）
一、 （20分）试论述系统仿真的目的、意义、分类及应用与发展概况。
解：系统仿真的目的：在分析系统各要素性质及其相互关系的基础上，建立能描述系统结构或行为过程的、且具有一定逻辑关系或数量关系的仿真模型，据此进行试验或定量分析，以获得正确决策所需的各种信息。
系统仿真的意义：CAD不是简单的使用计算机代替人工计算、制图等“传统的设计方法”，而是通过CAD系统与设计者之间强有力的“信息交互”作用，从本质上增强设计人员的想象力与创造力，从而有效地提高设计者的能力与设计结果的水平，因此，CAD技术中所涉及的“设计”应该是以提高社会生产力的水平、加快社会进步为目的的创造性的劳动。
系统仿真的分类：按模型分类分为：物理仿真和数学仿真，物理仿真又分为实物仿真、实时仿真、半实物仿真、在线仿真；数学仿真又分为数字仿真、非实时仿真、模拟仿真、离线仿真。
系统仿真的应用：现代仿真技术经过近50年的发展与完善，已经在各行业做出卓越贡献，同时也充分体现出其在科技发展与社会进步中的重要作用。仿真技术广泛应用在航空与航天工业、电力工业、原子能工业、石油、化工及冶金工业中。仿真技术还广泛应用在医学、社会学、宏观经济与商业策略的研究等非工程领域中。
系统仿真的发展概况：（1）在硬件方面，基于多CPU并行处理技术的全数字仿真系统将有效提高系统仿真的速度，从而使仿真系统“实时性”得到进一步的加强。（2）随着网络技术的不断完善与提高，分布式数字仿真系统将为人们广泛采用，从而达到“投资少、效果好”的目的。（3）在应用软件方面，直接面向用户的高效能的数字仿真软件不断推陈出新，各种专家系统与智能化技术奖更深入的应用于仿真软件开发中，使得在人—机界面、结果输出、综合评判等方面达到更理想的境界。（4）虚拟现实技术的不断完善，为控制系统数字仿真与CAD开辟了一个新时代。（5）随着FMS与CIMS技术的应用于发展，“离散事件系统”越来越多的为仿真领域所重视，离散事件仿真从理论到实现给我们带来许多新的问题。随着管理科学、柔性制造系统、计算机集成制造系统的不断发展，“离散事件系统仿真”问题越来越显示出它的重要性。

二、 （20分）用欧拉法和二阶龙格库塔法求下面系统的输出响应y(t)在0≤t≤1上，h=0.1时的数值。要求保留4位小数，并将结果与真解比较。
解：欧拉法（前向欧拉法,可以自启动）其几何意义：把f(t,y)在[]区间内的曲边面积用矩形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示

（1） m文件程序为 h=0.1;
disp('函数的数值解为'); %显示 ‘’中间的文字%
disp('y='); %同上%
y=1;
for t=0:h:1
m=y;
disp(y); %显示y的当前值%
y=m-m*h;
end
保存文件q2.m
在matalb命令行中键入>> q2
得到结果 函数的数值解为
y= 1 0.9000 0.8100 0.7290 0.6561 0.5905 0.5314 0.4783 0.4305 0.3874 0.3487

（2）另建一个m 文件求解在t[0,1]的数值 （ %是的真解%）
程序为h=0.1;
disp('函数的离散时刻解为');
disp('y=');
for t=0:h:1
y=exp(-t);
disp(y);
end 保存文件q3.m
在matlab命令行中键入>> q3
函数的离散时刻解为
y= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679

比较欧拉方法求解与真值的差别

欧拉
1
0.9000
0.8100
0.7290
0.6561
0.5905
0.5314
0.4783
0.4305
0.3874
0.3487
真值
1
0.9048
0.8187
0.7408
0.6703
0.6065
0.5488
0.4966
0.4493
0.4066
0.3679
误差
0
-0.0048
-0.007
–0.0118
–0.0142
–0.0160
–0.0174
–0.0183
–0.0188
-0.0192
-0.0192
显然误差与为同阶无穷小，欧拉法具有一阶计算精度，精度较低，但算法简单。
我们经常用到 预报-校正法 的二阶龙-格库塔法， 此方法可以自启动，具有二阶计算精度，几何意义：把f(t,y)在[]区间内的曲边面积用上下底为和、高为h的梯形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示
（1）m文件程序为 h=0.1;
disp('函数的数值解为');
disp('y=');
y=1;
for t=0:h:1
disp(y);
k1=-y;
k2=-(y+k1*h);
y=y+(k1+k2)*h/2;
end
保存文件q4.m
在matlab的命令行中键入 >> q4 显示结果为
函数的数值解为
y= 1 0.9050 0.8190 0.7412 0.6708 0.6071 0.5494 0.4972 0.4500 0.4072 0.3685

（1） 比较欧拉法与二阶龙格-库塔法求解.（误差为绝对值）

真值
1
0.9048
0.8187
0.7408
0.6703
0.6065
0.5488
0.4966
0.4493
0.4066
0.3679
龙库
1
0.9050
0.8190
0.7412
0.6708
0.6071
0.5494
0.4972
0.4500
0.4072
0.3685
误差
0
0.0002
0.0003
0.0004
0.0005
0.0006
0.0006
0.0006
0.0007
0.0006
0.0006

明显误差为的同阶无穷小，具有二阶计算精度，而欧拉法具有一阶计算精度，二阶龙格-库塔法比欧拉法计算精度高。

三、 （20分）分别使用解微分方程方法、控制工具箱、simulink求解具有如下闭环传递函数的系统的阶跃响应。
解：（1）用解微分方程方法：将转化为状态方程，利用matlab语句
>> num=[10];
>> den=[1 8 36 40 10];
>> [A B C D]=tf2ss(num,den)
得到结果：

A = -8 -36 -40 -10
1 0 0 0
0 1 0 0
0 0 1 0
B = 1
0
0
0

C = 0 0 0 10

D =0
得到状态方程
编写m文件求解微分方程组
function dx=wffc(t,x)
u=1; %阶跃响应，输入为1%
dx=[-8*x(1)-36*x(2)-40*x(3)-10*x(4)+u;x(1);x(2);x(3)];

保存文件 wffc.m %注意：保存文件的名字与函数名一致！%

在命令行键入>> [t,x]=ode45('wffc',[0,8],[0;0;0;0]);
>> y=10*x(:,4);
>> plot(t,y);
>> grid

得到结果为下图所示：

（2）控制工具箱:在matlab命令行中键入>> num=[10];
>> den=[1 8 36 40 10];
>> sys=tf(num,den);
>> step(sys);
>> grid
得到阶跃响应结果如图所示：

（3）simulink求解:在simulink模型窗口中建立如下模型，键入该题的传递函数。

start后，观察scope中的仿真波形如下：

四、 （20分）单位反馈系统的开环传递函数已知如下：

用matlab语句 、函数求取系统闭环零极点，并求取系统闭环状态方程的可控标准型实现。用单变量系统四阶龙格-库塔法求解输入阶跃函数r(t)=1(t);T=1s时的输出量y（t）的动态响应数值解。
解：已知开环传递函数，求得闭环传递函数为
在matlab命令行里键入>> a=[1 0];
>> b=[1 4.6];
>> c=[1 3.4 16.35];
>> d=conv(a,b);
>> e=conv(d,c)
e = 1.0000 8.0000 31.9900 75.2100 0
>> f=[0 0 0 5 100];
>> g=e+f
g = 1.0000 8.0000 31.9900 80.2100 100.0000
%以上是计算闭环传递函数的特征多项式%
>> p=roots(g) %计算特征多项式的根，就是闭环传递函数的极点%
p =
-0.9987 + 3.0091i
-0.9987 - 3.0091i
-3.0013 + 0.9697i
-3.0013 - 0.9697i
>> m=[5 100];
>> z=roots(m)
z = -20 %计算零点%
综上：当闭环传函形如时，可控标准型为：
；
所以可控标准型是

解：m文件为：function y=hs(A,B,C,D,R,T,h) %T为观测时间,h为计算步长，R为输入信号幅值%
disp('数值解为');
y=0;
r=R;
x=[0;0;0;0];
N=T/h;
for t=1:N;
k1=A*x+B*R;
k2=A*(x+h*k1/2)+B*R;
k3=A*(x+h*k2/2)+B*R;
k4=A*(x+h*k3)+B*R;
x=x+h*(k1+2*k2+2*k3+k4)/6;
y(t)=C*x+D*R;
end
在命令行里键入A= B= C= D= R= T= h=
y=hs(A,B,C,D,R,T,h) 得到结果。

五、 （20分）系统结构图如图所示，
（1） 写出该系统的联结矩阵和，并写出联结矩阵非零元素阵

（2）上图中，若各环节的传递函数已知为：

但；重新列写联接矩阵和非零元素矩阵，将程序sp4_2.m完善后，应用sp4_2.m求输出的响应曲线。
解： （1） 根据图中,拓扑连结关系，可写出每个环节输入受哪些环节输出的影响，
现列如入下:

把环节之间的关系和环节与参考输入的关系分别用矩阵表示出来，

即=，=，
（2）

程序为：
%输入数据%
%系统中不能出现纯比例、纯微分环节，含有微分项系数的环节不直接与外加参考输入联接
P=input('请按各环节输入参数矩阵P：');
Wij=input('请输入非零元素连接阵Wij:');
n=input('请输入环节个数（系统阶次）n=');
Y0=input('请输入阶跃输入幅值Y0=');
Yt0=input('请输入各环节初值Yt0=');
h=input('请输入计算步长h=');
L1=input('请输入打印间隔点数（每隔点l1输出一次）：');
T0=input('请输入起始时间T0=');
Tf=input('请输入终止时间Tf=');
nout=input('请输入输出环节编号：');
%形成闭环各系数阵%
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);%建立初始W（n*n型方阵）、W0阵（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;Qn=inv(Q);%求Q和Q的逆
R=C*W-A;V1=C*W0;%求R和V1阵
Ab=Qn*R;b1=Qn*V1;%形成闭环系数阵
%数值积分求解%
Y=Yt0';y=Y(nout);t=T0;%置初值，做好求解准备
N=round((Tf-T0)/(h*L1));
for i=1:N;%每循环一次，输出一点数据
for j=1:L1;%每输出点之间计算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+2*K2+2*K3+K4)/6;
end
y=[y,Y(nout)];
t=[t,t(i)+h*L1];
end
shuchu=[t',y'];
plot(t,y)
执行后在命令窗口输入：
请按各环节输入参数矩阵P：[1 0.01 1 0;0 0.085 1 0.17;1 0.01 1 0;0 0.051 1 0.15;1 0.0067 70 0;1 0.15 0.21 0;0 1 130 0;1 0.01 0.1 0;1 0.01 0.0044 0]
请输入非零元素连接阵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.1
请输入打印间隔点数（每隔点l1输出一次）：3
请输入起始时间T0=0
请输入终止时间Tf=10
请输入输出环节编号：7
执行结果为：