田宏遠(yuǎn)，譚火彬，梅雅淇

（北京航空航天大學(xué) 軟件學(xué)院，北京 100191）

0 引言

隨著新能源汽車的發(fā)展，仿真軟件逐漸成為新能源汽車研發(fā)與生產(chǎn)過程中不可缺少的工具，并逐漸主導(dǎo)著研發(fā)流程［1］。新能源汽車仿真涉及整車、電控邏輯、空氣動力學(xué)和熱力學(xué)等多個領(lǐng)域，包含了物理性質(zhì)、行為表現(xiàn)等多個維度。汽車電控單元（Electronic Control Unit，ECU）是指汽車內(nèi)部各個系統(tǒng)的控制模塊，比較典型的有新能源汽車的電池充放電控制單元。本文將聚焦于汽車電控單元系統(tǒng)行為仿真，重點(diǎn)討論如何從零開始構(gòu)建一個領(lǐng)域系統(tǒng)的仿真軟件。

1 相關(guān)工作

在電控單元仿真領(lǐng)域，既有成熟的商用軟件Simulink，又有開源的如SciLab 內(nèi)的Xcos 軟件，這些軟件都采用圖形界面建模，并內(nèi)置了多種求解器進(jìn)行求解仿真。功能更為強(qiáng)大的Simulink 還支持代碼生成，將模型轉(zhuǎn)化為實(shí)際的代碼，以便作進(jìn)一步驗(yàn)證［2］。Simulink 通常用于系統(tǒng)級別的設(shè)計(jì)和仿真，尤其是控制系統(tǒng)，可以對物理硬件、嵌入式軟件、算法和系統(tǒng)運(yùn)行環(huán)境組成的復(fù)雜虛擬系統(tǒng)進(jìn)行仿真［3］。此外，對于需要進(jìn)行上位機(jī)驗(yàn)證的系統(tǒng)，還可以生成嵌入式代碼進(jìn)行后續(xù)的測試與驗(yàn)證［4］。在建模語言方面，Modelica 是一個基于方程的面向?qū)ο蟮挠?jì)算機(jī)語言，支持跨領(lǐng)域復(fù)雜物理系統(tǒng)的建模，例如機(jī)械、電子以及面向過程的子模型系統(tǒng)［5］。Modelica 采用數(shù)學(xué)方程式的形式描述各個領(lǐng)域的物理現(xiàn)象和規(guī)律，并根據(jù)系統(tǒng)拓?fù)浣Y(jié)構(gòu)和語言內(nèi)在關(guān)聯(lián)對模型進(jìn)行集成，并通過對微分方程或代數(shù)方程的求解實(shí)現(xiàn)仿真運(yùn)算［6］。

用戶在使用上述軟件時會面臨3 個主要問題：①軟件安裝過程繁瑣，僅使用軟件的部分功能卻要對軟件進(jìn)行全量安裝；②軟件對于安裝機(jī)器的性能配置要求高，增加了用戶的硬件成本；③仿真流程不透明，很難針對用戶的特定需求進(jìn)行適應(yīng)性二次開發(fā)和模型庫擴(kuò)展。

為此，本文提出基于B/S 架構(gòu)、系統(tǒng)建模語言（Systems Modeling Language，SysML）與因果框圖（Causal Block Diagram，CBD）的建模仿真設(shè)計(jì)方案以解決上述問題，并實(shí)現(xiàn)了一套易于拓展的仿真流程，方便用戶進(jìn)行二次開發(fā)和拓展。

2 相關(guān)理論

2.1 因果框圖

因果框圖（Causal Block Diagram，CBD）是一種通過不同形狀的圖形和連接線來表示模塊間依賴關(guān)系的數(shù)據(jù)結(jié)構(gòu)［7］，線條所表示的數(shù)據(jù)依賴關(guān)系也可以用來表示數(shù)據(jù)的流動方向。圖1 表示一個簡單的CBD，由A、B、C、D 4 個基礎(chǔ)模塊組成。從依賴關(guān)系的角度來看，依賴關(guān)系與連線箭頭的方向相反，A 不依賴于任何模塊，B 依賴于A，C 依賴于A 和B，D 依賴于B 和C；從數(shù)據(jù)流動的角度來看，數(shù)據(jù)從A流出，到達(dá)B 和C，又從B 流向C，最后從B 和C 一起匯入D。

Fig.1 Example of CBD圖1 CBD示例

與數(shù)據(jù)結(jié)構(gòu)中圖的定義類似，指向模塊的線的數(shù)量稱為該模塊的入度，從該模塊指出的線的數(shù)量稱為該模塊的出度。入度為0 的模塊即為整個模型的輸入，稱為輸入模塊；出度為0 的模塊即為整個模型的輸出，稱為輸出模塊。入度和出度均不為0的模塊稱為功能模塊。

CBD 最早來自于電路領(lǐng)域，不同的形狀框可以代表不同的電路元件，有著不同功能。線條則代表導(dǎo)線，象征著電流的流動方向。電路圖本身也是一種對實(shí)體電路的建模，當(dāng)電源施加不同參數(shù)時，將相應(yīng)的電學(xué)元件表現(xiàn)出的輸出用來表征電路系統(tǒng)的行為，其實(shí)就是一個典型的仿真建模過程［8］。

