艾鳳明 馮首鴻 梁興壯 李偉林 林秦州

(1. 中國(guó)航空工業(yè)集團(tuán)公司沈陽(yáng)飛機(jī)設(shè)計(jì)研究所 沈陽(yáng) 110035； 2. 西北工業(yè)大學(xué)自動(dòng)化學(xué)院 西安 710129)

1 引言

復(fù)雜系統(tǒng)的設(shè)計(jì)已經(jīng)證實(shí)了仿真是從概念設(shè)計(jì)到樣品實(shí)現(xiàn)過(guò)程中不可或缺的步驟。事物仿真需要考慮初步估計(jì)、擇優(yōu)甚至是現(xiàn)有方案在實(shí)施前的重新設(shè)計(jì)，因此會(huì)降低風(fēng)險(xiǎn)。建立元器件本身及其之間交互作用的高保真模型，能夠使仿真的結(jié)果更具有確信度。

聯(lián)合仿真是現(xiàn)階段最有效的耦合系統(tǒng)仿真方案。從數(shù)學(xué)方面分析，聯(lián)合仿真由一組常微分方程和微分代數(shù)方程構(gòu)成，而仿真時(shí)間包括子系統(tǒng)之間的通信節(jié)點(diǎn)時(shí)間和每一個(gè)子系統(tǒng)內(nèi)部的積分時(shí)間。分析一個(gè)聯(lián)合仿真系統(tǒng)是否符合要求需要對(duì)其穩(wěn)定性、準(zhǔn)確性和快速性進(jìn)行分析。

聯(lián)合仿真系統(tǒng)中，每一個(gè)仿真器都可以看成是一個(gè)動(dòng)態(tài)系統(tǒng)，它們之間具有信號(hào)的傳遞，傳遞的信號(hào)與上一次信號(hào)傳遞的時(shí)間間隔為一個(gè)宏步長(zhǎng)，該步長(zhǎng)的大小反映了系統(tǒng)仿真計(jì)算時(shí)間的長(zhǎng)短。通常情況下，仿真過(guò)程結(jié)果都希望具有最小的通信時(shí)間，然而考慮到系統(tǒng)穩(wěn)定性，該計(jì)算時(shí)間不能無(wú)限變小，所以需要找出兼顧計(jì)算效率和計(jì)算穩(wěn)定性的方法來(lái)滿足系統(tǒng)在仿真需求。為了提高仿真速度，聯(lián)合仿真的并行化速度也需要提高，然而并行的子模型和相關(guān)求解器之間的松耦合會(huì)引起積分誤差，為了保證誤差在預(yù)先要求的范圍內(nèi)，同樣需要兼顧考慮仿真效率和仿真精確度的關(guān)系。

Dymola是一款基于Modelica語(yǔ)言的多領(lǐng)域仿真平臺(tái)。Modelica語(yǔ)言可以實(shí)現(xiàn)各模塊庫(kù)之間的連接，進(jìn)行多領(lǐng)域仿真。文獻(xiàn)[1]通過(guò)Modelica語(yǔ)言搭建了包括交流發(fā)電機(jī)、變壓整流器、電源系統(tǒng)接觸器和過(guò)欠壓、過(guò)欠頻檢測(cè)模塊等部件的數(shù)字仿真模型和供電系統(tǒng)模型，并完成系統(tǒng)邏輯設(shè)計(jì)。通過(guò)特殊的接口，Dymola可以與多款軟件實(shí)現(xiàn)聯(lián)合仿真，現(xiàn)在國(guó)內(nèi)外的主要研究方向?yàn)閷?duì)電機(jī)的控制和汽車系統(tǒng)建模。文獻(xiàn)[2]利用Simulink、CarSim、Dymola軟件實(shí)現(xiàn)對(duì)負(fù)荷傳感液壓轉(zhuǎn)向系統(tǒng)的聯(lián)合仿真；Dymola符合FMI接口協(xié)議，所以可以輸入輸出FMU模型；文獻(xiàn)[3]在Dymola中實(shí)現(xiàn)二極管模型的FMU模型輸出；文獻(xiàn)[4]在Dymola中實(shí)現(xiàn)熱流體模型的FMU模型輸出；文獻(xiàn)[5]是對(duì)自動(dòng)嵌入式物理模型的FMI研究；文獻(xiàn)[6]將Modelica模型通過(guò)FMI輸出到Excel表格中，通過(guò)FMI Toolbox for Matlab工具生成FMU模型，最后導(dǎo)入Dymola平臺(tái)中進(jìn)行聯(lián)合仿真，實(shí)現(xiàn)控制器算法模型和Dymola整車聯(lián)合仿真測(cè)試。

