(整理)控制系统数字仿真第二章习题答案

控制系统数字仿真与CAD第二章习题答案
2-1 思考题：
（1）数学模型的微分方程，状态方程，传递函数，零极点增益和部分分式五种形式，各有什么特点？
（2）数学模型各种形式之间为什么要互相转换？
（3）控制系统建模的基本方法有哪些？他们的区别和特点是什么？
（4）控制系统计算机仿真中的“实现问题”是什么含意？
（5）数值积分法的选用应遵循哪几条原则？
答：（1）微分方程是直接描述系统输入和输出量之间的制约关系，是连续控制系统其他数学模型表达式的基础。

状态方程能够反映系统内部各状态之间的相互关系，适用于多输入多输出系统。

传递函数是零极点形式和部分分式形式的基础。

零极点增益形式可用于分析系统的稳定性和快速性。

利用部分分式形式可直接分析系统的动态过程。

（2）不同的控制系统的分析和设计方法，只适用于特定的数学模型形式。

（3）控制系统的建模方法大体有三种：机理模型法，统计模型法和混合模型法。

机理模型法就是对已知结构，参数的物理系统运用相应的物理定律或定理，经过合理的分析简化建立起来的各物理量间的关系。

该方法需要对系统的内部结构和特性完全的了解，精度高。

统计模型法是采用归纳的方法，根据系统实测的数据，运用统计规律和系统辨识等理论建立的系统模型。

该方法建立的数学模型受数据量不充分，数据精度不一致，数据处理方法的不完善，很难在精度上达到更高的要求。

混合法是上述两种方法的结合。

（4）“实现问题”就是根据建立的数学模型和精度，采用某种数值计算方法，将模型方程转换为适合在计算机上运行的公式和方程，通过计算来使之正确的反映系统各变量动态性能，得到可靠的仿真结果。

（5）数值积分法应该遵循的原则是在满足系统精度的前提下，提高数值运算的速
度和并保证计算结果的稳定。

