汽車電控單元仿真軟件設計與實現

田宏遠，譚火彬，梅雅淇

（北京航空航天大學 軟件學院，北京 100191）

0 引言

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

1 相關工作

在電控單元仿真領域，既有成熟的商用軟件Simulink，又有開源的如SciLab 內的Xcos 軟件，這些軟件都采用圖形界面建模，并內置了多種求解器進行求解仿真。功能更為強大的Simulink 還支持代碼生成，將模型轉化為實際的代碼，以便作進一步驗證［2］。Simulink 通常用于系統級別的設計和仿真，尤其是控制系統，可以對物理硬件、嵌入式軟件、算法和系統運行環境組成的復雜虛擬系統進行仿真［3］。此外，對于需要進行上位機驗證的系統，還可以生成嵌入式代碼進行后續的測試與驗證［4］。在建模語言方面，Modelica 是一個基于方程的面向對象的計算機語言，支持跨領域復雜物理系統的建模，例如機械、電子以及面向過程的子模型系統［5］。Modelica 采用數學方程式的形式描述各個領域的物理現象和規律，并根據系統拓撲結構和語言內在關聯對模型進行集成，并通過對微分方程或代數方程的求解實現仿真運算［6］。

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

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

2 相關理論

2.1 因果框圖

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

Fig.1 Example of CBD圖1 CBD示例

與數據結構中圖的定義類似，指向模塊的線的數量稱為該模塊的入度，從該模塊指出的線的數量稱為該模塊的出度。入度為0 的模塊即為整個模型的輸入，稱為輸入模塊；出度為0 的模塊即為整個模型的輸出，稱為輸出模塊。入度和出度均不為0的模塊稱為功能模塊。

CBD 最早來自于電路領域，不同的形狀框可以代表不同的電路元件，有著不同功能。線條則代表導線，象征著電流的流動方向。電路圖本身也是一種對實體電路的建模，當電源施加不同參數時，將相應的電學元件表現出的輸出用來表征電路系統的行為，其實就是一個典型的仿真建模過程［8］。

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

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

利用CBD 可以完成系統的建模，在下一節中，本文將為靜態圖賦予語義，使其能夠實現動態的仿真功能。

2.2 離散事件系統仿真

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

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

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

模型圖通常可以經過處理后轉化為一個有向無環圖。對于含有數據依賴關系的有向無環圖［10］，可以采用拓撲排序的方法得到模塊的先后順序，該先后順序就是指稱語義的計算順序。拓撲排序的思想是從入度為零的節點開始（對應于CBD 中的輸入模塊），斷開其與節點的連線，并將這些入度為零的節點添加到結果列表中，標識為已排序，然后去尋找圖中入度為零的節點，以此類推，直至所有節點都被訪問過。由于圖中沒有圈的存在，因此算法執行到最后一定能訪問完所有節點。

如圖3 所示，圖的入度為零的模塊（即輸入模塊）是A，出度為零的模塊（即輸出模塊）是E。拓撲排序算法首先從A 開始，刪除其指向B 和C 的連線，并將A 加入已排序的結果列表sorted。之后，B 模塊為唯一入度為零的模塊，算法訪問B 模塊，刪除其指向C 和D 的連線，并將B 加入已排序的結果列表。以此類推，算法執行完畢后得到圖的拓撲排序結果為［A，B，C，D，E］，且已遍歷完所有模塊。

Fig.3 Example of topological sort圖3 拓撲排序示例

操作語義所要完成的工作就是對于拓撲排序后的結果列表，在每個時間步內逐一遍歷并執行相應模塊的指稱語義。下面將以偽代碼的形式描述操作語義指導下的指稱語義計算過程，如算法1所示。

算法1 離散事件系統仿真算法

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

輸出：仿真結果

2.3 連續事件系統仿真

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

2.4 代數環

代數環（algebraic loop）［13］是指CBD 中模塊之間的依賴關系構成一個環，即形成模塊間的循環依賴。具體地，在仿真模型中，如果一個模塊的值取決于另一個模塊的值，而該模塊的值又依賴于第一個模塊的值，則會形成一個代數環［14］。這會導致仿真系統無法求解變量的值，因為每個模塊都需要已知其它模塊的值才能計算出自己的值，類似于操作系統中的“死鎖”。

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

Fig.4 Example of algebraic loop圖4 代數環示例

如前文所述，CBD 表達了模塊間數據的依賴性，因此可以采用拓撲排序的方式對基礎模塊進行排序［15］，得到的順序就是正確的計算順序。然而，由于代數環的存在，用戶建立好的CBD 模型并不能直接用于拓撲排序，而是需要經過一定的特殊處理，這種處理被稱為代數環的消解。

關于代數環的消解，目前學界有多種解決辦法，最常用的就是將代數環轉化為數學方程的形式進行矩陣運算求解［16］。但并非所有形式的代數環均有解，因此在建模過程中要盡可能避免代數環的出現。

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

Fig.5 Breaking loop according to non-direct feedthrough module圖5 根據非直接饋通模塊斷開環路

3 仿真系統設計

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

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

3.1 軟件總體設計

軟件采取分層架構，根據功能自頂向下分為數據交互層、數據解析層和仿真運算層。軟件體系結構如圖7所示。

Fig.7 Software architecture圖7 軟件體系結構

數據交互層主要負責模型接收和遠程模型調用，以及仿真結果的返回。接收到的數據將交給數據解析層，將結果返回給建模端進行展示。

