控制系统数字仿真与CAD-课后习题答案

第一章习题

IT什么是仿真？它所遵循的基本原则是什么？

答：仿真是建立在控制理论，相似理论，信息处理技术和计算技术等理论基础之

上的，以计算机和其他专用物理效应设备为工具，利用系统模型对真实或假

想的系统进行试验，并借助专家经验知识，统计数据和信息资料对试验结果

进行分析和研究，进而做出决策的一门综合性的试验性科学。

它所遵循的基本原则是相似原理。

1-2在系统分析与设计中仿真法与解析法有何区别？各有什么特点？

答：解析法就是运用已掌握的理论知识对控制系统进行理论上的分析，计算。它

是一种纯物理意义上的实验分析方法，在对系统的认识过程中具有普遍意

义。由于受到理论的不完善性以及对事物认识的不全面性等因素的影响，其

应用往往有很大局限性。

仿真法基于相似原理，是在模型上所进行的系统性能分析与研究的实验方

法。

b3数字仿真包括那几个要素？其关系如何？

答：通常情况下，数字仿真实验包括三个基本要素，即实际系统，数学模型与计

算机。由图可见，将实际系统抽象为数学模型，称之为一次模型化，它还涉

及到系统辨识技术问题，统称为建模问题；将数学模型转化为可在计算机上

运行的仿真模型，称之为二次模型化，这涉及到仿真技术问题，统称为仿真

实验。

1-4为什么说模拟仿真较数字仿真精度低？其优点如何？。

答：由于受到电路元件精度的制约和容易受到外界的干扰，模拟仿真较数字仿真

精度低

但模拟仿真具有如下优点：

（1）描述连续的物理系统的动态过程比较自然和逼真。

（2）仿真速度极快，失真小，结果可信度高。

（3）能快速求解微分方程。模拟计算机运行时各运算器是并行工作的，模

拟机的解题速度与原系统的复杂程度无关。

（4）可以灵活设置仿真试验的时间标尺，既可以进行实时仿真，也可以进

行非实时仿真。

（5）易于和实物相连。

1-5什么是CAD技术？控制系统CAD可解决那些问题？

答：CAD技术，即计算机辅助设计（ComputerAidedDesign），是将计算机高速

而精确的计算能力，大容量存储和数据的能力与设计者的综合分析，逻辑判

断以及创造性思维结合起来，用以快速设计进程，缩短设计周期，提高设计

质量的技术。

控制系统CAD可以解决以频域法为主要内容的经典控制理论和以时域法为

主要内容的现代控制理论。止匕外，自适应控制，自校正控制以及最优控制等

现代控制测略都可利用CAD技术实现有效的分析与设计。

『6什么是虚拟现实技术？它与仿真技术的关系如何？

答：虚拟现实技术是一种综合了计算机图形技术，多媒体技术，传感器技术，显

示技术以及仿真技术等多种学科而发展起来的高新技术。

1-7什么是离散系统？什么是离散事件系统？如何用数学的方法描述它们？

答：本书所讲的“离散系统”指的是离散时间系统，即系统中状态变量的变化仅

发生在一组离散时刻上的系统。它一般采用差分方程，离散状态方程和脉冲

传递函数来描述。

离散事件系统是系统中状态变量的改变是由离散时刻上所发生的事件所驱

动的系统。这种系统的输入输出是随机发生的，一般采用概率模型来描述。

1-8如图1T6所示某卫星姿态控制仿真实验系统，试说明：

（1）若按模型分类，该系统属于那一类仿真系统？

（2）图中“混合计算机”部分在系统中起什么作用？

（3）与数字仿真相比该系统有什么优缺点？

答：（1）按模型分类，该系统属于物理仿真系统。

（2）混合计算机集中了模拟仿真和数字仿真的优点，它既可以与实物连接进

行实时仿真，计算一些复杂函数，又可以对控制系统进行反复迭代计算。其

