王 宇，孫永維，宋省身，史云龍

（1.空軍航空大學(xué) 軍事仿真技術(shù)研究所，吉林 長春 130022；2.空軍航空大學(xué) 飛訓(xùn)基地，吉林 長春 130062；3.空軍航空大學(xué) 航空航天情報系，吉林 長春 130022；4.中國空氣動力研究與發(fā)展中心 四川 綿陽 621000）

實時子步積分算法的設(shè)計與實現(xiàn)

王 宇1，孫永維2，宋省身3，史云龍4

（1.空軍航空大學(xué) 軍事仿真技術(shù)研究所，吉林 長春 130022；2.空軍航空大學(xué) 飛訓(xùn)基地，吉林 長春 130062；3.空軍航空大學(xué) 航空航天情報系，吉林 長春 130022；4.中國空氣動力研究與發(fā)展中心 四川 綿陽 621000）

積分算法是實時仿真的關(guān)鍵所在，傳統(tǒng)的單步法結(jié)構(gòu)簡單，但精度較差。而多步高階積分算法雖然精度高，卻由于結(jié)構(gòu)復(fù)雜實時性不強(qiáng)。文章在對實時仿真積分算法研究的基礎(chǔ)上，提出了實時子步積分算法，使不同的積分算法擁有相同的子步控制結(jié)構(gòu)。由于每個子步的運算量一樣，所以在步長和幀數(shù)一致的情況下，不同的積分算法運算時間是一致的。經(jīng)仿真實驗表明，該算法的子步結(jié)構(gòu)顯著提高了多步積分法的實時性，且誤差精度在可控范圍內(nèi)。

計算機(jī)仿真；實時性；子步積分算法；子步控制結(jié)構(gòu)

進(jìn)行實時仿真試驗時，仿真計算機(jī)與實物相連，仿真計算機(jī)必須在與真實系統(tǒng)同步的條件下獲取動態(tài)輸入信號，并實時地產(chǎn)生動態(tài)響應(yīng)。要實現(xiàn)實時仿真，必須采取相應(yīng)的實時仿真算法。而實時仿真中采用的算法需要有快速性、可取性、可靠性和相容性等特性[1]。傳統(tǒng)的單步積分法由于其步長易控性并有較好的實時性，在實時仿真中應(yīng)用較多，但是計算精度較差。多步積分法雖然計算精度高，但是計算步驟較復(fù)雜，實時性不強(qiáng)。

文中通過分析多種常用積分法結(jié)構(gòu)，首次提出了實時子步積分算法，該算法利用特殊的子步結(jié)構(gòu)改進(jìn)了幾種常用積分法，有效地整合了單步法與多步法的優(yōu)點，使新算法具有較好的實時性和精度[2]。

1 通用數(shù)字仿真算法

仿真算法在數(shù)字仿真的各個階段起著核心和支撐作用，一個好的仿真算法不僅能耗費較少的時間，而且也能得到比較高的精度。算法選擇是實時數(shù)字仿真首先碰到的問題，選擇合適的仿真算法可以使實時數(shù)字仿真以足夠的精度順利完成實驗任務(wù)。否則將會事倍功半，甚至在試驗過程中遭受種種挫折，嚴(yán)重影響進(jìn)度和質(zhì)量。

1.1 單步法

單步法中比較典型的有Runge-Kutta算法 （簡稱RK法）[3]。

一般的s級RK法有如下的形式：

最常用的Runge-Kutta法為四階RK法，亦稱為經(jīng)典RK法：

RK法的主要特點是計算tm+1點的值時，只需要用到tm點的值。其基本思想是在[tm，tm+1]內(nèi)多計算幾個點的斜率值Ki，然后將它們加權(quán)平均作為平均斜率K*，以構(gòu)造出具有更高精度的計算格式。顯式RK法在計算時，只用到K1，K2，…，Ki-1了的值。 而隱式 RK法在計算Ki時，不僅用到了 K1，K2，…，Ki-1的值，還可能用到Ki，Ki+1的值。顯然顯式RK法具有很好的實際意義。針對不同的需求，通過應(yīng)用Taylor級數(shù)匹配原理來確定公式中的參數(shù)值，可以構(gòu)造出各種RK法的公式。

1.2 線性多步法

Adams法是一種比較特殊的線性多步法。它計算ym+1的值，用到了 tm，ym的值和時刻 tm以前的導(dǎo)數(shù)值 fj=f（tj，yj），j＜m。Adams法也分為顯式和隱式2種。顯式Adams法是線性k+1步算法而隱式Adams法是線性k步算法。