本文對(duì)聯(lián)合仿真系統(tǒng)模型進(jìn)行了數(shù)學(xué)分析，建立了系統(tǒng)微分方程，在穩(wěn)定和精度要求下提出了線性隱性計(jì)算方法，并提出了基于FMI協(xié)議的聯(lián)合仿真計(jì)算流程，最后通過(guò)一個(gè)工程例子說(shuō)明了該方法的實(shí)際應(yīng)用效果。

2 聯(lián)合仿真系統(tǒng)模型分析

用非線性微分方程描述一個(gè)連續(xù)離散混合動(dòng)態(tài)系統(tǒng)Σ形式如下

圖1 被分割并行的系統(tǒng)

因此系統(tǒng)Σ可以重新寫成

式中，U1是子系統(tǒng) ∑1的輸入，U2是子系統(tǒng) ∑2的輸入。也就是說(shuō)，根據(jù)得到的解耦結(jié)果判斷是不是空集。同樣，1Y是子系統(tǒng) 1∑得到的輸出，Y2是子系統(tǒng) ∑2得到的輸出，即

為了更好數(shù)值運(yùn)算整個(gè)復(fù)雜系統(tǒng)，每一個(gè)仿真器都需要在通信節(jié)點(diǎn)進(jìn)行數(shù)據(jù)交換，來(lái)為其他子模塊提供所需要的數(shù)據(jù)。為了提高積分運(yùn)算速度，并行模塊要盡量獨(dú)立，所以子模型之間的同步時(shí)間P要遠(yuǎn)大于他們內(nèi)部的積分時(shí)間。因此在通信節(jié)點(diǎn)之間，每一個(gè)仿真器都有自己的積分速率(假設(shè)使用變步長(zhǎng)求解器)，并且在此期間，從其他模型處得到的數(shù)據(jù)將保持常數(shù)。

大的通信間隔很可能提高積分運(yùn)算速度，但是也可能會(huì)產(chǎn)生積分誤差累計(jì)從而降低最終結(jié)果的精確度。所以通過(guò)非嚴(yán)格同步來(lái)減小建模誤差是找到提高積分速度和精確度之間平衡的有效方式的第 一步。

圖2 進(jìn)行數(shù)據(jù)交換的子模型∑1和∑2

式中，ωmax為系統(tǒng)最大角頻率，滿足

所以系統(tǒng)的宏步長(zhǎng)滿足

然而這個(gè)宏步長(zhǎng)的邊界太大，由經(jīng)驗(yàn)求得，宏步長(zhǎng)的頻率要遠(yuǎn)大于Nyquist頻率，有時(shí)會(huì)達(dá)到100倍。如果宏步長(zhǎng)的頻率不夠大，則會(huì)導(dǎo)致系統(tǒng)的最終發(fā)散。為了解釋只有很小的宏步長(zhǎng)才能使耦合系統(tǒng)穩(wěn)定，需要對(duì)整個(gè)系統(tǒng)的穩(wěn)定性進(jìn)行分析，按照文獻(xiàn)[8]中所述的方法可以對(duì)這類環(huán)路采樣系統(tǒng)進(jìn)行穩(wěn)定性分析。此時(shí)，子系統(tǒng)通過(guò)零階保持器來(lái)耦合變量。

圖3為兩個(gè)子系統(tǒng)的聯(lián)合仿真方框圖，其主要在每一步離散化了耦合變量，這模擬了子系統(tǒng)可以分別使用數(shù)值求解器進(jìn)行內(nèi)部計(jì)算。對(duì)于兩個(gè)系統(tǒng)的耦合連接，通過(guò)線性化狀態(tài)空間所求出的離散時(shí)間傳遞函數(shù)來(lái)分析穩(wěn)定性[9]。

圖3 具有兩個(gè)子系統(tǒng)的聯(lián)合仿真方框圖

3 FMI結(jié)構(gòu)中的實(shí)現(xiàn)過(guò)程

若要保證系統(tǒng)穩(wěn)定，需要足夠小的宏步長(zhǎng)，這雖然保證了系統(tǒng)的精確度，但卻需要消耗過(guò)多的計(jì)算時(shí)間，降低了計(jì)算效率[10]。線性隱性穩(wěn)定方法利用子系統(tǒng)的Jacobian矩陣來(lái)線性化系統(tǒng)，這種方法允許使用相對(duì)較大的宏步長(zhǎng)來(lái)實(shí)現(xiàn)系統(tǒng)的穩(wěn)定[11]。

整個(gè)耦合系統(tǒng)的數(shù)學(xué)模型如下所示

對(duì)于整個(gè)系統(tǒng)，內(nèi)部子系統(tǒng)之間的連接是通過(guò)耦合方程來(lái)實(shí)現(xiàn)的。

