楊文強(qiáng)，吳文淵，陳經(jīng)緯，馮勇

（1.中國科學(xué)院重慶綠色智能技術(shù)研究院，重慶 400714；2.中國科學(xué)院大學(xué) 重慶學(xué)院，重慶 400714）

企業(yè)為了提高在全球化市場中的核心競爭力，提高產(chǎn)品的開發(fā)效率，降低產(chǎn)品的開發(fā)成本，優(yōu)化產(chǎn)品性能和結(jié)構(gòu)設(shè)計(jì)，計(jì)算機(jī)建模仿真技術(shù)將逐漸成為其必不可少的工具[1]。基于對仿真技術(shù)的需求，應(yīng)運(yùn)而生了數(shù)字化設(shè)計(jì)與制造技術(shù)，包含了基于云計(jì)算和互聯(lián)網(wǎng)的用戶需求挖掘與新產(chǎn)品開發(fā)策劃、計(jì)算機(jī)輔助設(shè)計(jì)（Computer Aided Design，CAD）、計(jì)算機(jī)輔助工藝設(shè)計(jì)（Computer Aided Process Planning，CAPP）、計(jì)算機(jī)輔助工程（Computer Aidedengineering，CAE）、計(jì)算機(jī)輔助制造（Computer Aided Manufacture，CAM）、數(shù)字樣機(jī)（Digital Mock-Up，DMU）等[2]，并結(jié)合產(chǎn)品“V”模式開發(fā)，能夠大幅縮短產(chǎn)品研發(fā)和測試的時(shí)間，有效避免不合理的設(shè)計(jì)造成的成本浪費(fèi)。

近年來，隨著計(jì)算機(jī)仿真技術(shù)的發(fā)展，通過建立數(shù)學(xué)模型和軟件算法相結(jié)合的形式來解決工程問題的方法得到了廣泛應(yīng)用。實(shí)現(xiàn)從設(shè)計(jì)到生產(chǎn)的過渡過程中軟、硬件結(jié)合的橋梁——硬件在環(huán)技術(shù)（Hardware in the Loop Simulation，HILs）逐漸變得至關(guān)重要。它利用計(jì)算機(jī)建模仿真技術(shù)對受控對象的運(yùn)行狀態(tài)進(jìn)行數(shù)字孿生，并依靠處理器實(shí)時(shí)輸出仿真結(jié)果，通過I/O 接口與控制器連接，對控制器進(jìn)行統(tǒng)的測試，有效地提高了系統(tǒng)的安全性、可行性，降低了長期設(shè)計(jì)成本[3]。新技術(shù)的提出，使仿真技術(shù)逐步發(fā)展出了離線（Off-line）仿真、在線（On-line）仿真及離線和在線相結(jié)合的仿真模式，極大地促進(jìn)和保證了虛擬還原現(xiàn)實(shí)的可靠性和準(zhǔn)確性。

1 仿真建模中的DAE 方程

當(dāng)前，伴隨著大數(shù)據(jù)的發(fā)展，作為世界制造工廠的我國，正面臨著數(shù)字化、智能化改革與創(chuàng)新的挑戰(zhàn)和機(jī)遇。新產(chǎn)品日新月異，更加廣泛地面向個性化定制的需求，在產(chǎn)品結(jié)構(gòu)和功能上更加復(fù)雜和豐富，諸如機(jī)器人、無人汽車等多智能體形式的現(xiàn)代高科技產(chǎn)品，由于其魯棒性、可靠性、高效性、可擴(kuò)展性等特性，在計(jì)算機(jī)網(wǎng)絡(luò)、機(jī)器人、制造業(yè)、交通控制、虛擬現(xiàn)實(shí)等方面得到廣泛應(yīng)用。與傳統(tǒng)的機(jī)電系統(tǒng)、電氣系統(tǒng)等不同，現(xiàn)代高科技產(chǎn)品是更高層面上的多個學(xué)科領(lǐng)域子系統(tǒng)集成于一體的復(fù)雜大系統(tǒng)，即多領(lǐng)域（multi-domain）耦合的復(fù)雜大系統(tǒng)。該系統(tǒng)具有多體系統(tǒng)、多物理場、多學(xué)科交叉融合的特征，是現(xiàn)階段制造技術(shù)發(fā)展的必由之路。要實(shí)現(xiàn)對現(xiàn)代高科技產(chǎn)品的設(shè)計(jì)，必然離不開利用計(jì)算機(jī)建模仿真技術(shù)對其進(jìn)行多領(lǐng)域的協(xié)同仿真，以實(shí)現(xiàn)產(chǎn)品整體性能的高精準(zhǔn)表述[1]。

1.1 仿真建模軟件的分類

目前，市場上仿真軟件主要分為專業(yè)的仿真軟件和通用的仿真軟件。專業(yè)的仿真軟件，如多體機(jī)械系統(tǒng)仿真的ADAMS、電磁場仿真的ANSOFT 和流體力學(xué)仿真的ANSYS 等[1]，由于其專注性太強(qiáng)，缺乏其他領(lǐng)域的可擴(kuò)展性，因此不適合應(yīng)用在多領(lǐng)域統(tǒng)一建模中。具有代表性的多領(lǐng)域統(tǒng)一建模通用仿真軟件主要有 MatlabSimulink，以及基于 Modelica 語言的SimulationX、MapleMaplesim、System Modeler、Dymola、Openmodelica、Mworks 等。其中，以MatlabSimulink軟件為代表的因果建模需要明確定義元件的輸入和輸出，而以開源的Modelica 語言軟件為代表的非因果建模則無需考慮元件之間的信號流。同時(shí)，MatlabSimulink 因其在矩陣計(jì)算中的突出特點(diǎn)，使用最為廣泛，同時(shí)兼具了計(jì)算能力強(qiáng)，計(jì)算效率高等特點(diǎn)，但建模過程復(fù)雜?；贛odelica 語言的非因果建模方便將模型投影到數(shù)學(xué)表達(dá)式，更加方便使用者直觀地理解，并且在符號計(jì)算層面上保障了計(jì)算的可靠性。特別地，還需要考慮到MatlabSimulink 并非開源軟件，Mathworks 公司于2020 年初受美國政府的干預(yù)對中國部分高校和科研機(jī)構(gòu)進(jìn)行了制裁[4]。

