《計算機仿真技術與CAD》習題答案

1、第 0 章 緒論0-1 什么是仿真？它所遵循的基本原則是什么？答：仿真是建立在控制理論、相似理論、信息處理技術和計算機技術等理論基礎之上的，以計 算機和其他專用物理效應設備為工具，利用系統(tǒng)模型對真實或假想的系統(tǒng)進行試驗， 并借助專 家經(jīng)驗知識、 統(tǒng)計數(shù)據(jù)和信息資料對試驗結果進行分析和研究，進而做出決策的一門綜合性的 試驗性科學。它所遵循的基本原則是相似原理。0-2 仿真的分類有幾種？為什么？答：依據(jù)相似原理來分：物理仿真、數(shù)學仿真和混合仿真。物理仿真： 就是應用幾何相似原理， 制作一個與實際系統(tǒng)相似但幾何尺寸較小或較大的物 理模型（例如飛機模型放在氣流場相似的風洞中）進行實驗研究。數(shù)學仿真：就

2、是應用數(shù)學相似原理，構成數(shù)學模型在計算機上進行研究。它由軟硬件仿真 環(huán)境、動畫、圖形顯示、輸出打印設備等組成?；旌戏抡嬗址Q數(shù)學物理仿真，它是為了提高仿真的可信度或者針對一些難以建模的實體， 在系統(tǒng)研究中往往把數(shù)學仿真、 物理仿真和實體結合起來組成一個復雜的仿真系統(tǒng)， 這種在仿 真環(huán)節(jié)中有部分實物介入的混合仿真也稱為半實物仿真或者半物理仿真。0-3 比較物理仿真和數(shù)學仿真的優(yōu)缺點。答：在仿真研究中，數(shù)學仿真只要有一臺數(shù)學仿真設備（如計算機等） ，就可以對不同的控制 系統(tǒng)進行仿真實驗和研究，而且，進行一次仿真實驗研究的準備工作也比較簡單，主要是受控 系統(tǒng)的建模、控制方式的確立和計算機編程。數(shù)學仿真

3、實驗所需的時間比物理仿真大大縮短， 實驗數(shù)據(jù)的處理也比物理仿真簡單的多。與數(shù)學仿真相比，物理仿真總是有實物介入，效果直觀逼真，精度高，可信度高，具有實 時性與在線性的特點；但其需要進行大量的設備制造、安裝、接線及調(diào)試工作，結構復雜，造 價較高，耗時過長，靈活性差，改變參數(shù)困難，模型難以重用，通用性不強。 0-4 簡述計算機仿真的過程。答：第一步：根據(jù)仿真目的確定仿真方案 根據(jù)仿真目的確定相應的仿真結構和方法，規(guī)定仿真的邊界條件與約束條件。 第二步：建立系統(tǒng)的數(shù)學模型 對于簡單的系統(tǒng)，可以通過某些基本定律來建立數(shù)學模型。而對于復雜的系統(tǒng)，則必須利 用實驗方法通過系統(tǒng)辯識技術來建立數(shù)學模型。數(shù)學模

4、型是系統(tǒng)仿真的依據(jù)，所以，數(shù)學模型 的準確性是十分重要。第三步：建立仿真模型即通過一定算法對原系統(tǒng)的數(shù)學模型進行離散化處理， 就連續(xù)系統(tǒng)言， 就是建立相應的差分方程。第四步：編制仿真程序?qū)τ诜菍崟r仿真，可用一般高級語言或仿真語言。對于快速的實時仿真，往往需要用匯編 語言。第五步：進行仿真實驗并輸出仿真結果通過實驗對仿真系統(tǒng)模型及程序進行校驗和修改，然后按系統(tǒng)仿真的要求輸出仿真結果。0-5 什么是 CAD 技術？控制系統(tǒng) CAD 可解決哪些問題？答：CAD 技術，即計算機輔助設計( Computer Aided Design )，是將計算機高速而精確的計算 能力、大容量存儲和數(shù)據(jù)處理能力與設計者