在耦合方程中，每一個(gè)子系統(tǒng)的輸入向量是其他子系統(tǒng)輸出向量的函數(shù)，K為具有以下特點(diǎn)的關(guān)聯(lián)矩陣：①K是一個(gè)方陣，認(rèn)為每一個(gè)子系統(tǒng)的輸出都對(duì)應(yīng)一個(gè)子系統(tǒng)的輸入；② 只包含1或者0；③ 矩陣中每一行和每一列中只有一個(gè)1。

對(duì)系統(tǒng)進(jìn)行線性化處理，規(guī)定Jacobian矩陣為

在使用隱性穩(wěn)定方法時(shí)，需要以下假設(shè)：① 在沒(méi)有考慮耦合變量的代數(shù)環(huán)情況下，假設(shè)DK的乘積為冪零；② 在聯(lián)合仿真宏步長(zhǎng)中，線性時(shí)不變系統(tǒng)由線性化(在點(diǎn)和系統(tǒng)常微分方程得到；③ 通過(guò)Backward Euler方法來(lái)離散線性化后的系統(tǒng)。

根據(jù)假設(shè)②，線性化后的系統(tǒng)數(shù)學(xué)模型如下 所示

式中，ξ、η和ω分別在線性化系統(tǒng)中對(duì)應(yīng)于x、y和u的值。因此，在一個(gè)宏步長(zhǎng)內(nèi)，重新描述線性化處理后的系統(tǒng)耦合方程為

式(9)為整個(gè)系統(tǒng)的微分代數(shù)方程，因?yàn)椴豢紤]代數(shù)環(huán)，微分代數(shù)方程可以通過(guò)下面的計(jì)算改寫成常微分方程。根據(jù)假設(shè)③，常微分方程首先進(jìn)行離散化處理

Backward Euler:

同樣可以得到輸出的關(guān)系式

Backward Euler:

為了考慮缺少代數(shù)環(huán)，式(10)可以通過(guò)假設(shè)②估計(jì)得到，即輸入ω是通過(guò)在t=tk處線性化系統(tǒng)的輸出y得到的。

Backward Euler:

將式(13)代入到式(12)中，并根據(jù)假設(shè)①設(shè)定DK的乘積為冪零，得

故式(10)可以寫成

根據(jù)式(16)，每一步結(jié)束后輸出值會(huì)在子系統(tǒng)中進(jìn)行求解和傳遞

4 FMI的實(shí)現(xiàn)過(guò)程

每個(gè)功能模型單元(Function mockup unit，F(xiàn)MU)只負(fù)責(zé)本地系統(tǒng)運(yùn)算，不包含任何耦合和仿真環(huán)境的相關(guān)信息[12]。協(xié)同仿真過(guò)程中，計(jì)算機(jī)只通過(guò)主求解器來(lái)聯(lián)合系統(tǒng)中不同從求解器的信息，并向從求解器提供方法來(lái)求解其他代表耦合系統(tǒng)的線性部分的常微分方程(Ordinary differential equation，ODE)。這些信息一部分來(lái)自于每一個(gè)求解器的模型描述結(jié)構(gòu)中，一部分來(lái)自于求解FMU所包含的系統(tǒng)方向?qū)?shù)的函數(shù)[13]。

FMU求解器需要具有計(jì)算所有狀態(tài)變量和所連的輸入輸出的方向?qū)?shù)的能力，在FMI協(xié)議中，通過(guò)標(biāo)記“provide Direction Derivative”來(lái)說(shuō)明，并通過(guò)方程“fmi Get Directional Derivative”來(lái)實(shí)現(xiàn)。穩(wěn)定方法“基于模型的外推法”不與FMI提供的“基于歷史的外推法”相兼容，所以屬性“can Interpolate Inputs”需要擴(kuò)展到其所指定的外推類型。

計(jì)算流程包括兩步，第一步在通信節(jié)點(diǎn)處進(jìn)行，第二步發(fā)生在求解在通信節(jié)點(diǎn)之間的連續(xù)時(shí)間DOE系統(tǒng)處。

第一步中狀態(tài)向量的倒數(shù)求得為

它的輸出為

Jacobian矩陣為

上式中當(dāng)k等于0時(shí)，表示為全局的初始狀態(tài)。

根據(jù)式(19)可知，出現(xiàn)在式(21)中的輸入變量包含了其他的狀態(tài)變量，并且輸入變量現(xiàn)在相當(dāng)是狀態(tài)變量的初始狀態(tài)。通過(guò)求解器內(nèi)的數(shù)值積分求解出的狀態(tài)向量是，這里us實(shí)際是，通過(guò)主求解器來(lái)連接線性系統(tǒng)的輸出和求解器的輸入，因此連接的結(jié)構(gòu)信息只用主求解器獲得。

