實時子步積分算法的設計與實現

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

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

實時子步積分算法的設計與實現

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

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

積分算法是實時仿真的關鍵所在，傳統的單步法結構簡單，但精度較差。而多步高階積分算法雖然精度高，卻由于結構復雜實時性不強。文章在對實時仿真積分算法研究的基礎上，提出了實時子步積分算法，使不同的積分算法擁有相同的子步控制結構。由于每個子步的運算量一樣，所以在步長和幀數一致的情況下，不同的積分算法運算時間是一致的。經仿真實驗表明，該算法的子步結構顯著提高了多步積分法的實時性，且誤差精度在可控范圍內。

計算機仿真；實時性；子步積分算法；子步控制結構

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

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

1 通用數字仿真算法

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

1.1 單步法

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

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

最常用的Runge-Kutta法為四階RK法，亦稱為經典RK法：

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

1.2 線性多步法

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

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

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

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

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

1.3 單步法與多步法的比較

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

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

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

③當進行變步長運算時，單步法使用靈活，高階多步法相對來說計算復雜，因此在實時仿真中多用到單步法。

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

2 子步積分算法

2.1 設計思想

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

表1 子步積分法存儲結構Tab.1 The storage structure of Substep integral method

其中$IR1列為每種算法的啟動項，當選擇了某種算法，該算法即按$IR1->$IR2->$IR3->$IR4->$IR5->$IR2…依次循環推進，直到仿真結束。積分子步每向前推進一步，則計算一次右函數。以下是幾種經改進算法的具體結構設計。

2.1.1 AB 法

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

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

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

推進導數值的更新，并由前三幀導數計算幀n+1的值。

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

2.1.2 ABM法

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

1）ABMST（算法啟動）

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

經過預估-校正，在ABMP4計算結果，另外三種不同階的ABM法設計思路也相同。

2.1.3 RK4法

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

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

2.1.4 小 結

通過循環積分子步的設計，不同的積分法有了相同計算結構，統一了計算量的大小。同時根據偽代碼的設計，每個子步都可以通過p->xl輸出計算結果，大大增強了算法的實時性。

2.2 設計流程

根據仿真系統的控制結構的需求分析，編寫算法的流程，如圖1，主要設置了7個子模塊，分別是：

1）#input（ ）

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

2）#ici（ ）

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

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

3）#der（ ）

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

4）#intiri（ ）

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

5）#int（ ）

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

6）#output（ ）

#output的主要功能是把輸出變量的積分數值到out.dat文件。

7）#init（ ）

#init的主要功能是初始化分配內存。程序的數據結構應該包括：xl（長字）、xi（初值）、xp（導數）、xpa（前一幀導數）、xpb（前二幀導數）、xpc（前三幀導數）、method（積分方法），所以我們可以把數據結構定義為：

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

3 測試對比

MATLAB求解常微分方程組主要應用ODE函數，利用具有非線性擾動的非線性Van Der Pol方程來進行測試，

van der pol微分方程[9]：

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

3.1 仿時真間

由于子步積分算法結構統一，所以在步長和幀數一致的情況下，不同積分算法的仿真時間是一致的，這是子步積分算法最大的特點。

3.2 誤差分析

我們實驗仿真步長為0.01，幀數為100，對比3種子步積分算法的輸出結果與對應的MATLAB算法輸出結果，分析數據可知，3組對比結果存在一定的誤差，且誤差都在穩定誤差區間內，根據算法結構分析，誤差產生的原因與算法精度的設置有關。

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

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

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

圖2 AB4法輸出數據對比曲線Fig.2 Output data contrast curve of AB4

圖3 ABM4法輸出數據對比曲線Fig.3 Output data contrast curve of ABM4

圖4 RK4法輸出數據對比曲線Fig.4 Output data contrast curve of RK4

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

4 結 論

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

[1]劉懷，費樹岷.控制系統中實時任務的動態優化調度算法[J].控制與決策，2005，20（3）：246-250.

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

[2]王永炎，王強，王宏安，等.基于優先級表的實時調度算法及其實現[J].軟件學報，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，Fohler 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].系統仿真學報，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]黃振全.實時數字仿真算法的研究[M].南京：東南大學2006.

[6]金士堯，黨崗，凌云翔，等.銀河高性能分布仿真系統的設計和實現[J].計算機研究與發展，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].計算機學報，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]王震，孫衛.多自由度Vanderpol振子極限環計算[J].計算機工程與應用，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—），男，吉林吉林人，碩士研究生，助理工程師。研究方向：虛擬仿真。