控制系统数字仿真与CAD-全习题答案.docx

《控制系统数字仿真与CAD-全习题答案.docx》由会员分享，可在线阅读，更多相关《控制系统数字仿真与CAD-全习题答案.docx（59页珍藏版）》请在课桌文档上搜索。

1、第一章习题1-1什么是仿真？它所遵循的基本原则是什么？答：仿真是建立在控制理论，相似理论，信息处理技术和计算技术等理论基础之上的，以计算机和其他专用物理效应设备为工具，利用系统模型对真实或假想的系统进行试验，并借助专家经验知识，统计数据和信息资料对试验结果进行分析和研究，进而做出决策的一门综合性的试验性科学。它所遵循的基本原则是相似原理。1-2在系统分析与设计中仿真法与解析法有何区别？各有什么特点？答：解析法就是运用已掌握的理论知识对控制系统进行理论上的分析，计算。它是一种纯物理意义上的实验分析方法，在对系统的认识过程中具有普遍意义。由于受到理论的不完善性以及对事物认识的不全面性等因素的影响，

2、其应用往往有很大局限性。仿真法基于相似原理，是在模型上所进行的系统性能分析与研究的实验方法。1-3数字仿真包括那几个要素？其关系如何？答：通常情况下，数字仿真实验包括三个基本要素，即实际系统，数学模型与计算机。由图可见，将实际系统抽象为数学模型，称之为一次模型化，它还涉及到系统辨识技术问题，统称为建模问题；将数学模型转化为可在计算机上运行的仿真模型，称之为二次模型化，这涉及到仿真技术问题，统称为仿真实验。1-4为什么说模拟仿真较数字仿真精度低？其优点如何？。答：由于受到电路元件精度的制约和容易受到外界的干扰，模拟仿真较数字仿真精度低但模拟仿真具有如下优点：（1）描述连续的物理系统的动态过程比较

3、自然和逼真。（2）仿真速度极快，失真小，结果可信度高。（3）能快速求解微分方程。模拟计算机运行时各运算器是并行工作的，模拟机的解题速度与原系统的复杂程度无关。（4）可以灵活设置仿真试验的时间标尺，既可以进行实时仿真，也可以进行非实时仿真。(5)易于和实物相连。1-5什么是CAD技术？控制系统CAD可解决那些问题？答：CAD技术，即计算机辅助设计(ComPUterAidedDesign),是将计算机高速而精确的计算能力，大容量存储和数据的能力与设计者的综合分析，逻辑判断以及创造性思维结合起来，用以快速设计进程，缩短设计周期，提高设计质量的技术。控制系统CAD可以解决以频域法为主要内容的经典控制理

4、论和以时域法为主要内容的现代控制理论。此外，自适应控制，自校正控制以及最优控制等现代控制测略都可利用CAD技术实现有效的分析与设计。1-6什么是虚拟现实技术？它与仿真技术的关系如何？答：虚拟现实技术是一种综合了计算机图形技术，多媒体技术，传感器技术，显示技术以及仿真技术等多种学科而发展起来的高新技术。卜7什么是离散系统？什么是离散事件系统？如何用数学的方法描述它们？答：本书所讲的“离散系统”指的是离散时间系统，即系统中状态变量的变化仅发生在一组离散时刻上的系统。它一般采用差分方程，离散状态方程和脉冲传递函数来描述。离散事件系统是系统中状态变量的改变是由离散时刻上所发生的事件所驱动的系统。这种系

5、统的输入输出是随机发生的，一般采用概率模型来描述。1-8如图1T6所示某卫星姿态控制仿真实验系统，试说明：(1)若按模型分类，该系统属于那一类仿真系统？(2)图中“混合计算机”部分在系统中起什么作用？(3)与数字仿真相比该系统有什么优缺点？答：(1)按模型分类，该系统属于物理仿真系统。(2)混合计算机集中了模拟仿真和数字仿真的优点，它既可以与实物连接进行实时仿真，计算一些复杂函数，又可以对控制系统进行反复迭代计算。其数字部分用来模拟系统中的控制器，而模拟部分用于模拟控制对象。(4)与数字仿真相比，物理仿真总是有实物介入，效果逼真，精度高，具有实时性与在线性的特点，但其构成复杂，造价较高，耗时过