1.2 仿真建模求解流程

仿真建模求解流程基于Modelica 語言仿真器的編譯求解過程，大致分為編譯、分析優(yōu)化和仿真求解3 個階段[5]，見圖1。

圖1 基于Modelica 語言的仿真器求解過程Fig.1 Solving process based on Modelica language simulation

其中，編譯階段可進(jìn)一步分為詞法分析、語法分析、語義分析和平坦化4 個部分，主要完成將Modelica模型的程序語言源代碼進(jìn)行計(jì)算機(jī)數(shù)學(xué)表達(dá)方程化處理的過程，利用詞法規(guī)則、語法規(guī)則進(jìn)行程序驗(yàn)證和分析，提取出其中的方程信息特征，平坦化后得到一組方程及其對應(yīng)的常量、參數(shù)和變量。

分析優(yōu)化階段是利用數(shù)學(xué)方法完成對方程求解前的相容性分析和模型簡化的預(yù)處理工作。其中，相容性分析是利用二部圖等方法檢查是否完美匹配，即方程和變量關(guān)系是否相等，也就是系統(tǒng)是否為恰約束系統(tǒng)。同時(shí)，對過約束和欠約束問題，能夠協(xié)助用戶高效地實(shí)現(xiàn)問題分析和解決方案擬定，幫助實(shí)現(xiàn)模型的校正。模型簡化有利于方程的快速求解，涉及方程表達(dá)式的符號約簡、DAE 指標(biāo)約簡、方程系統(tǒng)的分塊化處理，以及轉(zhuǎn)化為狀態(tài)空間表達(dá)式等技術(shù)，得到一個DAE 系統(tǒng)的可順序求解的方程子集序列。

仿真求解階段與ODE 求解不同，DAE 求解中不可或缺的重要環(huán)節(jié)是對代數(shù)約束方程進(jìn)行一致性初值（Consistent Initial Value）的求解。求解器根據(jù)DAE方程子集求解序列的類型篩選出數(shù)值求解包中最優(yōu)的求解函數(shù)，并按照求解函數(shù)格式將方程子集求解序列自動生成基于C 語言的求解算法的可執(zhí)行代碼，通過編譯器進(jìn)行執(zhí)行并輸出結(jié)果。

Simulink 仿真求解過程與基于Modelica 語言的仿真器的求解過程類似，也可分為仿真建模、編譯、分析優(yōu)化和求解這4 個過程。Simulink 仿真的區(qū)別在于：仿真建模階段，Simulink 不僅需要確定傳送寬度、數(shù)據(jù)類型、計(jì)算塊參數(shù)以及分配內(nèi)存，還需要明確信號流方向和采樣時(shí)間；分析優(yōu)化過程中，為了更好地發(fā)揮矩陣計(jì)算的效率，需要用戶根據(jù)情況選擇合適的求解器對其進(jìn)行數(shù)值求解；在求解的過程中，計(jì)算是離散分塊的數(shù)據(jù)流迭代求解的過程。

目前，編譯階段及分析優(yōu)化階段的相容性分析和模型簡化都屬于計(jì)算機(jī)程序語言的設(shè)計(jì)工作，都已經(jīng)有了很好的解決方案。在指標(biāo)約簡及求解階段，雖然也有比較成熟的解決方案，但仍然存在理論上的漏洞和巨大的可發(fā)展和優(yōu)化的空間。

目前，在軟件的底層算法上，都需要將多領(lǐng)域物理模型平坦化為數(shù)學(xué)模型的微分代數(shù)方程，進(jìn)而對DAE 進(jìn)行求解。因此，在數(shù)學(xué)上，多領(lǐng)域統(tǒng)一建模的仿真問題就是對應(yīng)DAE 方程的求解問題[5]。

DAE 方程的一般形式可以描述為：

式中：t為自變量，t∈I；I為非空區(qū)間，I?R；x為n維因變量，x=x(t) = [x1(t), …,x n(t)]；x(k)表示x(t)的k階導(dǎo)數(shù)，1≤k≤l，k∈Z。

與ODE 不同的是，DAE 由含有x的最高次導(dǎo)數(shù)x(l)的微分方程和不含x的最高次導(dǎo)數(shù)x(l)的代數(shù)方程2 個部分構(gòu)成，故而DAE 的解取決于因變量x和它的導(dǎo)數(shù)，而不像ODE 的解那樣僅取決于x自身。

DAE 方程的求解實(shí)質(zhì)上是將DAE 方程中關(guān)于因變量x導(dǎo)數(shù)的隱含約束方程找出來，即進(jìn)行指標(biāo)約簡使其轉(zhuǎn)化為ODE，然后再對ODE 進(jìn)行數(shù)值求解。其中，涉及指標(biāo)約簡、一致性初始值、常微分方程數(shù)值精確求解等幾個關(guān)鍵技術(shù)。