5、的綜合分析、邏輯判斷以及創(chuàng)造性思維結合起來， 以加快設計進程、縮短設計周期、提高設計質(zhì)量的技術。控制系統(tǒng) CAD 可以解決以頻域法為主要內(nèi)容的經(jīng)典控制理論和以時域法為主要內(nèi)容的現(xiàn) 代控制理論。此外，自適應控制、自校正控制以及最優(yōu)控制等現(xiàn)代控制策略都可利用 CAD 技 術實現(xiàn)有效的分析和設計。第 1 章 仿真軟件MATLAB1- 1 對于矩陣 A=1 2;3 4 ， MATLAB 以下四條命令：A.a(0.5) ; AA(0.5) ; sqrt(A) ; sqrtm(A) 所得結果相同嗎？它們中哪個結果是復數(shù)矩陣，為什么？ 答：A.A(0.5)=1.00001.4142;1.73212.0000

6、 ；AA(0.5)= 0.5537 + 0.4644i 0.8070 - 0.2124i;1.2104 - 0.3186i1.7641 + 0.1458i ；sqrt(A)= 1.00001.4142;1.73212.0000 ；sqrtm(A) = 0.5537 + 0.4644i 0.8070 - 0.2124i;1.2104 - 0.3186i1.7641 + 0.1458i ；其中，“ A.a(0.5) ”表示向量的乘方，“ Aa(0.5) ”表示矩陣的乘方，“ sqrt(A) ”只定義在矩陣的單個元素上，即分別對矩陣的每個元素進行運算，“ sqrtm(A) ”表示對矩陣(方陣)的超越