数字部分用来模拟系统中的控制器，而模拟部分用于模拟控制对象。

（4）与数字仿真相比，物理仿真总是有实物介入，效果逼真，精度高，具有实

时性与在线性的特点，但其构成复杂，造价较高，耗时过长，通用性不强。

太阳模拟器

星光模拟器地球模拟器

射

频混合计算机

模

拟数字部分接口模拟部分

器卫星动力学

角度读转台电子

三轴机械转台力矩器

出装置驱动器

指令与控制台

题1-8卫星姿态控制仿真试验系统

第二章习题

2-1思考题：

(1)数学模型的微分方程，状态方程，传递函数，零极点增益和部分分式五种

形式，各有什么特点？

(2)数学模型各种形式之间为什么要互相转换？

(3)控制系统建模的基本方法有哪些？他们的区别和特点是什么？

(4)控制系统计算机仿真中的“实现问题”是什么含意？

(5)数值积分法的选用应遵循哪几条原则？

答：(1)微分方程是直接描述系统输入和输出量之间的制约关系，是连续控制系

统其他数学模型表达式的基础。状态方程能够反映系统内部各状态之间的相互关

系，适用于多输入多输出系统。传递函数是零极点形式和部分分式形式的基础。

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

析系统的动态过程。

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

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

型法。机理模型法就是对已知结构，参数的物理系统运用相应的物理定律或定理，

经过合理的分析简化建立起来的各物理量间的关系。该方法需要对系统的内部结

构和特性完全的了解，精度高。统计模型法是采用归纳的方法，根据系统实测的

数据，运用统计规律和系统辨识等理论建立的系统模型。该方法建立的数学模型

受数据量不充分，数据精度不一致，数据处理方法的不完善，很难在精度上达到

更高的要求。混合法是上述两种方法的结合。

(4)“实现问题”就是根据建立的数学模型和精度，采用某种数值计算方法，

将模型方程转换为适合在计算机上运行的公式和方程，通过计算来使之正确的反

映系统各变量动态性能，得到可靠的仿真结果。

(5)数值积分法应该遵循的原则是在满足系统精度的前提下，提高数值运

算的速度和并保证计算结果的稳定。

2-2.用matlab语言求下列系统的状态方程、传递函数、零极点增益、和部分分

式形式的模型参数，并分别写出其相应的数学模型表达式：

§3+7/+24s+24

(1)

+1053+35sN+50s+24

-2.25-5-1.25-0.5-彳

2.25-4.25-1.25-0.252

(2)x=u

0.25-0.5-1.25-12

1.25-1.75-0.25-0.75_o_

y=[0202]X

(1)解：(1)状态方程模型参数:编写matlab程序如下

»num=[l72424];

»den=[l10355024];

»[ABCD]=tf2ss(num,den)

--10-35-50-24'~r

10000

得到结果：A=3=,C=[172424],D=[0]

01000

00100

--10-35-50-24r

10000

所以模型为：x=X+u,y=[l724241X

01000，LJ

00100

⑵零极点增益:编写程序»num=[l72424];

»den=[l10355024];

»[ZPK]=tf2zp(num,den)

得到结果Z=-2.7306+2.8531,-2.7306-2.8531i,-1.5388

P=-4,-3,-2,-l

K=1

(3)部分分式形式:编写程序>>num=[l72424];

»den=[l10355024];

»[RPH]=residue(num,den)

得到结果R=4.0000,-6.0000,2.0000,1.0000

P=-4.0000,-3.0000,-2.0000,-1.0000

H=[]

4-621

G(s)=+++——

5+4s+35+2s+1

(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=[4220]';

»C=[0202];

»D=[0];

»[numden]=ss2tf(A,B,C,D)

得到结果

num=04.000014.000022.000015.0000

den=1.00004.00006.25005.25002.2500

4s3+14s2+22s+15

G(s)=;；

s4+4S3+6.25S2+5.25s+2.25

⑵零极点增益模型参数:编写程序》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=[4220]);

»C=[0202];

»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

表达式G⑻i年+―)卜+1+12停)