2 DAE 指標(biāo)約簡

DAE 系統(tǒng)的指標(biāo)按照不同的定義和用途等分為微分指標(biāo)、結(jié)構(gòu)指標(biāo)、Kronecker 指標(biāo)、Perturbation指標(biāo)等[6]。雖然它們都在一定程度上代表了DAE 的求解困難和求解精度，指標(biāo)越高求解越困難，求解精度越差，但究其本質(zhì)最終目的都是求微分指標(biāo)。微分指標(biāo)是用來度量它與ODE 之間的“距離”，也就是將其轉(zhuǎn)換為ODE 系統(tǒng)而必須對其中部分或全部方程求微分的最小次數(shù)[7]。理論上，經(jīng)過恰當(dāng)次數(shù)的微分操作后，DAE 中所有次數(shù)導(dǎo)數(shù)都可以用新變量替換，得到等價(jià)表達(dá)的多項(xiàng)式系統(tǒng)[8]。高指標(biāo)（指標(biāo)≥2）的DAE 系統(tǒng)求解不可避免會用到數(shù)值微分運(yùn)算，不同于數(shù)值積分的情況，步長作為數(shù)值微分的分母項(xiàng)，減小步長會導(dǎo)致數(shù)值微分誤差的急劇增大，經(jīng)常會引起數(shù)值求解的不穩(wěn)定，造成結(jié)果錯誤，因此無法通過減小步長來實(shí)現(xiàn)誤差控制[9]?，F(xiàn)常用的DAE 求解器通常只能直接求解低指標(biāo)的問題或特殊形式的高指標(biāo)問題。

2.1 常見DAE 指標(biāo)約簡方法

對高指標(biāo)DAE 的數(shù)值求解問題，一般有2 種思路：一種是直接對DAE 進(jìn)行求解，但容易出現(xiàn)降階、漂移和不穩(wěn)定的現(xiàn)象；另一種是先將高指標(biāo)DAE 約簡為指標(biāo)為1 的DAE 或ODE，然后再對約簡后的系統(tǒng)方程進(jìn)行求解[6]。Gear 方法[10]、Pantelides 方法[8]、啞變量方法[11]、負(fù)權(quán)二部圖法[1]、加權(quán)二部圖法[12]、Pryce[13]方法是目前具有代表性的幾種指標(biāo)約簡方法，可對一般形式的DAE 系統(tǒng)實(shí)現(xiàn)指標(biāo)約簡。

Gear 方法是一種純符號的計(jì)算方法，反復(fù)地對代數(shù)方程進(jìn)行微分及部分低階微分變量替換，最終將DAE 系統(tǒng)轉(zhuǎn)化為ODE 系統(tǒng)。雖然該方法能完全實(shí)現(xiàn)指標(biāo)約簡，但符號算法復(fù)雜、計(jì)算效率低，且容易出現(xiàn)一些不必要的微分。

Pantelides 方法是一種結(jié)構(gòu)化分析方法，它通過分析方程系統(tǒng)結(jié)構(gòu)對應(yīng)的雅可比矩陣，尋找出其中奇異的最小結(jié)構(gòu)子集，然后對最小結(jié)構(gòu)子集微分，如此往復(fù)直到雅可比矩陣非奇異。該方法由于采用數(shù)值方法實(shí)現(xiàn)，與Gear 方法的符號方法相比，其計(jì)算量大幅降低，且避免了重復(fù)性微分，時(shí)間復(fù)雜度低，目前應(yīng)用在Dymola 和Maplesim 仿真軟件中。

啞變量方法通過對函數(shù)變量引進(jìn)新的啞變量來解決Pantelides 方法中動態(tài)變量的選擇問題，本質(zhì)上是Pantelides 方法的一種改進(jìn)版。

負(fù)權(quán)二部圖法和加權(quán)二部圖法是Pantelides 方法的補(bǔ)充方法，2 種方法相似，基于賦權(quán)二部圖的直接實(shí)現(xiàn)高指標(biāo)DAE 的約簡，且與Pantelides 方法具有相同的時(shí)間復(fù)雜度。

Pryce 方法是一種Pantelides 方法的等效方法，基于DAE 系統(tǒng)符號矩陣（Signature Matrix）的指標(biāo)約簡方法，通過求解對偶問題的整數(shù)規(guī)劃問題（Integer Linear Programming，IPL），能夠直接地獲取系統(tǒng)的最優(yōu)偏移量和篩選出所需指標(biāo)約簡的方程。

其中：DAE 系統(tǒng)F對應(yīng)的n×n維的符號矩陣σ中的元素構(gòu)造如下[13]：

對偶問題的整數(shù)規(guī)劃問題的形式如下：

式中：c為方程偏移量，c=[c1,…,cn]；d為導(dǎo)數(shù)偏移量，d=[d1,…,dn]，δ(F)為最優(yōu)解問題的最優(yōu)值。

2.2 DAE 指標(biāo)約簡修正方法

DAE 系統(tǒng)指標(biāo)約簡的成功與否取決于指標(biāo)約簡后的雅可比（Jacobian）矩陣是否奇異，如果不奇異，則指標(biāo)約簡成功；其中，雅可比矩陣為DAE 系統(tǒng)指標(biāo)約簡后因變量最高次導(dǎo)數(shù)x（l）的系數(shù)矩陣。值得注意的是，將DAE 系統(tǒng)方程經(jīng)過非奇異的線性矩陣A進(jìn)行線性重組處理，即A·F，或?qū)ζ渲心承┓匠踢M(jìn)行平方處理，即=0，形成新的DAE 系統(tǒng)，采用上述基于結(jié)構(gòu)的指標(biāo)約簡方法，如Pantelides 方法、Pryce 方法等，很可能會出現(xiàn)指標(biāo)約簡不成功的情況，究其本質(zhì)，也是由雅可比矩陣奇異導(dǎo)致的。對帶平方的DAE 進(jìn)行指標(biāo)約簡后，生成的ODE 系統(tǒng)中會存在冗余方程，影響數(shù)值求解的進(jìn)行。

