PAGEPAGE47第二章經典的連續系統仿真建模方法學本章討論經典的連續系統數字仿真的原理與方法，內容包括連續系統數字仿真的基本概念、經典的數值積分法、經典的線性多步法等。在數字計算機上進行連續系統仿真，首先要將連續模型離散化，因此，2.1節首先討論離散化原理及要求，這是連續系統仿真的基礎。然后，2.2節對經典的數值積分法龍格-庫塔法及其它典型的數值積分法仿真建模原理進行詳細分析，并通過實例說明其應用要點；而2.3節對經典的線性多步法進行了介紹．2.1離散化原理及要求在數字計算機上對連續系統進行仿真時，首先遇到的問題是如何解決數字計算機在數值及時間上的離散性與被仿真系統數值及時間上的連續性這一基本問題。從根本意義上講，數字計算機所進行的數值計算僅僅是“數字”計算，它表示數值的精度受限于字長，這將引入舍入誤差；另一方面，這種計算是按指令一步一步進行的，因而，還必須將時間離散化，這樣就只能得到離散時間點上系統性能。用數字仿真的方法對微分方程的數值積分是通過某種數值計算方法來實現的。任何一種計算方法都只能是原積分的一種近似。因此，連續系統仿真，從本質上是對原連續系統從時間、數值兩個方面對原系統進行離散化，并選擇合適的數值計算方法來近似積分運算，由此得到的離散模型來近似原連續模型。如何保證離散模型的計算結果從原理上確能代表原系統的行為，這是連續系統數字仿真首先必須解決的問題。設系統模型為：，其中u(t)為輸入變量，y(t)為系統變量；令仿真時間間隔為h，離散化后的輸入變量為，系統變量為，其中表示t=kh。如果，，即，（對所有k=0,1,2,…），則可認為兩模型等價，這稱為相似原理（參見圖2.1）。原連續模型原連續模型仿真模型u(t)hy(t)-+圖2.1相似原理實際上，要完全保證是很困難的。進一步分析離散化引入的誤差，隨著計算機技術的發展，由計算機字長引入的舍入誤差可以忽略，關鍵是數值積分算法，也稱為仿真建模方法。相似原理用于仿真時，對仿真建模方法有三個基本要求：（1）穩定性：若原連續系統是穩定的，則離散化后得到的仿真模型也應是穩定的。關于穩定性的詳細討論將在2.4節中進行。（2）準確性：有不同的準確性評價準則，最基本的準則是：絕對誤差準則：相對誤差準則：其中規定精度的誤差量。（3）快速性：如前所述，數字仿真是一步一步推進的，即由某一初始值出發，逐步計算，得到，每一步計算所需時間決定了仿真速度。若第k步計算對應的系統時間間隔為計算機由計算需要的時間為Tk，則，若Tk=hk稱為實時仿真，Tkhk稱為超實時仿真，而大多數情況是Tkhk，對應于離線仿真。圖2.2數值積分法原理連續系統數字仿真中離散化最基本的算法是數值積分算法。對于形如的系統，已知系統變量的初始條件，現在要求隨時間變化的過程計算過程可以這樣考慮（參見圖2.2）：首先求出初始點的，微分方程可以寫作：圖2.2數值積分法原理（2.1）圖2.2所示曲線下的面積就是，由于難以得到f(y,u,t)積分的數值表達式，人們對數值積分方法進行了長期探索，其中歐拉法是最經典的近似方法。歐拉法用矩形面積近似表示積分結果，也就是當t=t1時，的近似值為：(2.2)重復上述作法，當時所以，對任意時刻tk+1，有：（2.3）令稱為第步的計算步距。若積分過程中步距不變，可以證明，歐拉法的截斷誤差正比于。為進一步提高計算精度，人們提出了“梯形法”。梯形法近似積分形式如式(2.4)所示，令：已知：時的近似值,那么：(2.4)可見，梯形法是隱函數形式。采用這種積分方法最簡單的預報校正方法是用歐拉法估計初值，用梯形法校正，即：(2.5)(2.6)式(2.6)稱作預報公式，采用歐拉法，式(2.5)為校正公式，采用梯形法。用歐拉法估計一次的值，代入校正公式得到的校正值。設是規定的足夠小正整數，稱作允許誤差，若i=0,i+1=1稱作第一次校正；稱作第二次校正；通過反復迭代，直到滿足，這時是滿足誤差要求的校正值。上述方法是針對(2.3)式所示的微分方程在已知初值情況下進行求解，因此也稱為微分方程初值問題數值計算法，為統一起見，本書中稱為數值積分法。連續系統數字仿真的離散化方法有兩類，它們是數值積分方法和離散相似方法，本章討論數值積分法。數值積分方法采用遞推方式進行運算，而采用不同的積分方法會引進不同的計算誤差，為了提高計算精度，往往會增加運算量。就同一種積分算法而言，為提高計算精度，減小積分步距,計算量增大，影響系統運算速度。因此，計算精度和速度是連續系統仿真中常迂到的一對矛盾，也是數字仿真中要求解決的問題之一。也就是說，選擇合適的算法、合適的軟、硬件環境，在保證計算精度的前提下，考慮怎樣提高仿真的速度。經典的數值積分法可分為兩類：單步法與多步法，下面我們將分別來介紹這兩類方法中的最常用的算法。龍格庫塔法2.2.1龍格-庫塔法基本原理在上一節中我們已經講過，在連續系統的仿真中，主要的數值計算工作是對的一階微分方程進行求解。因為若令：,則有(2.7)因此主要問題就是如何對進行數值求解，即如何對進行積分。通常稱作“右端函數”計算問題。已知，假設我們從跨出一步，，時刻為，可以在附近展開臺勞級數，只保留項，則有：(2.8)在式(2.8)中括號后的下標0表示括號中的函數將用均同。我們假設這個解可以寫成如下形式：(2.9)其中，對式右端的函數在處展成臺勞級數，保留項，可得：將代入(2.9)式，則有：將上式與(2.8)式進行比較，可得：可見有四個未知數但只有三個方程，因此有無窮多個解，若限定，則可得其中一個解：將它們代入(2.9)式可得一組計算公式：(2.10)其中，若寫成一般遞推形式，即為：(2.11)其中，由于(2.8)式只取了兩項，而將以上的高階項略去了，所以這種遞推公式的截斷誤差正比于的三次方，又由于計算時只取了及項，故這種方法被稱為二階龍格庫塔法（簡稱RK-2）。根據上述原理，若在展成臺勞級數時保留及項，那么可得一套截斷誤差正比于的四階龍格--庫塔法（簡稱RK-4）公式：(2.12)其中：由于這組計算公式有較高的精度，所以在數字仿真中應用較為普遍。為了幫助讀者更好地掌握這種方法的使用，下面我們對龍格--庫塔法的特點作一些介紹。由于在解時，可以得到許多種龍格--庫塔公式，經常使用的除(2.11)及(2.12)式給出的兩組公式外，還有：(2.13)其中(2.14)其中13)式為二階龍格--庫塔公式，(2.14)式是Shampine提出的，稱為RKS3-4公式。對上述公式可以這樣解釋：因為滿足這個微分方程，已知處的導數為。根據外推原理，假設在--這個區間中的導數不變，則可得在的估計值0。令為這個估計值處的導數，因此(2.11)式就是按這個導數線性外推而得的。而(2.12)式則是4點導數的加權平均和，由于第2及第3點是以為步長，精度較高，所以加權平均時各取2份，第1、4點各取1份。不論幾階龍格--庫塔法，它們的計算公式總是由兩部分組成的，即是上一步的結果及步長乘以--中各點導數的加權平均和。因此各種龍格庫塔法可以寫成如下一般形式：(2.15)其中式中各系數滿足以下關系(2.16)一般來講，為了減少計算量總是希望減少每步計算右端函數的次數（即減少計算的次數）。可以證明，1階龍格--庫塔公式至少要計算一次右端函數，即(2.11)式中的s(稱為級數)至少等于1，而2階公式；3階公式；4階公式；依此類推。有時為了某種特殊需要，可以選擇的計算公式。表2.1為常見的各階龍格庫塔公式，其中第二項是4階4級RK-4公式1，而第三項是4階5級RK-4計算公式(簡稱RKM45)。其中，后三項有誤差估計公式，所以可以用誤差來控制積分步長，即每積分一步，利用誤差估計公式估計出這一步的截斷誤差。若它大于允許的誤差，則縮小步長，反之可以放大步長。表2.1常見的龍格--庫塔公式名稱公式特點二階龍格庫塔法四階龍格庫塔法四階龍格庫塔--默森法有誤差估計公式：二階龍格庫塔--費爾別格法有誤差估計公式：四階龍格--庫塔--夏普法有誤差估計公式：所有龍格庫塔公式都有以下特點：在計算時只用到，而不直接用等項。換句話說，在后一步的計算中，僅僅利用前一步的計算結果，所以稱為單步法。顯然它不僅能使存儲量減小，而且此法可以自啟動，即已知初值后，不必用別的方法來幫助，就能由初值逐步計算得到后續各時間點上的仿真值。步長在整個計算中并不要求固定，可以根據精度要求改變，但是在一步中，為計算若干個系數(俗稱：龍格庫塔系數)，則必須用同一個步長。龍格--庫塔法的精度取決于步長的大小及方法的階次。許多計算實例表明：為達到相同的精度，四階方法的可以比二階方法的大10倍，而四階方法的每步計算量僅比二階方法大一倍，所以總的計算量仍比二階方法小。正是由于上述原因，一般系統進行數字仿真常用四階龍格--庫塔公式。值得指出的是：高于四階的方法由于每步計算量將增加較多，而精度提高不快，因此使用得也比較少。若在展開臺勞級數時，只取這一項，而將以及以上的項都略去，則可得：這就是歐拉公式，因此歐拉公式也可看作是一階龍格--庫塔公式。它的截斷誤差正比于，是精度最低的一種數值積分公式。讀者也許會問：如果我們將步長取得十分小，那么歐拉公式也可以獲得很高的計算精度嗎？從理論上講是這樣，但是實際上由于計算機字長有限，在計算中存在著舍入誤差，它與計算次數成正比。計算步長很小，則一個系統的過渡過程的計算就會增加許多，舍入誤差就會十分明顯地表現出來，這樣一來，很難保證有較高的精度。2.2.2龍格--庫塔法的誤差估計及步長控制一個高精度的仿真方法必須將步長的控制作為必要的手段。實現步長控制涉及局部誤差估計和步長控制策略兩方面的問題。本節將對數值積分中單步法的誤差估計和步長控制作一介紹。龍格--庫塔數值積分方法的“誤差估計”和“步長控制”的基本想法如圖2.2所示,也就是，每積分一步都設法估計出本步的計算誤差k，然后判斷是否滿足允許誤差E，據此，選擇相應的步長控制策略，調整步長，再作下一步積分運算。允許允許誤差改變下一步計算步長-k+步長控制策略積分算法誤差估計本步計算圖2.2龍格庫塔法步長控制示意(1)龍格--庫塔法的誤差估計龍格--庫塔法的誤差估計，通常采用的方法是設法找到另一個低階(一般是低一階)的龍格--庫塔公式，則兩個公式計算結果之差可以被看作是誤差。例如，龍格--庫塔--默森法(Runge-Kutta-Merson)法[6]計算公式為：(2.17)其中它是一個4階5級公式。另外還可推出一個3階4級公式(2.18)用作為誤差估計，即(2.19)式(2.17)、(2.18)、及(2.19)簡稱為RKM3-4法。又比如，E.Fehlberg推導出的5階龍格--庫塔--費爾別格法[4]，它的計算公式為一個5階6級方法，另外用了一個四階5級方法求，也是用來估計誤差，這一套計算公式被公認是對非病態系統進行仿真最為有效的方法之一。由于它是5階精度，4階的誤差估計，因此被稱為RKF4-5階公式對，簡稱RK4-5法。具體系數如表2.2。表2.2Fehlberg公式中的各系數()(2.20)表中表示4階公式中的系數，因此誤差估計為(2.21)1016/13525/21621/41/40033/83/329/326656/128251408/2565412/131932/2197-7200/21977296/219728561/564302197/410451439/216-83680/513-845/4104-9/50-1/561/2-8/272-3544/25661859/4104-11/402/550當精度要求不太高時，RKF1-2法被廣采用。它的計算公式為(2.22)其中用另一個一階公式(2.23)來估計誤差(2.24)RKM3-4法與RKF4-5法是兩種應用得極為廣泛的數值積分算法。它們的共同特點是計算量較大，比如RKM3-4法，每次要計算5次,而只獲得四階精度及三階誤差估計，RKF4-5法每步要計算6次,能獲得5階精度及4階誤差估計。對于一般仿真問題，中等精度即可滿足要求，因此用RKM3-4公式對就夠了。1978年，夏普勒(Shampine)提出了RKS3-4公式對[7]，它每次只計算4次卻能獲得4階精度與3階誤差估計，如(2.25)式所示。(2.25)其中另外引入了一個3階公式(2.26)其中(2.27)它正好是下一次計算時的，因此只是在第一步要多算一次，以后仍每步計算四次。RKS3-4的誤差估計為(2.28)(2).步長控制有以下幾種步長控制的方法：加倍－減半法設定一個最小誤差限和最大誤差限，當估計的局部誤差大于最大誤差時將步長減半，并重新計算這一步；當誤差在最大誤差與最小誤差之間時，步長不變；當誤差小于最小誤差時將步長加倍。每步的局部誤差通常取以下形式：(2.29)其中為利用本節的公式計算出來的誤差估計。由(2.29)式可知，當較大時，是相對誤差，而當的絕對值很小時，就成了絕對誤差。這樣做的目的是避免當的值很小時，變得過大。上述步長控制的策略可以表示為下式：如果則選擇(2.30)其中max和min為最大、最小誤差限。這種步長控制的方法簡便易行，每步附加計算量小，但是這種方法不能實現每步都是最優步長。 最優步長法為了使每個積分步在保證精度的前提下能取最大步長（或稱最優步長），可以設法根據本步誤差的估計，近似確定下一步可能的最大步長。。這種方法可以做到在規定的精度下取得最大步長，因此減少了計算量。具體策略如下：給定相對誤差限0，設本步步長為,本步相對誤差估計值為(如2.29式所示)。假定所采用的積分算法為階，則(2.29)式中的可表示為(2.31)其中是在積分區間中某處一些偏導數的組合，通常可取，因此有(2.32)若則本步積分成功。現在來確定下一步的最大步長，假定足夠小，即可認為，故下一步的誤差可能為為使，則有將(2.32)代入上式，得(2.33)若，則本步失敗。此時仍可采用(2.33)式，但它表示重新積分的本步步長。由于我們假定了足夠小，因此基本不變，故必須限制步長的無限放大和縮小，一般可限制的最大放、縮系數為10，即要求(2.34)當函數中含有間斷特性時，采用上述兩這控制方法在間斷點附近會出現步長頻繁放大、縮小的振蕩現象。由于最優步長控制方法是以本步誤差外推下一步步長，因此振蕩現象更為嚴重。解決這一問題的方法將在第7章討論。實時龍格－庫塔法仿真模型的運行速度往往與實際系統運行的速度不同。然而，當有實際的裝置或被訓練的人介入仿真過程時，就要求仿真模型的運行速度往往與實際系統運行的速度保持一致，這稱為實時仿真。一般的數值積分法難以滿足實時仿真的要求，這不僅僅是因為由這些方法所得到的模型的執行速度較慢，而且這些方法的機理不符合實時仿真的特點。下面讓我們以二階龍格-庫塔為例來分析。假設對如下一般形式的系統進行仿真：RK-2公式如下：即在一個計算步內分兩子步：首先在tk時刻利用當前的uk，yk計算K1。假定在的時間內計算機正好計算一次右端函數，然后，在tk+h/2時刻計算K2，盡管此時yk+1/2已經得到，但uk+1則無法得到。實時仿真除了對模型執行速度滿足要求外，還要求實時接收外部輸入，并實時產生輸出。為此，要么對uk+1也進行預報（這將加大仿真誤差），要么將仿真執行延遲h/2。后者的RK-2公式的計算流程可表示成圖2.4。這就是說，輸出要遲后半個圖2.4RK-2圖2.4RK-2的計算流程與此相類似，(2.12)式所示的RK-4公式也不適合于實時仿真，對此，讀者可自行分析。為了克服這個缺陷，人們提出了如下形式的實時2階龍格－庫塔法：圖2.5實時RK-2公式計算流程它的計算流程如圖2.5所示。由圖不難看出，首先在tk時刻利用當前的uk，yk計算K1。然后，在tk+h/2時刻計算K2，此時yk+1/2已經得到，由于計算一次右端函數需要的時間，uk+1/2也可得到，K2的計算就不會引入新的誤差，并能實時輸出yk+1。圖2.5實時RK-2公式計算流程文獻[1]介紹了幾種高階的實時RK公式，例如，4階5級的實時計算公式如下：其中它的計算流程如圖2.6所示：圖圖2.6實時RK4-5公式計算流程該公式的特點是：將一個仿真步分為5個子步，每個子步均能對外部輸入實時采樣，并保證了實時輸出。2.2.4面向方程的龍格－庫塔法仿真舉例為了加深對數值積分法的理解，以便能使用仿真程序來解決一兩個實際問題，本小節將給出面向方程的龍格－庫塔法仿真實例。“面向方程”是指仿真對象是以一階微分方程組形式描述的系統，如(2.35)式所示：(2.35)其中，是維狀態向量；是k維參數向量；是維控制向量；是維函數向量；是給定的初始狀態。本例所用仿真程序CSS1.C采用定步長RK-4公式作數值積分。CSS1.C是由C語言編制的面向方程的仿真程序(原程序見附錄)，其程序結構如圖2.7所示。程序中主要變量說明：T1仿真總時間T2積分步長T仿真時間Y(1)Y(20)系統狀態變量Y(21)Y(30)用戶自定義變量G(1)G(20)系統狀態變量的導數A(1)A(5)用戶輸出變量P(1)P(20)系統的參數變量N1方程階次N2參數個數N3輸出點數（在0TT1區間內，設輸出時間間隔為T3，則N3＝T1/T3+1）下面對仿真對象作一說明。設衛星在空中運行的運動方程為：(2.36)其中是重力系數(=401408km3/s)。衛星軌道采用極坐標表示，通過仿真，研究發射速度對衛星軌道的影響。這是一個二階微分方程組，為此，首先要將其轉換成一階微分方程組。若設，則有：(2.37)這是四個一階微分方程，有4個狀態變量及。今希望用直角坐標輸出，故還引入兩個定義變量，這樣又得到兩個代數方程：(2.38)根據衛星的發射速度，可以建立起方程組(2.37)的初始條件：km(衛星到地心的距離，即為地球之半徑)，，仿真研究的目的是要計算出當發射速度為8km/s，10km/s及12km/s時的衛星軌道。已知y(4)=y/y(1)，由y及可求出，它們分別是：0.00125，0.0015625，0.001875(1/s)，按CSS1.C的要求，為了仿真該系統，用戶要輸入系統的運動方程，即要將(2.37)及(2.38)式寫到程序中去。CSS01.C規定這部分寫在系統模型輸入(modsub)程序塊中；比如例中要求輸出則A(1)，A(2)，A(3)，A(4)四個數組都要用到，CSS1.C規定這部分寫在輸出轉換(repsub)程序塊中；狀態變量的初始數據及要改變的參量數值，CSS1.C規定它放在主程序中。此系統沒有參數變量。下面就是為仿真衛星軌道系統由用戶寫的一段仿真源程序。在系統模型輸入程序塊(modsub)中寫入以下程序：Y(21)=Y(1)cos(Y(2))Y(22)=Y(1)sin(Y(2))G(1)=Y(3)G(2)=Y(4)G(3)=-401408/(Y(1)Y(1))+Y(1)Y(4)Y(4))G(4)=-2Y(3)Y(4)/Y(1)在輸出轉換(repsub)程序塊中寫入以下程序：A(1)=Y(21)A(2)=Y(22)A(3)=Y(1)A(4)=Y(2)如果需要，可以在輸出打印(outsut)程序塊中加入Y(4)0=Y(4)的打印語句。程序寫好后則可按以下步驟進行操作：將CSS1.C程序及用戶寫好的仿真源程序輸入到計算機中。運行CSS1.C按計算機提問輸入必要的數據。計算機執行CSS1.C并在終端上或打印機上輸出運行結果。下面給出用CSS1對發射速度為8Km/s時衛星軌道仿真的運行過程和結果。C:/>CSS1.CInputtotalsimulationtimeT1,stept2=10000,200InputtheorderofsystemN1=4InputthenumberofunknownparametersN2=0InputthenumberofoutputpointsN3=51symbol191\f"Symbol"\s10.5InputthenumberofrunsJ1=3Inputinitialvaluesofstastevariables=6400，0，0，00125Inputthevaluesofparameters=0OUTPUTDEVICE2(PRINTER),5CRT?然后機器開始運算，下面是運算一次后的結果：Y(I)Y(I)maxY(I)min16666.64600.00000211.97503.000003.15935-.160054.00125.00000216400.00000-6665.12300226484.59500-6475.95700SIMULATIONSTEPT2=200.00000TA(1)A(2)A(3)A(4)200.000006205.06800-2625.199006403.97800.24990400.000005632.80400460.996306415.65600.49918600.000004719.816005857.741006434.28800.74728800.000003523.71300-5904.655006458.69700.993721000.000002118.523005751.833006487.375001.238141200.00000589.16300-6475.957006518.586001.480291400.00000-974.40360-3196.890006550.498001.720101600.00000-2482.78100-5245.305006581.283001.957621800.00000-3852.12900-3377.111006609.230002.193022000.00000-5008.32900-3907.278006632.829002.426572200.00000-5890.24000-1619.592006650.829002.658672400.00000-6452.096004470.647006662.297002.889732600.00000-6665.123006484.595006666.646003.120242800.00000-6518.49300-2042.907006663.655003.350703000.00000-6019.65600-2370.215006653.475003.581623200.00000-5194.113005789.948006636.627003.813483400.00000-4084.58200-3236.519006613.982004.046743600.00000-2749.507003792.250006586.734004.281773800.00000-1260.826005677.029006556.350004.518884000.00000299.14620-4177.597006524.511004.758254200.000001841.686004245.863006493.026004.999984400.000003276.52700990.172706463.740005.243984600.000004517.38300-1451.660006438.412005.490074800.000005487.769003588.771006418.604005.737905000.0000061