6、长，通用性不强。星光模拟器地球模拟器太阳模拟器射频模拟器混合计算机数字部分接口模拟部分卫星动力学转台电子驱动器指令与控制台题1-8卫星姿态控制仿真试验系统第二章习题2-1思考题：(1)数学模型的微分方程，状态方程，传递函数，零极点增益和部分分式五种形式，各有什么特点？(2)数学模型各种形式之间为什么要互相转换？(3)控制系统建模的基本方法有哪些？他们的区别和特点是什么？(4)控制系统计算机仿真中的“实现问题”是什么含意？(5)数值积分法的选用应遵循哪几条原则？答：(1)微分方程是直接描述系统输入和输出量之间的制约关系，是连续控制系统其他数学模型表达式的基础。状态方程能够反映系统内部各状态之间的

7、相互关系，适用于多输入多输出系统。传递函数是零极点形式和部分分式形式的基础。零极点增益形式可用于分析系统的稳定性和快速性。利用部分分式形式可直接分析系统的动态过程。(2)不同的控制系统的分析和设计方法，只适用于特定的数学模型形式。(3)控制系统的建模方法大体有三种：机理模型法，统计模型法和混合模型法。机理模型法就是对已知结构，参数的物理系统运用相应的物理定律或定理,经过合理的分析简化建立起来的各物理量间的关系。该方法需要对系统的内部结构和特性完全的了解，精度高。统计模型法是采用归纳的方法，根据系统实测的数据，运用统计规律和系统辨识等理论建立的系统模型。该方法建立的数学模型受数据量不充分，数据精

8、度不一致，数据处理方法的不完善，很难在精度上达到更高的要求。混合法是上述两种方法的结合。(4)“实现问题”就是根据建立的数学模型和精度，采用某种数值计算方法,将模型方程转换为适合在计算机上运行的公式和方程，通过计算来使之正确的反映系统各变量动态性能，得到可靠的仿真结果。(5)数值积分法应该遵循的原则是在满足系统精度的前提下，提高数值运算的速度和并保证计算结果的稳定。2-2,用InatIab语言求下列系统的状态方程、传递函数、零极点增益、和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：Cz、53+7j2+241s+24(1)G(s)=-7754+10?+3552+505+24-4.2

9、5-0.5-1.75-1.25-1.25-1.25-0.25y=0 2 0 2 X错误！未找到引用源。2.252.250.251.25(1)解：(1)状态方程模型参数:编写matlab程序如下num=l72424;den=ll10355024;ABCD=tf2ss(num,den)得到结果：A=-10100-35-50-24,B=000,C=172424,D=0010001000-10-35-50-24所以模型为：X=10010000X+00u,y=172424X00100零极点增益:编写程序num=ll724241;den=l10355024;ZPK=tf2zp(num,den)得到结果Z=

10、-2.7306+2.8531,-2.7306-2.853Ii,-1.5388P=-4,-3,-2,-1K=I(3)部分分式形式:编写程序num=72424;den=l10355024;RPH=residue(num,den)得到结果R=4.0000,-6.0000,2.0000,1.0000P=-4.0000,-3.0000,-2.0000,-1.0000H=5+4S+35+25+1(2)解：(1)传递函数模型参数：编写程序A=2.25-5-1.25-0.52.25-4.25-1.25-0.250.25-0.5-1.25-11.25-1.75-0.25-0.75；B=4220J,;C=IO20

11、2J;D=IO;numden=ss2tf(A,B,C,D)得到结果num=04.000014.000022.000015.0000den=1.00004.00006.25005.25002.2500、4s3+14s2+22s+15G(三)=：s4+4S3+6.25S2+5.25s+2.25(2)零极点增益模型参数:编写程序A=2.25-5-1.25-0.52.25-4.25-1.25-0.250.25-0.5-1.25-11.25-1.75-0.25-0.75；B=4220;C=IO202J;D=IO;Z,P,K=ss2zp(A,B,C,D)得至U结果Z=LOOOO+L2247i-LOOOO-

12、1.2247i-1.5000P=-0.5000+0.8660i-0.5000-0.8660i-1.5000-1.5000K=4.0000表达式G(.)=74(sl-L2247i)(slL2247i)(s+0.5-0.866i)(s+0.5+0.866i)(s+1.5)(3)部分分式形式的模型参数：编写程序A=2.25-5-1.25-0.52.25-4.25-1.25-0.250.25-0.5-1.25-11.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)得到结果R=4.0000-0.