為了進(jìn)行層次化的系統(tǒng)建模，CBD 中也引入了分層的概念。CBD 中的一個模塊可以是由若干其他模塊和連線組成的子圖，子圖中的子模塊間也存在著依賴關(guān)系，整個子圖又在外層中與其他模塊形成依賴關(guān)系。圖2 給出了一個含有層次關(guān)系的CBD。其中，B 模塊由一個子圖組成，子圖中B3依賴于B2，B2依賴于B1。而在整個CBD 中，C 依賴于B，B 依賴于A。因此，B1、B2 和B3 均依賴于A，C依賴于B1、B2 和B3。這種層次關(guān)系建?？梢载S富模型的表達(dá)能力，支持將系統(tǒng)分解為子系統(tǒng)，為復(fù)雜系統(tǒng)建模提供了解耦和抽象功能。

Fig.2 CBD with layers圖2 含有層次的CBD

利用CBD 可以完成系統(tǒng)的建模，在下一節(jié)中，本文將為靜態(tài)圖賦予語義，使其能夠?qū)崿F(xiàn)動態(tài)的仿真功能。

2.2 離散事件系統(tǒng)仿真

離散事件系統(tǒng)（Discrete Event System，DEV）是指一種離散狀態(tài)、事件驅(qū)動的系統(tǒng)，其狀態(tài)演變完全取決于時間上異步離散事件的發(fā)生。DEV 僅由離散狀態(tài)空間和事件驅(qū)動狀態(tài)轉(zhuǎn)換機(jī)制組成［9］。離散事件系統(tǒng)嘗試用一個個時間點(diǎn)和狀態(tài)的變化去描述一個系統(tǒng)，在一個特定的時刻，每個基本模塊都有自己所處的狀態(tài)，并有相應(yīng)的狀態(tài)變化函數(shù)，用于計(jì)算該模塊下一時刻所處的狀態(tài)。

仿真模型由一個個基礎(chǔ)仿真模塊通過連接組成。使用CBD 和SysML 表示離散事件系統(tǒng)，需要定義好基礎(chǔ)模塊圖的指稱語義（Denotational Semantics）和操作語義（Operational Semantics）。指稱語義即為每個基礎(chǔ)模塊定義計(jì)算公式；操作語義即指定基礎(chǔ)模塊的計(jì)算順序，計(jì)算序列按順序計(jì)算指稱語義中定義的各個公式的值，得到最終的仿真結(jié)果。

指稱語義的定義因不同的領(lǐng)域需求而有所不同。以電控仿真領(lǐng)域?yàn)槔?，通常需要為加法器、積分器、增益器等常見基礎(chǔ)仿真模塊定義指稱語義。每個基礎(chǔ)模塊含有輸入和輸出端口，模塊從輸入端口獲取數(shù)據(jù)后，經(jīng)指稱語義計(jì)算后，由輸出端口輸出。連線用來連接兩個基礎(chǔ)仿真模塊的輸入和輸出，連線兩端的數(shù)據(jù)是相同的。特別的，不含輸入端口的模塊稱為輸入模塊，是整個模型的輸入；不含輸出端口的模塊稱為輸出模塊，是整個模型的輸出。

模型圖通?？梢越?jīng)過處理后轉(zhuǎn)化為一個有向無環(huán)圖。對于含有數(shù)據(jù)依賴關(guān)系的有向無環(huán)圖［10］，可以采用拓?fù)渑判虻姆椒ǖ玫侥K的先后順序，該先后順序就是指稱語義的計(jì)算順序。拓?fù)渑判虻乃枷胧菑娜攵葹榱愕墓?jié)點(diǎn)開始（對應(yīng)于CBD 中的輸入模塊），斷開其與節(jié)點(diǎn)的連線，并將這些入度為零的節(jié)點(diǎn)添加到結(jié)果列表中，標(biāo)識為已排序，然后去尋找圖中入度為零的節(jié)點(diǎn)，以此類推，直至所有節(jié)點(diǎn)都被訪問過。由于圖中沒有圈的存在，因此算法執(zhí)行到最后一定能訪問完所有節(jié)點(diǎn)。

如圖3 所示，圖的入度為零的模塊（即輸入模塊）是A，出度為零的模塊（即輸出模塊）是E。拓?fù)渑判蛩惴ㄊ紫葟腁 開始，刪除其指向B 和C 的連線，并將A 加入已排序的結(jié)果列表sorted。之后，B 模塊為唯一入度為零的模塊，算法訪問B 模塊，刪除其指向C 和D 的連線，并將B 加入已排序的結(jié)果列表。以此類推，算法執(zhí)行完畢后得到圖的拓?fù)渑判蚪Y(jié)果為［A，B，C，D，E］，且已遍歷完所有模塊。

Fig.3 Example of topological sort圖3 拓?fù)渑判蚴纠?/p>

操作語義所要完成的工作就是對于拓?fù)渑判蚝蟮慕Y(jié)果列表，在每個時間步內(nèi)逐一遍歷并執(zhí)行相應(yīng)模塊的指稱語義。下面將以偽代碼的形式描述操作語義指導(dǎo)下的指稱語義計(jì)算過程，如算法1所示。

算法1 離散事件系統(tǒng)仿真算法

輸入：表示模型的有向無環(huán)圖G，仿真時長T，仿真步長step

輸出：仿真結(jié)果

2.3 連續(xù)事件系統(tǒng)仿真

