(完整)系統建模與仿真習題答案(forstudents)(完整)系統建模與仿真習題答案(forstudents)(完整)系統建模與仿真習題答案(forstudents)第一章習題1-1什么是仿真？它所遵循的基本原則是什么？答：仿真是建立在控制理論，相似理論，信息處理技術和計算技術等理論基礎之上的,以計算機和其他專用物理效應設備為工具，利用系統模型對真實或假想的系統進行試驗，并借助專家經驗知識，統計數據和信息資料對試驗結果進行分析和研究，進而做出決策的一門綜合性的試驗性科學。它所遵循的基本原則是相似原理。1-2在系統分析與設計中仿真法與解析法有何區別?各有什么特點?答：解析法就是運用已掌握的理論知識對控制系統進行理論上的分析，計算。它是一種純物理意義上的實驗分析方法,在對系統的認識過程中具有普遍意義。由于受到理論的不完善性以及對事物認識的不全面性等因素的影響，其應用往往有很大局限性.仿真法基于相似原理，是在模型上所進行的系統性能分析與研究的實驗方法.1-3數字仿真包括那幾個要素？其關系如何？答:通常情況下，數字仿真實驗包括三個基本要素，即實際系統，數學模型與計算機。由圖可見，將實際系統抽象為數學模型，稱之為一次模型化，它還涉及到系統辨識技術問題,統稱為建模問題；將數學模型轉化為可在計算機上運行的仿真模型，稱之為二次模型化，這涉及到仿真技術問題，統稱為仿真實驗.1—4為什么說模擬仿真較數字仿真精度低？其優點如何?.答：由于受到電路元件精度的制約和容易受到外界的干擾,模擬仿真較數字仿真精度低但模擬仿真具有如下優點：描述連續的物理系統的動態過程比較自然和逼真。仿真速度極快，失真小，結果可信度高。能快速求解微分方程.模擬計算機運行時各運算器是并行工作的，模擬機的解題速度與原系統的復雜程度無關.可以靈活設置仿真試驗的時間標尺，既可以進行實時仿真，也可以進行非實時仿真.易于和實物相連。1-5什么是CAD技術？控制系統CAD可解決那些問題？答：CAD技術，即計算機輔助設計（ComputerAidedDesign)，是將計算機高速而精確的計算能力,大容量存儲和處理數據的能力與設計者的綜合分析，邏輯判斷以及創造性思維結合起來，用以加快設計進程，縮短設計周期，提高設計質量的技術.控制系統CAD可以解決以頻域法為主要內容的經典控制理論和以時域法為主要內容的現代控制理論。此外，自適應控制，自校正控制以及最優控制等現代控制測策略都可利用CAD技術實現有效的分析與設計.1—6什么是虛擬現實技術？它與仿真技術的關系如何?答：虛擬現實技術是一種綜合了計算機圖形技術,多媒體技術，傳感器技術，顯示技術以及仿真技術等多種學科而發展起來的高新技術。虛擬現實技術不斷完善，為控制系統數字仿真與CAD開辟了一個新時代。1—7什么是離散系統？什么是離散事件系統？如何用數學的方法描述它們？答：本書所講的“離散系統”指的是離散時間系統，即系統中狀態變量的變化僅發生在一組離散時刻上的系統。它一般采用差分方程，離散狀態方程和脈沖傳遞函數來描述。離散事件系統是系統中狀態變量的改變是由離散時刻上所發生的事件所驅動的系統。這種系統的輸入輸出是隨機發生的，一般采用概率模型來描述.1-8如圖1-16所示某衛星姿態控制仿真實驗系統，試說明：若按模型分類,該系統屬于那一類仿真系統？圖中“混合計算機”部分在系統中起什么作用？與數字仿真相比該系統有什么優缺點？答:(1)按模型分類，該系統屬于物理仿真系統.(2）混合計算機集中了模擬仿真和數字仿真的優點，它既可以與實物連接進行實時仿真，計算一些復雜函數，又可以對控制系統進行反復迭代計算。其數字部分用來模擬系統中的控制器,而模擬部分用于模擬控制對象.與數字仿真相比,物理仿真總是有實物介入，效果逼真,精度高，具有實時性與在線性的特點,但其構成復雜，造價較高，耗時過長，通用性不強。第二章習題2—1思考題：(1)數學模型的微分方程，狀態方程,傳遞函數,零極點增益和部分分式五種形式，各有什么特點？答：微分方程是直接描述系統輸入和輸出量之間的制約關系,是連續控制系統其他數學模型表達式的基礎.狀態方程能夠反映系統內部各狀態之間的相互關系，適用于多輸入多輸出系統.傳遞函數是零極點形式和部分分式形式的基礎。零極點增益形式可用于分析系統的穩定性和快速性。利用部分分式形式可直接分析系統的動態過程。(2)數學模型各種形式之間為什么要互相轉換？答:不同的控制系統的分析和設計方法，只適用于特定的數學模型形式。（3）控制系統建模的基本方法有哪些？他們的區別和特點是什么？答：控制系統的建模方法大體有三種：機理建模法,實驗建模法和綜合建模法。機理建模法就是對已知結構,參數的物理系統運用相應的物理定律或定理，經過合理的分析簡化建立起來的各物理量間的關系。該方法需要對系統的內部結構和特性完全的了解，精度高。實驗建模法是采用歸納的方法,根據系統實測的數據，運用統計規律和系統辨識等理論建立的系統模型。該方法建立的數學模型受數據量不充分，數據精度不一致，數據處理方法的不完善,很難在精度上達到更高的要求。綜合建模法是上述兩種方法的結合.（4）控制系統計算機仿真中的“實現問題”是什么含意？答：“實現問題”就是根據建立的數學模型和精度，采用某種數值計算方法，將模型方程轉換為適合在計算機上運行的公式和方程，通過計算來使之正確的反映系統各變量動態性能，得到可靠的仿真結果。（5）數值積分法的選用應遵循哪幾條原則？答：數值積分法應該遵循的原則是在滿足系統精度的前提下，提高數值運算的速度和并保證計算結果的穩定.2-2。用matlab語言求下列系統的狀態方程、傳遞函數、零極點增益、和部分分式形式的模型參數，并分別寫出其相應的數學模型表達式:（1)G(s）= （2）MACROBUTTONMTEditEquationSection2SEQMTEqn\r\h\＊MERGEFORMATSEQMTSec\h=y=[0202］X解：(1)狀態方程模型參數：編寫matlab程序如下〉〉num=［172424］；〉〉den=[110355024］;〉>[ABCD]=tf2ss(num,den）得到結果:A=，B=,C=,D=[0]所以模型為：=X+u,y=X（2）零極點增益：編寫程序>〉num=[172424]；〉〉den=［110355024]；>〉［ZPK]=tf2zp(num，den）得到結果Z=—2。7306+2.8531i，-2。7306—2。8531i，—1.5388P=-4，-3,—2，—1K=1（3)部分分式形式:編寫程序〉>num=[172424]；>〉den=[110355024]；〉〉［RPH］=residue(num,den)得到結果R=4.0000，-6。0000,2.0000，1。0000P=-4.0000,-3。0000,-2。0000，-1.0000H=［]G（s）=（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=［4220］'；〉>C=[0202］；>〉D=［0];>>［numden］=ss2tf(A,B,C，D)得到結果num=04.000014.000022.000015。0000den=1。00004.00006。25005。25002。2500(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=[0202］;〉>D=[0];>〉［Z,P,K］=ss2zp（A,B，C,D）得到結果Z=-1。0000+1.2247i—1。0000—1.2247i-1.5000P=—0.5000+0.8660i—0。5000-0。8660i-1。5000—1。5000K=4。0000表達式(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.00000.0000—2.3094i0.0000+2.3094iP=—1。5000—1.5000—0.5000+0.8660i-0。5000-0.8660iH=[]2—3.用歐拉法求下面系統的輸出響應y（t)在0≤t≤1上，h=0.1時的數值.要求保留4位小數，并將結果與真解比較。解：歐拉法(前向歐拉法，可以自啟動）其幾何意義：把f（t，y)在[]區間內的曲邊面積用矩形面積近似代替。利用matlab提供的m文件編程，得到算法公式。如下所示(1）m文件程序為h=0。1;disp(’函數的數值解為’）；%顯示‘'中間的文字％disp（’y=’)；％同上%y=1;fort=0：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文件求解在t［0，1］的數值（％是的真解%）程序為h=0.1;disp('函數的離散時刻解為’);disp（'y='）；fort=0：h：1y=exp（-t）；disp(y）;end保存文件q3.m在matalb命令行中鍵入>〉q3函數的離散時刻解為y=10。90480.81870.74080.67030。60650。54880。49660。44930。40660。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的數值解，并于歐拉法求得的結果比較。解：我們經常用到預報—校正法的二階龍-格庫塔法，此方法可以自啟動，具有二階計算精度，幾何意義：把f（t,y)在［］區間內的曲邊面積用上下底為和、高為h的梯形面積近似代替。利用matlab提供的m文件編程，得到算法公式。如下所示（1）m文件程序為h=0.1;disp（'函數的數值解為')；disp(’y='）;y=1;fort=0：h:1disp（y）；k1=—y；k2=-(y+k1*h);y=y+（k1+k2)＊h/2；end保存文件q4.m在matlab的命令行中鍵入〉〉q4顯示結果為函數的數值解為y=10。90500。81900.74120.67080。60710.54940.49720。45000.40720。3685比較歐拉法與二階龍格—庫塔法求解.（誤差為絕對值）真值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步距取得較小時，誤差是很小的.編輯m文件程序h=0.1；disp('四階龍格-庫塔方法求解函數數值解為’）;disp（'y=’）;y=1;fort=0:h:1disp(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=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。40660.3679誤差00000000000顯然四階龍格—庫塔法求解精度很高，基本接近真值。三種方法比較可以得到精度（四階）〉精度(二階）>精度（歐拉）2—6．已知二階系統狀態方程為寫出取計算步長為h時，該系統狀態變量X=[]的四階龍格-庫塔法遞推關系式。解：四階龍格—庫塔法表達式所以狀態變量的遞推公式可以寫作：A=,B=，可以寫成則遞推形式2—7單位反饋系統的開環傳遞函數已知如下 用matlab語句、函數求取系統閉環零極點，并求取系統閉環狀態方程的可控標準型實現。解:已知開環傳遞函數，求得閉環傳遞函數為在matlab命令行里鍵入>〉a=[10］；>〉b=[14.6］;〉〉c=[13。416.35］；〉>d=conv（a，b)；>〉e=conv（d,c)e=1.00008.000031.990075.21000〉〉f=[0005100]；〉〉g=e+fg=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%計算零點％綜上：當閉環傳函形如時,可控標準型為：；所以可控標準型是2-8用matlab語言編制單變量系統三階龍格—庫塔法求解程序，程序入口要求能接收狀態方程各系數陣（A，B，C，D）,和輸入階躍函數r（t)=R＊1（t）；程序出口應給出輸出量y（t）的動態響應數值解序列。解：m文件為:functiony=hs(A,B，C,D,R，T，h）％T為觀測時間，h為計算步長,R為輸入信號幅值%disp(’數值解為’）；y=0；r=R；x=[0;0；0;0]；N=T/h；fort=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=,B=,C=,D=［0]在命令行里鍵入〉>A=[010000100001—100—80。21-31.99—8];〉〉B=［0001］’；>〉C=［-100500];>>D=[0];〉>T=1;>〉R=1；>〉h=0。01;>〉y=hs(A,B,C，D，R,T，h)數值解為08。3333e-0075。8659e-0061。8115e—0053。9384e-0057.0346e-005。。。。%僅取一部分%2-10。用式（2—34）梯形法求解試驗方程，分析對計算步長h有何限制，說明h對數值穩定性的影響.解：編寫梯形法程序為得到穩定系統最終漸進收斂。系統穩定則計算得。h的選取不能超出上述范圍，否則系統不穩定。2-11如圖2-27所示斜梁滾球系統,若要研究滾球在梁上的位置可控性,需首先建立其數學模型，已知力矩電機的輸出轉矩M與其電流i成正比，橫梁為均勻可自平衡梁（即當電機不通電且無滾球時，橫梁可處于=0的水平狀態），是建立系統的數學模型,并給出簡化后系統的動態結構圖。解：設球的質心到桿的距離為0，該系統為特殊情況下的球棒系統。另令分別表示棒的慣量、球的質量和球的慣量。則球質心的位置和速度為其中,。因而動能的移動部分為因而動能的移動部分為球棒系統的旋轉動能為因而，系統總的動能等于其中為常數。此系統的拉格朗日方程組為綜合以上公式的系統的方程組為設系統在平衡點附近，，，則系統方程可化為對上式進行拉普拉斯變換并化簡后可得到。參考文獻：［1]Hauser,S。Sestry，andP。Kokotovic.“Nonlinearcontrolviaapproximateinput-outputlinearization"。IEEETrans。onAutomaticControl，vol。37：pp。392-398，1992。［2］R.Sepulchre.“Slowpeakingandlow-gaindesignsforglobalstabilizationofnonlinearsystems”.submittedforIEEETAC1999.［3]R.Sepulchre,M。Jankovic，andP。KokotovicConstructiveNonlinearControl。Springer-Verlag，1997。[4]R.Teel。“UsingSaturationtostabilizeaclassofsingle-inputpartiallylinearcompositesystems”.IFACNOLCOS'92Symposium,pages369-374，June1992.2-12如圖2-28所示雙水箱系統中,為流入水箱1的液體流量，為流出水箱2的液體流量,試依據液容與液阻的概念,建立的系統動態結構圖.解：根據液容和液阻的概念，可分別列出兩個水箱的數學模型對上式進行在零初始條件下進行拉普拉斯變換得化簡后可得第三章習題3—2設典型閉環結構控制系統如圖4-47所示，當階躍輸入幅值時,用sp4_1.m求取輸出的響應。解:用sp4_1。m求解過程如下：在MATLAB語言環境下，輸入以下命令語句>〉a=［0.0160.8643。273。421];〉〉b=[3025]；>>X0=[0000]；％系統狀態向量初值為零>>V=2；%反饋系數>〉n=4；>〉T0=0；Tf=10；>〉h=0.01；R=20;%仿真步長h=0.01，階躍輸入幅值〉>sp4_1%調用sp4_1.m函數>>plot(t，y)運行結果為:附：sp4_1.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(1，n—1），1]’；m1=length（b)；C=[fliplr(b),zeros(1，n—m1）］；Ab=A—B＊C＊V;X=X0'；y=0；t=T0；N=round（(Tf-T0)/h)；fori=1:NK1=Ab*X+B*R;K2=Ab*(X+h＊K1/2)+B*R;K3=Ab*（X+h*K2/2)+B＊R；K4=Ab*(X+h*K3)+B*R；X=X+h*(K1+2*K2+2*K3+K4)/6;y=[y，C＊X］；t=［t,t（i）+h］;end3—4系統結構圖如圖3-55，寫出該系統的聯結矩陣和，并寫出聯結矩陣非零元素陣。解：根據圖3-55中,拓撲連結關系，可寫出每個環節輸入受哪些環節輸出的影響，現列如入下：把環節之間的關系和環節與參考輸入的關系分別用矩陣表示出來,即=，=，3-6若系統為圖4-5b雙輸入-雙輸出結構，試寫出該系統的聯接矩陣，，說明應注意什么?解：根據圖4—5b中,拓撲連結關系，可列寫如下關系式:轉換成矩陣形式為所以聯接矩陣=，=此時應注意輸入聯接矩陣變為型。3—-8求圖3-56非線性系統的輸出響應y(t），并與無非線性環節情況進行比較。解:（1）不考慮非線性環節影響時，求解過程如下：先將環節編號標入圖中.2）在MATLAB命令窗口下，按編號依次將環節參數輸入P陣；>〉P=［0。110。51；01200；2110;10110］；3）按各環節相對位置和聯接關系，有聯接矩陣如下：，，所以非零元素矩陣>>WIJ=[101；14-1;211;321；431］；4）由于不考慮非線性影響,則非線性標志向量和參數向量均應賦零值；〉〉Z=［0000];S=［0000］;5）輸入運行參數：開環截至頻率約為1，故計算步長h取經驗公式值，即，取h=0。01；每0。25秒輸出一點。故取=25。〉>h=0.01;〉>L1=25;〉〉n=4；>>T0=0>>Tf=20;〉>nout=4；〉>Y0=10;〉>sp4_4;>〉plot（t，y，'r')〉>holdon運行結果如圖中紅色實線所示。（2）考慮非線性環節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，1);W=zeros（n,n）;fork=1:mif（WIJ(k,2）==0)；W0（WIJ(k，1))=WIJ(k，3）;elseW(WIJ(k,1）,WIJ(k,2）)=WIJ（k,3）;end；end;fori=1:nif（A（i）==0）；FI（i)=1；FIM（i)=h*C(i）/B（i）;FIJ(i）=h*h＊（C（i）/B（i））/2;FIC(i）=1;FID(i）=0；if（D(i)～=0)；FID（i)=D（i)/B（i）;elseendelseFI（i)=exp(—h＊A(i）/B（i）)；FIM（i）=（1-FI（i)）*C(i)/A(i）;FIJ（i）=h＊C（i）/A(i）-FIM（i）＊B(i）/A（i);FIC(i）=1；FID（i)=0;if（D(i）~=0）；FIC（i）=C(i)/D（i)—A(i）/B(i）;FID(i)=D（i)/B(i）；elseendendendY=zeros(n,1);X=Y;y=0；Uk=zeros（n,1)；Ubb=Uk;t=T0:h＊L1:Tf；N=length(t);fork=1:N-1fori=1：L1Ub=Uk;Uk=W＊Y+W0＊Y0;fori=1:nif(Z（i)~=0）if(Z（i）==1)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）/h;Uf=2*Uk—Ub；X=FI'。＊X+FIM'