例如四階Adams顯式（Adams-Bashforth，簡稱AB法）：

線性多步法比較典型的算法為Adams法[4]：

為了進(jìn)一步提高算法精度，又提出了基于AB法和AM法的Adams預(yù)估-校正法，利用AB法作預(yù)估，再用AM法作校正，以四階為例：

四階 Adams隱式（Adams-Monlton，簡稱AM法）：

1.3 單步法與多步法的比較

由于單步法和線性多步法采取了不同的策略來提高仿真算法的精度，因而可以對它們進(jìn)行對比研究[5]。

①單步法在計算ym+1的值時只用到y(tǒng)m的值，不需要ym-1，ym-2，…，各項的值，即單步法最大的特點是可以自啟動，而多步法計算 ym+1需要用到 ym，ym-1， …，ym-k+1以及對應(yīng)的 fm，fm-1，…，fm-k+1的值，所以要占用額外的儲存空間來保存前k個yj和 fj的 值[6-7]。

②多步法計算ym+1，需要計算多個右函數(shù)來給出ki值。

③當(dāng)進(jìn)行變步長運算時，單步法使用靈活，高階多步法相對來說計算復(fù)雜，因此在實時仿真中多用到單步法。

在基于控制系統(tǒng)的數(shù)字計算機(jī)中，由于對實時性要求強(qiáng)，多采用單步法來計算，這大大限制了算法的選擇性，本文針對這個問題，提出了實時子步積分算法[8]。

2 子步積分算法

2.1 設(shè)計思想

根據(jù)常用的幾種積分算法，我們利用子步積分結(jié)構(gòu)改進(jìn)四階以內(nèi)的 Adams顯式（AB）法，Adams預(yù)估-校正（ABM）法和Runge-Kutta（RK）法。首先，把每種積分方法都劃分成五個積分子步$IR1—$IR5，如表1所示。

表1 子步積分法存儲結(jié)構(gòu)Tab.1 The storage structure of Substep integral method

其中$IR1列為每種算法的啟動項，當(dāng)選擇了某種算法，該算法即按$IR1->$IR2->$IR3->$IR4->$IR5->$IR2…依次循環(huán)推進(jìn)，直到仿真結(jié)束。積分子步每向前推進(jìn)一步，則計算一次右函數(shù)。以下是幾種經(jīng)改進(jìn)算法的具體結(jié)構(gòu)設(shè)計。

2.1.1 AB 法

AB法的4個積分子步都是相同的，以四階的AB4法為例：i、ABST（算法啟動）

啟動了AB算法，設(shè)置了幀n-1的導(dǎo)數(shù)值xpa、幀n-2的導(dǎo)數(shù)值xpb和幀n-3的導(dǎo)數(shù)值xpc，使之全部等于右函數(shù)xp的初值。

xl為長字，用來儲存計算結(jié)果，step_time為積分步長。

推進(jìn)導(dǎo)數(shù)值的更新，并由前三幀導(dǎo)數(shù)計算幀n+1的值。

另外三種不同階的AB法設(shè)計思路相同，只是右函數(shù)計算項p->xl有所不同。

2.1.2 ABM法

由于ABM為預(yù)估-校正算法，需要兩幀一輸出，以四階ABM4法為例：

1）ABMST（算法啟動）

啟動了ABM算法，與AB初始條件設(shè)置一樣，由于為兩幀一輸出，故步長為AB法的兩倍。

經(jīng)過預(yù)估-校正，在ABMP4計算結(jié)果，另外三種不同階的ABM法設(shè)計思路也相同。

2.1.3 RK4法

RK法需要計算的值，四幀一輸出，以4階RK法為例：其步長為ABM法的2倍，即為AB法的4倍。

4個子步中的p->xpa分別計算了，然后在RK44中p->xl計算結(jié)果。

2.1.4 小 結(jié)

通過循環(huán)積分子步的設(shè)計，不同的積分法有了相同計算結(jié)構(gòu)，統(tǒng)一了計算量的大小。同時根據(jù)偽代碼的設(shè)計，每個子步都可以通過p->xl輸出計算結(jié)果，大大增強(qiáng)了算法的實時性。

2.2 設(shè)計流程

根據(jù)仿真系統(tǒng)的控制結(jié)構(gòu)的需求分析，編寫算法的流程，如圖1，主要設(shè)置了7個子模塊，分別是：

1）#input（ ）

#input的主要功能是控制程序的輸入，從in.dat文件讀入變量的初值，和步長、幀數(shù)、維數(shù)以及算法的選擇。

2）#ici（ ）

