熊 濤，丁建完，陳立平

（華中科技大學(xué)國家CAD支撐軟件工程技術(shù)研究中心，湖北 武漢 430074）

面向半物理仿真的陳述式模型求解方法研究*

熊 濤，丁建完，陳立平

（華中科技大學(xué)國家CAD支撐軟件工程技術(shù)研究中心，湖北 武漢 430074）

通過符號操作和數(shù)值計算相結(jié)合，提出了一種求解半物理仿真模型的新方法。為了滿足半物理仿真對實(shí)時性的要求，在模型編譯階段將代表數(shù)值積分的隱式離散公式插入到仿真模型中，增廣后的方程系統(tǒng)伴隨著非線性方程的出現(xiàn)，需要在積分的每一步對這些非線性方程進(jìn)行迭代求解，而求解非線性方程的時間復(fù)雜度隨維度的變大成指數(shù)增加，因此引入代數(shù)環(huán)撕裂減小代數(shù)方程塊耦合變量數(shù)，以滿足實(shí)時求解對粒度的要求。最后通過實(shí)例對文中提出的方法進(jìn)行了驗(yàn)證。

半物理仿真；微分代數(shù)離散；代數(shù)環(huán)撕裂；實(shí)時積分

1 引言

仿真技術(shù)是以相似原理、信息技術(shù)、系統(tǒng)技術(shù)及其應(yīng)用領(lǐng)域有關(guān)的專業(yè)技術(shù)為基礎(chǔ)，以計算機(jī)和各種物理設(shè)備為工具，利用系統(tǒng)模型對實(shí)際的或設(shè)想的系統(tǒng)進(jìn)行試驗(yàn)研究的一門綜合性技術(shù)。仿真技術(shù)應(yīng)用廣泛，從早期的航空、航天、火力發(fā)電和核動力發(fā)電部門擴(kuò)展到今天的軍事、電子、通信、交通、冶金、建筑、氣象、機(jī)械、輕工等多種行業(yè)和部門，其應(yīng)用已滲透到系統(tǒng)生命周期的全過程［1］。半物理仿真作為工程領(lǐng)域內(nèi)一種應(yīng)用廣泛的仿真技術(shù)，涉及機(jī)電液控多個領(lǐng)域和學(xué)科?！鞍胛锢矸抡妗庇址Q硬件在環(huán)仿真或物理—數(shù)學(xué)仿真［2］，是指針對仿真研究內(nèi)容，將控制系統(tǒng)、被控制對象或系統(tǒng)部件用數(shù)學(xué)模型來描述，并轉(zhuǎn)化為高速運(yùn)行的實(shí)時仿真模型，然后與實(shí)物連接成為一個集計算機(jī)技術(shù)、數(shù)學(xué)建模技術(shù)及數(shù)據(jù)采集技術(shù)于一體的半物理仿真系統(tǒng)。然后借助物理效應(yīng)模型，進(jìn)行數(shù)學(xué)與物理的實(shí)時聯(lián)合仿真，通過仿真試驗(yàn)對控制系統(tǒng)的控制策略、控制功能以及系統(tǒng)可靠性等進(jìn)行測試和評估［3］。目前，半物理仿真已被廣泛地應(yīng)用于軍事、汽車制造、工程機(jī)械和高校實(shí)驗(yàn)室建設(shè)領(lǐng)域。發(fā)展勢頭迅猛的多領(lǐng)域物理系統(tǒng)實(shí)時仿真主要應(yīng)用方式就為半物理仿真，如減速箱［4］、柴油機(jī)控制［5］、汽車［6］、無 人汽 車［7］、鋰 電 池［8］和 航 電模擬［9］等 。

本文對基于Modelica語言面向半物理仿真的陳述式模型實(shí)時求解算法進(jìn)行研究，組織結(jié)構(gòu)如下，首先對半物理仿真及其在工程實(shí)際中的應(yīng)用進(jìn)行簡單介紹；其次對半物理模型的數(shù)學(xué)表達(dá)—微分代數(shù)方程DAEs（Differential Algebraic Equations）及其常用的求解方法進(jìn)行分析，給出半物理仿真的難點(diǎn)；然后針對求解中的困難提出在模型編譯階段對DAEs進(jìn)行離散，并對離散后能夠縮減問題規(guī)模的代數(shù)環(huán)撕裂進(jìn)行討論；最后給出該方法的優(yōu)勢和展望。

