線性多步量化狀態系統算法及在溢流閥上應用

傅東金，樊志華，沈漢武，李志華

(杭州電子科技大學機械工程學院，浙江 杭州 310018)

0 引 言

常微分方程(Ordinary Differential Equation, ODE)系統廣泛應用于機械、電子、控制、航空航天等工程領域，是建模仿真、數值計算、高性能計算等領域的重要研究方向之一。近年來，國家十分重視工業軟件的開發，其中求解算法是工業軟件的底層核心[1]。實際工程系統中的常微分方程常表現出剛性問題，即系統中同時存在快速與慢速變化量，剛性問題一直是傳統數值積分方法難以解決的問題[2]。溢流閥的仿真計算中，由于閥門容腔壓力巨變、油液壓縮和多變量耦合等問題，導致溢流閥的動力學模型十分復雜，呈現高度的非線性和剛性，仿真求解十分困難[3]，而高效精確的求解溢流閥動力學問題是其設計的基礎[4]。幾十年來，盡管提出了不少求解剛性問題的方法，如隱式龍格庫塔法、Gear法及其他一些特色方法[5-7]，但這些都是基于時間離散的方法，其算法性能的提升受到一定的限制。量化狀態系統(Quantized State System, QSS)[8]是一種在計算過程中認為狀態變量是不連續的算法，其邏輯與傳統時間離散的積分方法不同，在QSS中可以通過量化函數將狀態變量以“量子”大小進行分割，并行計算每個狀態變量狀態改變所需要的時間，這一過程稱之為“躍遷”，再以最小躍遷時間為統一步長，采取個體異步、整體同步的時間推進思想，從而達到推進積分的目的，是一種并行式事件驅動型算法[9]。在求解一般ODE方程時，由于在求解過程中不進行迭代，故而其求解速度得到大幅提升，并且因為其本身的求解特性導致它天然具有強穩定、高精度等優點。

在QSS基礎上，又提出了后向量化狀態系統(Back Ward Quantized State System, BQSS)、中心量化狀態系統(Centered Quantized State System, CQSS)等算法[10-12]，在不同工程領域取得了一定效果，但這些算法依然存在一些問題。首先這些算法沒有考慮多變量耦合情況下多狀態變量并行計算躍遷條件與單變量計算躍遷條件的差異，導致積分難以推進；其次這些算法在每一積分步均以量子大小進行躍遷，求解精度不高；最后已有的QSS系列算法均是使用單步法，沒有充分利用已有的信息來提高算法的求解性能。目前，國內學者對QSS研究還較少，比如文獻[13-15]對QSS方法進行了初步應用。本文以1個典型剛性算例和1個工程實例為研究對象，借鑒QSS方法和Adams多步法，提出一種線性多步量化狀態系統算法(Linear Multi-Step Quantized State System, LMSQSS)，以提高剛性方程數值求解的精度和效率。

1 量化狀態系統

一組由常微分方程表示的狀態方程系統如下：

(1)

式中，x(t)為狀態向量，x(t)∈R，u(t)為輸入向量，QSS方法將式(1)量化為：

(2)

式中，q(t)為狀態變量的量化向量，每個量化變量qj(t)通過量化函數得到：

(3)

式中，ΔQj為量子，qj(t-)為前一次的量化變量值。

量子由初始設定，量子的大小與算法最終求解誤差呈正相關，文獻[16]給提出了選擇量子的方法。

圖1 狀態變量xj(t)和量化變量qj(t)的軌跡圖

QSS的狀態變量xj(t)和量化變量qj(t)的軌跡如圖1所示，從圖1可以看出，qj(t)軌跡為分段常數，當xj(t)的變化超過了1個“量子”時，引起qj(t)的變化，這個行為稱之為“躍遷”，這一特性保證了算法的穩定性和全局誤差界限，因為qj(t)與xj(t)的差值不會超過1個量子。

從圖1還可以看出，任意軌跡的狀態變量與相對應的量化變量共存在9種狀態，而QSS未考慮到多變量耦合情況下單獨狀態變量躍遷條件與非耦合情況下的不同，只考慮圖1中的1，2，5，6，8這5種狀態[17]，本文算法完善后的躍遷條件表示為如下9種狀態：

2 線性多步量化狀態系統

本文所提的LMSQSS算法通過顯式形式表達，避免了迭代求解，提高了求解效率；并利用求解過程中已知信息進行修正，得到更加精確的結果；LMSQSS算法本身包含隱式求解思想，將量化變量值作為系統狀態變量預測值，進一步提升了算法的穩定性和精確性。

為了抑制仿真求解過程中的振蕩，需要保證積分推進過程中xj的取值一直向著qj取值的方向靠攏，故qj的量化函數為：

(4)

(5)

根據不同躍遷速度確定躍遷時間為：

(6)

從而求出每推進一次仿真所需的時間為：

Δt=min{Δtj}

(7)

此時，狀態變量的仿真推進時間為：

ti+1=ti+Δt

(8)

躍遷狀態變量為：

(9)

非躍遷狀態變量為：

(10)

式中，n表示多步法的階數。

(11)