雅可比矩陣奇異的情形可以分為2 類：一類是變量在求解區(qū)間的某些孤立的點(diǎn)上的奇異，但整個區(qū)間上是連續(xù)可積的，如例1 所示，這一類的奇異不會影響DAE 系統(tǒng)求解的質(zhì)量，通常在求解過程中可以通過坐標(biāo)系變換等方法來擺脫這一類奇異的影響；另一類則是退化，即在整個求解區(qū)間上都是奇異的。進(jìn)一步，退化問題又可以按照能否通過符號計(jì)算的方式約簡為0，分為符號退化（如例2 所示）和數(shù)值退化（如例3所示）2 類，其中，符號退化問題通常表現(xiàn)為結(jié)構(gòu)上可約導(dǎo)致的雅可比矩陣奇異，包含結(jié)構(gòu)奇異（Structurally Singularity）和完全奇異（Identically Singularity）。

例1，活塞缸的曲柄滑塊機(jī)構(gòu)（見圖2）如圖2所示，根據(jù)歐拉-拉格朗日方程分析得到該曲柄滑塊機(jī)構(gòu)的DAE 方程如下。

圖2 曲柄滑塊機(jī)構(gòu)Fig.2 Slider-crank mechanism

假設(shè)轉(zhuǎn)動慣量J1= 1、J2= 2，桿長l1= 1、l2= 2；當(dāng)燃料燃燒對活塞（滑塊）做功，推動活塞運(yùn)動，其運(yùn)動規(guī)律可控，滿足δ˙˙ - sint=0，此時(shí)有該DAE 的雅可比矩陣行列式值為 4cosθ1·sinθ2·cosθ2+ 4sinθ1·(cosθ2)2，當(dāng)2 根桿重疊時(shí)，曲柄滑塊機(jī)構(gòu)處于死點(diǎn)位置，雅可比矩陣奇異。此位置如果為初始啟動位置，勢必造成機(jī)構(gòu)無法啟動，仿真求解失敗，可以采取額外的措施增加方程數(shù)量修正雅可比矩陣的奇異；當(dāng)運(yùn)動過程中經(jīng)過該位置，則由慣性會繼續(xù)運(yùn)動，避免卡死，對位置的仿真通常會由計(jì)算誤差或步長而忽略不計(jì)。

例2，符號退化問題可以分為以下3 類。

1） 結(jié)構(gòu)奇異（不完美匹配）[14]

2） 結(jié)構(gòu)奇異（完美匹配）[15]

3） 完全奇異[14]——非線性單擺模型

對二部圖不具有完美匹配的結(jié)構(gòu)奇異問題，如例2 中1）所述，已經(jīng)不是一個單獨(dú)的DAE 系統(tǒng)了，多個系統(tǒng)之間相互獨(dú)立，應(yīng)對其單獨(dú)進(jìn)行仿真建模，而無共同仿真的實(shí)際意義。對二部圖具有完美匹配的結(jié)構(gòu)奇異問題，如例2 中2）所述，可以通過符號約簡的方法對其進(jìn)行預(yù)處理，消除其中的冗余部分[16]。

對于完全奇異問題，如例2 中3）所述，雖然雅可比矩陣J奇異，但仍然存在指標(biāo)約簡方法未能找出全部的隱含約束條件的情況，因此需要進(jìn)一步對其進(jìn)行修正。最早的修正方法是Murota[17]提出的通過組合松弛算法及其改良方法[18-19]，對線性的DAE 系統(tǒng)進(jìn)行重組表達(dá)。隨后，在其修正框架上發(fā)展出改進(jìn)的修正方法：線性組合方法（LinearCombination，LC）、表達(dá)式替換法（ExpressionSubstitution，ES）[15]、替換法（Substitution Method）和增廣法（Augmentation Method）[14]。

基于組合松弛方法中提出的修正框架概括為以下3 步。

1）計(jì)算等效整數(shù)規(guī)劃問題的最優(yōu)解(c,d)，如果該問題沒有解，則DAE 將無法實(shí)現(xiàn)完美匹配，算法將以失敗告終。

2）判斷雅可比矩陣是否奇異，如果不是，則返回ODE 系統(tǒng)F(c)。

大多數(shù)的退化問題都是符號退化的情形，只有少數(shù)，例如在含參數(shù)的微分代數(shù)方程中，參數(shù)在取值范圍變化會引起雅可比矩陣的變化，當(dāng)取到某些特殊參數(shù)時(shí)，難以避免會出現(xiàn)數(shù)值退化問題，此時(shí)，現(xiàn)有的修正方法將不能勝任。

圖3 共線彎矩作用下梁的疊加變形Fig.3 Superposition deformation of beams under collinear bending moments

2.3 一致性初值問題

將DAE 系統(tǒng)進(jìn)行指標(biāo)約簡后，只是完成了對隱含約束方程部分的處理，但還缺乏對初值問題或邊界值問題相關(guān)約束條件的描述，因此還需要解決一致性初值問題，才能將DAE 系統(tǒng)完全地轉(zhuǎn)換成為一個ODE 系統(tǒng)，同時(shí)，結(jié)合例3，可知初值問題的求解有利于數(shù)值退化問題的分析。