#ici的主要功能是儲存變量的初值xi，并把子步積分控制順序的變量index初始化為1，好進(jìn)行積分子步的推進(jìn)。

圖1 算法流程圖Fig.1 Algorithm flow chart

3）#der（ ）

#der的主要功能有2個，一是儲存方程，二是在程序運行中求導(dǎo)函數(shù)，所以統(tǒng)一的方程格式是xp=···（xp代表變量x的導(dǎo)數(shù)），變量的導(dǎo)數(shù)就等于等號右邊。由于方程不局限于一維，所以用數(shù)組形式 xp[0]=···，xp[1]=···，···，xp[n]=···來儲存方程組。這樣，方程的儲存和讀入，以及變量右函數(shù)求導(dǎo)問題就解決了。

4）#intiri（ ）

#intiri的主要功能是進(jìn)行積分子步的控制，由1，2，3，4，5，2，3，4，5，2，3，4，5…的順序向前推進(jìn)。

5）#int（ ）

#int的主要功能是儲存10種積分算法，在程序需要調(diào)用計算積分的時候，按照需求，提供相應(yīng)的積分子步，完成程序的計算。

6）#output（ ）

#output的主要功能是把輸出變量的積分?jǐn)?shù)值到out.dat文件。

7）#init（ ）

#init的主要功能是初始化分配內(nèi)存。程序的數(shù)據(jù)結(jié)構(gòu)應(yīng)該包括：xl（長字）、xi（初值）、xp（導(dǎo)數(shù)）、xpa（前一幀導(dǎo)數(shù)）、xpb（前二幀導(dǎo)數(shù)）、xpc（前三幀導(dǎo)數(shù)）、method（積分方法），所以我們可以把數(shù)據(jù)結(jié)構(gòu)定義為：

根據(jù)要求，把前6個參數(shù)定義為雙精度，積分算法為單精度就可以。

3 測試對比

MATLAB求解常微分方程組主要應(yīng)用ODE函數(shù)，利用具有非線性擾動的非線性Van Der Pol方程來進(jìn)行測試，

van der pol微分方程[9]：

在步長一定的情況下，把實時仿真得到的結(jié)果與MATLAB的結(jié)果進(jìn)行對比，這里選擇3種代表性的算法AB4、ABM4和RK4法測試。

3.1 仿時真間

由于子步積分算法結(jié)構(gòu)統(tǒng)一，所以在步長和幀數(shù)一致的情況下，不同積分算法的仿真時間是一致的，這是子步積分算法最大的特點。

3.2 誤差分析

我們實驗仿真步長為0.01，幀數(shù)為100，對比3種子步積分算法的輸出結(jié)果與對應(yīng)的MATLAB算法輸出結(jié)果，分析數(shù)據(jù)可知，3組對比結(jié)果存在一定的誤差，且誤差都在穩(wěn)定誤差區(qū)間內(nèi)，根據(jù)算法結(jié)構(gòu)分析，誤差產(chǎn)生的原因與算法精度的設(shè)置有關(guān)。

通過對比數(shù)據(jù)曲線圖，我們發(fā)現(xiàn)圖2中4條數(shù)據(jù)曲線基本擬合；圖3中數(shù)據(jù)曲線每隔一點分離、重合一次，可知分離點為ABM4法的預(yù)估部分，重合點為ABM4法的校正部分；圖4中數(shù)據(jù)曲線每隔三點分離、重合一次，因為RK4法每四步輸出一次結(jié)果。

對輸出數(shù)據(jù)進(jìn)行差值對比，對比方式為AB4法每幀對比，ABM4法隔一幀取值，RK4隔三幀取值，計算3組差值的方差，如表2所示。

可以看出，AB4->ABM4->RK4差值的波動越來越小，趨于穩(wěn)定，誤差逐漸減小，精度變高。實驗結(jié)果證明實時子步積分方法是可行的，得到的結(jié)果與預(yù)期一致，算法的結(jié)構(gòu)設(shè)計比較成功。

圖2 AB4法輸出數(shù)據(jù)對比曲線Fig.2 Output data contrast curve of AB4

圖3 ABM4法輸出數(shù)據(jù)對比曲線Fig.3 Output data contrast curve of ABM4

圖4 RK4法輸出數(shù)據(jù)對比曲線Fig.4 Output data contrast curve of RK4

表2 輸出結(jié)果差值方差Tab.2 Difference values variance of output

4 結(jié) 論