2 微分方程求解方法

實(shí)時仿真在仿真軟件應(yīng)用領(lǐng)域發(fā)展的目標(biāo)之一是能夠以快速的采樣率實(shí)時模擬越來越復(fù)雜的模型，滿足半物理仿真的特性，也就是在很短的時間間隔內(nèi)對物理模型的數(shù)學(xué)描述進(jìn)行求解。工程中通常采用DAEs描述物理系統(tǒng)的動態(tài)特性和行為，DAEs通常由相互耦合的微分方程和代數(shù)方程混合而成，由于代數(shù)約束的存在，使DAEs的求解與常微分方程ODEs（Ordinary Differential Equations）產(chǎn)生了極大的不同，隨著一系列特殊求解方法的提出，DAEs在工程實(shí)踐中得到了廣泛應(yīng)用。其中，指標(biāo)［10］是DAEs的一個重要屬性，文獻(xiàn)［11］對低指標(biāo)問題進(jìn)行了詳細(xì)研究，高指標(biāo)問題可以通過指標(biāo)縮減轉(zhuǎn)化為ODEs，指標(biāo)縮減代表性的方法有：Gear方法［12，13］、Pantelides方法［14］和啞導(dǎo)方法［15，16］等。

微分方 程 的 求 解 方 法 分 為 變 步 長［17］和 定步長［18］兩種，其中變步長方法需要在每個積分步對收斂性進(jìn)行判斷，且必須考慮求解精度對于步長的影響，效率較低?！坝布诃h(huán)”是半物理仿真的顯著特點(diǎn)，由于回路中接入實(shí)物，因此仿真是實(shí)時進(jìn)行的，即仿真模型的時間比例尺和自然時間比例尺相同，同時受交互硬件固有頻率的影響，使變步長算法在半物理仿真中的應(yīng)用很困難。實(shí)際操作中半物理仿真往往使用定步長的方法，定步長算法又分為顯式和隱式兩種，對于常用的顯式方法，當(dāng)模型由機(jī)械、電氣、液壓或者熱力組件耦合而成時，往往兼有快速和慢速子系統(tǒng)，導(dǎo)致仿真變成剛性問題，也就是說，模型中的時間常數(shù)具有較大的跨度，較快的時間常數(shù)決定仿真的計算量（步長），這將帶來過多的計算量，降低了求解效率。使用隱式的求解算法雖然解決了數(shù)值穩(wěn)定性問題，并且允許使用較大的步長，但求解精度制約著步長，而且使用隱式方法意味著由微分方程離散而來的非線性方程需要在每一積分步迭代求解，實(shí)時求解規(guī)模為n的非線性方程系統(tǒng)存在一定問題：操作的時間復(fù)雜度是O（n3），并且問題達(dá)到收斂所需的迭代次數(shù)可能隨不同的步長而變動。綜上，如何保證半物理仿真的實(shí)時性便成為一個重要的課題，除去仿真計算機(jī)自身的速度和精度，仿真模型的求解算法是影響仿真速度的一大重要因素，本文給出一種符號操作與數(shù)值計算相結(jié)合的半物理模型求解方法。該方法的基本思想為，使用代表求解方法的隱式公式在求解之前對DAEs進(jìn)行離散，并借助相關(guān)符號操作實(shí)現(xiàn)求解的規(guī)劃與分解，達(dá)到問題分而治之的目的。

3 DAEs的離散

在模型編譯階段將代表數(shù)值積分運(yùn)算的隱式歐拉離散公式插入到方程系統(tǒng)中，實(shí)現(xiàn)指標(biāo)約減后DAEs的增廣，從而將連續(xù)問題預(yù)先轉(zhuǎn)化為離散問題。完整約束的半仿真模型在指標(biāo)約減后具有下面的表現(xiàn)形式：

其中，E為微分方程，Φ為代數(shù)約束方程，x為狀態(tài)變量，y為代數(shù)變量，假設(shè)在一致性初始條件x（t0）= x0，˙x（t0）=x0，y（t0）=y0下，具有唯一且光滑的解［17］。