(s+0.5-0.866i)(s+0.5+0.866i)(s+1.5)

(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=[4220]';

»C=[0202];

»D=[0];

>>[numden]=ss2tf(A,B,C,D)

>>[R,P,H]=residue(num,den)

得至U结果R=4.0000-0.00000.0000-2.3094i0.0000+

2.3094i

P=-1.5000-1.5000-0.5000+0.8660i-0.5000-

0.8660i

H=[]

―、42.3094z2.3094z

G(s)=1

s+1.55+0.5-0.866/5+0.5+0.866z

2-3.用欧拉法求下面系统的输出响应y(t)在04W1上，h=0.1时的数值。

y'=-y/(o)=i

要求保留4位小数，并将结果与真解了(/)=丁比较。

%+i=%+力*/&，%)

解：欧拉法y=/(如力)(前向欧拉法，可以自启动)其几何意义：把f(t,y)

次o)=%

在%，力］区间内的曲边面积用矩形面积近似代替。利用matlab提供的m文件编

程，得到算法公式。如下所示

(1)m文件程序为h=0.1;

dispC函数的数值解为5％显示“中间的文字％

disp('y=');%同上％

y=i;

fort=O:h:l

m=y;

disp(y);%显示y的当前值％

y=m-m*h;

end

保存文件q2.m

在matalb命令行中键入>>q2

得到结果函数的数值解为

y=10.90000.81000.72900.65610.59050.53140.47830.43050.3874

0.3487

(2)另建一个m文件求解了=/'在te[0,1]的数值(％了=/是

/=-y,y(0)=l的真解%)

程序为h=0.1;

dispC函数的离散时刻解为)；

disp('y=');

fort=0:h:l

y=exp(-t);

disp(y)；

end保存文件q3.m

在matalb命令行中键入>>q3

函数的离散时刻解为

y=10.90480.81870.74080.67030.60650.54880.49660.44930.4066

0.3679

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

欧10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487

拉

真10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679

值

误0-0.0048-0.0007-0.0118-0.0142-0.0160-0.0174-0.0183-0.0188-0.0192-0.0192

差

显然误差与"为同阶无穷小，欧拉法具有一阶计算精度，精度较低，但算法简

单。

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

hZ77、

"+1+~(^1+左2)

解:我们经常用到同肛校正正的二阶龙-格库塔法,<%=//,%)此

&碣)

、/(")=y

方法可以自启动，具有二度，几何意义：把f(t,y)在心,%］区间内

的曲边面积用上下底为《和A”、高为h的梯形面积近似代替。利用matlab

提供的m文件编程，得到算法公式。如下所示

(Dm文件程序为h=0.1;

dispC函数的数值解为)；

disp('y=');

y=i;

fort=O:h:l

disp(y);

kl=-y;

k2=-(y+kl*h);

y=y+(kl+k2)*h/2;

end

保存文件q4.m

在matlab的命令行中键入>>q4显示结果为

函数的数值解为

y=10.90500.81900.74120.67080.60710.54940.49720.45000.4072

0.3685

(2)比较欧拉法与二阶龙格-库塔法求解.(误差为绝对值)

真10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679

值

龙10.90500.81900.74120.67080.60710.54940.49720.45000.40720.3685

库

误00.00020.00030.00040.00050.00060.00060.00060.00070.00060.0006

差

明显误差为“得同阶无穷小，具有二阶计算精度，而欧拉法具有以阶计算精度,

二阶龙格-库塔法比欧拉法计算精度高。

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

h

〕4+i=+乙(/+2k2+2k3+k4)

0

k\=f(tQk)

5

解：四阶龙格-库塔法表达式■心=//+%以+g左)，其截断误差为h

7r/h〃7、

左3=/&+,，力+万人2)