在仿真建模中，系統(tǒng)的初始狀態(tài)通?？梢圆煌?，如曲柄滑塊機(jī)構(gòu)在正反方向上都可以進(jìn)行啟動，但其運(yùn)動的狀態(tài)——拉伸或壓縮，卻是截然不同的。如果只是對一個初始方向進(jìn)行分析和優(yōu)化，則所得的結(jié)果都是局部的而非全局的。

例如，對代數(shù)方程形如雙曲線方程x2-y2-1 =0 這一類，如果采用牛頓迭代等數(shù)值方法對其進(jìn)行求解，勢必會造成部分一致性初值信息的丟失。而這些信息有可能會關(guān)系雅可比矩陣的奇異性，如果單純的放任不管，則可能導(dǎo)致DAE 系統(tǒng)與一致性初值的不匹配，造成求解失敗。

針對類似的多項(xiàng)式方程求解問題，常采用Sommese、Wampler、Versheclde 等提出的數(shù)值代數(shù)幾何學(xué)[20-21]的同倫延拓求解方法（Homotopy Continuation Mothed）來進(jìn)行處理。對于非線性較強(qiáng)（如含三角函數(shù)）的情況，該方法仍然無法實(shí)現(xiàn)求解。

3 ODE 求解器

在機(jī)械動力學(xué)和控制工程等領(lǐng)域的模擬仿真中，系統(tǒng)在某時(shí)刻后的狀態(tài)變化，通常被描述為線性常微分方程的初值問題。如果常微分方程在時(shí)間t上連續(xù)，在x上Lipschitz 連續(xù)，那么它在邊界條件的領(lǐng)域范圍內(nèi)一定存在唯一的精確解*x。目前，在線性常微分方程的初值問題和[22]邊界值問題[23-24]上，已有非常成熟的數(shù)值求解方法和p階誤差的分析方法，并且開發(fā)了很多面向應(yīng)用的求解器，例如Matlab[25]和Maple[26]等。據(jù)了解，關(guān)于如何嚴(yán)格給出由求解器數(shù)值解插值得到的近似解與精確解的距離上界，以及如何在數(shù)值解的基礎(chǔ)上提高近似解的質(zhì)量等方面的工作還比較少。

3.1 ODE 誤差估計(jì)方法

對誤差估計(jì)而言，通過殘差對其進(jìn)行界定具有許多優(yōu)點(diǎn)，最大的優(yōu)點(diǎn)就是它能夠用來評估整個求解區(qū)間內(nèi)的誤差，而不只是每個網(wǎng)格格點(diǎn)上的局部誤差[27]。研究人員通常喜歡將求得的數(shù)值解代回原方程求出殘差[28-30]，再進(jìn)行誤差估計(jì)。Constantinescu[31]在文章中提出利用時(shí)間步進(jìn)策略結(jié)合殘差給出一個p階“全局誤差”的估計(jì)方法，能夠更好地提高誤差估計(jì)的準(zhǔn)確性，但本質(zhì)上只是在格點(diǎn)上的誤差估計(jì)，而且舍去了高階項(xiàng)，不能說是準(zhǔn)確的誤差估計(jì)。Enright[32]在文章中提出了合理的估計(jì)殘差范數(shù)方法，但沒有考慮每個點(diǎn)導(dǎo)數(shù)的數(shù)值誤差，這些導(dǎo)數(shù)的不準(zhǔn)確進(jìn)一步導(dǎo)致殘差定義的不準(zhǔn)確，從而導(dǎo)致誤差估計(jì)的不準(zhǔn)確。另一些文獻(xiàn)[28,30,33]在考慮導(dǎo)數(shù)問題的誤差估計(jì)方面做了大量工作，但其中大部分都在努力解決局部誤差問題。

一般情況下，人們是無法求出常微分方程的精確解的，但幸運(yùn)的是，可以通過很多數(shù)值方法求解器求得數(shù)值解，如Matlab 等，可以進(jìn)一步減小步長來逼近精確解，基于此，很多人在這方面做了重要的研究工作[34-36]。無論步長取得多小，局部誤差控制得多好，人們?nèi)匀粺o法回答數(shù)值解和精確解之間到底相差多少這個問題。全局誤差估計(jì)方法[37-40]能夠部分地回答這個問題，但是這些方法計(jì)算量太大且條件過于苛刻。文獻(xiàn)[41]中提出利用Lipschitz 常數(shù)來協(xié)助誤差上界的估計(jì)，給出了誤差估計(jì)上界：