13、00000.0000-2.3094i0.0000+2.3094iP=-1.5000-1.5000-0.5000+0.8660i-0.5000-0.8660iH=42.3094z2.3094/G(三)=15+1.55+0.5-0.866/5+0.5+0.866/23用欧拉法求下面系统的输出响应y(t)在0tl上，h=0.1时的数值。y=-y,y()=要求保留4位小数，并将结果与真解NQ=/比较。=*M-)解:欧拉法y=九)(前向欧拉法,可以自启动)其几何意义：把f(t,y)M)=%在片,北1区间内的曲边面积用矩形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示(1)m文件程

14、序为h=0.1;dispC函数的数值解为)；显示“中间的文字disp(y=);%同上y=i;fort=O:h:1m=y;disp(y);%显示y的当前值y=m-m*h;end保存文件q2.m在matalb命令行中键入q2得到结果函数的数值解为y=10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487(2)另建一个m文件求解y=e-在t0,l的数值(y=/是V=-y,y(0)=l的真解)程序为h=0.1;dispC函数的离散时刻解为)；disp(y=);fort=O:h:ly=exp(-t);disp(y);end保存文件q3.

15、m在matalb命令行中键入q3函数的离散时刻解为y=10.90480.81870.74080.67030.60650.54880.49660.44930.406603679比较欧拉方法求解与真值的差别欧拉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

16、显然误差与为同阶无穷小，欧拉法具有一阶计算精度，精度较低，但算法简单。2-4用二阶龙格库塔法求解2-3的数值解，并于欧拉法求得的结果比较。hz.，、1=-(1+2)解:我们经常用到预报校正法的二阶龙格库塔法，k、=f(tk,yk)此%2=/&+%+g).f(hy)=y方法可以自启动，具有二阶计算精度,几何意义：把f(t,y)在匕,为区间内的曲边面积用上下底为人和A讨、高为h的梯形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示(Dm文件程序为h=0.1;dispC函数的数值解为)；disp(y=);y=l;fort=O:h:1disp(y);kl=-y;k2=-(y+kl

17、*h);y=y+(k1+k2)*h2;end保存文件q4.m在matlab的命令行中键入q4显示结果为函数的数值解为y=10.90500.81900.74120.67080.60710.54940.49720.45000.40720.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.000

18、50.00060.00060.00060.00070.00060.0006明显误差为3得同阶无穷小，具有二阶计算精度，而欧拉法具有以阶计算精度,二阶龙格-库塔法比欧拉法计算精度高。2-5.用四阶龙格库塔法求解题23数值解，并与前两题结果相比较。解：四阶龙格库塔法表达式h+=+(K+2k2+2勺+%)O与=/+*+)其截断误差为h与=/4+*”+知%=/&+%，+秘3)同阶无穷小，当h步距取得较小时，误差是很小的.（1）编辑01文件程序h=0.1;dispC四阶龙格-库塔方法求解函数数值解为上dispCy三,）;y=l;fort=O:h:1disp(y);kl=-y;k2=-(y+kl*h2);

19、k3=-(y+k2*h2);k4=-(y+k3*h);y=y+(k1+2*k2+2*k3+k4)*h6;end保存文件q5.m在matlab命令行里键入q5得到结果四阶龙格-库塔方法求解函数数值解为y=10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679(2)比较这几种方法：对于四阶龙格-库塔方法真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679龙库10.90480.81870.74080.67030.60650.54880.49660.44930.

20、40660.3679误差00000000000显然四阶龙格库塔法求解精度很高，基本接近真值。三种方法比较可以得到精度（四阶）精度（二阶）精度（欧拉）2-6.已知二阶系统状态方程为+伪小b2LX2(。)XlOx20写出取计算步长为h时，该系统状态变量X=xpx2的四阶龙格库塔法递推关系式。h+=yfk+公(4+2&+2%+%)O4=/4，)解：四阶龙格库塔法表达式2=+p+).hh.x无3=+产)k4=f(tlc+h,yk+hk3)所以状态变量的递推公式可以写作:an an“21 22,B=4,X=I_”可以写成X=AX+81.JLX2XR=XJk+9&+2&+2%+%)Ol=AXk+Bu则递推