在連續(xù)事件系統(tǒng)仿真中，需要對連續(xù)時間和狀態(tài)進(jìn)行建模與仿真，此時則需要引入新的描述形式——微分方程系統(tǒng)（Differential Equation System，DES）［12］。使用微分方程進(jìn)行系統(tǒng)建模時，并不會通過狀態(tài)轉(zhuǎn)移函數(shù)來計(jì)算下一時刻模塊的狀態(tài)變化，而是通過導(dǎo)函數(shù)的形式給出狀態(tài)的變化率來描述狀態(tài)如何變化。在給定狀態(tài)值和時間段的情況下，任意時刻的狀態(tài)值需要根據(jù)導(dǎo)函數(shù)來計(jì)算?？梢圆扇?shù)值積分的方法對連續(xù)事件系統(tǒng)的狀態(tài)進(jìn)行求解，學(xué)界對此已有很多研究，例如歐拉法、龍格庫塔（Runge-Kutta）法（根據(jù)不同情況分為多種階數(shù)，常見的有二階和四階）［11］。

2.4 代數(shù)環(huán)

代數(shù)環(huán)（algebraic loop）［13］是指CBD 中模塊之間的依賴關(guān)系構(gòu)成一個環(huán)，即形成模塊間的循環(huán)依賴。具體地，在仿真模型中，如果一個模塊的值取決于另一個模塊的值，而該模塊的值又依賴于第一個模塊的值，則會形成一個代數(shù)環(huán)［14］。這會導(dǎo)致仿真系統(tǒng)無法求解變量的值，因?yàn)槊總€模塊都需要已知其它模塊的值才能計(jì)算出自己的值，類似于操作系統(tǒng)中的“死鎖”。

圖4 展示的即是建模過程中產(chǎn)生的一個代數(shù)環(huán)。其中，G 是一個增益模塊，目的是將輸入放大為原來的2 倍后輸出。圖中明顯形成了一個圈回路。要想計(jì)算C 的值，就需要先計(jì)算C 的增益與2 相加后的結(jié)果，從而造成了循環(huán)依賴。從數(shù)學(xué)形式上，可以通過解方程得到C 的值為-2。但是如果放入到仿真流程中，經(jīng)過一圈的計(jì)算，還是回到了C 的輸出上，僅憑語句無法直接得出C 的值，除非經(jīng)過進(jìn)一步的數(shù)學(xué)推理（即方程求解）。

Fig.4 Example of algebraic loop圖4 代數(shù)環(huán)示例

如前文所述，CBD 表達(dá)了模塊間數(shù)據(jù)的依賴性，因此可以采用拓?fù)渑判虻姆绞綄A(chǔ)模塊進(jìn)行排序［15］，得到的順序就是正確的計(jì)算順序。然而，由于代數(shù)環(huán)的存在，用戶建立好的CBD 模型并不能直接用于拓?fù)渑判?，而是需要?jīng)過一定的特殊處理，這種處理被稱為代數(shù)環(huán)的消解。

關(guān)于代數(shù)環(huán)的消解，目前學(xué)界有多種解決辦法，最常用的就是將代數(shù)環(huán)轉(zhuǎn)化為數(shù)學(xué)方程的形式進(jìn)行矩陣運(yùn)算求解［16］。但并非所有形式的代數(shù)環(huán)均有解，因此在建模過程中要盡可能避免代數(shù)環(huán)的出現(xiàn)。

然而，并非所有在模型中出現(xiàn)的圈都是代數(shù)環(huán)。在CBD 模型中，當(dāng)存在回路并且回路中不存在非直接饋通模塊（non-direct feedthrough block）時，才會出現(xiàn)代數(shù)環(huán)。非直接饋通模塊是指模塊無需知道輸入信號即可計(jì)算當(dāng)前時間步的輸出，通常維護(hù)一個state 變量來存儲過往的計(jì)算值直接輸出，典型的非直接饋通模塊有積分模塊（Integrator）和單位時延模塊（Unit Delay）。因此，在環(huán)路中可以斷開所有連接到非直接饋通模塊的連線，使得環(huán)路斷開，此時則不再有代數(shù)環(huán)的問題，如圖5所示。

Fig.5 Breaking loop according to non-direct feedthrough module圖5 根據(jù)非直接饋通模塊斷開環(huán)路

3 仿真系統(tǒng)設(shè)計(jì)

用戶使用建模仿真軟件進(jìn)行系統(tǒng)建模仿真的一個完整流程如圖6 所示。首先用戶在建模端進(jìn)行相應(yīng)的模型搭建，并設(shè)定好仿真參數(shù)（如仿真時長、仿真步長等）；然后建模端會以JSON（JavaScript Object Notation）文件的形式將模型保存，并發(fā)送到仿真端；仿真端接收到文件后，則進(jìn)行相應(yīng)的仿真運(yùn)行操作。仿真器將根據(jù)用戶設(shè)置的參數(shù)運(yùn)行模型。在運(yùn)行時，會判斷當(dāng)前時間步是否到達(dá)仿真時長，如果沒有到達(dá)則繼續(xù)進(jìn)行仿真運(yùn)算；如果到達(dá)，則視為仿真完成。運(yùn)行完成后將結(jié)果返回給建模端進(jìn)行可視化展示和數(shù)據(jù)下載。從提交仿真模型直至仿真結(jié)束所涉及的全部過程稱為仿真階段。本文軟件關(guān)注的就是仿真階段的全部功能。

Fig.6 Flow of modeling and simulating圖6 建模仿真流程

3.1 軟件總體設(shè)計(jì)

軟件采取分層架構(gòu)，根據(jù)功能自頂向下分為數(shù)據(jù)交互層、數(shù)據(jù)解析層和仿真運(yùn)算層。軟件體系結(jié)構(gòu)如圖7所示。

Fig.7 Software architecture圖7 軟件體系結(jié)構(gòu)

數(shù)據(jù)交互層主要負(fù)責(zé)模型接收和遠(yuǎn)程模型調(diào)用，以及仿真結(jié)果的返回。接收到的數(shù)據(jù)將交給數(shù)據(jù)解析層，將結(jié)果返回給建模端進(jìn)行展示。