DAEs離散步驟：

步驟1 根據(jù)是否具有微分形式將變量劃分為集合x和集合y，滿足dim（W）=dim（x）+dim （y），W=［x，y］。

步驟2 為保證方程系統(tǒng)的恰定性，取x為已知量，˙x和y為未知量，對DAEs進(jìn)行強(qiáng)連通分量SCC（Strongly Connected Components）凝聚［19］，以及拓?fù)渑判颍瑥亩梢缘玫巾樞蚯蠼獾姆匠套蛹蛄小?/p>

步驟3 就xi而言，若˙xi可通過匹配的方程顯式地求解出來，則向DAEs中添加如下的離散方程：

其中res（xi）表示將xi當(dāng)作撕裂變量，以迭代的方式求解xi。

若˙xi無法顯式求解，向DAEs中添加：

步驟4 若與˙xj或者yk匹配的賦值形式方程en處于代數(shù)環(huán)中，即結(jié)構(gòu)關(guān)聯(lián)矩陣中所處的對角塊維度超過1，則在en的表達(dá)式中添加一個新項(xiàng)res（˙xi）或者res（yk），標(biāo)記為撕裂變量，即：

步驟5 如果res（xi）、res（˙xj）或是res（yk）所處的代數(shù)環(huán)不能被撕裂，則移除相應(yīng)的res（）操作符。

步驟6 以x、˙x和y為未知量，重新對離散后的DAEs進(jìn)行SCC凝聚和拓?fù)渑判?，從而可以得到順序求解的方程子集序列?/p>

下面通過一個簡單的例子進(jìn)行說明。

其中，ω、D和k為常數(shù)，f和g為函數(shù)，x、x1和x2為 微 分 變 量，u和y為 代 數(shù) 變 量。若 使 用DASSL［17］直接求解，則需對全部的5個變量進(jìn)行迭代。DAEs離散方法處理如下：由于微分變量˙x、˙x1和˙x2均具有顯式表達(dá)，所以使用下面的離散公式（8）～公式（10）對式（7）進(jìn)行增廣：

增廣后的DAEs如下：

就式（11）而言，為保證方程的恰定性需同時將˙x、˙x1、˙x2、x、x1和x2視為未知量，消去增廣式中的微分變量˙x、˙x1和˙x2，得：

式（12）的二部圖最大匹配及SCC凝聚如圖1和圖2所示，圖1中雙向邊代表最大匹配，其中需要注意的是式f4中的x1，可以通過式f3轉(zhuǎn)換為關(guān)于x2的表達(dá)。

Figure 1 Direct constraint graphˉG of equ.（12）圖1 式（12）的約束有向圖ˉG

Figure 2 Cohesion graphˉC of SCCˉG圖2 ˉG強(qiáng)連通分量凝聚圖ˉC

圖2中，強(qiáng)連通分量C1=｛f1|x｝、C2=｛f2|y｝、C3=｛f3|x1｝、C4=｛f4|x2｝和C5=｛f5| u｝，這5個強(qiáng)連通分量構(gòu)成了代數(shù)環(huán)，當(dāng)選擇變量x作為撕裂變量，反映在圖中為切斷強(qiáng)連通分量C1和C2之間的有向弧，可得順序求解的方程序列C2→C4→C3→C5→C1，方程求解順序與圖ˉC中有向弧的方向相反。

觀察式（13）可知，只需要對變量x進(jìn)行迭代，即可實(shí)現(xiàn)方程的求解，所以移除標(biāo)記res（x1）和res（x2）。同理，可選擇變量x1或x2作為撕裂變量。下面介紹對于縮減問題規(guī)模具有重要意義的代數(shù)環(huán)撕裂。

4 代數(shù)環(huán)撕裂

為了減小DAEs離散后非線性求解塊的粒度，使用代數(shù)環(huán)撕裂作為后處理，以提高求解效率。代數(shù)環(huán)撕裂的基本思想：利用代數(shù)環(huán)結(jié)構(gòu)上的稀疏性，通過切斷某些耦合關(guān)系，將某些方程轉(zhuǎn)化為賦值等式的形式，從而縮減需要聯(lián)立求解部分的規(guī)模，提高求解效率?？紤]一般形式的方程組，如式（14）所示：