《秘3)

同阶无穷小，当h步距取得较小时，误差是很小的.

⑴编辑m文件程序h=0.1;

dispC四阶龙格-库塔方法求解函数数值解为上

disp('y=');

y=i;

fort=O:h:l

disp(y)；

kl=-y;

k2=-(y+kl*h/2);

k3=-(y+k2*h/2);

k4=-(y+k3*h);

产y+(kl+2*k2+2*k3+k4)*h/6;

end保存文件q5.m

在matlab命令行里键入》q5

得到结果四阶龙格-库塔方法求解函数数值解为

y=10.90480.81870.74080.67030.60650.54880.49660.44930.4066

0.3679

(2)比较这几种方法：

对于四阶龙格-库塔方法

真10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679

值

龙10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679

库

误00000000000

差

显然四阶龙格-库塔法求解精度很高，基本接近真值。三种方法比较可以得到

精度(四阶)〉精度(二阶)〉精度(欧拉)

X1

2-6.已知二阶系统状态方程为

写出取计算步长为h时，该系统状态变量X=[X],%]的四阶龙格-库塔法递推关

系式。

加=%+沙+24+2%+h)

左=/(%%)

解：四阶龙格-库塔法表达式卜2=+g，力

7r/卜〃7、

左3=/(。+]，%+万人2)

左4=/&+〃，”+跟3)

所以状态变量的递推公式可以写作:

ananA

A=,B=，X=可以写成X=4R+3M

^^21^^22b2

h

X=X-\k2k+2k.+kJ

Mk+o[+2

k]-AX卜+Bu

贝1J递推形式1左2=4(X^+勺*///2)+5"

k3=A(Xk+k2^h/2)+Bu

左4=A(Xk+&*〃)+5M

2-7单位反馈系统的开环传递函数已知如下

〜、55+100

G(s)=1

s(s+4.6)(s+3.4s+16.35)

用matlab语句、函数求取系统闭环零极点，并求取系统闭环状态方程的

可控标准型实现。

解：已知开环传递函数，求得闭环传递函数为

E、55+100

G(s)=5

s(s+4.6)(1+3.4s+16.35)+5s+100

在matlab命令行里键入》a=[l0];

»b=[l4.6];

»c=[l3.416.35];

»d=conv(a,b)；

»e=conv(d,c)

e=1.00008.000031.990075.21000

»f^[0005100];

»g=e+f

g=1.00008.000031.990080.2100100.0000

%以上是计算闭环传递函数的特征多项式％

»p=roots(g)%计算特征多项式的根，就是闭环

传递函数的极点％

p=

-0.9987+3.0091i

-0.9987-3.0091i

-3.0013+0.9697i

-3.0013-0.9697i

»m=[5100];

»z=roots(m)

z=-20%计算零点％

综上：当闭环传函形如G(s)=…+*s+”时,可控标准型为:

nn

s+axs+...+*-s+an

010...o'0

001...00

；C=[b,

A=\＞；B=b吁i

0010

-an一%,I」」

Xi

010

001

000

-100-80.21-31.99

所以可控标准型是

r=[-ioo500]+[0]M

2-8用matlab语言编制单变量系统三阶龙格-库塔法求解程序，程序入口要求能

接收状态方程各系数阵(A,B,C,D)，和输入阶跃函数r(t户R*l(t);程序出口应给

出输出量y(t)的动态响应数值解序列

为，九……。

解：m文件为：functiony=hs(A,B,C,D,R,T,h)%T为观测时间,h为计算步长，R

为输入信号幅值％

dispC数值解为)；

y=o；

r=R;

x=[0;0;0;0];

N=T/h;

fort=l:N;

kl=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).

0100

0010

解:,C=[-100500],D=[0]

0001

-100-80.21-31.99-8

在命令行里键入》A=[0100

0010

0001

-100-80.21-31.99-8];

»B=[0001]';

»C=[-100500];

»D=[0];

»T=1;

»R=l;

»h=0.01;

»y=hs(A,B,C,D,R,T,h)

数值解为

0

8.3333e-007

5.8659e-006

1.8115e-005

3.9384e-005

7.0346e-005

o%仅取一部分％

2-10用式(2-34)梯形法求解试验方程j=--v,分析对计算步长h有何限制,

T

说明h对数值稳定性的影响。

h7、

%+1="+5(Z74+k2)

,1

解：编写梯形法程序为

左=一yk

T

71/17、

左2=一一(力一一V，)

TT

hh2

得到加="1二+彳)稳定系统最终渐进收敛。

系统稳定则T+营

<1计算得0<7z<2ro

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

2-11如图2-27所示斜梁滚球系统，若要研究滚球在梁上的位置可控性，需首先

建立其数学模型，已知力矩电机的输出转矩M与其电流i成正比，横梁为均匀

可自平衡梁(即当电机不通电且无滚球时，横梁可处于6=0的水平状态)，是建

立系统的数学模型，并给出简化后系统的动态结构图。

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

人,//分别表示棒的惯量、球的质量和球的惯量。则球质心的位置和速度为

xc=(xcos0,xsin0)

vc=(vcos^-xtysin^,vsin^+x^cos0)

其中文=v,0=CDO因而动能的移动部分为

因而动能的移动部分为「=刎”]("%2)

球棒系统的旋转动能为K0,=L02+_l/2(与

22r

因而，系统总的动能K=K.s+Ke等于

K=—(/.+mx2)a)2+—Amv2

22

其中4=1+与＞1为常数。

此系统的拉格朗日方程组为

fd.dT.dT.a

dtdxa

d,dTdT.

——(--)=ki-mgcosf)