數(shù)據(jù)解析層主要負(fù)責(zé)解析數(shù)據(jù)交互層接收到的數(shù)據(jù)，包括將模型文件解析為內(nèi)存中的對象、狀態(tài)機(jī)語句值計(jì)算以及遠(yuǎn)程模型文件解析。解析后的模型信息將傳給下一層仿真運(yùn)算層進(jìn)行仿真計(jì)算。

最底層是仿真運(yùn)算層，也是整個仿真工具的核心。在這一層中，首先需要有基礎(chǔ)模塊庫支撐整個模型的搭建和仿真運(yùn)行。其次，需要求解器支持連續(xù)模型的方程求解，即數(shù)值積分計(jì)算。當(dāng)模型被解析完成后，會通過順序表構(gòu)建、模型自檢進(jìn)行仿真運(yùn)行前的準(zhǔn)備，之后即可運(yùn)行仿真。如果涉及聯(lián)合仿真，將會與其他模型進(jìn)行數(shù)據(jù)交互，實(shí)現(xiàn)模型之間的同步運(yùn)行。后文將重點(diǎn)關(guān)注仿真運(yùn)算層的設(shè)計(jì)與實(shí)現(xiàn)。

3.2 基礎(chǔ)模塊庫設(shè)計(jì)

基礎(chǔ)模塊庫是支撐建模與仿真的基礎(chǔ)?；A(chǔ)模塊是指參與建模的最基本的計(jì)算單元（對應(yīng)CBD 中的模塊），基礎(chǔ)模塊間通過有方向的連接線的連接構(gòu)成模型圖［17］。

所有基礎(chǔ)模塊都具有相似的一部分功能和屬性，也有其自身的特點(diǎn)，因此可以使用接口的形式將公共的功能抽象出來。

對于模塊共有的屬性，可以采取繼承父類子類的形式將共有屬性傳遞到每個子類中，或者選擇組合的方式——即持有該公共屬性實(shí)例，作為字段出現(xiàn)在其他類中來實(shí)現(xiàn)公共屬性的統(tǒng)一管理。圖8 是所有模塊公共屬性抽象出的共有結(jié)構(gòu)體BasicProperty，擁有編號id、組件名name 和是否啟用enable 3個屬性。

Fig.8 Class diagram of public basic property圖8 公共基本屬性類圖

所有基本操作方法組成的特征稱為BasicOperation，如圖9 所示。具體方法主要包括：①獲取模塊類型：用來判斷模塊屬于哪一類；②獲取模塊的輸出值：類型是Value，可能是矩陣，也可能是一個單值；③模塊自檢：用來檢查模塊是否滿足仿真規(guī)范；④模塊更新：用來更新模塊狀態(tài)值或輸出值；⑤模塊重置：用來重制并清空模塊歷史狀態(tài)信息；⑥獲取連接數(shù)：用來獲取有多少其他模塊連接到當(dāng)前模塊，即有多少輸入端口；⑦獲取第i 個輸入端值：用來獲取連接到當(dāng)前模塊的第i 個模塊的輸出值（當(dāng)前狀態(tài)值）；⑧連接：用來將一個模塊連接到當(dāng)前模塊。此外，還有用來判斷模塊是否是離散模塊以及直接饋通模塊的方法。

Fig.9 Class diagram of public basic operation圖9 公共方法組成的特征類圖

對于每一個具體模塊，可以定義其結(jié)構(gòu)體，并為該結(jié)構(gòu)體實(shí)現(xiàn)BasicOperation 接口，同時將BasicProperty 作為該結(jié)構(gòu)體的屬性字段。對于模塊自身特有的行為和屬性，也可以為該模塊的結(jié)構(gòu)體增添相應(yīng)的屬性字段和方法。

對于模型中的每一個基礎(chǔ)模塊，完成實(shí)例化之后，需要存儲在一個集合中供仿真器訪問運(yùn)行。然而，每一個基礎(chǔ)模塊都是一個不同的結(jié)構(gòu)體，如何使用一個集合類型存儲不同的結(jié)構(gòu)體類型是一個重要問題。在面向?qū)ο蟮恼Z言中，可以考慮采用多態(tài)的形式——父類類型（或接口類型）指向子類實(shí)例的方式，以使用一種類型代表多種不同類型［18］。

3.3 仿真流程框架

離散與連續(xù)事件的仿真計(jì)算是指根據(jù)拓?fù)渑判虻玫降幕A(chǔ)模塊排序序列依次遍歷每個基礎(chǔ)模塊，根據(jù)模塊的性質(zhì)進(jìn)行離散狀態(tài)更新或連續(xù)數(shù)值積分計(jì)算的過程。后文將針對這兩種計(jì)算的特點(diǎn)進(jìn)行詳細(xì)介紹。

仿真計(jì)算階段的計(jì)算流程如圖10 所示。在仿真計(jì)算階段，仿真器將根據(jù)模型提供的信息，從用戶指定的開始時間到結(jié)束時間內(nèi)，以固定間隔連續(xù)計(jì)算系統(tǒng)的狀態(tài)和輸出。計(jì)算模型狀態(tài)和輸出的這些連續(xù)時間點(diǎn)即稱為時間步，時間步之間的間隔稱為步長。步長可以取決于很多因素，如求解器類型、過零檢測［19］和采樣時間等。

Fig.10 Flow of simulation calculation圖10 仿真計(jì)算流程