式（14）為三個相互耦合的方程，約束表示圖如圖3a所示，這三個方程組成了一個強(qiáng)連通分量。撕裂該代數(shù)環(huán)的符號操作為：打斷組成強(qiáng)連通分量中的某條邊，檢查圖中是否存在新的強(qiáng)連通分量，若存在，繼續(xù)上述的操作；否則，結(jié)束。需要注意的是，每打斷約束表示圖中的一條邊，為保證方程系統(tǒng)的恰定性，需要引入一個輔助方程和輔助變量，如圖3b所示。

Figure 3 Direct graph of equ.（14）圖3 式（14）的有向表示圖

依次打斷邊（f1，x1）和（f2，x2），得到圖3b所示的形式，方程f'1和變量x'1是打斷邊（f1，x1）時引入的輔助方程和輔助變量，f'2和x'2是打斷邊（f2，x2）時引入的輔助方程和輔助變量，符號操作后的式（14）可以通過下面的計算過程求得：

BEGIN

LOOP

UNTIL|x'1—x1|＜εAND|x'2—x2|＜ε

END

fi_calc_xi表示通過方程fi計算變量xi，ε表示迭代收斂所允許的誤差。

下面給出代數(shù)環(huán)撕裂的一般方法，假設(shè)代數(shù)環(huán)具有式（15）指定的形式：

撕裂就是從代數(shù)環(huán)中分離出一部分方程和變量，即｛g1|g1∈g｝和｛z1|z1∈z｝。假定z1為已知量，可由剩余的方程｛g2|g2∪g1=g｝，依次求解剩余的變量｛z2|z2∪z1=z｝?？梢酝ㄟ^牛頓法求解分裂后的方程系統(tǒng)，該過程描述如下：

步驟1 選擇合適的撕裂變量z1，并賦初值z1（0）；

步驟2 通過剩余的方程g2求出z2，即z2= g2（z1）；

步驟3 計算殘差δ=g1（z1，z2）；

步驟4 若δ小于收斂誤差ε，結(jié)束，否則，轉(zhuǎn)步驟2。

通過代數(shù)環(huán)撕裂，可以將耦合方程的維數(shù)從dim（h1）+dim（h2）減小到dim（h1）。然而就撕裂過程而言，選出合適的z1是十分困難的，不同的選取方式會帶來計算量和誤差上的差異，同時在實(shí)際操作中必須考慮初值z1（0）以及h1（z1，z2）迭代的收斂性。為此，采用的方法是取DAEs離散后標(biāo)記為res的變量作為撕裂變量，來驗(yàn)證能否達(dá)到代數(shù)環(huán)撕裂的目的。

5 實(shí)例分析

本文提出的方法在多領(lǐng)域建模與仿真平臺MWorks［20］中進(jìn)行實(shí)驗(yàn)，下面通過笛卡爾坐標(biāo)系下三自由度挖掘機(jī)的轉(zhuǎn)臂模型進(jìn)行驗(yàn)證。

轉(zhuǎn)臂由包括R1、R2、R3、R6、R7、R8、R9、R11、R12、R13、R14、R15和P4、P5、P10在內(nèi)的12個轉(zhuǎn)動副和3個移動副組成，為使存在環(huán)路的模型可解，取R6、R8、R11、R14為切割鉸。記r表示R. phi，用以描述轉(zhuǎn)動副的相位轉(zhuǎn)角；p表示P.s，描述移動副的伸張長度。指標(biāo)約減后選擇r7、r12、r15作為狀態(tài)變量的DAEs具有下面的表現(xiàn)形式：

Figure 4 Excavator jib model with three degrees of freedom圖4 三自由度挖掘機(jī)的轉(zhuǎn)臂模型

Figure 5 Directed constraint graph of augmented DAEs圖5 轉(zhuǎn)臂模型增廣DAEs的有向約束圖h

Figure 6 Cohesion graphˉChof SCCh圖6的強(qiáng)連通分量凝聚圖

由圖6可得各強(qiáng)連通分量中需要聯(lián)立求解的變量數(shù)分別為5、3和3，求解順序?yàn)椤?，將各?qiáng)連通分量中標(biāo)記為res的變量作為撕裂變量執(zhí)行代數(shù)環(huán)撕裂，均可以獲得迭代變量數(shù)為1的求解序列，如下所示：