dtQQdO

综合以上公式的系统的方程组为

mXx—mxO2+mgsin(6)=0

(，+mx2)O+2mxx0+mgxcos(^)=ki

设系统在平衡点附近方土0,cosd^l,sm0^0,则系统方程可化为

mAx+mgO=0

(4+mx2)3+mgx=ki

对上式进行拉普拉斯变换并化简后可得到支应。

/(s)

参考文献：

[1]Hauser,S.Sestry,andP.Kokotovic."Nonlinearcontrolviaapproximate

input-outputlinearization”.IEEETrans,onAutomaticControl,vol.37:pp.392-398,

1992.

[2]R.Sepulchre."Slowpeakingandlow-gaindesignsforglobalstabilizationof

nonlinearsystems”,submittedforIEEETAC1999.

[3]R.Sepulchre,M.Jankovic,andP.KokotovicConstructiveNonlinearControl.

Springer-Verlag,1997.

[4]R.Teel."UsingSaturationtostabilizeaclassofsingle-inputpartiallylinear

compositesystems".IFACNOLCOS'92Symposium,pages369-374,June1992.

2-12如图2-28所示双水箱系统中，孙“为流入水箱1的液体流量，/“，为流出水

箱2的液体流量，试依据液容与液阻的概念，建立

Qout⑻8口加⑸,凡(S),0⑻,9⑸］的系统动态结构图。

解：根据液容和液阻的概念，可分别列出两个水箱的数学模型

q=&

J-Out

、火2

对上式进行在零初始条件下进行拉普拉斯变换得

'。皿⑻&⑻/⑻

Gs凡⑻=&(s)-。。“⑻

%

以心)=等

、火2

化简后可得

&“(s)_________________]

2

。加(s)&3R2c2s+(R£+R2C2+&G)s+1

Q“,(s)1

。⑻R2C2S+1

Q”⑻1

H](s)&R2c2s+7?2+1

Q”,(s)_i

凡(s)R2

第三章习题

4-2设典型闭环结构控制系统如图4-47所示，当阶跃输入幅值&=20时，用sp4」.m求取