仿真開始時，模型會設(shè)定系統(tǒng)的初始狀態(tài)和輸出，然后在每個時間步中計(jì)算系統(tǒng)輸入、狀態(tài)和輸出的新值，并對模型進(jìn)行更新（不涉及拓?fù)浣Y(jié)構(gòu)的更改）以應(yīng)用這些計(jì)算結(jié)果。當(dāng)仿真結(jié)束時，模型將輸出用戶所關(guān)心的模塊狀態(tài)及模塊的最終值。

在仿真計(jì)算流程的一開始，執(zhí)行引擎會調(diào)用計(jì)算輸出模塊的輸出方法使仿真器開始工作，仿真器將根據(jù)模塊狀態(tài)更新順序表指定的計(jì)算順序依次調(diào)用模塊的更新方法。對于基礎(chǔ)模塊的狀態(tài)更新，如果模型中僅有離散狀態(tài)模塊，仿真器會根據(jù)采樣時間計(jì)算仿真步長（可以是固定步長，也可以是可變步長），并調(diào)用模型的更新方法。模型會按順序調(diào)用每個基礎(chǔ)模塊的更新方法。

對于包含連續(xù)狀態(tài)的模型，仿真器將使用求解器對模型進(jìn)行求解。具體地，根據(jù)求解器的數(shù)值積分算法，求解器會進(jìn)入子時間步循環(huán)。在子時間步中，模塊會以更小的仿真步長進(jìn)行計(jì)算。在該循環(huán)中，求解器會重復(fù)調(diào)用模型的更新方法，使得在一個主時間步中可以按照連續(xù)時間間隔計(jì)算模型的輸出。使用該方法可以提高數(shù)值積分計(jì)算的準(zhǔn)確度。

最后，仿真器會判斷是否到達(dá)仿真結(jié)束時間，如果沒到，則循環(huán)上述過程；如果到達(dá)，則結(jié)束仿真并輸出結(jié)果。

3.4 離散狀態(tài)模塊計(jì)算

離散模塊的典型代表就是單位時延模塊。單位時延模塊在每個時間步的輸出結(jié)果是上一步結(jié)果的值，在0 時刻的初值是一個指定值。以上定義是在采樣時間與仿真步長相同的情況下，在實(shí)際情況中，采樣時間與仿真步長可能不一致。仿真步長是指模型在仿真時運(yùn)算的間隔時間，例如設(shè)置為0.1 的步長意味著每隔0.1 s，整個模型系統(tǒng)就要進(jìn)行一次計(jì)算。而模塊的采樣時間是一個參數(shù)，其指示在仿真過程中模塊何時生成輸出，并在適當(dāng)時更新其內(nèi)部狀態(tài)［20］。內(nèi)部狀態(tài)包括但不限于記錄的連續(xù)狀態(tài)和離散狀態(tài)。模塊的采樣時間可以由用戶指定，如果無需關(guān)注每個仿真步長的結(jié)果，通?？梢詫⒉蓸訒r間設(shè)置為仿真步長的整數(shù)倍。因此，修改后的時延模塊定義如下：單位時延模塊在每個采樣時間的輸出結(jié)果是上一個采樣周期的值，在0 時刻的初值是一個指定值。例如，若仿真步長為0.1，采樣時間為1，則在0～0.9 的仿真時間內(nèi)，模塊的輸出均為x(0)；在1.0～1.9 的仿真時間內(nèi)，模塊的輸出為x(1)，以此類推。

3.5 連續(xù)狀態(tài)模塊計(jì)算

連續(xù)狀態(tài)系統(tǒng)通常采用微分方程的形式進(jìn)行系統(tǒng)建模，因此連續(xù)狀態(tài)模型的求解依賴于數(shù)學(xué)上的數(shù)值積分算法。連續(xù)狀態(tài)求解并不是在真正的連續(xù)時間下進(jìn)行計(jì)算，而是將時間分割成足夠密集的時間點(diǎn)去模擬連續(xù)計(jì)算。

該軟件實(shí)現(xiàn)了歐拉法和四階龍格庫塔法兩種數(shù)值積分求解器，從而支持連續(xù)狀態(tài)模型的求解。歐拉法的計(jì)算相對簡單，此處不再贅述，本節(jié)將著重介紹四階龍格庫塔法。四階龍格庫塔法簡稱RK4 法，采用高階的非因果數(shù)值積分求解方法能在保證高精度的同時，具有良好的通用性、穩(wěn)定性和靈活性［21］。

四階龍格庫塔法采用離散化的時間步長來解決連續(xù)問題，在每個時間步內(nèi)計(jì)算若干個點(diǎn)的函數(shù)斜率，并使用這些點(diǎn)斜率的平均值去估計(jì)下一個時間步的函數(shù)值。其大致流程如下：

（1）在每個時間步內(nèi)，計(jì)算函數(shù)在當(dāng)前仿真時間點(diǎn)的導(dǎo)數(shù)以求得斜率，并使用若干個仿真時間點(diǎn)與下一個時間點(diǎn)的中間點(diǎn)計(jì)算斜率，之后將上述斜率的加權(quán)平均值作為當(dāng)前仿真步的斜率。

（2）使用上一步得到的斜率更新下一個時間步函數(shù)值。

使用數(shù)學(xué)公式描述上述流程如下：

微分方程的初值表示為：

則使用RK4對該方程的求解為：

其中：

其中，h 代表仿真步長。通過上述計(jì)算，函數(shù)在下一時刻的值由當(dāng)前值與當(dāng)前值附近若干點(diǎn)的斜率和仿真步長所決定。

