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

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

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

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

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

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

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

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

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

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

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

1.1 單步法

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

一般的s級(jí)RK法有如下的形式：

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

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

1.2 線性多步法

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

例如四階Adams顯式（Adams-Bashforth，簡(jiǎn)稱AB法）：

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

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

四階 Adams隱式（Adams-Monlton，簡(jiǎn)稱AM法）：

1.3 單步法與多步法的比較

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

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

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

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

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

2 子步積分算法

2.1 設(shè)計(jì)思想

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

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

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

2.1.1 AB 法

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

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

xl為長(zhǎng)字，用來儲(chǔ)存計(jì)算結(jié)果，step_time為積分步長(zhǎng)。

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

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

2.1.2 ABM法

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

1）ABMST（算法啟動(dòng)）

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

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

2.1.3 RK4法

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

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

2.1.4 小 結(jié)

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

2.2 設(shè)計(jì)流程

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

1）#input（ ）

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

2）#ici（ ）

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

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

3）#der（ ）

#der的主要功能有2個(gè)，一是儲(chǔ)存方程，二是在程序運(yùn)行中求導(dǎo)函數(shù)，所以統(tǒng)一的方程格式是xp=···（xp代表變量x的導(dǎo)數(shù)），變量的導(dǎo)數(shù)就等于等號(hào)右邊。由于方程不局限于一維，所以用數(shù)組形式 xp[0]=···，xp[1]=···，···，xp[n]=···來儲(chǔ)存方程組。這樣，方程的儲(chǔ)存和讀入，以及變量右函數(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的主要功能是儲(chǔ)存10種積分算法，在程序需要調(diào)用計(jì)算積分的時(shí)候，按照需求，提供相應(yīng)的積分子步，完成程序的計(jì)算。

6）#output（ ）

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

7）#init（ ）

#init的主要功能是初始化分配內(nèi)存。程序的數(shù)據(jù)結(jié)構(gòu)應(yīng)該包括：xl（長(zhǎng)字）、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個(gè)參數(shù)定義為雙精度，積分算法為單精度就可以。

3 測(cè)試對(duì)比

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

van der pol微分方程[9]：

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

3.1 仿時(shí)真間

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

3.2 誤差分析

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

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

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

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

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

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

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

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

4 結(jié) 論

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

[1]劉懷，費(fèi)樹岷.控制系統(tǒng)中實(shí)時(shí)任務(wù)的動(dòng)態(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)先級(jí)表的實(shí)時(shí)調(diào)度算法及其實(shí)現(xiàn)[J].軟件學(xué)報(bào)，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]劉德貴，朱珍民.實(shí)時(shí)仿真算法的研究進(jìn)展[J].系統(tǒng)仿真學(xué)報(bào)，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í)時(shí)數(shù)字仿真算法的研究[M].南京：東南大學(xué)2006.

[6]金士堯，黨崗，凌云翔，等.銀河高性能分布仿真系統(tǒng)的設(shè)計(jì)和實(shí)現(xiàn)[J].計(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ì)算機(jī)控制中的模糊調(diào)度設(shè)計(jì)[J].計(jì)算機(jī)學(xué)報(bào)，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].計(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稿件編號(hào)201212059

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