《计算机仿真技术与CAD》习题答案

.

第0章

绪论

0-1什么是仿真？它所遵循的基本原则是什么？

答：

仿真是建立在控制理论、相似理论、信息处理技术和计算机技术等理论基础之上的，以计算机和其他专用物理效应设备为工具，利用系统模型对真实或假想的系统进行试验，并借助专家经验知识、统计数据和信息资料对试验结果进行分析和研究，进而做出决策的一门综合性的试验性科学。

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

0-2仿真的分类有几种？为什么？

答：

依据相似原理来分：物理仿真、数学仿真和混合仿真。

物理仿真：就是应用几何相似原理，制作一个与实际系统相似但几何尺寸较小或较大的物理模型（例如飞机模型放在气流场相似的风洞中）进行实验研究。

数学仿真：就是应用数学相似原理，构成数学模型在计算机上进行研究。它由软硬件仿真环境、动画、图形显示、输出打印设备等组成。

混合仿真又称数学物理仿真，它是为了提高仿真的可信度或者针对一些难以建模的实体，在系统研究中往往把数学仿真、物理仿真和实体结合起来组成一个复杂的仿真系统，这种在仿真环节中有部分实物介入的混合仿真也称为半实物仿真或者半物理仿真。

0-3比较物理仿真和数学仿真的优缺点。

答：

在仿真研究中，数学仿真只要有一台数学仿真设备（如计算机等），就可以对不同的控制系统进行仿真实验和研究，而且，进行一次仿真实验研究的准备工作也比较简单，主要是受控系统的建模、控制方式的确立和计算机编程。数学仿真实验所需的时间比物理仿真大大缩短，实验数据的处理也比物理仿真简单的多。

与数学仿真相比，物理仿真总是有实物介入，效果直观逼真，精度高，可信度高，具有实时性与在线性的特点；但其需要进行大量的设备制造、安装、接线及调试工作，结构复杂，造价较高，耗时过长，灵活性差，改变参数困难，模型难以重用，通用性不强。

0-4简述计算机仿真的过程。

答：

第一步：根据仿真目的确定仿真方案

根据仿真目的确定相应的仿真结构和方法，规定仿真的边界条件与约束条件。

第二步：建立系统的数学模型

对于简单的系统，可以通过某些基本定律来建立数学模型。而对于复杂的系统，则必须利用实验方法通过系统辩识技术来建立数学模型。数学模型是系统仿真的依据，所以，数学模型的准确性是十分重要。

.

.

第三步：建立仿真模型

即通过一定算法对原系统的数学模型进行离散化处理，就连续系统言，就是建立相应的差分方程。

第四步：编制仿真程序

对于非实时仿真，可用一般高级语言或仿真语言。对于快速的实时仿真，往往需要用汇编语言。

第五步：进行仿真实验并输出仿真结果

通过实验对仿真系统模型及程序进行校验和修改，然后按系统仿真的要求输出仿真结果。

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

答：

CAD技术，即计算机辅助设计（ComputerAidedDesign），是将计算机高速而精确的计算能力、大容量存储和数据处理能力与设计者的综合分析、逻辑判断以及创造性思维结合起来，以加快设计进程、缩短设计周期、提高设计质量的技术。

控制系统CAD可以解决以频域法为主要内容的经典控制理论和以时域法为主要内容的现代控制理论。此外，自适应控制、自校正控制以及最优控制等现代控制策略都可利用CAD技术实现有效的分析和设计。

.

第1章

.

仿真软件——MATLAB

1-1对于矩阵A=[12;34]，MATLAB以下四条命令：

A.^(0.5)；A^(0.5)；sqrt(A)；sqrtm(A)

所得结果相同吗？它们中哪个结果是复数矩阵，为什么？

答：

A.^(0.5)=[1.00001.4142;1.73212.0000]；

A^(0.5)=[0.5537+0.4644i0.8070-0.2124i;1.2104-0.3186i1.7641+0.1458i]；

sqrt(A)=[1.00001.4142;1.73212.0000]；

sqrtm(A)=[0.5537+0.4644i0.8070-0.2124i;1.2104-0.3186i1.7641+0.1458i]；

其中，“A.^(0.5)”表示向量的乘方，“A^(0.5)”表示矩阵的乘方，“sqrt(A)”只定义在矩阵的单个元素上，即分别对矩阵的每个元素进行运算，“sqrtm(A)”表示对矩阵（方阵）的超越函数进行运算。

1-4求二元函数方程组：sin(x-y)=0,cos(x+y)=0的解。

答：