使用隱式歐拉法對該模型進(jìn)行求解，對比一般方法（需對11個變量同時進(jìn)行迭代）和文中提出的方法，仿真時間均為10 s，表1給出了不同方式下積分時間的對比。

Table 1 Comparison of different methods in solving time表1 不同方式求解時間對比

采用DAEs離散和代數(shù)環(huán)撕裂的方法能夠滿足實(shí)時性的要求（7.35 s＜10 s），同時為了保證提出方法的正確性，圖7給出r2的仿真結(jié)果對比。

Figure 7 Simulation results of r2in two ways圖7 兩種方式下r2仿真結(jié)果對比

就轉(zhuǎn)臂模型而言，借助DAEs離散和代數(shù)環(huán)撕裂將需要聯(lián)立求解的模型（迭代變量數(shù)為11）轉(zhuǎn)化為三個代數(shù)方程塊的順序求解，且每個方程塊的迭代變量數(shù)為1，并且數(shù)值積分時間小于模型仿真時間，能夠滿足半物理仿真的要求。

6 結(jié)束語

本文主張在模型編譯階段采用隱式歐拉法對DAEs進(jìn)行離散，并通過代數(shù)環(huán)撕裂來減小增廣后的方程系統(tǒng)規(guī)模。絕大多數(shù)情況下本文提出的方法能夠?qū)⒄麄€方程系統(tǒng)分解為多個可單獨(dú)迭代求解的子系統(tǒng)，從而提高求解效率，滿足半物理仿真的要求，并通過實(shí)例分析驗(yàn)證了該方法的正確性?？偟膩碚f，只有在離散后的代數(shù)環(huán)無法撕裂的情況下，所需迭代的變量數(shù)才與通用積分方法相等。在未來的研究中希望能夠?qū)Ψ抡婺Ｐ椭械目焖贍顟B(tài)與慢速狀態(tài)進(jìn)行區(qū)分，從而能夠在DAEs的離散中使用顯式的離散公式。

［1］ Liu Rui-ye.Elements of computer simulation［M］.Beijing：Publishing House of Electronics Industry，2011.（in Chinese）

［2］ Aranya C，Xin Yu-feng.Hardware-in-the-loop simulations and verifications of smart power systems over an exo-GENI testbed［C］∥Proc of the 2nd GENI Research and Educational Experiment Workshop，2013：16-19.

［3］ Kang Feng-ju，Yang Hui-zhen，Gao Li-e.Technology and application of modern simulation［M］.Beijing：National Defense Industry Press，2006.（in Chinese）

［4］ Casella F，Donida F，Lovera M.Beyond simulation：Computer aided control system design using equation-based object oriented modeling for the next decade［C］∥Proc of the 2nd International Workshop on Equation-Based Object-Oriented Languages and Tools，2008：35-45.

［5］ Song Bai-ling.Semi-physical simulation technology research of diesel electronic control system［D］.Harbin：Harbin Engineering University，2009.（in Chinese）

［6］ Mattsson S E，Elmqvist H，Otter M.Physical system modeling with Modelica［J］.Control Engineering Practice，1998，6 （4）：501-510.

［7］ Tang Guo-ming.Design of hardware-in-the-loop simulation system for driverless vehicle［D］.Hefei：University of Science and Technology of China，2012.（in Chinese）

［8］ Luciano S，Ines C，Manuela G.A design methodology for semi-physical fuzzy models applied to the dynamic characterization of LiFePO4 batteries［J］.Applied Soft Computing，2014，14（B）：269-288.

［9］ Jiang Zhan-si.Parameter optimization and inference solving for models based on Modelica［D］.Wuhan：Huazhong University of Science&Technology，2008.（in Chinese）

［10］ Unger J，Kr?ner A，Marquardt W.Structural analysis of differential-algebraic equation systems-theory and applications ［J］.Computers&Chemical Engineering，1995，19（8）：867-882.

［11］ Negrut D，Jay L O，Khude N.A discussion of low order numerical integration formulas for rigid and flexible multibody dynamics［J］.ASME Journal of Computational and Nonlinear Dynamics，2007，4（2）：1-11.