输出y«）的响应。

解：用sp4_Lm求解过程如下：

在MATLAB语言环境下，输入以下命令语句

»a=[0.0160.8643.273.421];

»b=[3025];

»X0=[0000];%系统状态向量初值为零

»V=2;%反馈系数v=2

»n=4;

»T0=0;Tf^l0;

»h=0.01;R=20;%仿真步长h=0.01,阶跃输入幅值

7?=20

»sp4_l%调用sp4_l.m函数

»plot(t,y)

运行结果为：

附：sp4」.m函数为

b=b/a(1);a=a/a(1);A=a(2:n+1);

A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];

B=[zeros(l,n-1)」T

ml=length(b);

C=[fliplr(b),zeros(1,n-m1)];

Ab=A-B*C*V;

X=X01;

y=0;t=T0;

N=round((Tf-T0)/h);

fori=l:N

Kl=Ab*X+B*R;

K2=Ab*(X+h*Kl/2)+B*R;

K3=Ab*(X+h*K2/2)+B*R;

K4=Ab*(X+h*K3)+B*R;

X=X+h*(Kl+2*K2+2*K3+K4)/6;

y=[y,c*x]；

t=[t,t(i)+h];

end

4-4系统结构图如图4-48,写出该系统的联结矩阵少和叫,并写出联结矩阵非零元素

阵皈。

解：根据图4-48中％,尤拓扑连结关系，可写出每个环节输入%受哪些环节输出入的

影响，

现列如入下：

%=%

u2=yl-y9

W3=>2

以=%一%

<

4=”一％0

%=y6

%=以

w9=y^

、%o=y-：

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

U=WY+W^0

-o0000000000-

〃2100000000-10%0

u010000000000

3y3

"400100000-1000

y4

%000100000000

—*+*为

U0000100000-10

6y6

"700000100000%0

的00000100000以0

u900000010000%0

“1000000010000_%。0

101

211

-

00000000000「「

29-1

100000000-100

321

010000000000

431

00100000-1000

48-1

000100000000

即用=，%=，w=541

0000100000-10u

651

000001000000

610-1

000001000000

761

000000100000

861

000000100000

-971

1071

4-6若系统为图4-5b双输入-双输出结构，试写出该系统的联接矩阵少，Wo,说明应注意

什么？

%必

1►2

%yn6

—►44►6

解：根据图4-5b中/，天拓扑连结关系，可列写如下关系式:

%

忧2=%

%-%

<

%=%+%

牝=%

4=y4

转换成矩阵形式为

001

000

000即

100乂2

000

010

0000

1000

0100

所以联接矩阵少=

0010

0000

0001

此时应注意输入联接矩阵w0变为6x2型。

4-8求图4-49非线性系统的输出响应y（t）,并与无非线性环节情况进行比较。

解：（1）不考虑非线性环节影响时，求解过程如下：

1）先将环节编号标入图中。

2）在MATLAB命令窗口下，按编号依次将环节参数输入P阵;

»P=[0,110.51;01200;2110;10110）;

3）按各环节相对位置和联接关系，有联接矩阵如下：

101

「000-f丁

14-1

10000

w=，%=所以非零元素矩阵wu=211

01000

321

0010_0,

431

»WIJ=[101;14-1;211;321;431];

4）由于不考虑非线性影响，则非线性标志向量和参数向量均应赋零值;

»Z=[0000];S=[0000];

5）输入运行参数：开环截至频率4/约为1,故计算步长h取经验公式值，即

/z<—,―=0.02,取h=0.01；每0.25秒输出一点。故取，=25。

50g

»h=0.01;

»L1=25;

»n=4;

»T0=0

»Tf^20;

»nout=4;

»Y0=10;

»sp4_4;

»plot(t,y,T')

»holdon

运行结果如图中红色实线所示。

⑵考虑非线性环节N影响时，只需将非线