>>[x,y]=solve('sin(x-y)=0','cos(x+y)=0','x','y')

x=

-1/4*pi

1/4*pi

y=

-1/4*pi

1/4*pi

1-5求函数y(t)=exp(-t)*|sin[cost]|的最大值（0<=t<inf）。

答：

>>f='(-1)*exp(-(abs(x)))*abs(sin(cos(abs(x))))';

>>x=fminsearch(f,0),ymax=exp(-(abs(x)))*abs(sin(cos(abs(x))))

x=

0

ymax=

0.8415

1-6设D2y-3Dy+2y=x,y(0)=1,Dy(0)=0，求y(0.5)的值。

答：

>>f='D2y-3*Dy+2*y=x';g=dsolve(f,'y(0)=1,Dy(0)=0','x');x=0.5;y=eval(g)

y=

0.6100

1-7求方程cos(t)^2*exp(-0.1t)=0.5t的解。

.

.

答：

>>t1=solve('cos(t)^2*exp(-0.1*t)=0.5*t','t');t=eval(t1)

t=

0.8329

1-8求方程组：x^2+y^2=1,xy=2的解。

答：

>>[x,y]=solve('x^2+y^2=1','x*y=2','x','y')

x=

-1/2*(1/2*5^(1/2)+1/2*i*3^(1/2))^3+1/4*5^(1/2)+1/4*i*3^(1/2)

-1/2*(1/2*5^(1/2)-1/2*i*3^(1/2))^3+1/4*5^(1/2)-1/4*i*3^(1/2)

-1/2*(-1/2*5^(1/2)+1/2*i*3^(1/2))^3-1/4*5^(1/2)+1/4*i*3^(1/2)

-1/2*(-1/2*5^(1/2)-1/2*i*3^(1/2))^3-1/4*5^(1/2)-1/4*i*3^(1/2)

y=

1/2*5^(1/2)+1/2*i*3^(1/2)

1/2*5^(1/2)-1/2*i*3^(1/2)

-1/2*5^(1/2)+1/2*i*3^(1/2)

-1/2*5^(1/2)-1/2*i*3^(1/2)

1-9求f(kT)=kexp(-akT)的Z变换表达式。

答：

>>symsktz;f=k*exp(-a*t);F=ztrans(f,t,z)

f=

k*z/exp(-a)/(z/exp(-a)-1)

1-10求一阶微分方程Dx=ax+by(t),x(0)=x0的解。

答：

>>f='Dx=a*x+b*y';x=dsolve(f,'x(0)=x0','t')

x=

-b*y/a+exp(a*t)*(b*y+x0*a)/a

1-12求以下方程组边值问题的解。

Df=3f+4g,Dg=-4f+3g,f(0)=0,g(0)=1

答：

>>f='Dx1=3*x1+4*x2,Dx2=-4*x1+3*x2';[x1,x2]=dsolve(f,'x1(0)=0,x2(0)=1','t')

x1=

exp(3*t)*sin(4*t)

x2=

exp(3*t)*cos(4*t)

.

.

.

xt

xt

ut

1

1

1

-2-3

11

1

11

xt

.

第2章控制系统的数学模型及其转换

2-1已知系统的传递函数为

Gs

s

3

s2

6s

s1211s6

试用MATLAB建立其状态空间表达式。

答：

>>num=[111];den=[16116];[A,B,C,D]=tf2ss(num,den)

A=

-6-11-6

100

010

B=

1

0

0

C=

111

D=

0

2-2已知系统的状态空间表达式为

0110

x2tx2tu2t

10xt

yt

2

试用MATLAB求其传递函数阵。

答：

>>A=[01;-2-3];B=[10;11];C=[10;11];D=zeros(2,2);

>>[num1,den1]=ss2tf(A,B,C,D,1),[num2,den2]=ss2tf(A,B,C,D,2)

num1=

00004.0000

00002.0000

den1=

132

num2=

00.00001.0000

01.00001.0000

.

xt

x

t

1

1

1

1-1

.

den2=

132

2-3已知两子系统的传递函数分别为

Gs1

1

s1s2

，

Gs

2

1

ss3

试利用MATLAB求两子系统串联和并联时系统的传递函数。答：

>>num1=1;den1=[132];num2=1;den2=[130];>>[num,den]=series(num1,den1,num2,den2)

num=

00001

den=

161160

>>num1=1;den1=[132];num2=1;den2=[130];>>[num,den]=parallel(num1,den1,num2,den2)

num=

00262

den=

161160

2-4设系统的状态空间表达式为

011

x2t-2-3x2t2utxt

yt30

x2t

若取线性变换阵