［12］ Gear C W，Petzold L R.ODE methods for the solution of differential/algebraic systems［J］.SIAM Journal on Numerical Analysis，1982，21（4）：716-728.

［13］ Gear C W.Differential-algebraic equation index transformations［J］.SIAM Journal on Scientific and Statistical Computing，1988，9（1）：39-47.

［14］ Pantelides C C.The consistent initialization of differentialalgebraic systems［J］.SIAM Journal on Scientific and Statistical Computing，1988，9（2）：213-231.

［15］ Mattsson S E，S?derlind G.Index reduction in differentialalgebraic equations using dummy derivatives［J］.SIAM Journal on Scientific Computing，1993，14（3）：677-692.

［16］ Feehery W F，Barton P I.A differentiation-based approach to dynamic simulation and optimization with high-index differential-algebraic equations［C］∥Proc of the 2nd International Workshop on Computational Differentiation，1996：239-252.

［17］ Brenan K E，Campbell S L，Petzold L R.Numerical solution of initial-value problems in differential-algebraic equations ［M］.Philadelphia：SIAM，1996.

［18］ Mathews J H，F(xiàn)ink K K.Numerical methods using Matlab ［M］.London：Prentice Hall，2004.

［19］ Xavier A.On the complexity of strongly connected components in directed hypergraphs［J］.Algorithmica，2014，69 （2）：335-369.

［20］ http：∥www.tongyuan.cc/.

附中文參考文獻(xiàn)：

［1］ 劉瑞葉.計算機(jī)仿真技術(shù)基礎(chǔ)［M］.北京：電子工業(yè)出版社，2011.

［3］ 康鳳舉，楊惠珍，高立娥.現(xiàn)代仿真技術(shù)與應(yīng)用［M］.北京：國防工業(yè)出版社，2006.

［5］ 宋百齡.柴油機(jī)控制系統(tǒng)半物理仿真技術(shù)研究［D］.哈爾濱：哈爾濱工業(yè)大學(xué)，2009.

［7］ 唐國明.無人駕駛汽車半物理仿真系統(tǒng)設(shè)計［D］.合肥：中國科學(xué)技術(shù)大學(xué)，2012.

［9］ 蔣占四.Modelica仿真模型的參數(shù)優(yōu)化及推理求解研究［D］.武漢：華中科技大學(xué)，2008.

熊濤（1990），男，河南確山人，博士生，研究方向?yàn)槎囿w系統(tǒng)動力學(xué)、多領(lǐng)域建模與求解。E-mail：xiongtao39@163.com

XIONG Tao，born in 1990，Ph D candidate，his research interests include dynamics of multibody systems，multi-domain modeling and solving.

A method for solving semi-physical simulation declarative models

XIONG Tao，DING Jian-wan，CHEN Li-ping
（National CAD Support Software Engineering Research Center，Huazhong University of Science and Technology，Wuhan 430074，China）

We propose a new method for solving semi-physical simulation models using a mixed symbolic and numeric approach.In order to meet the real-time requirement of semi-physical simulations，implicit discretization formulae representing the numerical integration algorithm are inserted into the DAEs symbolically at compile stage.Then nonlinear equations will appear in the augmented equation system and all these nonlinear equations should be solved together at each integration step.In order to meet the fine granularity required for solving the models in real-time，tearing algebraic loop is introduced.After that the dimensions of nonlinear equation blocks can be reduced as the time complexity of calculating nonlinear equations increases exponentially with the growth of dimensions.Finally，an example is given to show that the proposed method is not only easy to implement but also efficient.

semi-physical simulation；discrete DAEs；tearing algebraic loops；real-time integration

TP301

A

10.3969/j.issn.1007-130X.2015.08.018

1007-130X（2015）08-1540-06

2014-06-27；

2014-10-11

國家科技支撐計劃資助項(xiàng)目（2012BAF16G02）

通信地址：430074湖北省武漢市華中科技大學(xué)國家CAD支撐軟件工程技術(shù)研究中心

Address：National CAD Support Software Engineering Research Center，Huazhong University of Science and Technology，Wuhan

430074，Hubei，P.R.China