如果主求解器不支持穩(wěn)定性，式(22)的右邊會(huì)降為0，這樣從求解器的實(shí)際輸入將不變：如果沒(méi)有設(shè)置標(biāo)記“can Interpolate Inputs”，那么零階保持器外推因主求解器不支持穩(wěn)定性而在FMI協(xié)議中定義。這種情況下，通過(guò)主求解器在每一個(gè)通信節(jié)點(diǎn)調(diào)用方程“fmi Get Real”和“fmi Set Real”，每一個(gè)從求解器的輸入是通過(guò)式(20)直接給出的，是從求解器的輸出，這些從求解器會(huì)求出第s個(gè)求解器的輸入。

主求解器方面需要設(shè)置回調(diào)方程(Callback function)，這個(gè)方程用于求解相連的從求解器的動(dòng)態(tài)輸入。在聯(lián)合仿真之前，從求解器的模型結(jié)構(gòu)信息將體現(xiàn)在準(zhǔn)備矩陣和式(23)的計(jì)算中。隨后，在聯(lián)合仿真過(guò)程中，這些元素將會(huì)調(diào)用方程“fmi Get Continuous States”“fmi Get Derivatives”和“fmi Get Directional Derivatives”在每一次的信息交換節(jié)點(diǎn)上更新。當(dāng)對(duì)式(21)、(23)數(shù)值積分時(shí)，系統(tǒng)通過(guò)從數(shù)值求解器調(diào)用這些方程。

根據(jù)式(23)可知，回調(diào)方程需要返回穩(wěn)定輸入的微分向量，這需要通過(guò)主仿真環(huán)境在從設(shè)備實(shí)現(xiàn)方程屬性“fmi Callback Functions”來(lái)實(shí)現(xiàn)。圖4簡(jiǎn)略地描述了主求解器中的算法流程[14]。

圖4 主求解器中的算法流程

利用FMI2.0 for Co-simulation標(biāo)準(zhǔn)的線性隱性穩(wěn)定方法可以很好對(duì)系統(tǒng)加速，其改善了聯(lián)合仿真在平衡穩(wěn)定或精確度和計(jì)算效率上的問(wèn)題。隨著FMI2.0的問(wèn)世，一個(gè)巨大的提升是出現(xiàn)了方向?qū)?shù)矩陣的接口，這使得聯(lián)合仿真系統(tǒng)的強(qiáng)耦合性更有效[15-16]。

5 測(cè)試模型與測(cè)試結(jié)果

基于二自由度液壓模型如圖5所示，其中包括了兩個(gè)液壓子系統(tǒng)，可將這個(gè)系統(tǒng)看為兩個(gè)子系統(tǒng)實(shí)現(xiàn)聯(lián)合仿真。

圖5 在AMESim中搭建的二自由度液壓模型

所有測(cè)試都采用同一套參數(shù)，測(cè)試結(jié)果如表1所示。

表1 在不同的方案下的測(cè)試結(jié)果對(duì)比

聯(lián)合仿真的精度是通過(guò)方均根誤差來(lái)表示的。在所有的聯(lián)合仿真測(cè)試中，CPU時(shí)間通過(guò)設(shè)定一定的間隔時(shí)間為50 μs。從表1中可以看出，混合使用C/Python標(biāo)準(zhǔn)的#5方案相比#3方案在CPU耗時(shí)上更差，這是由于在兩個(gè)仿真環(huán)境中帶有很多的轉(zhuǎn)換接口；相反，對(duì)比#2和#6，同樣使用C標(biāo)準(zhǔn)的#3和#7(或#8)，系統(tǒng)明顯實(shí)現(xiàn)了加速。說(shuō)明利用FMI2.0標(biāo)準(zhǔn)的線性隱性穩(wěn)定方法可以很好地改善聯(lián)合仿真在平衡穩(wěn)定或精確度和計(jì)算效率上的問(wèn)題。

6 結(jié)論

本文的研究目的是提高聯(lián)合仿真系統(tǒng)數(shù)值積分運(yùn)算速度，同時(shí)兼顧系統(tǒng)穩(wěn)定運(yùn)行并保證積分誤差在預(yù)期范圍內(nèi)。論文研究了一種基于FMI通用接口標(biāo)準(zhǔn)的聯(lián)合仿真方法，得到如下結(jié)論。

(1) 其基本方法是將整個(gè)系統(tǒng)分為若干個(gè)子系統(tǒng)，各部分之間在通信節(jié)點(diǎn)處交換數(shù)據(jù)，數(shù)值分析時(shí)需要將系統(tǒng)合理解耦分析，在保證運(yùn)算速度的同時(shí)使系統(tǒng)具有良好的穩(wěn)定性和準(zhǔn)確性。

(2) 線性隱性穩(wěn)定方法在保證系統(tǒng)穩(wěn)定的前提下對(duì)仿真運(yùn)算時(shí)間具有一定的加速，其基于聯(lián)合仿真的FMI標(biāo)準(zhǔn)下的工程實(shí)例也證明了該方法使系統(tǒng)具有良好的穩(wěn)定性和準(zhǔn)確性。