设新的状态变量为

xP

1

x

11

P

，则利用MATLAB求在新状态变量下，系统状态空间表达式。

答：

>>A=[01;-2-3];B=[1;2];C=[30];D=[0];P=[11;1-1];>>[A1,B1,C1,D1]=ss2ss(A,B,C,D,P)

A1=

-20

3-1

B1=

3

-1

C1=

1.50001.5000

.

1

1

13

1

1

.

D1=

0

2-5已知离散系统状态空间表达式

xk101xkx2k1x2kxk

yk11

x2k

试用MATLAB求其系统的脉冲传递函数。

答：

0

ut

>>A=[01;13];B=[0;1];C=[11];D=0;T=1;[A1,B1,C1,D1]=c2dm(A,B,C,D,T)

A1=

2.95987.3357

7.335724.9669

B1=

1.9598

7.3357

C1=

11

D1=

0

.

xt

x

t

1

1

1

.

第3章连续系统的数字仿真

3-1已知线性定常系统的状态空间表达式为

012

x2t-5-6x2t0ut

xt

yt12

x2t

且初始状态为零，试利用四阶-龙格库塔法求系统的单位阶跃响应。答：

%ex3_1.m

r=1;A=[01;-5-6];B=[2;0];C=[12];d=0;

Tf=5;h=0.1;

x=[zeros(length(A),1)];y=0;t=0;

fori=1:Tf/h

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+h*(K1+2*K2+2*K3+K4)/6;

=[y;C*x];t=[t;t(i)+h];

end

plot(t,y)

3-2设单位反馈系统的开环传递函数

Gs

4

ss2

试利用二阶-龙格库塔法求系统的单位阶跃响应。

.

.

答：

%ex3_2.m

r=1;numo=4;deno=[1,2,0];[num,den]=cloop(numo,deno);

[A,b,C,d]=tf2ss(num,den);

Tf=5;h=0.1;

x=[zeros(length(A),1)];y=0;t=0;

fori=1:Tf/h

K1=A*x+b*r;

K2=A*(x+h*K1)+b*r;

x=x+h*(K1+K2)/2;y=[y;C*x];t=[t;t(i)+h];

end

plot(t,y)

3-4利用input()函数修改例3-1所给程序ex3_1.m，将其中给定的参数r，numo，deno，numh和denh利用键盘输入，使其变为连续控制系统面向传递函数的通用数字仿真程序。

答：

.

.

3-5利用input()函数修改例3-2所给程序ex3_2.m，将其中给定的参数r，P，W，W0和Wc利用键盘输入，使其变为连续控制系统面向结构图的通用数字仿真程序。答：

.

.

第4章连续系统按环节离散化的数字仿真

4-1已知非线性习题如图题4-1所示，试利用连续系统按环节离散化的数字仿真方法，求输出量y的动态响应，并与无非线性环节进行比较。（图略）

答：

%ex4_1.m%主程序

R=10;

P=[0.110.5155;011000;212000;10110000];

W=[000-1;1000;0100;0010];

W0=[1;0;0;0];Wc=[0001];

Tf=25;T=0.02;

A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);

FZ=P(:,5);S=P(:,6);

n=length(A);

fori=1:n

if(A(i)~=0)

if(B(i)==0)

E(i)=0;F(i)=0;G(i)=0;H(i)=0;

L(i)=(C(i)+D(i)/T)/A(i);Q(i)=-D(i)/(A(i)*T);

else

(i)=exp(-A(i)*T/B(i));

(i)=(D(i)/B(i)-C(i)/A(i))*((1-E(i))*B(i)/(A(i)*T)-1);

(i)=(D(i)/B(i)-C(i)/A(i))*(1+(E(i)-1)*(1+B(i)/(A(i)*T)));

(i)=1;L(i)=D(i)/B(i);Q(i)=0;

end

else

if(B(i)~=0)

E(i)=1;F(i)=0.5*C(i)*T/B(i);G(i)=F(i);

H(i)=1;L(i)=D(i)/B(i);Q(i)=0;

else

disp('A(i)=B(i)=0');

end

end

end

x=[zeros(length(A),1)];x0=x;z=x;

u=[zeros(length(A),1)];u0=u;

y=[zeros(length(Wc(:,1)),1)];t=0;

forj=1:Tf/T

u1=u;u=W*x+W0*R;

fori=1:n

if(FZ(i)~=0)

if(FZ(i)==1)u(i)=saturation(u(i),S(i));end

if(FZ(i)==2)u(i)=deadzone(u(i),S(i));end

.

.

