連續(xù)系統(tǒng)仿真的基礎算法

1、2.1 2.1 數值積分算法數值積分算法 2.2 2.2 數值積分算法的基本分析數值積分算法的基本分析 2.3 2.3 連續(xù)系統(tǒng)仿真的離散相似算法連續(xù)系統(tǒng)仿真的離散相似算法2.4 2.4 常用快速數字仿真算法常用快速數字仿真算法2.5 2.5 實時數字仿真算法實時數字仿真算法 小結小結3.1.1 數值積分算法的基本原理3.1.2 歐拉法3.1.3 龍格-庫塔法3.1.4 微分方程數值積分的矩陣分析 上一頁下一頁返回連續(xù)系統(tǒng)仿真的數值積分算法就是利用數值積分法將常微分方程（組）描述的連續(xù)系統(tǒng)變換成離散形式的仿真模型差分方程（組）。為了能在計算機上進行求解，首先要把被仿真系統(tǒng)的數學模型表示為一階微

2、分方程組或狀態(tài)空間模型。上一頁下一頁返回 假設一階向量微分方程及初值問題為00)(),(yyyfytt可以寫成一階微分方程組和初值問題形式00200210012121222111)()()(),(),(),(nnnnnnnytyytyytyyyytfyyyytfyyyytfy上一頁下一頁返回 (2.1)(2.2)110,kkttttt110d),()(d),()()(01kkkttkttktttttttyfyyfyy在上的解析解為希望用公式 kkkqyy11kykykq)(1kty)(kty1d),(kktttt yf來代替解析解，其中 ， ， 分別是 ， ， 的近似解。 上一頁下一頁返回 (

3、2.3)(2.4) 所謂數值解法就是尋求初值問題的解在一系列離散時間點 121,kktttt,21yy1,kkyy上的近似解。相鄰兩個離散時間點的間距 kktth1稱為步長或步距 。數值積分法的主要問題歸結為如何對向量函數 ),(yf t進行數值積分，求出 ),(yf t在區(qū)間 ),(1kktt上定積分的近似解 kq ，并且數值方法的共同特點是步進式的（即所謂遞推式的）。 上一頁下一頁返回 常用的基本算法可以分為單步法、多步法和預估-校正法三大類，而每一類算法又可以分為顯式算法和隱式算法兩種類型。 先從標量形式開始討論。考慮一階微分方程及初值問題 00)(),(ytyytfy 解析解為 11d

4、),()(d),()()(001kkkttkttktytftytytftyty希望用公式 kkkqyy1來代替解析解，其中 分別是 , , 的近似值。kkkqyy,1)(1kty)(kty1d),(kktttytf 上一頁下一頁返回 (2.5) (2.6) (2.7)00)(),(ytyytfy ),(1kktt 歐拉（Euler）法是一種最簡單的數值積分算法，對于 在區(qū)間 上求積分，有 1d),()()(1kkttkktytftyty若區(qū)間 ),(1kktt足夠小，則 ),(1kktt上的 ),( ytf可以近似地看成是常數),(kkytf上一頁下一頁返回 (2.8)11),()(kkkkk

5、yythfyty即有),(kkkythfq 歐拉法為 kkkhfyy1式中),(kkkytff 上一頁下一頁返回 (2.9) (2.10) 歐拉法的圖形表示： 為圖2.1中的非曲面面積actk+1tk，歐拉法用矩形面積abtk+1tk代替。1d),(kktttytf顯然，歐拉法僅適用于步長h很小的場合。 上一頁下一頁返回圖2.1龍格-庫塔（Runge-Kutta）法是求解常微分方程初值問題(2.5)式的各種數值積分算法中應用得最廣泛的一種，包括許多不同的公式。思路：用若干個時間點上 f 的函數值的線性組合來代替 f 的各階導數項，然后按泰勒公式展開確定其中的系數。 上一頁下一頁返回 RK法包含

6、有顯式、隱式和半隱式等算法。僅介紹顯式RK法。 一般形式： rickahyhctfkkWhyyijjijkikiriiikk, 2 , 10),(11111 式中， Wi為待定的權因子；r為使用的k值的個數；ki為不同時間點上導數 f 的值；ci，aij為待定系數。上一頁下一頁返回 (2.13)當r =1時，只有一個k1，就得到了歐拉法 kkkhfyy1當r =2時，),(),()(12122122111khayhctfkytfkkWkWhyykkkkkk將 kfk 1代入 2k中，并將 2k在 ),(kkyt 附近用泰勒級數展開,可得)()(),(),(),(221222122122hOff

7、 haf hcfhOyfytfatfchytffhayhctfkkykkkttyykkkkkkkkk上一頁下一頁返回 (2.15)(2.14)于是，有22111khWkhWyykk)()(3212222221hOffaWhfcWhfWWhykykkkk 將其與泰勒級數式逐項進行比較，令h，h2項的對應系數相等，得到以下關系 212112122221aWcWWW取c2=1，有 2121WW121a上一頁下一頁返回 (2.16)(2.17)從而，有RK2法 ),(),()(2121211hkyhtfkytfkkkhyykkkkkk當r=3時，可得RK3 法)32,32()3,3(),()3(423

8、121311hkyhtfkkhyhtfkytfkkkhyykkkkkkkk上一頁下一頁返回(2.18) (2.19) 當r=4時，可得到如下著名的四階RK法（亦稱為經典RK法，簡記為RK4法）。 ),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyykkkkkkkkkk 該算法是數字仿真中最常用的一種算法。計算量較大，但計算精度較高。在實際仿真中得到了廣泛的應用。 上一頁下一頁返回 (2.20)幾種數值積分法都可以看成是 在 附近用泰勒級數展開而產生的，只不過是泰勒級數所取項數的多少不同而已。歐拉法只取前兩項， RK2