根據(jù)上述數(shù)值積分算法可以得出，在每個仿真步長內(nèi)，涉及到數(shù)值積分的計(jì)算僅發(fā)生在當(dāng)前仿真時刻tn、tn+和tn+h 3 個時間點(diǎn)。因此，在涉及連續(xù)狀態(tài)的模型求解時，需要將單個仿真步長的計(jì)算細(xì)化為3 個階段，在此3 個階段內(nèi)完成離散狀態(tài)的求解與更新以及連續(xù)狀態(tài)下方程的數(shù)值積分求解。具體的，在tn和tn+h 兩個仿真時刻都會觸發(fā)離散模塊的計(jì)算，而由于下一個仿真時刻tn+1=tn+h，因此在下一個仿真時刻會發(fā)生離散模塊的重復(fù)計(jì)算。此類重復(fù)計(jì)算在存在特定模塊（如計(jì)數(shù)器模塊）時會影響結(jié)果。為避免這種情況的發(fā)生，需要在tn+h時刻設(shè)置離散模塊的狀態(tài)為未啟用來阻止含有狀態(tài)的離散模塊更新。

三階段的單步仿真流程如圖11 所示。在每一個單次仿真時間步內(nèi)，首先會按照模塊更新順序表更新除積分器外的其他基礎(chǔ)模塊。具體地，首先會針對時延模塊調(diào)用更新輸出值函數(shù)，使得時延模塊的初始值或上一個時間步計(jì)算的狀態(tài)值更新到輸出。接著會調(diào)用所有時延模塊和輸出模塊的更新函數(shù)，以計(jì)算時延模塊當(dāng)前時間步的值，并從所有輸出模塊出發(fā)，按照模塊更新順序表對應(yīng)的順序更新與輸出相連的模塊。在此過程中，可以計(jì)算RK4 中的k1。當(dāng)其他全部模塊的狀態(tài)更新完成后，將禁用上文中提及的模塊以避免產(chǎn)生重復(fù)更新錯誤，然后將時間步設(shè)置為，計(jì)算k2和k3。最后，將時間步設(shè)置為tn+h，計(jì)算k4和yn+1。計(jì)算完成后，啟用剛才禁用的全部模塊，并重復(fù)上述過程進(jìn)行下一時間步的計(jì)算，直到仿真結(jié)束。

Fig.11 Three stages of single-step simulation圖11 單步仿真3個階段

通過上述過程即可實(shí)現(xiàn)離散和連續(xù)系統(tǒng)的混合仿真。當(dāng)需要增添新的解算器時，需根據(jù)實(shí)際的數(shù)值積分算法進(jìn)行階段的微調(diào)，但總體計(jì)算順序不變。

4 測試與驗(yàn)證

為了對仿真系統(tǒng)實(shí)例進(jìn)行測試，本節(jié)根據(jù)基礎(chǔ)模塊框架分別搭建了3 種不同類型的基礎(chǔ)模型實(shí)例，并在Simulink 中搭建同樣的3 個模型，通過對比來驗(yàn)證基礎(chǔ)模型仿真求解器仿真結(jié)果的可靠性與準(zhǔn)確性。

4.1 離散系統(tǒng)仿真

本文以汽車制動系統(tǒng)中線性化的滑移率—附著系數(shù)關(guān)系模型作為離散系統(tǒng)仿真實(shí)例。防抱死制動系統(tǒng)（Antilock Brake System，ABS）的主要功能是通過控制車輪滑移率來動態(tài)調(diào)整車輪產(chǎn)生的最大摩擦力，進(jìn)而保證車輛在制動過程中能夠保持較為穩(wěn)定的轉(zhuǎn)向控制，并在最短距離內(nèi)停車。車輪與地面之間的摩擦系數(shù)μ又稱為路面附著系數(shù)，其與車輪滑移率s之間存在著一定的非線性關(guān)系。在此本文將其用兩段直線進(jìn)行近似處理，得到兩段式的路面附著系數(shù)與車輪滑移率的變化關(guān)系。其變化規(guī)律可由公式（4）所示的分段差分方程表示，其中μh表示路面附著系數(shù)的峰值，μg表示車輪完全抱死時的路面附著系數(shù)，s0表示路面附著系數(shù)處于峰值時所對應(yīng)的輪滑率：

車輪輪滑率的取值范圍為0～1，為模擬路面附著系數(shù)隨車輪輪滑率的線性變化，本文搭建了如圖12 所示的模型。模型覆蓋了常數(shù)模塊（Constant）、單位延遲模塊（Unit-Delay）、加法模塊（Sum）、If 分支模塊（IfBlock、IfAction）、信號合并模塊（Merge）以及輸出模塊（Output）。其中，對于公式（4）所涉及的變量，設(shè)置μh=0.8，μg=0.6，s0=0.2。該模型的輸入是常量0.001，代表一個時間步的步長。首先使用單位時延模塊，延后一個時間步對變量u1進(jìn)行0.001（即一個時間步）的累加，使u1代表當(dāng)前時間步長。將變量u1輸入給一個邏輯判斷模塊，如果u1的值小于或等于0.2，則進(jìn)行公式（4）中第一行的計(jì)算，否則進(jìn)行公式（4）中第二行的計(jì)算。二者的計(jì)算結(jié)果通過merge 模塊復(fù)合在一起，并輸出到顯示模塊，最終輸出函數(shù)圖像。

Fig.12 Relationship model between linearized slip rate and adhesion coefficient圖12 線性化滑移率—附著系數(shù)關(guān)系模型