2-2.用matlab语言求下列系统的状态方程、传递函数、零极点增益、和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：
(1) G（s）=
32
432
72424
10355024
s s s
s s s s
+++
++++
(2)
.
X=
2.25 -5 -1.25 -0.54
2.25 -4.25 -1.25 -0.252
0.25 -0.5 -1.25 -12
1.25 -1.75 -0.25 -0.75 0
X
⎡⎤⎡⎤
⎢⎥⎢⎥
⎢⎥⎢⎥
+⎢⎥⎢⎥
⎢⎥⎢⎥
⎣⎦⎣⎦
u y=[0 2 0 2] X
（1）解：（1）状态方程模型参数:编写matlab程序如下>> num=[1 7 24 24];
>> den=[1 10 35 50 24];
>> [A B C D]=tf2ss(num,den)
得到结果：A=
-10 -35 -50 -24
1 0 0 0
0 1 0 0
0 0 1 0
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎣⎦
,B=
1
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎣⎦
,C=[]
1 7 24 24,D=[0]
所以模型为:
.
X=
-10 -35 -50 -24
1 0 0 0
0 1 0 0
0 0 1 0
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎣⎦
X+
1
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎣⎦
u,y=[]
1 7 24 24X
(2)零极点增益:编写程序>> num=[1 7 24 24];
>> den=[1 10 35 50 24];
>> [Z P K]=tf2zp(num,den)
得到结果Z= -2.7306 + 2.8531 , -2.7306 - 2.8531i ,-1.5388
P= -4, -3 ,-2 ,-1
K=1
(3) 部分分式形式:编写程序>> num=[1 7 24 24];
>> den=[1 10 35 50 24];
>> [R P H]=residue(num,den)
得到结果R= 4.0000 ,-6.0000, 2.0000, 1.0000 P= -4.0000, -3.0000 , -2.0000 ,-1.0000 H=[]
G(s)=4621
4321
s s s s -+++
++++
（2）解：（1）传递函数模型参数：编写程序>> 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)
得到结果
num = 0 4.0000 14.0000 22.0000 15.0000 den =1.0000 4.0000 6.2500 5.2500 2.2500
32432
4 s + 14 s + 22 s + 15
()s + 4 s + 6.25 s + 5.25 s + 2.25
G s =
(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]；
>> B=[4 2 2 0]'; >> C=[0 2 0 2];
>> 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 -1.5000 K = 4.0000
表达式 ()()
()()()
4s+1-1.2247i s+1+1.2247i ()s+0.5-0.866i s+0.5+0.866i s+1.5G s =
（3）部分分式形式的模型参数：编写程序>> 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)
>> [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 =[]
4 2.3094 2.3094() 1.50.50.8660.50.866i i
G s s s i s i
=
-+
++-++
2-3.用欧拉法求下面系统的输出响应y(t)在0≤t ≤1上，h=0.1时的数值。

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

解：欧拉法1'0
0*(,)(,)()k k k k k k y y h f t y y f t y y t y
+=+⎧⎪
=⎨⎪=⎩（前向欧拉法,可以自启动）其几何意义：把f(t,y)在[,k k 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
y e-
=在t∈[0,1]的数值（%t
y e-
=是',(0)1
y y y
=-=的真解%）
程序为h=0.1;
disp('函数的离散时刻解为');
disp('y=');
for t=0:h:1
y=exp(-t);
disp(y);
end 保存文件q3.m
在matalb命令行中键入>> q3
函数的离散时刻解为
y= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679 比较欧拉方法求解与真值的差别
显然误差与2h为同阶无穷小，欧拉法具有一阶计算精度，精度较低，但算法简单。

2-4用二阶龙格库塔法求解2-3的数值解，并于欧拉法求得的结果比较。

解：我们经常用到预报-校正法的二阶龙-格库塔法，
112
1
21
'
()
2
(,)
(,)
(,)
k k
k k
k k
h
y y k k
k f t y
k f t h y hk
f t y y
+
⎧
=++
⎪
⎪⎪
=
⎨
⎪=++
⎪
=
⎪⎩
此方法
可以自启动，具有二阶计算精度，几何意义：把f(t,y)在[,k k t y ]区间内的曲边面积用上下底为k f 和1k f 、高为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
（2） 比较欧拉法与二阶龙格-库塔法求解.（误差为绝对值）
明显误差为3h 得同阶无穷小，具有二阶计算精度，而欧拉法具有以阶计算精度，二阶龙格-库塔法比欧拉法计算精度高。

2-5．用四阶龙格-库塔法求解题2-3数值解，并与前两题结果相比较。

解：四阶龙格-库塔法表达式
11234
1
21
32
43
(22)
6
(,)
(,)
22
(,)
22
(,)
k k
k k
k k
k k
k k
h
y y k k k k
k f t y
h h
k f t y k
h h
k f t y k
k f t h y hk
+
⎧
=++++
⎪
⎪
=
⎪
⎪⎪
=++
⎨
⎪
⎪
=++
⎪
⎪
=++
⎪⎩
，其截断误差为5h同阶无
穷小，当h步距取得较小时，误差是很小的.
(1)编辑m文件程序h=0.1;
disp('四阶龙格-库塔方法求解函数数值解为');
disp('y=');
y=1;
for t=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 保存文件q5.m
在matlab命令行里键入>> q5
得到结果四阶龙格-库塔方法求解函数数值解为
y= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679 (2)比较这几种方法：
对于四阶龙格-库塔方法
显然四阶龙格-库塔法求解精度很高，基本接近真值。

三种方法比较可以得到
精度（四阶）〉精度（二阶）〉精度（欧拉）
2-6．已知二阶系统状态方程为.1111210111.2222021222(0);;(0)x a a x x b x u x b x x
a a x ⎡⎤
⎡⎤⎡⎤
⎡⎤⎡⎤⎡⎤⎢⎥=+=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎢⎥⎣⎦
写出取计算步长为h 时，该系统状态变量X=[12,x x ]的四阶龙格-库塔法递推关系式。

解：四阶龙格-库塔法表达式1
12341213
243(22)6
(,)(,)22(,)22
(,)k k k k k k k k k k h y y k k k k k f t y h h k f t y k h h k f t y k k f t h y hk +⎧
=++++⎪⎪
=⎪⎪⎪
=++⎨⎪
⎪
=++⎪⎪=++⎪⎩
所以状态变量的递推公式可以写作：
A=11122122a a a a ⎡⎤⎢⎥⎣⎦
,B=12b b ⎡⎤⎢⎥⎣⎦,12x x x ⎡⎤
=⎢⎥⎣⎦可以写成.X AX Bu =+
则递推形式1
12341213
2
43*(22)6
(*/2)(*/2)(*)k k k k k k h X X k k k k k AX Bu k A X k h Bu k A X k h Bu
k A X k h Bu +⎧
=++++⎪⎪
=+⎪⎪
=++⎨⎪=++⎪⎪=++⎪⎩
2-7单位反馈系统的开环传递函数已知如下 25100
()( 4.6)( 3.416.35)
s G s s s s s +=
+++
用matlab 语句 、函数求取系统闭环零极点，并求取系统闭环状态方程的可控标准型实现。

解：已知开环传递函数，求得闭环传递函数为2
5100
()( 4.6)( 3.416.35)5100
s G s s s s s s +=
+++++
在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 %计算零点%
综上：当闭环传函形如1111
11...()...n n n
n n n n b s b s b G s s a s a s a ----+++=++++时，可控标准型为： 1010...00001...00;00101n A B a a ⎧⎫⎡⎤⎪⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥==⎨⎬⎢⎥⎪⎪⎢⎥⎪⎪
⎢⎥--⎪⎪⎩⎭⎣
⎦；[]1
1;0n n C b b b D -⎡⎤==⎣⎦
所以可控标准型是.11.22.33.44
123401000
001000
001010080.2131.9981
[100500][0]
x x x x u x x x x x x Y u x x ⎡⎤⎢⎥⎡⎤⎡
⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥----⎢⎥⎣⎦⎣⎦
⎣⎦⎢⎥⎢⎥⎣⎦
⎡⎤
⎢⎥⎢⎥=-+⎢⎥⎢⎥⎢⎥⎣⎦
2-8用matlab 语言编制单变量系统三阶龙格-库塔法求解程序，程序入口要求能接收状态方程各系数阵（A,B,C,D ）,和输入阶跃函数r(t)=R*1(t);程序出口应给出输出量y （t ）的动态响应数值解序列
01,,......,n y y y 。

解：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/3)+B*R; k3=A*(x+2*h*k2/3)+B*R;
x=x+h*(k1+3*k3)/4; y(t)=C*x+D*R; end
在命令行里键入A= B= C= D= R= T= h= y=hs(A,B,C,D,R,T,h) 得到结果。

2-9．用题2-8仿真程序求解题2-7系统的闭环输出响应y(t).
解：A=01000010000110080.2131.998⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥----⎣⎦,B=0001⎡⎤
⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦
,C=[100500]-,D=[0] 在命令行里键入>> A=[0 1 0 0
0 0 1 0
0 0 0 1
-100 -80.21 -31.99 -8]; >> B=[0 0 0 1]'; >> C=[-100 5 0 0]; >> D=[0]; >> T=1; >> R=1; >> h=0.01;
>> y=hs(A,B,C,D,R,T,h)
数值解为
8.3333e-007
5.8659e-006
1.8115e-005
3.9384e-005
7.0346e-005。

%仅取一部分%
2-10.用式（2-34）梯形法求解试验方程'1
y y
τ
=-，分析对计算步长h有何限制，说明h对数值稳定性的影响。

解：编写梯形法程序为
112 1
2
()
2
1
11
() k k
k
k k
h
y y k k k y
k y y h
τ
ττ
+
⎧
=++
⎪
⎪
⎪
=-
⎨
⎪
⎪
=--
⎪⎩
得到
2
12
(1)
2
k k
h h
y y
ττ
+
=-+稳定系统最终渐进收敛。

系统稳定则
2
2
11
2
h h
ττ
-+<计算得02
hτ
<<。

h的选取不能超出上述范围，否则系统不稳定。

2-11如图2-27所示斜梁滚球系统，若要研究滚球在梁上的位置可控性，需首先建立其数学模型，已知力矩电机的输出转矩M与其电流i成正比，横梁为均匀可自平衡梁（即当电机不通电且无滚球时，横梁可处于θ=0的水平状态），是建立系统的数学模型，并给出简化后系统的动态结构图。

解：设球的质心到杆的距离为0，该系统为特殊情况下的球棒系统。

另令12,,I m I 分别表示棒的惯量、球的质量和球的惯量。

则球质心的位置和速度为
(cos ,sin )
(cos sin ,sin cos )
c c x x x v v x v x θθθωθθωθ==-+
其中x v =，θω=。

因而动能的移动部分为
因而动能的移动部分为 222
211
()22trans c K mv m v x ω==+
球棒系统的旋转动能为 221211()22rot
v K I I r
ω=+ 因而，系统总的动能trans rot K K K =+等于
222111
()22
K I mx mv ωλ=
++ 其中2
2
11I mr
λ=+
>为常数。

此系统的拉格朗日方程组为
()sin ()cos d T T
mg dt x x
d T T ki mg dt θθθθ∙∙∂∂⎧-=-⎪∂⎪∂⎨
∂∂⎪-=-⎪∂∂⎩
综合以上公式的系统的方程组为
22
1
sin()0()2cos()m x mx mg I mx mxx mgx ki λθθθθθ⎧-+=⎪
⎨+++=⎪⎩ 设系统在平衡点附近0θ∙
≈，cos 1θ≈，sin θθ≈，则系统方程可化为
2
1
()m x mg I mx mgx ki λθθ+=⎧⎨++=⎩ 对上式进行拉普拉斯变换并化简后可得到()
()
X s I s 。

参考文献：
[1] Hauser, S. Sestry, and P. Kokotovic. “Nonlinear control via approximate input -output linearization”. IEEE Trans. on Automatic Control, vol.37:pp.392-398, 1992.
[2] R. Sepulchre. “S low peaking and low-gain designs for global stabilization of nonlinear systems”. submitted for IEEE TAC 1999.
[3] R. Sepulchre, M. Jankovic, and P. Kokotovic Constructive Nonlinear Control. Springer-Verlag, 1997.
[4] R. Teel. “Using Saturation to stabilize a class of single-input partially linear composite systems”. IFAC NOLCOS'92 Symposium, pages 369-374, June 1992.
2-12如图2-28所示双水箱系统中，in q 为流入水箱1的液体流量，out q 为流出水箱2的液体流量，试依据液容与液阻的概念，建立112()[(),(),(),()]out in Q s Q s H s Q s H s ∝的系统动态结构图。

解：根据液容和液阻的概念，可分别列出两个水箱的数学模型
1
1122112
112
2in out out dh C q q dt dh C q q dt h h q R h q R ⎧=-⎪⎪
⎪=-⎪⎪
⎨-⎪=
⎪
⎪⎪=⎪⎩
对上式进行在零初始条件下进行拉普拉斯变换得
11122
1121
122()()()
()()()
()()()()
()in out out
C sH s Q s Q s C sH s Q s Q s H s H s Q s R H s Q s R =-⎧⎪=-⎪⎪-=⎨⎪
⎪=⎪⎩
化简后可得
21122112221()1
()()1
out in Q s Q s R C R C s R C R C R C s =++++ 122()1
()1out Q s Q s R C s =+ 11222()1
()1
out Q s H s R R C s R =++ 22
()1
()out Q s H s R =
1 ()s
2()s。