9、法取了前三項，RK4法取了前五項。從理論上講，可以構造任意高階的計算方法。但是，精度的階數與計算函數f值的次數之間的關系并非等量增加的。 )(tyktt 上一頁下一頁返回每步計算f 的次數234567r8算法精度階數234456r-2由此可見，RK4法有其優(yōu)越性。 表2.1 f 的計算次數與算法精度階數的關系 上一頁下一頁返回 對于一階向量微分方程及初值問題 00)(),(yyyfyttRK4法的向量形式為),()2,2()2,2(),()22(6342312143211KyfKKyfKKyfKyfKKKKKyyhhthhthhtthkkkkkkkkkk上一頁下一頁返回 (2.21) (2.2

10、2)具體寫出來就是上一頁下一頁返回),()2,2,2,2()2,2,2,2(),()22(63,23, 213, 142,22, 212, 131,21, 211, 12, 2, 114321,1,nknkkkiinknkkkiinknkkkiiknkkkiiiiiikikihkyhkyhkyhtfkkhykhykhyhtfkkhykhykhyhtfkyyytfkkkkkhyy(2.23) 在控制系統(tǒng)的仿真中，最常見的向量微分方程是線性定常系統(tǒng)的狀態(tài)方程 uBAxx即有ubbbxxxaaaaaaaaaxxxnnnnnnnnn212121222211121121RK4法的4個iK可表示為)()2