式中：h為步長；M為的估計(jì)上界；L為Lipschitz 常數(shù)，由于Lipschitz 常數(shù)為正，導(dǎo)致估計(jì)上界隨時(shí)間增長呈指數(shù)增長，且增長速率過大，不具有實(shí)用性。文獻(xiàn)[29]通過推導(dǎo)條件數(shù)和殘差的關(guān)系來估計(jì)全局誤差，但求解條件數(shù)需要對齊次常微分方程的基本解（Fundamental Solution）進(jìn)行求解，在多數(shù)情況下這一點(diǎn)幾乎無法實(shí)現(xiàn)。專門針對穩(wěn)定系統(tǒng)的常微分方程，文獻(xiàn)[42]提出利用Lyapunov 定理構(gòu)造全局誤差控制方程來限制誤差上限，能夠起到誤差估計(jì)和提高解的精度的作用，但構(gòu)造全局誤差控制方程也需要基本解的求解。文獻(xiàn)[43]中，提出了一種區(qū)間求解的方法，僅適用于初值問題求解，該方法基于泰勒級數(shù)展開，估計(jì)的范圍和計(jì)算精度取決于級數(shù)展開的項(xiàng)數(shù)，精度提高會造成計(jì)算量急劇增大；同時(shí)，經(jīng)實(shí)驗(yàn)驗(yàn)證，區(qū)間求解方法需要大量的預(yù)處理時(shí)間；且在相同的步長和階數(shù)的情況下，區(qū)間求解方法的精度要低于龍格庫塔等數(shù)值求解算法。文獻(xiàn)[44-45]利用殘差和精確解之間的關(guān)系，通過理論推導(dǎo)和證明構(gòu)造了誤差估計(jì)上界，能夠準(zhǔn)確地實(shí)現(xiàn)線性O(shè)DE 初值問題的誤差上界的估計(jì)，避免了Lipschitz 常數(shù)的指數(shù)增長問題和求解基礎(chǔ)解的問題，但針對非線性和邊界值問題還缺少進(jìn)一步的研究。

在數(shù)值求解中，解的精度通常取決于步長，但更小的步長意味著更大的計(jì)算量，如何在固定步長的情況下減小全局誤差是一個很有難度的問題，目前還沒有太多相關(guān)的研究?！叭毕菪拚保―efect Correction）方法[46-47]，類似于牛頓迭代方法求解非線性方程組，將殘差代回方程進(jìn)行迭代數(shù)值求解，可用于非線性常微分方程，但它們都只是考慮了在網(wǎng)格點(diǎn)上的局部誤差。同時(shí)，在計(jì)算效率上，雖然現(xiàn)有的數(shù)值求解方法都很快，但仍然可以考慮利用技術(shù)手段將計(jì)算效率進(jìn)一步提升。

3.2 ODE 數(shù)值求解方法

現(xiàn)有的ODE 系統(tǒng)的求解器有很多，大致可以分為3 類：數(shù)值求解器、符號求解器和區(qū)間求解器。數(shù)值求解器主要基于歐拉法及其改進(jìn)格式、龍格庫塔格式、亞當(dāng)斯格式等[48]，其優(yōu)勢在于計(jì)算效率高，目前應(yīng)用最為廣泛，常見的為Matlab 提供的求解器，針對不同類型的方程與不同問題，比較典型的有針對常微分方程的ode45、針對時(shí)滯微分方程的dde23、針對邊值問題的bvp4c 和bvp5c 等。符號求解器是基于微分方程通解公式，主要應(yīng)用在Malpe 的求解器中，比較典型的求解器有dsolve、linearsol 等。目前只能針對部分格式的常微分方程進(jìn)行求解，如齊次方程、線性方程、Bernoulli 方程、Riccati 方程等格式。

Nedialkov 等[43]的區(qū)間求解方法需要采用高階的泰勒級數(shù)展開對微分方程進(jìn)行預(yù)處理，然后才能夠得到網(wǎng)格點(diǎn)所在的一個時(shí)間區(qū)間上的解區(qū)間，該網(wǎng)格點(diǎn)上的精確解屬于該解區(qū)間，目前應(yīng)用在VNODE-LP區(qū)間求解器。由于沒有將數(shù)值解插值作為近似解，因而區(qū)間求解方法只能給出網(wǎng)格點(diǎn)所在區(qū)間上解得的誤差估計(jì)區(qū)間，無法給出全局的誤差估計(jì)。

有限差分方法是一種應(yīng)用范圍更廣的數(shù)值求解方法，該方法不僅能夠用來求解一般的ODE、DAE、偏微分方程（Partial Differential Equation，PDE），且不區(qū)分初值問題與邊界值問題，還能夠用來求解ODE或DAE 相關(guān)的優(yōu)化問題。其基本思想是通過對定義域進(jìn)行網(wǎng)格劃分，將原問題在網(wǎng)格點(diǎn)上進(jìn)行離散化，然后用代數(shù)差分方程組去逼近原來的微分方程，同時(shí)需滿足對應(yīng)原問題的在離散點(diǎn)處的邊界條件，典型的有限差分方法包括配置法和有限元法[49]。

文獻(xiàn)[45]利用殘差二范數(shù)和精確解插值相結(jié)合的符號推導(dǎo)方法，給出了一種線性O(shè)DE 和部分非線性O(shè)DE 的優(yōu)化求解方法，同樣能夠解決初值問題和邊界值問題的求解。該方法首先需要在符號層面上完成ODE 方程的預(yù)處理后才可賦值進(jìn)行計(jì)算。雖然采用離線和在線的方法能夠極大地加速特定類別方程的計(jì)算效率，但針對非線性強(qiáng)的ODE 方程求解還需要進(jìn)一步研究。

隨著深度學(xué)習(xí)技術(shù)的發(fā)展，利用神經(jīng)網(wǎng)絡(luò)架構(gòu)求解常微分方程[50]、偏微分方程[51]成為了可能，這也是目前比較熱的研究領(lǐng)域，是一個值得研究的方向。Chen 首先提出了將殘差網(wǎng)絡(luò)連續(xù)變化為微分方程的思想，能夠高精度地對常微分方程的初值問題進(jìn)行求解，以及在此基礎(chǔ)上演變而來的[52-53]，并非通用的求解模型，針對不同類型的ODE 都需要采用訓(xùn)練集對神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練，效率較低。同理，基于數(shù)據(jù)驅(qū)動的回歸分析方法[54-55]來求解偏微分方程也存在這樣的問題。

