經(jīng)典的連續(xù)系統(tǒng)仿真建模方法學(xué)

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