數據解析層主要負責解析數據交互層接收到的數據，包括將模型文件解析為內存中的對象、狀態機語句值計算以及遠程模型文件解析。解析后的模型信息將傳給下一層仿真運算層進行仿真計算。

最底層是仿真運算層，也是整個仿真工具的核心。在這一層中，首先需要有基礎模塊庫支撐整個模型的搭建和仿真運行。其次，需要求解器支持連續模型的方程求解，即數值積分計算。當模型被解析完成后，會通過順序表構建、模型自檢進行仿真運行前的準備，之后即可運行仿真。如果涉及聯合仿真，將會與其他模型進行數據交互，實現模型之間的同步運行。后文將重點關注仿真運算層的設計與實現。

3.2 基礎模塊庫設計

基礎模塊庫是支撐建模與仿真的基礎。基礎模塊是指參與建模的最基本的計算單元（對應CBD 中的模塊），基礎模塊間通過有方向的連接線的連接構成模型圖［17］。

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

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

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

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

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

對于每一個具體模塊，可以定義其結構體，并為該結構體實現BasicOperation 接口，同時將BasicProperty 作為該結構體的屬性字段。對于模塊自身特有的行為和屬性，也可以為該模塊的結構體增添相應的屬性字段和方法。

對于模型中的每一個基礎模塊，完成實例化之后，需要存儲在一個集合中供仿真器訪問運行。然而，每一個基礎模塊都是一個不同的結構體，如何使用一個集合類型存儲不同的結構體類型是一個重要問題。在面向對象的語言中，可以考慮采用多態的形式——父類類型（或接口類型）指向子類實例的方式，以使用一種類型代表多種不同類型［18］。

3.3 仿真流程框架

離散與連續事件的仿真計算是指根據拓撲排序得到的基礎模塊排序序列依次遍歷每個基礎模塊，根據模塊的性質進行離散狀態更新或連續數值積分計算的過程。后文將針對這兩種計算的特點進行詳細介紹。

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

Fig.10 Flow of simulation calculation圖10 仿真計算流程

仿真開始時，模型會設定系統的初始狀態和輸出，然后在每個時間步中計算系統輸入、狀態和輸出的新值，并對模型進行更新（不涉及拓撲結構的更改）以應用這些計算結果。當仿真結束時，模型將輸出用戶所關心的模塊狀態及模塊的最終值。

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

對于包含連續狀態的模型，仿真器將使用求解器對模型進行求解。具體地，根據求解器的數值積分算法，求解器會進入子時間步循環。在子時間步中，模塊會以更小的仿真步長進行計算。在該循環中，求解器會重復調用模型的更新方法，使得在一個主時間步中可以按照連續時間間隔計算模型的輸出。使用該方法可以提高數值積分計算的準確度。

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

3.4 離散狀態模塊計算

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

3.5 連續狀態模塊計算

連續狀態系統通常采用微分方程的形式進行系統建模，因此連續狀態模型的求解依賴于數學上的數值積分算法。連續狀態求解并不是在真正的連續時間下進行計算，而是將時間分割成足夠密集的時間點去模擬連續計算。

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

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

（1）在每個時間步內，計算函數在當前仿真時間點的導數以求得斜率，并使用若干個仿真時間點與下一個時間點的中間點計算斜率，之后將上述斜率的加權平均值作為當前仿真步的斜率。

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

使用數學公式描述上述流程如下：

微分方程的初值表示為：

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

其中：

其中，h 代表仿真步長。通過上述計算，函數在下一時刻的值由當前值與當前值附近若干點的斜率和仿真步長所決定。

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

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

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

通過上述過程即可實現離散和連續系統的混合仿真。當需要增添新的解算器時，需根據實際的數值積分算法進行階段的微調，但總體計算順序不變。

4 測試與驗證

為了對仿真系統實例進行測試，本節根據基礎模塊框架分別搭建了3 種不同類型的基礎模型實例，并在Simulink 中搭建同樣的3 個模型，通過對比來驗證基礎模型仿真求解器仿真結果的可靠性與準確性。

4.1 離散系統仿真

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

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

Fig.12 Relationship model between linearized slip rate and adhesion coefficient圖12 線性化滑移率—附著系數關系模型

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

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

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

4.2 離散系統仿真

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

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

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

Fig.15 Model of bouncing ball system圖15 彈球系統模型

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

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

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

4.3 混合系統仿真

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

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

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

對于汽車理論速度的計算，其數學模型可由公式（7）表示。其中，x表示當前時刻汽車速度操縱機構的位置，vstandard表示根據汽車速度操縱機構位置計算出的汽車的理論速度。

對于汽車實際速度的計算，其數學模型可由公式（8）表示。其中，F表示汽車牽引力，m表示汽車質量，b表示摩擦阻力因子，vreal表示汽車的實際速度。

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

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

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

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

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

5 結語

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

本文仿真方法的優點在于采用了接口式設計，代碼各個模塊之間的耦合度較低，可以按照接口規范快速實現對基礎模塊庫和解算器的擴充，以提升模型的仿真能力。本文方法的局限性在于：對模型的仿真范圍在很大程度上依賴于基礎模塊庫，基礎模塊庫限制著模型的表達能力，而沒有提供一套可編程的形式供用戶自己編程實現相應基礎模塊的計算邏輯。