4 DAE 系統(tǒng)中的參數(shù)設(shè)計(jì)

當(dāng)前的仿真建模軟件還處于發(fā)展的初期階段，大多數(shù)情況下只是考慮如何完善對一般形式的DAE 系統(tǒng)的求解功能，而對含參數(shù)的DAE 系統(tǒng)問題涉及較少。誠然，一般形式的DAE 系統(tǒng)的求解功能是含參數(shù)的DAE 系統(tǒng)問題求解的基石。對產(chǎn)品進(jìn)行優(yōu)化設(shè)計(jì)、對儀器設(shè)備進(jìn)行優(yōu)化控制、對原材料和制備過程進(jìn)行優(yōu)化調(diào)節(jié)，是提高行業(yè)競爭力的最直接方法。如汽車領(lǐng)域，在設(shè)計(jì)時(shí)可以利用仿真器優(yōu)化行星齒輪機(jī)構(gòu)特征參數(shù)，使之與動力部件的最高轉(zhuǎn)速和整車設(shè)計(jì)最高車速匹配，使發(fā)動機(jī)工況得到了改善，整車燃油經(jīng)濟(jì)性和排放都得到大幅提升[56]。究其本質(zhì)，就是對DAE 系統(tǒng)進(jìn)行參數(shù)設(shè)計(jì)，即參數(shù)優(yōu)化，最常見的還有對動力系統(tǒng)的PID 參數(shù)調(diào)節(jié)。而對產(chǎn)品參數(shù)化設(shè)計(jì)，如果采用傳統(tǒng)的“V”模式開發(fā)，勢必會造成時(shí)間和金錢的浪費(fèi)，因此采用參數(shù)優(yōu)化方法結(jié)合仿真建模的方式會是未來仿真建模領(lǐng)域發(fā)展的主要趨勢。

對含參數(shù)DAE 系統(tǒng)而言，目前應(yīng)用最多的優(yōu)化方法都是局部優(yōu)化法，具體包括：利用極大值原理將含參數(shù)的DAE 模型一般化，然后再利用DAE 求解器進(jìn)行求解的變分法；將優(yōu)化含參數(shù)的DAE 模型為非線性規(guī)劃模型，進(jìn)而進(jìn)行優(yōu)化求解的聯(lián)立法；有限差分方法[49]。采用局部優(yōu)化方法收斂得到的局部最優(yōu)很可能是一種假象，實(shí)際情況是人們可能從來沒有找到過“局部最優(yōu)”。

針對含參數(shù)DAE 系統(tǒng)的全局最優(yōu)化方法方面的研究很少，可以參考在含參數(shù)ODE 問題求解的確定性方法和隨機(jī)性方法。隨機(jī)性方法（如啟發(fā)式方法），通過不斷搜索的方式，逐個比較極值，力求遍歷全局，但計(jì)算量大、效率低下，典型的有蟻群算法[57-58]、退火算法[59-60]、粒子算法[61]等。這類方法沒有嚴(yán)格的理論支撐，能夠收斂到極值，但不能保證獲得的最優(yōu)解就是全局最優(yōu)解。確定性方法是基于函數(shù)的解析性質(zhì)（如凸性、稠密性、單調(diào)性和Lipschitz 等），然后利用這些性質(zhì)確定含參數(shù)DAE 系統(tǒng)的界，進(jìn)一步確定全局最優(yōu)解[62]，主要有分支定界法、整數(shù)全局尋優(yōu)、聚類和隧道方法等。分支定界算法是將原問題通過松弛技術(shù)轉(zhuǎn)化為簡單的凸優(yōu)化問題[63]，已經(jīng)實(shí)現(xiàn)部分含參數(shù)的ODE 系統(tǒng)求解中的應(yīng)用[62-64]，但其在計(jì)算過程中需要求解ODE 系統(tǒng)的基本解來確定凸函數(shù)的上界和下界，因此只能針對特定的ODE 優(yōu)化系統(tǒng)進(jìn)行求解。文獻(xiàn)[65]中提出的分支定界方法框架下新的凸松弛方法能夠處理多項(xiàng)式形式的廣義幾何規(guī)劃等問題，但其僅限于特定形式的多項(xiàng)式問題。目前，關(guān)于分支定界方法的研究主要集中在凸松弛和分支策略上，其目的是通過緊度來保證解的質(zhì)量，通過裁剪分支來保證解的效率；其中，最重要的環(huán)節(jié)就是凸松弛技術(shù)，包括：線性松弛[66]、凸二次松弛[67-68]、拉格朗日松弛[69]、二階錐松弛[70]，以及半正定松弛[71]等。

5 未來研究方向

多領(lǐng)域統(tǒng)一建模仿真軟件的目標(biāo)就是通過計(jì)算機(jī)技術(shù)實(shí)現(xiàn)虛擬現(xiàn)實(shí)模擬及數(shù)字孿生，協(xié)助開發(fā)人員便捷、可靠、有效地進(jìn)行新產(chǎn)品的開發(fā)設(shè)計(jì)。為了更好地實(shí)現(xiàn)這一目標(biāo)，必須要滿足DAE 求解的可靠，而根據(jù)之前所述內(nèi)容可知，現(xiàn)有DAE 求解技術(shù)存在諸多局限，對仿真軟件求解器進(jìn)行優(yōu)化，構(gòu)建一個更加精確和完備的DAE 求解架構(gòu)，顯得非常有必要。未來在這一領(lǐng)域的研究可以歸納為以下幾個方向。