if(FZ(i)==3)[u(i),u0(i)]=backlash(u0(i),u(i),u1(i),S(i));end

if(FZ(i)==4)u(i)=sign1(u(i),S(i));end

end

end

x1=x;

fori=1:n

z(i)=E(i)*z(i)+F(i)*u(i)+G(i)*u1(i);

x(i)=H(i)*z(i)+L(i)*u(i)+Q(i)*u1(i);

end

fori=1:n

if(FZ(i)~=0)

if(FZ(i)==5)x(i)=saturation(x(i),S(i));end

if(FZ(i)==6)x(i)=deadzone(x(i),S(i));end

if(FZ(i)==7)[x(i),x0(i)]=backlash(x0(i),x(i),x1(i),S(i));end

if(FZ(i)==8)x(i)=sign1(x(i),S(i));end

end

end

y=[y,Wc*x];t=[t,t(j)+T];

end

plot(t,y)

%saturation.m%子程序

functionx=saturation(u,s)

if(abs(u)>=s)

if(u>0)x=s;

elsex=-s;

end

else

x=u;

end

修改“P=[0.110.5100;011000;212000;10110000];”

.

.

>>ex4_1

4-2针对例3-2所给线性定常系统，试利用第4章所给程序，求系统的单位阶跃响应，并对其结果进行比较。

答：

>>ex3_2

>>ex4_1

.

.

4-3针对例4-1所给系统，去掉饱和非线性环节后求系统的单位阶跃响应，并与例4-1所得结果进行比较。

答：

>>ex4_1

修改“P=[0.110.5100;011000;212000;10110000];”>>ex4_1

.

.

4-4利用input()函数修改例4-1所给程序ex4_1.m，将其中给定的参数R，P，W，W0和Wc利用键盘输入，使其变为连续控制系统按环节离散化的通用数字仿真程序。答：

略

.

.

第5章采样控制系统的数字仿真

5-1已知采样控制系统的结构图如图题5-1所示（图略）。试利用采样控制系统的数字仿真方法，求当采样周期T=0.1s，且初始状态为零时，离散系统的单位阶跃响应。答：

%ex5_1.m

R=1;

Gr=[1];Fr=[0];

P=[111000;121000];

W=[00;10];W0=[1;0];Wc=[01];

Tf=25;Tm=0.1;

T=0.01;

A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);FZ=P(:,5);S=P(:,6);

n=length(A);n1=length(Fr);m1=length(Gr);

fori=1:n

if(A(i)~=0)

if(B(i)==0)

E(i)=0;F(i)=0;G(i)=0;H(i)=0;

L(i)=(C(i)+D(i)/T)/A(i);Q(i)=-D(i)/(A(i)*T);

else

(i)=exp(-A(i)*T/B(i));

(i)=(D(i)/B(i)-C(i)/A(i))*((1-E(i))*B(i)/(A(i)*T)-1);

(i)=(D(i)/B(i)-C(i)/A(i))*(1+(E(i)-1)*(1+B(i)/(A(i)*T)));

(i)=1;L(i)=D(i)/B(i);Q(i)=0;

end

else

if(B(i)~=0)

E(i)=1;F(i)=0.5*C(i)*T/B(i);G(i)=F(i);

H(i)=1;L(i)=D(i)/B(i);Q(i)=0;

else

disp('A(i)=B(i)=0');

end

end

end

x=[zeros(length(A),1)];x0=x;z=x;

u=[zeros(length(A),1)];u0=u;

y=[zeros(length(Wc(:,1)),1)];

t=0;Ur=[zeros(n1,1)];Er=[zeros(m1,1)];

forij=0:Tf/Tm;

e=R-x(n);Er=[e;Er(1:m1-1)];

ur=-Fr*Ur+Gr*Er;Ur=[ur;Ur(1:n1-1)];

forj=1:Tm/T

u1=u;u=W*x+W0*ur;

.

.

fori=1:n

if(FZ(i)~=0)

if(FZ(i)==1)u(i)=saturation(u(i),S(i));end

if(FZ(i)==2)u(i)=deadzone(u(i),S(i));end

if(FZ(i)==3)[u(i),u0(i)]=backlash(u0(i),u(i),u1(i),S(i));endif(FZ(i)==4)u(i)=sign1(u(i),S(i));end

end

end

x1=x;

fori=1:n

z(i)=E(i)*z(i)+F(i)*u(i)+G(i)*u1(i);

x(i)=H(i)*z(i)+L(i)*u(i)+Q(i)*u1(i);

end

fori=1:n

if(FZ(i)~=0)

if(FZ(i)==5)x(i)=saturatio