基于功率模擬的聯(lián)合仿真技術(shù)研究

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Backward Euler:

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

Backward Euler:

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

Backward Euler:

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

故式(10)可以寫成

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

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

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

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

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

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

它的輸出為

Jacobian矩陣為

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

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

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

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

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

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

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

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

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

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

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

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

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

6 結(jié)論

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

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

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