表1 各階數對應的系數αk

LMSQSS算法具有不同的階數，各階系數通過Adams-Bashforth法得到，各階數對應的系數如表1所示。

3 仿真算例實驗

仿真實驗在Windows7系統上的MATLAB2012a平臺上進行，處理器為Intel i5-4260。分別采用LMSQSS算法、ODE45，ODE15s，ODE23s，QSS算法進行實驗。ODE45是傳統求解的首選算法，是顯式變步長單步Runge_kutta法，當解的光滑性比較好時，其精度比較高；ODE15s和ODE23s都是針對剛性問題的優秀求解算法，ODE15s是隱式Gear變步長多步求解法，ODE23s是隱式Rosenbrock變步長單步求解法。為了研究不同階數對LMSQSS算法性能的影響，實驗選取2階、3階和4階階數進行仿真求解。對所有算法設置相同誤差容限。

系統各狀態變量相對誤差計算公式如下：

(12)

式中，k=1,2,…,n，um[k]為各個算法求得的解，u[k]為DASSL求解器在很小誤差設定(10-9)下求得的解，作為基準值用于計算相對誤差[18]。

一個線性二階ODE系統[19]如下：

(13)

圖2 LMSQSS算法對算例的仿真結果

表2 不同算法對算例1仿真結果

從圖2中可以看出，隨著仿真時間的推進，變量x1和x2并未出現振蕩或發散。

從表2可以看出，相比ODE45，2階LMSQSS算法的效率提高了約435倍，精度提高22倍；相比傳統剛性求解器中最好的ODE15s，效率提高了2倍，精度提高了69倍；相比ODE23s，效率提高了5倍，精度提高了82倍。在求解剛性二階系統問題時，QSS算法的結果產生發散。

由表2可知，在相同量子情況下，隨著階數的增加，LMSQSS算法的求解精度明顯增高，而CPU仿真時間變化并不大。

1—主閥左腔；2—主閥右腔彈簧；3—主閥右腔；4—過濾器；5—先導閥彈簧。圖3 溢流閥結構簡圖

4 工程實例

4.1 工程實例模型

一個先導式溢流閥[20]的結構簡圖如圖3所示。圖3中，油液從主閥左腔流入，經主閥右腔和過濾器到達先導閥，在油液壓力到達先導閥彈簧所設定壓力值后，先導閥開啟，油液流入，油液再經上部回油管路流回，致使主閥左右腔形成壓力差，直到壓力能夠克服主閥彈簧彈力，使主閥芯移動，油液再經主閥口流出。采用功率鍵合圖法對其進行建模[21-22]，得到數學模型的微分代數方程如下：

(16)

仿真初值：[x1,x2,x3,x4,x5,x6]=[V1,V2,P1,P2,Y1,Y2]=[0.1,6.37×10-3,0,0,0,0]，仿真時長為0.5 s。其中，V1為泵至主閥左腔管道容腔中用來補償油液壓縮量及管道受壓變形量的液壓油總容積，單位為cm3；V2為主閥右腔中用來補償油液壓縮量的液壓油總容積，單位為cm3；P1和P2分別為主閥芯和先導閥閥芯的動量，單位為N·s/cm2；Y1和Y2分別為主閥芯和先導閥閥芯的位移，單位為cm。其余參數在表3中進行說明。

表3 溢流閥參量值

4.2 工程實例仿真結果分析

采用LMSQSS算法對溢流閥數學模型進行求解，得到圖4閥芯位移軌跡圖和圖5入口壓力圖。從圖4可以直觀看到溢流閥的閥芯位移變化情況，其中主閥芯位移軌跡很平穩，并未因為系統耦合而發生振蕩；先導閥軌跡出現短暫波動后(0.02 s)，趨于平穩。從圖5中可以看出，閥的入口壓力在油液進入后開始急劇上升，其動態超調量在22%左右，上升時間在0.002 0 s左右，峰值出現的時間在0.003 5 s左右。在先導閥開啟后，壓力迅速下降，短暫波動后，在0.006 0 s后趨于穩定。由此可見，LMSQSS算法具有較好的穩定性和收斂性，能夠勝任溢流閥動力學的實際求解。

圖4 閥芯位移的仿真結果

圖5 入口壓力的仿真結果

不同算法對溢流閥模型的求解結果如表4所示。

表4 不同算法對溢流閥仿真結果

從表4中可以看出，對比傳統方法中結果最好的ODE15s，2階LMSQSS算法的仿真精度提高了42倍，仿真效率提高了12倍。ODE45和QSS算法在求解該問題時出現發散，不能得到有效解。還可以看出，在相同量子情況下，隨著階數的增加，LMSQSS算法的求解精度明顯增高，而CPU仿真時間變化并不大。

5 結束語

針對剛性ODE求解難題，本文提出一種有新的線性多步量化狀態系統算法LMSQSS。對QSS躍遷理論不足之處加以完善，擴大了算法的應用范圍。相比于傳統求解方法，無論在求解精度還是求解效率上，LMSQSS都更具優勢。下一步計劃研究如何將算法實現并行求解，進一步提高求解效率。