11、(2)2(2)(3423121htuhhtuhhtuhtukkkkkkkkBKxAKBKxAKBKxAKBAxK上一頁下一頁返回 (2.24)(2.25) (2.26)3.2.1 單步法和多步法3.2.2 顯式算法和隱式算法3.2.3 截斷誤差和舍入誤差3.2.4 數值積分算法的計算穩(wěn)定性3.2.5 數值積分算法的計算精度、速度、穩(wěn)定性與步長的關系3.2.6 數值積分算法的選擇原則3.2.7 誤差估計與步長控制3.2.8 數值積分算法仿真實例上一頁下一頁返回上節(jié)中介紹的幾種數值積分算法都有一個共同特點：在計算yk+1時只用到了yk ,而不直接用yk-1, yk-2 ， 等項的值，即在本次計算中

12、僅僅用到了前一次的計算結果。這類算法稱為單步法。單步法運算有如下優(yōu)點：需要存儲的數據量少，占用的存儲空間少；只需要知道初值，就可以啟動遞推公式進行運算，即可以自啟動；容易實現變步長運算。 上一頁下一頁返回與單步法相對應的還有一類數值積分算法，稱為多步法。在它的計算公式中，本次計算不僅要用到前一次的計算結果，還要用到更前面的若干次結果。例如AB4法:)9375955(243211kkkkkkffffhyy式中 3 , 2 , 1 , 0),(iyihtffikkikl 多步法與單步法相比，欲達到相同的精度，計算工作量要少得多。因此，在相同條件下多步法比單步法要快。多步法的主要缺點是，它不是自啟動

13、方法，必須用其它方法求初始幾步的值。另外，多步法不容易實現變步長運算。 上一頁下一頁返回 (2.28)如果數值積分算法在計算yk+1時所用到數據均已求出，則稱其為顯式算法。 如果數值積分算法的計算公式右端中隱含有未知量 yk+1，則稱其為隱式算法。例如AM4法就是隱式算法 : )5199(242111kkkkkkffffhyy上一頁下一頁返回(2.29)隱式算法也不是自啟動的算法，需要用另一個顯式算法估計一個初值，然后再用隱式算法進行迭代運算，這就構成了預估-校正算法。 四階Adams預估-校正算法： )5199(24)9375955(2421)0(11321)0(1kkkkkkkkkkkkf

14、fffhyyffffhyy上一頁下一頁返回 (2.30) 在分析數值積分算法的精度時，通常用泰勒級數作為工具。假設前一步求得結果是準確的，即有 )(kktyy 則用泰勒級數求得在1kt處的解析解為 )()(!1)(! 21)()()(1)(2rkrrkkkkhOtyhrtyhtyhtyhty 不同的數值積分算法相當于在上式中取了不同項數之和而得到的近似解。 上一頁下一頁返回 (2.31)(2.32)這種由數值積分算法單獨一步引進的附加誤差稱為局部截斷誤差。它是數值積分算法給出的解與微分方程的解析解之間的差，故又稱為局部離散誤差。步長h越小，局部截斷誤差就越小。 若數值積分算法的局部截斷誤差為

15、，則方法為r階的。)(1rhO 算法的階數可以作為衡量算法精確度的一個重要標志。上一頁下一頁返回 以上的分析是在假設前一步所得結果是準確的前提下得出的，即 )(kktyy 成立時，有 )()(111rkkhOyty 但在求解過程中，實際上只有當k=0 時， 才成立。而當k=1,2,3,時， 并不成立。因而r階算法求得的解的實際誤差要大一些的。 )(kktyy 設yk是在無舍入誤差情況下由r階方法算出的微分方程式的近似解， y(t)為微分方程的解析解， 為算法的整體截斷誤差，則可以證明， 。即整體截斷誤差比局部截斷誤差低一階。 kkkyty)()(rkhO上一頁下一頁返回(2.33)1. 計算穩(wěn)

16、定性的概念 2. 歐拉法的計算穩(wěn)定性 3. 龍格-庫塔法的計算穩(wěn)定性 4. 多步法的計算穩(wěn)定性 上一頁下一頁返回連續(xù)系統(tǒng)的數字仿真，實質上就是將給定的微分方程（組）變換成差分方程（組），然后從初值開始，逐步進行遞推計算。 【例2.1】 已知微分方程及初值 5 . 1031)0(30tyyy 試比較在取不同步長時，其精確解與數值積分算法解之間的差異。 上一頁下一頁返回 【解】該初值問題的解析解為 tety3031)(上一頁下一頁返回對應的結果曲線如圖2.2（a）所示。圖2.2 (a) 解析解5 . 0 ,075. 0 , 1 . 0h5 . 1t)(ty取，分別用歐拉法和RK4法計算處的，所得結

17、果曲線如圖2.2 （b）（g）所示。 圖2.2 (b) 歐拉法（h=0.1） 圖2.2 (c) RK4法（h=0.1） 上一頁下一頁返回圖2.2 (d) 歐拉法（h=0.075） 圖2.2 (e) RK4法（h=0.075） 圖2.2 (f) 歐拉法（h=0.05） 圖2.2 (g) RK4法（h=0.05） 上一頁下一頁返回 從圖中可以看出，解析解單調下降并迅速收斂到0。 當 時，歐拉法和RK4法的解曲線均發(fā)散，數值積分算法的解是錯誤的。 1 . 0h 當 時，歐拉法的解曲線仍然發(fā)散，對應的解是錯誤的；RK4法的解曲線單調下降并收斂到0，對應的解是正確的。 075. 0h 當 時，歐拉法和R

18、K4法的解均收斂到0（雖然歐拉法的解曲線是振蕩收斂的）。如果只要求得到 處的 的解，則兩種數值積分算法的解都可以認為是正確的。事實上，此時有 05. 0h5 . 1t)(ty 上一頁下一頁返回歐拉法 RK4法 10101044085. 3)5 . 1 ( 4)5 . 1 (y解析解 21105417286. 9)5 . 1 (y為什么會出現上述情況呢？這是因為數值積分算法只是一種近似積分方法，在反復的遞推計算中會引進誤差（這種誤差通常是由初始數據的誤差及計算過程中的舍入誤差產生）。如果誤差的累積越來越大（見【例2.1】）,將使計算出現不穩(wěn)定，從而得到錯誤的結果。系統(tǒng)的

19、穩(wěn)定性與計算穩(wěn)定性是兩個不同的概念。由于選用的數值積分算法不同，即使對于同一系統(tǒng)，差分方程也各不相同，計算穩(wěn)定性也就各不一樣。上一頁下一頁返回通常用一個簡單的一階微分方程來考查數值積分算法的計算穩(wěn)定性。 微分方程及初值問題0)0(yyyy稱為測試方程（Test Equation），其中， 為方程的特征根。根據穩(wěn)定性理論，當特征根位于左半s平面，即 時，原系統(tǒng)穩(wěn)定。 j0Re上一頁下一頁返回 (2.34) 用歐拉法對測試方程進行計算，有 , 1 , 0)1 (1kyhhyyykkkk 假定僅僅在初始值y0處引進了初始誤差 ，而在遞推計算過程沒有引進任何其它誤差。顯然，此時 的誤差01ky 1k僅

20、由 ky的誤差 k引起的，所以有 )(1 ()(11kkkkkkkkyhyhyy用（2.36）式減去（2.35）式，得到誤差方程 （2.35）（2.36）上一頁下一頁返回kkh)1 (1以此類推，有 0121)1 ()1 ()1 (kkkkhhh 當 時， 。表明若在遞推計算過程中的某步引入了誤差，則隨著遞推步數的增加，這個誤差將逐漸擴大，最終導致差分方程的解完全失真。 11 hkk1 反之，當 時， 。表明隨著遞推步數的增加，引入的誤差會被逐漸減小或保持有界。在這種情況下，稱差分方程是計算穩(wěn)定的。 11 hkk1上一頁下一頁返回 (2.37) (2.38) 顯然，對于歐拉法，合理地選擇步長h

21、使其滿足 是保持其計算穩(wěn)定性的充要條件。 11 h 令 ，由 ，可知歐拉法在實軸上的穩(wěn)定區(qū)域為（-2，0），即有 hh11 h2hh因此， 越大，步長h就應該取得越小。 Re 為了保證計算穩(wěn)定性而對步長h有所限制的數值積分算法稱為條件穩(wěn)定算法。在任何步長h0的情況下計算都穩(wěn)定的數值積分算法稱為絕對穩(wěn)定算法（亦稱為無條件穩(wěn)定算法或恒穩(wěn)算法）。歐拉法是一種條件穩(wěn)定的計算方法。 上一頁下一頁返回 將RK法應用于測試方程，容易得到下列誤差方程 )()(!1)(! 211121rkrkhOhrhh令 ，代入上式，可得該式穩(wěn)定的條件為 hh1!1! 21121rhrhh可得如表2.2所示的各階RK法在實軸

22、上的穩(wěn)定區(qū)域。 上一頁下一頁返回 （2.40）(2.39)顯式RK法都是條件穩(wěn)定算法，步長h的大小除了與所選用的算法有關外，還與方程本身的性質有關。對于實際系統(tǒng)，由于其特征值不一定為實數，因此滿足（2.40）式的也是復數。一般而言，步長h必須滿足不等式 (2.41) 1hh12/12hh 6/2/132hhh24/6/2/1432hhhhr 所在的實負穩(wěn)定區(qū)域1（-2，0）2（-2，0）3（-2.51，0）4（-2.785，0） 表2.2 RK法在實軸上的穩(wěn)定區(qū)域 32h上一頁下一頁返回將多步法應用于測試方程，可以得到類似的誤差方程。結論：在同階的多步法中，隱式算法的穩(wěn)定區(qū)域比顯式算法的大得多

23、；隨著階數的增加，多步法的穩(wěn)定區(qū)域逐漸減小。 上一頁下一頁返回方法的階數1234顯式（-2，0）（-1，0）（-6/11，0）（-3/10，0）隱式（-，0）（-，0）（-6，0）（-3，0） 表2.3 Adams法在實軸上的穩(wěn)定區(qū)域 除了AM1法和AM2法（隱式Adams法）為絕對穩(wěn)定算法外，其它算法都是條件穩(wěn)定算法。這就是說，步長h必須滿足不等式 式中，為由積分算法確定的常數。上一頁下一頁返回|h|plot(tout,yout,k);grid; 上一頁下一頁返回圖2.8 (exam2_4.mdl) 得到如圖2.9 (a)所示的階躍響應曲線。顯然，該曲線不是很理想，超調量較大。為此，可以將外

24、環(huán)的PI控制器參數調整為 ，并分別選擇=0.17，0.5，1，1.5，可以得到如圖2.9 (b)所示的仿真結果?？梢钥闯觯绻x擇PI控制器為 ，就能夠得到較為滿意的控制效果。sas85. 0) 1(ss85. 015 . 1圖2.9 （a）圖2.9 （b）上一頁下一頁返回 【例2.5】考慮著名的Van de Pol方程 0) 1(2yyyy 25. 0)0(, 0)0(yy試繪制其相軌跡（ ）?！窘狻窟x擇狀態(tài)變量 yxyx21則有如下非線性狀態(tài)方程 1221221) 1(xxxxxx上一頁下一頁返回200 t 第1個方程可以看成是將x2(t)信號作為一個積分器的輸入端，積分器的輸出將成為x1

25、(t)信號。同樣，x2(t)信號本身也可以看成是一個積分器的輸出，而積分器的輸入端信號應該為 。于是，利用Simulink提供的各種模塊可以得到如圖2.10所示的仿真結構。 1212) 1(xxx上一頁下一頁返回 圖2.10(exam2_5.mdl) Simulink模型中有些模塊需要將輸入端和輸出端（通常用于反饋路徑）掉換一下方向。為此，可以用鼠標單擊需要掉換方向的模塊選中它，則選中的模塊的四個角出現黑點，表明它處于選中的狀態(tài)。然后， 打開Simulink的Format菜單(見圖2.11)，選擇其中的翻轉子菜單（Flip block）即可。 上一頁下一頁返回圖2.11 啟動仿真后，仿真結果將

26、賦給MATLAB工作空間內的變量t，x1和x2。在MATLAB命令窗口中運行指令 plot(x1,x2,k);grid; 得到相軌跡，如圖2.12所示。 上一頁下一頁返回圖2.12Simulink模型中的參數可以是實際數值，也可以是用字母表示的變量名。用字母表示的變量可以在MATLAB的工作空間中賦值，或用M文件賦值，然后進行仿真計算。 上一頁下一頁返回 【例2.6】 含有磁滯回環(huán)非線性環(huán)節(jié)的控制系統(tǒng)如圖2.13 (a)所示，其中，磁滯回環(huán)非線性環(huán)節(jié)的特性如圖2.13（b）所示。試研究該非線性環(huán)節(jié)對系統(tǒng)性能的影響。 上一頁下一頁返回圖2.13 【解】 構建系統(tǒng)的Simulink模型，如圖2.1

27、4所示。為了便于研究問題的方便起見，將Backlash模塊的參數deadband width設置為C1。 在MATLAB命令窗口依次運行C1的不同值(C1=0.5,1,2)的指令后啟動仿真，并在仿真結束后在MATLAB命令窗口中運行指令plot(tout,yout,k); 上一頁下一頁返回圖2.14 (exam2_6.mdl)得到如圖2.15所示的仿真曲線。顯然，不同磁滯寬度的階躍響應曲線的形狀不相同。需要注意的是，為了將不同C1值對應的階躍響應曲線繪制在同一坐標軸下，在第一次繪制圖形后，應該在MATLAB命令窗口中運行指令hold on 上一頁下一頁返回圖2.15 (r(t)=21(t) 此

28、外，如果輸入的幅值增大或減小，則原來系統(tǒng)響應曲線的形狀將可能出現不同。圖2.16為C1=1時，輸入幅值等于3和0.6的仿真結果。顯然，這一點與線性系統(tǒng)不相同。 上一頁下一頁返回圖2.16如果Simulink模型中存在復雜的非線性環(huán)節(jié)或復雜的邏輯運算，而在MATLAB提供的所有工具箱中都找不到相應的模塊，可以自己編制一個M函數，嵌入Simulink模型中。上一頁下一頁返回 【例2.7】 某非線性系統(tǒng)如圖2.17所示，試求r(t)=2.1(t) 時系統(tǒng)的動態(tài)響應。 【解】構建系統(tǒng)的Simulink模型，如圖2.18所示。為了便于研究問題的方便起見，不采用Discontinuities 子庫中的Sa

29、turation模塊，而選擇User-Defined-Functions子庫中的MATLAB Fcn模塊, 并將參數MATLAB Function設置為：saturation _zone。 上一頁下一頁返回圖2.17返回上一頁下一頁返回圖2.18 (exam2_7.mdl)然后編制M函數內容如下：% saturation_ zone functionfunction uo = saturation_ zone (ui)if ui=1 uo=1;elseif uia=0 0;1 -1; b=1;0; c=0 1; d=0; T=0.05; a1,b1,c1,d1= c2dm(a,b,c,d,T,

30、 zoh); 并分別鍵入K=1，10，30。選擇定步長的discrete（no continuous states）算法，步長為0.05。啟動仿真后，得到如圖2.27所示的仿真結果。 )(T)(T上一頁下一頁返回 圖2.26 (exam2_8.mdl)圖2.27(a) K=1圖2.27(b) K=10圖2.27(c) K=30上一頁下一頁返回 在許多實際應用場合下，希望能盡可能快地獲取仿真結果。這就構成了實際的快速需求和仿真速度之間的矛盾。當采用普通微型計算機對某些控制系統(tǒng)實施仿真時，這種矛盾往往變得非常突出，從而提出了快速仿真的需求。 快速數字仿真算法的種類很多，它們各有優(yōu)缺點，目前還沒有一

31、種算法能夠適用于所有問題。本節(jié)將介紹常用快速數字仿真算法。上一頁下一頁返回3.4.1 仿真中對快速性的需求3.4.2 相匹配原理 3.4.3 替換法3.4.4 根匹配法 上一頁下一頁返回 通常需要采用快速仿真的問題可以分為下列幾個方面： 利用仿真技術進行控制系統(tǒng)的參數優(yōu)化設計時，需要反 復地對系統(tǒng)進行仿真計算，以獲得滿意的工程參數； 在數學-物理混合仿真中，由于有實物系統(tǒng)介入仿真模型，要求仿真模型的時間比例尺與真實系統(tǒng)的時間比例尺完全相同（即所謂實時仿真），如果系統(tǒng)比較復雜或者方程個數很多，要在規(guī)定時間內完成一次系統(tǒng)中每一個方程的計算，就要求仿真的計算速度比較快；上一頁下一頁返回 在復雜系統(tǒng)的

32、控制中，常常需要在線用仿真方法對被控系統(tǒng)的狀態(tài)進行預測，以確定系統(tǒng)的控制策略，這就要求仿真模型的時間比例尺小于系統(tǒng)的運行時間比例尺（即所謂超實時仿真），從而對仿真速度提出了更高的要求。 一般而言，對于第一種情況，仿真速度慢意味著不能及時獲得仿真結果，并且占用大量的機時，費用較高，效率低下；而對于后兩種情況，仿真速度慢意味著仿真的目的達不到或者失去實際應用價值。上一頁下一頁返回 算法的計算精度和速度是一對矛盾，在研究快速數字仿真算法時，計算速度成為矛盾的主要方面。通常在滿足計算穩(wěn)定性及工程精度要求的前提下，應當盡可能提高仿真計算速度。 一般快速數字仿真算法有如下兩點基本要求： 每步計算量要小；

33、算法具有良好的穩(wěn)定性，允許采用較大的計算步長，同時又能保證必要的計算精度。 應當根據實際情況，盡可能地減小仿真中附加信息的計算和處理。 上一頁下一頁返回相匹配的含義是，如果被仿真系統(tǒng)的數學模型是穩(wěn)定的，則其仿真模型也應該是穩(wěn)定的，并且二者的瞬態(tài)、穩(wěn)態(tài)特性一致。如果對于同一輸入信號，二者的輸出具有相一致的時域特性，或者二者具有相一致的頻率特性，則稱仿真模型與原系統(tǒng)模型相匹配。將高階系統(tǒng)的數學模型G(s)轉換成與之相匹配、每步計算量較小、允許采用較大步長且具有合理精度的仿真模型G(z) ，就可以利用對應的差分方程進行快速仿真。從G(s)直接推導出G(z)的方法有兩種。一種稱為替換法，另一種稱為根匹

34、配法。 上一頁下一頁返回1替換法的基本思想2雙線性變換法3雙線性變換法步驟4雙線性變換法特性分析5雙線性變換法仿真實例上一頁下一頁返回根據控制理論，s域和z域之間的準確映射關系為 （2.67）或 （2.68）式中，T為采樣周期，亦為仿真計算時的步長。（2.68）式是超越方程，很難由 得到一個關于變量u和y的線性差分方程。 Tsez zTsln1)()()(zUzGzY上一頁下一頁返回ln(z)可以展開成下列的無窮級數 （2.70） 取該級數的第一項，則可以得到s域與z域之間一種近似映射關系 （2.71）或 （2.72）利用（2.71）式可以把G(s)轉換為G(z)。數學上稱這種變換方法為雙線性

35、變換法（或Tustin法）。 121233) 1() 1() 12(1) 1() 1(31112lnmmzzmzzzzz112zzTs2/12/1TsTsz上一頁下一頁返回 雙線性變換法是否滿足相匹配原理呢？ 首先，根據相匹配原理，要求當G(s)穩(wěn)定時，G(z)也應該是穩(wěn)定的。 將 代入（2.72）式，得 即有 （2.73） 由（2.73）式可知,雙線性變換法將左半s平面映射到z平面的單位圓內。也就是說，如果原來系統(tǒng)的G(s)是穩(wěn)定的，則通過雙線性變換法得到的G(z)必然是穩(wěn)定的。js)(21)(21jTjTz22222)2()21 ()2()21 (TTTTz上一頁下一頁返回 其次，雙線性變

36、換后的穩(wěn)態(tài)增益不變。 設連續(xù)系統(tǒng)的傳遞函數為 其穩(wěn)態(tài)增益為bn/an 。若將G(s)進行雙線性變換，有顯然，G(z)的穩(wěn)態(tài)增益仍然是bn/an。 nnnnmmnasasbsbsG11)(nnnnmmnazzTazzTbzzTbzG11)112()112()112()(上一頁下一頁返回uy ssUsYsG1)()()(11112)()()(zzTzUzYzG) 1()(2) 1()(kukuTkyky) 1()(2) 1(kykyTky上一頁下一頁返回 最后，雙線性變換法具有一定的精度。 設有方程 （2.74）即 （2.75）將（2.71）式代入上式，得 （2.76）由此可得差分方程 （2.77

37、）這相當于應用AM2法（也稱為梯形法）于（2.74）式 。 綜上所述，雙線性變換法符合相匹配原理，是一種具有一定計算精度的絕對穩(wěn)定算法。 雙線性變換法適合于以分子分母多項式形式描述的G(s)。 上一頁下一頁返回【例2.9】 已知連續(xù)系統(tǒng)的傳遞函數為 （2.78）試采用雙線性變換法求出對應的脈沖傳遞函數和差分方程，并對所得結果進行分析?！窘狻?將（2.71）式代入上式，得脈沖傳遞函數 （2.79） 12)(2ssssG22)2()2() 1)(1(21)112(2)112(112)()()(TzTzzTzzTzzTzzTzUzYzG22122)22()22(211)2(2zTTzTTzTT上一頁

38、下一頁返回)2()()2(2)2()22() 1()22(2)(22kukuTTkyTTkyTTky122TT 于是,差分方程為 由（2.79）式可知，因為 所以，G(z)是穩(wěn)定的。 G(s)的分子多項式為1階，分母多項式為2階，而G(z)的分子、分母多項式的階次相同，均為2階。 G(s)的穩(wěn)態(tài)增益為0， G(z)的穩(wěn)態(tài)增益也為0。 上一頁下一頁返回為了考慮雙線性變換法所得仿真模型的精度，除了可以在時域中進行討論外，也可以從頻域的角度進行分析。事實上將，分別代入（2.78）式和（2.79）式，可得 （2.80） （2.81）取T=1s，分別求出（2.80）式和（2.81）式的幅頻特性和相頻特性

39、，如表2.5所示。 js Tjez1)(2)()(2jjjjG2232)1 ()(2j22)22() 1)(1()2(2)(TTeeeTTeGTjTjTjTj上一頁下一頁返回表2.5 G(s)與G(z)頻率特性比較(rad/sec)0.81.01.21.5幅頻特性連續(xù)系統(tǒng)模型0.099010.27520.44120.48780.500.49180.4615離散化模型0.09910.2770.4470.4930.4980.4760.417相頻特性連續(xù)系統(tǒng)模型78.5856.6028.0712.680-10.39-22.62離散化模型78.5756.3626.519.57-5.0

40、6-17.67-33.56上一頁下一頁返回 設連續(xù)系統(tǒng)的傳遞函數為 （2.82）采用雙線性變換法求取仿真模型的步驟如下： 將G(s)中的s按（2.71）式替換得到G(z)，并整理成有理分式形式 （2.83） 根據得到的G(z)，由Y(z)=G(z)U(z)，兩邊取Z反變換，得到便于計算機遞推計算的差分方程。 （2.84） mnasasbsbsUsYsGnnnnmmn11)()()(nnnnzzzzzG111101)()() 1()()() 1()(101nkukukunkykykynn上一頁下一頁返回 雙線性變換法是一種絕對穩(wěn)定的算法，可采用大步長。 雙線性變換法是由（2.70）式取一次項近似

41、得到的，所以對應的差分方程雖具有一定的精度，但當T取得太大時精度會降低。因此在選取T時，主要應考慮仿真精度和仿真速度的需求。經驗表明，若取 ( 為被仿真系統(tǒng)的自然頻率)，可以保證較好的仿真精度。 雙線性變換法得到的差分方程的每項系數，都是預先一次計算確定的，在仿真遞推過程中不需要重新計算，因此仿真遞推的計算量較小。對于n階連續(xù)系統(tǒng)，每步計算量約為(2n+1)次乘法。所以，該法適合于快速仿真。nT152n上一頁下一頁返回雙線性變換法中由G(s)到G(z)的轉換， 以及由G(z)到差分方程的轉換，都可很容易地編程實現。雙線性變換法具有串聯性。如圖2.29所示，G1(s)和G2(s)相串聯。若令 G

42、(s)= G1(s). G2(s),則有 G(z)= G1(z). G2(z) 圖2.29式中G(z) ,G1(z),G2(z) 分別是G(s),G1(s) ,G2(s)的雙線性變換結果。這就是說，雙線性變換沒有破壞拉氏變換的串聯性，可以用簡單的低階環(huán)節(jié)組成復雜的高階系統(tǒng)。上一頁下一頁返回雙線性變換公式不僅可以方便地應用于傳遞函數，也可以應用于狀態(tài)空間模型。雙線性變換后系統(tǒng)的階次不變，且分子、分母多項式具有相同的階次。將（2.71）式代入（2.82）式，有 （2.85）可見，G(z)的分子、分母多項式均為n次，在分子上增添了nm個z=-1的零點。由此得到的差分方程是多步關系式（當n1時）。由于

43、被仿真系統(tǒng)一般僅給出在初始時刻處nnnnmmnazzTazzTbzzTbzG11)112()112()112()(niimiimnnnpzzzzab11)()() 1(上一頁下一頁返回的系統(tǒng)的初始值，而在以前各采樣點處的信息需要經過變換求得，這會給雙線性變換法的應用帶來一定的麻煩?？梢宰C明，雙線性變換后得到的脈沖傳遞函數頻率特性與原連續(xù)系統(tǒng)的頻率特性在低頻段相差較小，在高頻段差別較大。這一點也可以從表2.5看出。所以，雙線性變換法主要用于有限帶寬的系統(tǒng)。 上一頁下一頁返回 利用MATLAB的控制系統(tǒng)工具箱中的c2dm函數，很容易實現連續(xù)系統(tǒng)傳遞函數到離散系統(tǒng)脈沖傳遞函數之間的轉換?！纠?.10

44、】考慮如圖2.30所示的單位反饋系統(tǒng)，試用雙線性變換法對其進行仿真 (r(t)=1(t) 。上一頁下一頁返回圖2.30 【解】采用如圖2.31所示的參數化Simulink模型，其中Discrete Filter模塊中參數Numerator設置為：nd1，Denominator為：dd1，Sample time為：0.05。其余兩個Discrete Filter模塊中的參數分別為nd2，dd2；nd3，dd3；Sample time均為0.05。 在MATLAB命令窗口運行下列指令T=0.05; n1=1 3.33; d1=1 24; n2=10; d2=1 0; n3=1; d3=1 10 1

45、6; nd1,dd1=c2dm(n1,d1,T,tustin);nd2,dd2=c2dm(n2,d2,T, tustin); nd3,dd3= c2dm(n3,d3,T,tustin);上一頁下一頁返回圖2.31 (exam2_10.mdl)選擇定步長的discrete（no continuous states）算法，步長為0.05，啟動仿真，可以得到如圖2.32所示的仿真結果。上一頁下一頁返回圖2.321根匹配法原理 2根匹配法步驟 3根匹配法特性分析 上一頁下一頁返回根匹配法的基本思想是要使離散化模型的瞬態(tài)特性和穩(wěn)態(tài)特性與原連續(xù)系統(tǒng)保持一致。更明確地說，就是要使離散化后所得脈沖傳遞函數的零

46、點和極點與原連續(xù)系統(tǒng)傳遞函數的零點和極點相匹配。所以，這種算法又稱為零極點匹配法。先根據s域和z域之間的準確映射關系 ,直接由G(s)的零點和極點確定對應G(z)的零點和極點；再根據一定的原則確定G(z)的增益，則可以很方便地由G(s)求出與之對應的G(z)，使二者的瞬態(tài)特性和穩(wěn)態(tài)特性保持一致。顯然，根匹配法適合于以零極點形式描述的G(s) 。 Tsez 上一頁下一頁返回Tsez 求出連續(xù)系統(tǒng)的零點和極點，即將系統(tǒng)的傳遞函數寫成下列的零極點形式 （2.87） 利用關系 在z平面上一一對應地確定零極點的位置。對于（2.87）式表示的系統(tǒng)，與z1,z2,zm相對應的零點為eTz1,eTz2,eTz

47、m;與p1,p2,pm相對應的極點為eTp1,eTp2,eTpm 。于是，有 （2.88） mnpspspszszszsKsUsYsGnm)()()()()()()(2121)()()()()(2121nmTpTpTpTzTzTzzezezezezezezKzG上一頁下一頁返回s0Tez)()()()()(2121nmTpTpTpmnTzTzTzzezezezzezezezKzG上一頁下一頁返回 在G(z)的分子上配上n-m個附加零點，使G(z)的分子多項式和分母多項式的階次相同。因為當nm，即連續(xù)系統(tǒng)的分母多項式階次高于分子多項式階次時，在s平面的無窮遠處實際上還存在n-m個零點。因此，在z

48、平面上必須再配上n-m個相應的零點。如果認為n-m個零點位于s平面上負實軸的無窮遠處，即 ,那么z平面上相應的零點由可知，應配在原點處。于是，對應的脈沖傳遞函數為 （2.89） 根據G(s)的特征確定Kz，一般有兩種方法：一種是選擇適當的輸入信號，應用終值定理，分別確定連續(xù)系統(tǒng)模型G(s)和離散化仿真模型G(z)的終值，然后根據終值（不等于零或無窮大）相等的原則確定G(z)的增益Kz。另一種是使由G(s)和G(z)分別得到的頻率特性在某一臨界頻率處相同，從而確定Kz。一般對于低通濾波器應使其低頻增益相同，而對于帶通濾波器則應使中頻增益相同。 根據得到的G(z)，由Y(z)=G(z)U(z)，兩

49、邊取Z反變換，得到便于計算機遞推計算的差分方程。 上一頁下一頁返回)()()(1trtytyT)1(111)(111TsTsTsG1/)(TTzezzKzG上一頁下一頁返回 【例2.11】 設某連續(xù)系統(tǒng)的微分方程為 試用根匹配法確定其離散化模型。 【解】 首先寫出系統(tǒng)的傳遞函數 它有一個極點p1=-1/T1和一個無窮遠處的零點。 由（2.89）式可得脈沖傳遞函數 式中，T為采樣周期。 根據拉氏變換的終值定理，有 而根據Z變換的終值定理，可得 令兩終值相等，有 增益Kz為1111lim1)(lim)()(lim)(1000ssTssssGsRssGysss11/111111lim1)(1lim)

50、()(1lim)(TTzTTzzzzeKzzezzKzzzzzGzzzRzGzzy111/TTzeK1/1TTzeK上一頁下一頁返回 于是，求得的離散化模型為 根據可以進一步求出差分方程為 )()1 () 1()(11/krekyekyTTTT1/111111)1 ()(zeeezzezGTTTTTTTT上一頁下一頁返回102)()()(222nnnsssUsYsG)()(212pspssGn22, 11nnjp22211211TjTTpTjTTpnnnneeezeeez上一頁下一頁返回 【例2.12】 二階系統(tǒng)的傳遞函數為試用根匹配法求其離散化模型。 【解】首先，將原系統(tǒng)的傳遞函數寫成零極點

51、形式式中根據（2.86）式，得到z平面上的兩個極點 原系統(tǒng)有兩個無窮遠處的零點，相應地應在z平面上配上兩個原點處的零點，則對應的脈沖傳遞函數為根據終值定理分別求出單位階躍輸入作用下的終值為式中 TnTzTTjTjTzTjTTjTznnnnnnnnnnezTezzKezeeezzKeezeezzKzG222221122112)1cos(2)()()(2222baKzzzGzzysssGyzzs11)(1lim)(11)(lim)(10TnTnnebTea22)1cos(2上一頁下一頁返回11baKzbaKz1212211)1 ()(bzazbabazzzbazG)()1 ()2() 1()(ku

52、bakbykayky令兩終值相等，得于是 最后，得到離散化模型為對應的差分方程為上一頁下一頁返回2) 1()()()(sssUsYsG2)() 1()(TzezzzKzG22) 1()(,1)(zTzzUssU上一頁下一頁返回Tsez 【例2.13】 連續(xù)系統(tǒng)的傳遞函數為 試用根匹配法求其離散化模型。 【解】根據 以及G(s)的零極點，可得 在利用終值定理確定Kz時，由于在s平面上存在一個原點處的零點，加上躍階信號后，終值為零。為了使終值取得不為零的有限值，應當加上單位斜坡信號，即于是，有 從而得到最后，離散化模型為11) 1(lim1)(lim)(22020sssssssGyss222121

53、)1 () 1()() 1(1lim) 1()(1lim)(TzTzzzeTKzTzezzzKzzzTzzGzzyTeKTz2)1 (211222)1 (1)1 ()() 1()1 ()(zezTeezzzTezGTTTT上一頁下一頁返回由根匹配法所采用映射關系 可以看出，如果原系統(tǒng)模型G(s)是穩(wěn)定的，則它的全部極點都位于左半s平面。根匹配法具有一定的精度。 以【例2.13】為例，比較G(s)和G(z)的頻率特性。 取采樣周期T=1s，由將G(s)和G(z)的幅頻特性和相頻特性列于表2.6。Tsez 上一頁下一頁返回2121)() 1()1 ()(ezzzezG2)3679. 0() 1(4

54、 . 0zzz （rad/sec）0.81.01.22.0幅頻特性 連續(xù)系統(tǒng)模型0.09900.27520.44120.48780.50.49180.4 離散化模型0.09540.26590.43070.48120.50.50010.449相頻特性 連續(xù)系統(tǒng)模型75.5856.628.0712.680-10.39-36.87 離散化模型80.5062.3739.627.919.1012.450.35表2.6 G(s)與G(z)頻率特性比較上一頁下一頁返回 比較表2.6中G(s)和G(z)的頻率特性可知，二者的幅頻特性是比較接近的，而離散化模型G(z)在相位上超前較大。當=1時

55、，G(s)的相位為0，而G(z)的相位為19.1。一般而言，如果T 取得比較大，計算精度還會降低。 上一頁下一頁返回實時數字仿真通常是指把一個數字仿真過程嵌入到一個具有實物模型的實際系統(tǒng)或仿真系統(tǒng)的運行過程中，這就要求仿真模型的時間比例尺必須完全等于原始模型的時間比例尺。當仿真系統(tǒng)有實際裝置或操作人員介入試驗回路時，為了真實反映實際系統(tǒng)的工作情況，要求仿真系統(tǒng)必須按照實際系統(tǒng)運行的時序要求來完成數字仿真過程的各個步驟。 上一頁下一頁返回3.5.1 實時數字仿真的概念 3.5.2 實時數字仿真算法的特性 3.5.3 常用的實時數字仿真算法 上一頁下一頁返回【例3.14】 飛行器控制系統(tǒng)的半實物仿

56、真 一個飛行器控制系統(tǒng)可以分為被控對象和控制器兩個子系統(tǒng)。在導彈、運載火箭的設計和飛行試驗前的仿真中，控制對象的運動通常為彈體運動，包括質心運動，繞質心的旋轉運動和彈體振動等。顯然，在實驗室里難以直接實現被控對象的運動。所以，通常采用一定的數值方法，將描述彈體運動的數學模型進行離散化，建立數字仿真模型，然后編程運行來代替彈體運動，并將其稱為數字處理過程。而控制器及其執(zhí)行機構，包括舵式控制發(fā)動機等則用實物表示，并將其稱為實物系統(tǒng)過程。于是，由實物系統(tǒng)過程和數字處理過程等組成仿真回路處理過程等組成仿真上一頁下一頁返回回路，即半實物數字仿真模型。由于在這樣的飛行器控制系統(tǒng)的仿真中介入了實物，因而在彈

57、體運動的數字處理過程中信息的采樣輸入、仿真模型的計算、數據的輸出都必須在規(guī)定的時間內按照實際運行的時序完成，進行實時數字仿真。 可以將控制系統(tǒng)的實時數字仿真描述如下。 假設一個實際的控制系統(tǒng)是由實物系統(tǒng)過程A和實物系統(tǒng)過程B兩部分組成，如圖2.33所示。而在實時仿真過程中，常常將系統(tǒng)的一部分（例如，實物系統(tǒng)過程A）用數學模型代替。在實際實現時，用計算機的數字處理過程A來代替，如圖2.34所示。 上一頁下一頁返回上一頁下一頁返回圖2.33 圖2.34圖2.34中計算機數字處理過程A的輸入是z(t)，輸出是y(t)。對于任意一組輸入值，經過采樣控制和A/D轉換，由計算機數字處理過程A的計算得到相應

58、的輸出值ym。將ym稱為數字處理過程A對于輸入zm的響應。響應ym對于輸入zm的時間延遲稱為響應時間。過程B為實物系統(tǒng)，通過D/A轉換和輸出控制接受過程A的輸出量ym，并輸出z(t)。計算機數字處理過程A必須在與實物系統(tǒng)同步的條件下獲取動態(tài)輸入信號，并實時地產生動態(tài)輸出響應。這時，仿真模型的輸入和輸出都是具有固定采樣周期的數值序列。在數字處理過程A滿足系統(tǒng)各項功能要求的情況下，對于任意特定的輸入z(t)，響應時間都能滿足系統(tǒng)所要求的時間限制，則這種數字仿真過程通常為實時數字仿真過程。 上一頁下一頁返回 實時仿真模型中包含有實物或人，它具有如下特性。 1實時性 實時性主要體現在如下幾個方面。（1

59、）實物模型的運行時間 這是實物模型在輸入作用下的運行時間。若用計算機來代替，則響應時間的要求就是實時要求，其具體數值必須滿足隨機尖峰負載時的處理要求。（2）數據轉換（如數模D/A、模數A/D轉換）的時間 從發(fā)出轉換命令到完成的時間。D/A、A/D轉換本身需要一定的時間。在仿真中，希望越快越好，但轉換速度越快，系統(tǒng)成本越高。因此要根據實際需要確定。上一頁下一頁返回 （3）輸入/輸出響應時間 計算機接收動態(tài)輸入到產生動態(tài)輸出的響應時間，有固定的最小的要求。它是計算過程的延時，可以通過采用不同的算法和計算步長來達到要求。如果無法達到要求，則必須更換仿真系統(tǒng)所用的計算機。 2周期性 在實時仿真模型中，

60、模型和各個子模型都有規(guī)律性和周期要求。它們以設定的周期時間接受輸入信息并一幀一幀地產生輸出信息。每一個子系統(tǒng)都必須按照規(guī)定的順序，在分配給它們的時間片內完成信息采樣、計算、恢復和輸出等任務。上一頁下一頁返回 3可靠性 由于有實物包含在內，與一般計算機仿真系統(tǒng)相比，實時仿真系統(tǒng)要復雜得多，并且一般均為專用的仿真系統(tǒng)。因此，仿真系統(tǒng)的可靠性是第一位的。在【例2.14】的飛行器半實物仿真中，系統(tǒng)運行的任何差錯都可能導致仿真的失敗。而引起差錯的原因，一是數學模型或仿真模型或相應的計算機軟件有錯；二是外界干擾破壞了程序的正常運行；三是硬件系統(tǒng)故障。第一類錯誤可以通過修改模型、選擇恰當的仿真算法或調整軟件