文中分析了幾種常用的積分算法，提出的實時子步積分算法結(jié)構(gòu)，改進(jìn)了幾種常用積分算法，很好地解決了實時仿真中的變步長計算和運算量不一致問題。使得積分算法的選擇不再局限，增強(qiáng)了在工程計算上的擴(kuò)展性。實現(xiàn)了算法的程序編寫，經(jīng)與MATLAB仿真結(jié)果對比，驗證了子步積分算法的精確性，分析了誤差的產(chǎn)生原因。該算法結(jié)構(gòu)為控制領(lǐng)域的工程實踐提供了一個很好的設(shè)計思想，而如何結(jié)合實時仿真控制系統(tǒng)將是下一步研究方向。

[1]劉懷，費樹岷.控制系統(tǒng)中實時任務(wù)的動態(tài)優(yōu)化調(diào)度算法[J].控制與決策，2005，20（3）：246-250.

LIU Huai，F(xiàn)EI Shu-min.Optimal dynamic scheduling algorithm for real-time tasks in digital control systems[J].Control and Decision，2005，20（3）：246-250.

[2]王永炎，王強(qiáng)，王宏安，等.基于優(yōu)先級表的實時調(diào)度算法及其實現(xiàn)[J].軟件學(xué)報，2004，15（3）:360-370.

WANG Yong-yan，WANG Qiang，WANG Hong-an，et al.A real-time scheduling algorithm based on priority table and its implementation[J].Journal of Software，2004，15（3）:360-370.

[3]Marti P，F(xiàn)ohler G.Jitter compensation for real-time control system[J].Proceedings of the 22nd IEEE Real-time System Symposium，2001，19（2）:139-148.

[4]劉德貴，朱珍民.實時仿真算法的研究進(jìn)展[J].系統(tǒng)仿真學(xué)報，2003，24（2）:134-136.

LIU De-gui，ZHU Zhen-min.Research advance for real-time simulation algorithms[J].Acta Simulata Systematica Sinica，2003，24（2）:134-136.

[5]黃振全.實時數(shù)字仿真算法的研究[M].南京：東南大學(xué)2006.

[6]金士堯，黨崗，凌云翔，等.銀河高性能分布仿真系統(tǒng)的設(shè)計和實現(xiàn)[J].計算機(jī)研究與發(fā)展，2001，29（2）:457-459.

JIN Shi-yao，DANG Gang，LING Yun-xing，et al.Design and implementation of YH high performance distribute simulation system[J].Journal of Computer Research and Development，2001，29（2）:457-459.

[7]金宏，王宏安，唐雪梅，等.計算機(jī)控制中的模糊調(diào)度設(shè)計[J].計算機(jī)學(xué)報，2006，29（3）:414-422.

JIN Hong，WANG Hong-an，ANG Xue-mei，etal.Fuzzy scheduling design in computer control[J].Chinese Journal of Computers，2006，29（3）:414-422.

[8]Stankovic J，Ramamritham K.The spring kernel:a new paradigm for real-time systems[J].IEEE Software，1991，8（3）:62-72.

[9]王震，孫衛(wèi).多自由度Vanderpol振子極限環(huán)計算[J].計算機(jī)工程與應(yīng)用，2011，45（2）:85-89.

WANG Zhen，SUN Wei.Computingforlimitcycleof Vanderpoloscillatorwith multidegree offreedom[J].Computer Engineering and Applications，2011，45（2）:85-89.

Design and realization of real-time substep integral algorithm

WANG Yu1， SUN Yong-wei2， SONG Xing-shen3， SHI Yun-long4
（1.Dept.of Military Simulation Technology Institute， Aviation University of AirForce， Changchun 130022， China；2.Dept.of Base of Flight Training， Aviation University of AirForce， Changchun 130062， China；3.Dept.of Aeronautics and Astronautics Intelligence， Aviation University of AirForce， Changchun 130022， China；4.China Aerodynamic Research and Development Center， Mianyang 621000， China）

The integral algorithm is the key point in real-time simulation， and traditional single-step algorithms is simple， but the accuracy is poor.The multi-step high order integral algorithms have high accuracy，but the complicated structure make poor real-time nature.Based on integral algorithms for real-time simulation，this paper put forward the real-time substep integral algorithms，to make different integral algorithms has the same substep control structure.For every substep has the same calculations，different integral algorithms have the same runtime under the same step length and frame number.Experiments show that the substep control structure significantly enhances the timeliness of Multi-step integral algorithms，and the error precision is under control.

computer simulation； timeliness； substep integral algorithms； substep control structure

TP391.9

A

1674－6236（2013）08-0053-05

2012-12-10稿件編號201212059

王 宇（1988—），男，吉林吉林人，碩士研究生，助理工程師。研究方向：虛擬仿真。