21、形式卜2=4X人+K*力/2)+Bvk3=A(Xk+k2h2)+Bu&=4乂+右次)+32-7单位反馈系统的开环传递函数已知如下、5s+100G(三)=；5(5+4.6)(?+3.45+16.35)用matlab语句、函数求取系统闭环零极点，并求取系统闭环状态方程的可控标准型实现。解：已知开环传递函数，求得闭环传递函数为、55+100G(三)=S(S+4.6)(52+3.45+16.35)+5s+100在matlab命令行里键入a=lO;b=l4.6;c=l3.416.35;d=conv(a,b);e=conv(d,c)e=1.00008.000031.990075.21000f=000510

22、0;g=e+fg=1.00008.000031.990080.2100100.0000%以上是计算闭环传递函数的特征多项式传递函数的极点p=roots(g)%计算特征多项式的根，就是闭环P=-0.9987+3.0091i-0.9987-3.0091i-3.00130.9697i-3.0013-0.9697irn=5100;z=roots(m)z=-20%计算零点综上：当闭环传函形如G(三)=，二+%s+时,可控标准型为:S+S+.+Cln_ySUn所以可控标准型是000-100100-80.21010-31.99r = -100 50 0|XX2+ 10128用matlab语言编制单变量系统三

23、阶龙格库塔法求解程序，程序入口要求能接收状态方程各系数阵(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=Th;fort=l:N;kl=A*x+B*R;k2=A*(x+h*kl3)+B*R;k3=A*(x+2*h*k23)+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

24、,D,R,T,h)得到结果。2-9.用题2-8仿真程序求解题2-7系统的闭环输出响应y(t).0100000100解：A=0001,B=0,C=-100500,D=IOJ-100-80.21-31.99-81在命令行里键入A=010000100001-100-80.21-31.99-8;B=o00iy；C=I-100500;D=IO;T=I;R=I;h=0.01;y=hs(A,B,C,D,R,T,h)数值解为08.3333e-0075.8659e-0061.8115e-0053.9384e-0057.0346e-005O%仅取一部分210用式(234)梯形法求解试验方程y,分析对计算步长h有何

25、限制,说明h对数值稳定性的影响。i=+(+2)解：编写梯形法程序为,1K=二然,1/1,、k2=_(_然力l2h2得到m二M(I-4+3)稳定系统最终渐进收敛。2hIr系统稳定则1-：+今1计算得0为常数。mr此系统的拉格朗日方程组为d月T、T.八=-mgsndtx&d,TTl.n一()=Ki-mgcosdtd综合以上公式的系统的方程组为mx-mx2+mgsin(6)=0(1+mx)+2nxx+Mgxcos(O)=ki设系统在平衡点附近方方0,cos6*l,Sinea凡则系统方程可化为mx+m=0(1+Wtr2)+mgx=ki对上式进行拉普拉斯变换并化简后可得到3O/(三)参考文献：1Haus

26、er,S.Sestry,andP.Kokotovic.ttNonlinearcontrolviaapproximateinput-outputlinearization.IEEETrans,onAutomaticControl,vol.37:pp.392-398,1992.R.Sepulchre.Slowpeakingandlow-gaindesignsforglobalstabilizationofnonlinearsystems”,submittedforIEEETAC1999.3R.Sepulchre,M.Jankovic,andP.KokotovicConstructiveNonlin

27、earControl.Springer-Verlag,1997.41R.Teel.UsingSaturationtostabilizeaclassofsingle-inputpartiallylinearcompositesystems,IFACNOLCOS,92Symposium,pages369-374,June1992.212如图228所示双水箱系统中，外1为流入水箱1的液体流量，%,为流出水箱2的液体流量，试依据液容与液阻的概念，建立QoutG)Qin(三),Hi(三),Qi,%(三)I的系统动态结构图。解：根据液容和液阻的概念，可分别列出两个水箱的数学模型对上式进行在零初始条件下进行