在使用本仿真求解器進(jìn)行仿真運(yùn)算時，設(shè)置仿真時長為1 s，仿真步長為0.001 s，求解得到如圖13 所示的滑移率—附著系數(shù)關(guān)系曲線。同時在MATLAB（R2022a）/Simulink 中選擇同樣的仿真模型，設(shè)置步長類型為定步長，求解器為離散（無連續(xù)狀態(tài)）類型，通過仿真運(yùn)行得到如圖14所示的滑移率—附著系數(shù)關(guān)系曲線。

Fig.13 Curve depicting the relationship between slip rate and adhesion coefficient（simulation results of this software）圖13 滑移率—附著系數(shù)關(guān)系曲線（本文軟件仿真結(jié)果）

Fig.14 Curve depicting the relationship between slip rate and adhesion coefficient（simulation results of MATLAB（R2022a）/Simulink）圖14 滑移率—附著系數(shù)關(guān)系曲線（MATLAB（R2022a）/Simulink仿真結(jié)果）

4.2 離散系統(tǒng)仿真

本文以彈球模型作為連續(xù)系統(tǒng)仿真的實(shí)例。彈球模型是一個經(jīng)典的物理學(xué)問題，也是一個典型的基礎(chǔ)連續(xù)系統(tǒng)模型。當(dāng)一個球以一定的初速度和角度從一定高度自由落下時，其將受到重力的作用加速下落。當(dāng)球碰到地面時，其能量將轉(zhuǎn)換成彈性形變的勢能，并隨即釋放出來，使球反彈起來。根據(jù)牛頓第二定律，球的速度變化和位置變化都是連續(xù)的，如公式（5）所示。其中，g表示重力產(chǎn)生的加速度，x(t)表示球在t時刻的位置，v(t)表示球在t時刻的速度。

當(dāng)小球落到地面并與地面發(fā)生彈性碰撞時，其碰撞前后的速度變化可由公式（6）表示。其中，k是小球的恢復(fù)系數(shù)，v-是小球碰撞前的速度，v+是小球碰撞后的速度。

整個彈球系統(tǒng)模型如圖15 所示，模型覆蓋了常數(shù)模塊（Constant）、增益模塊（Gain）、積分模塊（Integrator）、關(guān)系運(yùn)算模塊（RelationalOperator）以及輸出模塊（Output）。其中，對于公式（5）和公式（6）所涉及的變量，設(shè)置模型參數(shù)g=9.81，k=0.8，彈球初始速度v0=15，初始位置x0=10。該模型的輸入是常量-9.81，代表重力系數(shù)，輸出是速度積分器和位置積分器在每一個時間步的結(jié)果（即小球的速度和位置變化）。對于速度積分器，常量-9.81 作為積分參數(shù)輸入，常量15 代表初速度，是積分器的初始值（僅在仿真開始時生效），-0.8 代表反彈速度轉(zhuǎn)換系數(shù)，是小球每次落地后經(jīng)過反彈得到的反向速度值（即若以v0的速度落地，則彈起的瞬時速度為-0.8v0）。對于位置積分器，常量10為小球的初始位置值，也作為位置積分器的初值（僅在仿真開始時生效）。將速度積分器的結(jié)果作為距離積分器的輸入進(jìn)行積分計(jì)算，得到小球在該步長時間內(nèi)的距離值。此外，兩個積分器會對距離是否大于或等于0 進(jìn)行判斷。當(dāng)距離結(jié)果由小于零轉(zhuǎn)換到大于或等于零時，會觸發(fā)速度積分器和距離積分器的上升沿信號進(jìn)行重置運(yùn)算，將距離結(jié)果重置為0，速度結(jié)果重置為當(dāng)前結(jié)果乘以-0.8。

Fig.15 Model of bouncing ball system圖15 彈球系統(tǒng)模型

在使用本仿真求解器進(jìn)行仿真運(yùn)算時，設(shè)置參數(shù)如下：設(shè)置仿真時長為20s，仿真步長為0.01s，得到如圖16 所示的有關(guān)小球速度和位置變化的仿真曲線。同時在MATLAB（R2022a）/Simulink 中選擇同樣的仿真模型，設(shè)置步長類型為定步長，求解器類型為ode4（Runge-Kutta），通過仿真運(yùn)行得到如圖17 所示的關(guān)于小球速度和位置變化的仿真曲線。

Fig.16 Velocity and position change curve（simulation results of this software）圖16 速度與位置變化曲線（本文軟件仿真結(jié)果）

Fig.17 Velocity and position change curve（simulation results of MATLAB（R2022a）/Simulink）圖17 速度與位置變化曲線（MATLAB（R2022a）/Simulink仿真結(jié)果）

4.3 混合系統(tǒng)仿真

混合系統(tǒng)中既包含連續(xù)狀態(tài)的求解，又包含離散狀態(tài)的求解，在實(shí)際仿真中應(yīng)用廣泛。例如在控制系統(tǒng)中，控制器的狀態(tài)通常是離散的，而被控制的系統(tǒng)狀態(tài)是連續(xù)的，因此需要使用混合系統(tǒng)模型來描述整個控制系統(tǒng)。另外，混合系統(tǒng)模型也被廣泛應(yīng)用于嵌入式系統(tǒng)中，例如汽車電控單元、飛機(jī)控制系統(tǒng)等。

本文以汽車行駛速度控制系統(tǒng)作為混合系統(tǒng)的仿真實(shí)例，如圖18 所示。汽車行駛控制系統(tǒng)控制著汽車在行駛過程中的速度變化，其工作原理是先測量汽車的實(shí)際速度，然后根據(jù)汽車速度操縱機(jī)構(gòu)的位置計(jì)算汽車的理論速度。通過這兩個速度的差值計(jì)算應(yīng)該施加的牽引力，并由此牽引力改變汽車速度，直到其穩(wěn)定在指定的速度為止。其中，對實(shí)際速度的計(jì)算是連續(xù)的，其余的數(shù)值計(jì)算是離散的。