7、函數(shù)進行運算。1- 4 求二元函數(shù)方程組： sin(x-y)=0,cos(x+y)=0 的解。 答：x,y=solve(sin(x-y)=0,cos(x+y)=0,x,y)x =-1/4*pi1/4*piy =-1/4*pi1/4*pi1- 5 求函數(shù) y(t)=exp(-t)*|sincost| 的最大值( 0=tf=(-1)*exp(-(abs(x)*abs(sin(cos(abs(x);x=fminsearch(f,0),ymax=exp(-(abs(x)*abs(sin(cos(abs(x)x =0ymax =0.84151-6 設 D2y-3Dy+2y=x,y(0)=1,Dy(0)=

8、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.61001-7 求方程 cos(t)A2*exp(-0.1t)=0.5t 的解。 答：t1=solve(cos(t)A2*exp(-0.1*t)=0.5*t,t);t=eval(t1)0.83291-8 求方程組：xA2+yA2=1,xy=2的解。答：x,y=solve(xA2+yA2=1,x*y=2,x,y)x =-1/2*(1/2*5A(1/2)+1/2*i*3A(1/2)A3+1/4*5A(1/2)+1/4*i*3A(1

9、/2) -1/2*(1/2*5A(1/2)-1/2*i*3A(1/2)A3+1/4*5A(1/2)-1/4*i*3A(1/2) -1/2*(-1/2*5A(1/2)+1/2*i*3A(1/2)A3-1/4*5A(1/2)+1/4*i*3A(1/2) -1/2*(-1/2*5A(1/2)-1/2*i*3A(1/2)A3-1/4*5A(1/2)-1/4*i*3A(1/2)y =1/2*5A(1/2)+1/2*i*3A(1/2)1/2*5A(1/2)-1/2*i*3A(1/2)-1/2*5A(1/2)+1/2*i*3A(1/2)-1/2*5A(1/2)-1/2*i*3A(1/2)1-9 求 f(kT

10、)=kexp(-akT) 的 Z 變換表達式。答：syms k t z;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)/a1-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(

11、0)=0,x2(0)=1,t)x1 =exp(3*t)*sin(4*t)x2 =exp(3*t)*cos(4*t)第2章控制系統(tǒng)的數(shù)學模型及其轉(zhuǎn)換2- 1已知系統(tǒng)的傳遞函數(shù)為試用MATLAB 建立其狀態(tài)空間表達式。答：n um=1 1 1;de n=1 6 11 6;A,B,C,D=tf2ss( num,de n) A =-6-11-6100010B =100C =1 1 1D =02- 2已知系統(tǒng)的狀態(tài)空間表達式為 試用MATLAB求其傳遞函數(shù)陣。答： A=0 1;-2 -3;B=1 0;1 1;C=1 0;1 1;D=zeros(2,2); nu m1,de n1=ss2tf(A,B,C,

12、D,1), nu m2,de n2=ss2tf(A,B,C,D,2)num1 =001.00002.00004.00002.0000den1 =132num2 =00.00001.000001.00001.0000den2 =1322- 3已知兩子系統(tǒng)的傳遞函數(shù)分別為1G1 s，G2 ss 1 s 2試利用MATLAB 求兩子系統(tǒng)串聯(lián)和并聯(lián)時系統(tǒng)的傳遞函數(shù)。答： num1=1;den1=1 3 2;num2=1;den2=1 3 0;num,den=series(num1,den1,num2,den2)num =00001den =161160x P 1x ，則利用 MA TLAB 求在新狀態(tài)

13、變量下，系統(tǒng)狀態(tài)空間表達式。 num1=1;den1=1 3 2;num2=1;den2=1 3 0; num,den=parallel(num1,den1,num2,den2)00 262den =16 11602-4 設系統(tǒng)的狀態(tài)空間表達式為若取線性變換陣P1 11 -1num =設新的狀態(tài)變量為答：A=0 1;-2 -3;B=1;2;C=3 0;D=0;P=1 1;1 -1; A1,B1,C1,D1=ss2ss(A,B,C,D,P)A1 =-2 03 -1B1 =3-1C1 =1.5000 1.5000D1 =02-5 已知離散系統(tǒng)狀態(tài)空間表達式試用 MATLAB 求其系統(tǒng)的脈沖傳遞函數(shù)

14、。 答：A=0 1;1 3;B=0;1;C=1 1;D=0;T=1;A1,B1,C1,D1=c2dm(A,B,C,D,T) A1 =2.9598 7.33577.3357 24.9669B1 =1.95987.3357C1 =1 1D1 =0第 3 章 連續(xù)系統(tǒng)的數(shù)字仿真3- 1 已知線性定常系統(tǒng)的狀態(tài)空間表達式為 且初始狀態(tài)為零，試利用四階- 龍格庫塔法求系統(tǒng)的單位階躍響應。答：%ex3_1.mr=1; A=0 1;-5 -6; B=2;0; C=1 2; d=0;Tf=5; h=0.1;x=zeros(length(A),1); y=0; t=0;for i=1:Tf/hK1=A* x+B

15、*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=y;C*x; t=t;t(i)+h;endplot(t,y)3- 2 設單位反饋系統(tǒng)的開環(huán)傳遞函數(shù)試利用二階 - 龍格庫塔法求系統(tǒng)的單位階躍響應。答：%ex3_2.mr=1; numo=4; deno=1,2,0; num,den=cloop(numo,deno); A,b,C,d=tf2ss(num,de n);Tf=5; h=0.1;x=zeros(length(A),1); y=0; t=0;for i=1:Tf

16、/hK1=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;endplot(t,y)3-4 利用 input( ) 函數(shù)修改例3-1 所給程序 ex3_1.m ，將其中給定的參數(shù)r，numo ，deno ，numh 和 denh 利用鍵盤輸入， 序。使其變?yōu)檫B續(xù)控制系統(tǒng)面向傳遞函數(shù)的通用數(shù)字仿真程答：3- 5 利用 input( ) 函數(shù)修改例 3-2 所給程序 ex3_2.m ，將其中給定的參數(shù)r，P，W，W0和 Wc 利用鍵盤輸入，使其變?yōu)檫B續(xù)控制系統(tǒng)面向結構圖的通用數(shù)字仿真程序。答:第 4 章 連續(xù)系統(tǒng)按環(huán)節(jié)離

17、散化的數(shù)字仿真4- 1 已知非線性習題如圖題 4-1 所示，試利用連續(xù)系統(tǒng)按環(huán)節(jié)離散化的數(shù)字仿真方法， 求輸出量 y 的動態(tài)響應，并與無非線性環(huán)節(jié)進行比較。 (圖略) 答：%ex4_1.m % 主程序R=10;P=0.1 1 0.5 1 5 5; 0 1 1 0 0 0; 2 1 2 0 0 0;10 1 10 0 0 0;W=0 0 0 -1; 1 0 0 0; 0 1 0 0; 0 0 1 0;W0=1;0;0;0;Wc=0 0 0 1;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

18、(A);for i=1:nif (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);elseE(i)=exp(-A(i)*T/B(i);F(i)=(D(i)/B(i)-C(i)/A(i)*(1-E(i)*B(i)/(A(i)*T)-1); G(i)=(D(i)/B(i)-C(i)/A(i)*(1+(E(i)-1)*(1+B(i)/(A(i)*T);H(i)=1;L(i)=D(i)/B(i);Q(i)=0;endelseif (B(i)=0)E(i)=1;F(i)=0

19、.5*C(i)*T/B(i);G(i)=F(i);H(i)=1;L(i)=D(i)/B(i);Q(i)=0;elsedisp( A(i)=B(i)=0 );endendend x=zeros(length(A),1);x0=x;z=x; u=zeros(length(A),1);u0=u;y=zeros(length(Wc(:,1),1);t=0;for j=1:Tf/Tu1=u; u=W*x+W0*R;for i=1:nif (FZ(i)=0)if (FZ(i)=1) u(i)=saturation(u(i),S(i); endif (FZ(i)=2) u(i)=deadzone(u(i),

20、S(i); endif (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); endendendx1=x;for i=1:nz(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);endfor i=1:nif (FZ(i)=0)if (FZ(i)=5) x(i)=saturation(x(i),S(i); endif (FZ(i)=6) x(i)=deadzone(x(i),S(i

21、); endif (FZ(i)=7) x(i),x0(i)=backlash(x0(i),x(i),x1(i),S(i); endif (FZ(i)=8) x(i)=sign1(x(i),S(i); endendendy=y,Wc*x;t=t,t(j)+T;endplot(t,y)%saturation.m % 子程序function x=saturation(u,s)if (abs(u)=s)if (u0) x= s;else x=-s;endelsex= u;end修改“ P=0.1 1 0.5 1 0 0; 0 1 1 0 0 0; 2 1 2 0 0 0;10 1 10 0 0 0;”

22、ex4_14- 2 針對例 3-2 所給線性定常系統(tǒng)， 試利用第 4 章所給程序， 求系統(tǒng)的單位階躍響應， 并對其結果進行比較。答：ex3_2ex4_14- 3 針對例 4-1 所給系統(tǒng)， 去掉飽和非線性環(huán)節(jié)后求系統(tǒng)的單位階躍響應， 并與例 4-1 所得結果進行比較。答：ex4_1修改“ P=0.1 1 0.5 1 0 0; 0 1 1 0 0 0; 2 1 2 0 0 0;10 1 10 0 0 0;”ex4_14-4 利用 input( )函數(shù)修改例 4-1 所給程序 ex4_1.m ，將其中給定的參數(shù)R， P， W，W0 和 Wc 利用鍵盤輸入， 使其變?yōu)檫B續(xù)控制系統(tǒng)按環(huán)節(jié)離散化的通用數(shù)

23、字仿真程序。 答：略第 5 章 采樣控制系統(tǒng)的數(shù)字仿真5- 1 已知采樣控制系統(tǒng)的結構圖如圖題 5-1 所示(圖略) 。試利用采樣控制系統(tǒng)的數(shù)字 仿真方法，求當采樣周期 T=0.1s ，且初始狀態(tài)為零時，離散系統(tǒng)的單位階躍響應。 答：%ex5_1.mR=1;Gr=1;Fr=0;P=1 1 1 0 0 0;1 2 1 0 0 0;W=0 0;1 0;W0=1;0; Wc=0 1;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=lengt

24、h(Gr);for i=1:nif (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);elseE(i)=exp(-A(i)*T/ B(i);F(i)=(D(i)/B(i)- C(i)/ A(i)*(1- E(i)* B(i)/( A(i)*T)-1);G(i)=(D(i)/B(i)- C(i)/ A(i)*(1+( E(i)-1)*(1+ B(i)/( A(i)*T);H(i)=1; L(i)=D(i)/ B(i); Q(i)=0;endelseif (B

25、(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;elsedisp( A(i)= B(i)=0 );endendendx=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);for ij=0:Tf/Tm;e=R-x(n);Er=e;Er(1:m1-1);ur=-Fr*Ur+ Gr*Er;Ur= ur;Ur(1:n1-1);for j=1:Tm/Tu1= u; u = W*x+W0*ur;for i=1:nif (FZ(i)=0)if (FZ(i)=1) u(i)=saturation(u(i), S(i);endif (FZ(i)=2) u(i)=deadzone(u(i), S(i);endif (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);endendendx1= x;for i=1:nz(i)=E(i)*z(i)+F(i)* u(i)+G(i)*u