28、拉普拉斯变换得GSHl(三)=20G)-QI(三)C2SH2(三)=Ql(三)-Qtfllt(三)G(三)=HGHH(三)段Qa“=陪K2化简后可得2“(三)=!Qin(三)舄。/?。?/+(RlG+&2。2+&G)5+12mG)1QI(三)R2C2s+(三)1HI(S)+&+1QW(三)=1H2(三)R2第三章习题4-2设典型闭环结构控制系统如图447所示，当阶跃输入幅值R=20时，用SP4.m求取输出y(t)的响应。r(t)解：用SP4.m求解过程如下:在MATLAB语百环境下，输入以下命令语句 a=0.016 0.864 3.27 3.42 1; b=30 25; X0=l0 0 0 0

29、1; V=2; n=4; TO=OjTf= 10; h=0.01;R=20 ;R = 20 sp4_l plot(t)运行结果为：%系统状态向量初值为零%反馈系数u = 2%仿真步长h=0.01,阶跃输入幅值%调用sp4_l.m函数附：SP4_Lm函数为b=ba(1);a=a/a(I);A=a(2:n+1);A=rot90(rot90(eye(n-1,n);-fliplr(A);B=zeros(l,n-l),l;ml=length(b);C=fliplr(b),zeros(l,n-m1);Ab=A-B*C*V;X=X0,;y=O;t=TO;N=round(Tf-TO)h);fori=l:NKI

30、=Ab*X+B*R;K2=Ab*(X+h*Kl+B*R;K3=Ab*(X+h*K2+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;end44系统结构图如图448,写出该系统的联结矩阵W和W0,并写出联结矩阵非零元素阵心。解:根据图4-48中,M拓扑连结关系，可写出每个环节输入/受哪些环节输出M的影响，现列如入下：%=%w2=J-y93=%4=Hf%=Mh=0.01;L1=25;n=4;T0=0Tf=20;nout=4;YO=10;sp4_4;plot（t,y,r,）holdon运行结果如图中红色实线所示。（2

31、）考虑非线性环节N影响时，只需将非线性标志向量Z和参数向量S的相应分量正确输入即可。在MATLAB命令窗口中输入下列语句：Z=4000;S=|5000;%第一个线性环节后有饱和非线性，参数值为5。sp4_4;plot(t,y,-)运行结果如图中蓝色虚线所示。从图中可以清楚的地看出，饱和非线性环节对线性系统输出响应的影响。附：sp4_4函数为：A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1);W0=zeros(n,l);W=zeros(n,n);fork=l:mif(WIJ(k,2)=0);WO(WIJ(k,1)=WIJ(k,3);els

32、eW(WU(ktl),WIJ(k,2)=WU(k,3);end;end;fori=l:nif(A(i)=0);FI(i)=l;FIM(i)=h*C(i)B(i);FlJ=h*h*(CB(i)2;FIC(i)=l;HD(i)=0;if(D(i)-=0);FID(i)=D(i)B(i);elseendelseFl(i)=exp(-h*A(i)B(i);FIM(i)=(l-FI(i)*C(i)A(i);FIJ(i)=h*C(i)A(i)-FIM(i)*B(i)A(i);11C(i)=kFID(i)=O;if(D(i)-=O);FIC(i)=C(i)D(i)-A(i)B(i);FID(i)=D(i)B

33、(i);elseendendendY=zeros(n,l);X=Y;y=O;Uk=zeros(n,I);Ubb=Uk;t=TO:h*Ll:Tf;N=length(t);fork=l:N-lfori=l:LlUb=Uk;Uk=W*Y+W0*Y0;fori=l:nif(Z(i)-=O)if(Z()=l)Uk(i)=satu(Uk(i),S(i);endif(Z(i)=2)Uk(i)=dead(Uk(i),S(i);endif(Z(i)=3)Uk(i),Ubb(i)=backlash(Ubb(i),Uk(i),Ub(i),S(i);endendendUdot=(Uk-Ub)Zh;Uf=2*Uk-Ub;X=H.*X+FIM,.*Uk+FIJ,.*Udot;Yb=Y;Y=FIC.*X+FID.*Uf;fori=l:nif(Z(ib=O)if(Z(i)=4)Y(i)=satu(Y(i),S(i);endif(Z(i)=5)Y(i)=dead(Y(iS(i);endif(Z(i)=6)Y,Ubb(i)=backlash(Ubb,Y(i),Yb(i),S(i);endendendendy=y,Y(nout);end附：饱和非线性函数salu.m为：functionUc=satu(Ur,S1)if(abs(Ur)=Sl)if(Ur0)Uc=Sl;el