Fig.18 Model of automobile speed control system圖18 汽車行駛速度控制系統(tǒng)模型

對于汽車?yán)碚撍俣鹊挠?jì)算，其數(shù)學(xué)模型可由公式（7）表示。其中，x表示當(dāng)前時刻汽車速度操縱機(jī)構(gòu)的位置，vstandard表示根據(jù)汽車速度操縱機(jī)構(gòu)位置計(jì)算出的汽車的理論速度。

對于汽車實(shí)際速度的計(jì)算，其數(shù)學(xué)模型可由公式（8）表示。其中，F(xiàn)表示汽車牽引力，m表示汽車質(zhì)量，b表示摩擦阻力因子，vreal表示汽車的實(shí)際速度。

汽車牽引力則由行駛控制器PID 根據(jù)汽車當(dāng)前速度與指定速度的差值計(jì)算得到，其數(shù)學(xué)模型可由公式（9）表示。其中，u(n)表示系統(tǒng)輸入的理論速度，y(n)表示系統(tǒng)輸出的牽引力；s(n)為一個中間量，表示對速度進(jìn)行時延累加；KP、KI、KD分別為PID 控制器的比例、積分與微分控制參數(shù)。

汽車行駛速度控制系統(tǒng)模型如圖18 所示，模型覆蓋了常數(shù)模塊（Constant）、增益模塊（Gain）、加法模塊（Sum）、減法模塊（Subtract）、單位延遲模塊（UnitDelay）、乘法模塊（Product）、除法模塊（Divide）以及輸出模塊（Output）。其中，對于公式（7）—公式（9）所涉及的變量，本文設(shè)置模型參數(shù)KP=1，KI=0.01，KD=0.01，m=1500，b=20，并設(shè)置延時模塊采樣時間為0.02s。在該模型中，輸入為操縱機(jī)構(gòu)的位置，用常數(shù)0.5 表示，輸出是汽車的實(shí)際速度和牽引力。首先根據(jù)公式（7）計(jì)算出汽車的理論速度作為公式（9）中的u(n)輸入系統(tǒng)，此時模型的計(jì)算分為3 路：最上路將u(n)直接輸出，對應(yīng)公式（9）中的KPu(n)；中路將u(n)與一個單位時延前的結(jié)果相加，即s(t)=s(t-1)+u(t)；下路將u(n)與一個單位時延前的結(jié)果相減，即d(n)=u(n)-u(n-1)。以上3 路結(jié)果按公式（9）中的牽引力輸出公式進(jìn)行計(jì)算，得到當(dāng)前時刻的牽引力。該牽引力參與公式（8）的積分計(jì)算，得到實(shí)際速度。

在使用本仿真求解器進(jìn)行仿真運(yùn)算時，設(shè)置參數(shù)如下：仿真時長為1 000 s，仿真步長為0.01 s，求解得到如圖19 所示的汽車牽引力與行駛速度變化曲線。同時，在MATLAB（R2022a）/Simulink 中建立同樣的仿真模型，設(shè)置步長類型為定步長，求解器類型為ode4（Runge-Kutta），仿真運(yùn)行得到如圖20 所示的汽車牽引力與行駛速度變化曲線。

Fig.19 Curve depicting the relationship between automobile traction force and velocity change（simulation results of this software）圖19 汽車牽引力與速度變化曲線（本文軟件仿真結(jié)果）

Fig.20 Curve depicting the relationship between automobile traction force and velocity change（simulation results of MATLAB（R2022a）/Simulink）圖20 汽車牽引力與速度變化曲線（MATLAB（R2022a）/Simulink仿真結(jié)果）

5 結(jié)語

本文主要圍繞如何從零開始搭建一個領(lǐng)域系統(tǒng)仿真軟件的問題，從仿真軟件的相關(guān)基礎(chǔ)理論開始，闡述了基于CBD 的仿真系統(tǒng)的仿真計(jì)算流程，并介紹了相應(yīng)算法。在軟件實(shí)現(xiàn)上，從軟件整體架構(gòu)出發(fā)，介紹了如何搭建仿真基礎(chǔ)模塊庫，并闡明了如何設(shè)計(jì)一套仿真流程以支持離散事件系統(tǒng)仿真和連續(xù)事件系統(tǒng)仿真兩大仿真功能。對于微分方程的求解，采用業(yè)界常用的四階龍格庫塔法，并提出在仿真流程中使用四階龍格庫塔法進(jìn)行微分方程求解的方法。最后針對多個領(lǐng)域的仿真典型案例，使用本文設(shè)計(jì)的軟件與常用仿真軟件MATLAB（R2022a）/Simulink進(jìn)行了對比測試，本文設(shè)計(jì)的軟件在準(zhǔn)確度和效率上符合預(yù)期。

本文仿真方法的優(yōu)點(diǎn)在于采用了接口式設(shè)計(jì)，代碼各個模塊之間的耦合度較低，可以按照接口規(guī)范快速實(shí)現(xiàn)對基礎(chǔ)模塊庫和解算器的擴(kuò)充，以提升模型的仿真能力。本文方法的局限性在于：對模型的仿真范圍在很大程度上依賴于基礎(chǔ)模塊庫，基礎(chǔ)模塊庫限制著模型的表達(dá)能力，而沒有提供一套可編程的形式供用戶自己編程實(shí)現(xiàn)相應(yīng)基礎(chǔ)模塊的計(jì)算邏輯。