1）在現(xiàn)有的仿真模型、編譯器和分析優(yōu)化器等環(huán)節(jié)引入模型的校驗(yàn)、驗(yàn)證和確認(rèn)（Verification，Validation and Accreditation，VV&A）技術(shù)，提高DAE仿真模型的正確性，解決其中的過約束和欠約束方程問題，這樣能夠提高DAE 求解的可行性、真實(shí)性和可靠性。值得考慮的是：可以采用硬件在環(huán)技術(shù)、機(jī)器學(xué)習(xí)技術(shù)（如對抗生成學(xué)習(xí)網(wǎng)絡(luò)、強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)等）來實(shí)現(xiàn)；特別是VV&A 技術(shù)在針對大型復(fù)雜的DAE仿真系統(tǒng)時(shí)，在精度和置信度的提高上，可以做進(jìn)一步的研究。

2）雖然高指標(biāo)DAE 系統(tǒng)的指標(biāo)約簡技術(shù)已經(jīng)較完善，但在DAE 系統(tǒng)出現(xiàn)退化（尤其是數(shù)值退化）的情形時(shí)，以及在DAE 中存在方程的線性重組或平方的情形時(shí)，現(xiàn)有的指標(biāo)約簡及其修正方法可能會失效，亟需對現(xiàn)有技術(shù)進(jìn)行修正或補(bǔ)充。

3）一致性初值問題是求解DAE 系統(tǒng)的關(guān)鍵問題之一，對于單分支的情況，采用數(shù)值求解約束條件方程是完全能夠勝任的，但對于多分支的情況，有必要遍歷每個分支，考慮DAE 多個狀態(tài)的初值，數(shù)值求解方法只能得到個別的解。目前可以采用的解決方法有同倫連續(xù)方法等，但缺乏其在DAE 一致性初值上的相關(guān)理論推導(dǎo)和證明，以及在非多項(xiàng)式約束條件方程會失效。

4）針對ODE 數(shù)值求解，現(xiàn)有仿真軟件中的數(shù)值求解器已經(jīng)非常成熟，但對非線性程度高的大規(guī)模復(fù)雜模型上還可以進(jìn)一步研究，可以考慮結(jié)合神經(jīng)網(wǎng)絡(luò)技術(shù)；全局誤差可控計(jì)算在高精度過程控制中有非常廣的應(yīng)用前景，但對其的研究還很少，目前的研究局限在近似線性O(shè)DE 系統(tǒng)，且符號計(jì)算不利于大規(guī)模求解，短期來看，可考慮利用泰勒級數(shù)展開等技術(shù)推廣到非線性的復(fù)雜模型上。

5）在ODE 解的誤差估計(jì)上，目前在局部誤差和全局誤差的估計(jì)上存在很多有效可行的方法，但全局誤差估計(jì)方法的研究還局限在線性的初值問題上，可考慮利用線性化技術(shù)或初值轉(zhuǎn)換技術(shù)推廣應(yīng)用到非線性和邊界值問題的情形，這也是一個不錯的研究方向。

6）含參數(shù)DAE 系統(tǒng)的優(yōu)化問題任重而道遠(yuǎn)，特別是全局參數(shù)優(yōu)化問題。短期來看可以利用指標(biāo)約簡技術(shù)，將DAE 方程轉(zhuǎn)化為ODE 方程，再利用現(xiàn)有的優(yōu)化方法（如分支定界法）實(shí)現(xiàn)特殊形式的微分代數(shù)方程的參數(shù)優(yōu)化（如低指標(biāo)、線性、常系數(shù)等），這類問題在控制工程中有著廣泛的應(yīng)用前景。長期來看，需要在理論上進(jìn)行創(chuàng)新，如在Singer[64]研究的基礎(chǔ)上，避開基礎(chǔ)解的求解，從根源上解決這一類問題。

6 結(jié)語

仿真建模的關(guān)鍵就是建立DAE 系統(tǒng)和求解DAE系統(tǒng)，前者是產(chǎn)品設(shè)計(jì)開發(fā)人員根據(jù)產(chǎn)品的物理、數(shù)學(xué)和經(jīng)驗(yàn)等構(gòu)建的有機(jī)結(jié)合，隨著模型的改變而改變；后者是仿真建模軟件的核心，其中涉及指標(biāo)約簡技術(shù)、數(shù)值求解技術(shù)、誤差分析技術(shù)及優(yōu)化求解技術(shù)等?，F(xiàn)有的指標(biāo)技術(shù)和基于此的修正方法能夠滿足絕大多數(shù)DAE 系統(tǒng)轉(zhuǎn)化為ODE 系統(tǒng)的需求，但針對特殊形式（如數(shù)值退化問題和多分支一致性初值問題）還需要進(jìn)一步研究。同時(shí)，現(xiàn)有的數(shù)值求解技術(shù)也相當(dāng)成熟，用戶可以根據(jù)需求去合理地選取，但在連續(xù)區(qū)間誤差函數(shù)可積的情況下，全局誤差最小化也是一個可以研究的方向。其次，現(xiàn)有的局部誤差估計(jì)方法和全局誤差估計(jì)方法都比較準(zhǔn)確，如何進(jìn)一步控制誤差估計(jì)上限，以及針對邊界值、非線性問題進(jìn)行準(zhǔn)確誤差估計(jì)，同樣是一個不錯的研究方向。最后，含參數(shù)的DAE 問題或含參數(shù)的ODE 問題的求解仍然存在諸多的困難和挑戰(zhàn)，關(guān)系到產(chǎn)品個性化發(fā)展的未來。