線性多步量化狀態(tài)系統(tǒng)算法及在溢流閥上應(yīng)用

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

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

0 引 言

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

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

1 量化狀態(tài)系統(tǒng)

一組由常微分方程表示的狀態(tài)方程系統(tǒng)如下：

(1)

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

(2)

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

(3)

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

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

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

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

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

2 線性多步量化狀態(tài)系統(tǒng)

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

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

(4)

(5)

根據(jù)不同躍遷速度確定躍遷時間為：

(6)

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

Δt=min{Δtj}

(7)

此時，狀態(tài)變量的仿真推進時間為：

ti+1=ti+Δt

(8)

躍遷狀態(tài)變量為：

(9)

非躍遷狀態(tài)變量為：

(10)

式中，n表示多步法的階數(shù)。

(11)

表1 各階數(shù)對應(yīng)的系數(shù)αk

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

3 仿真算例實驗

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

系統(tǒng)各狀態(tài)變量相對誤差計算公式如下：

(12)

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

一個線性二階ODE系統(tǒng)[19]如下：

(13)

圖2 LMSQSS算法對算例的仿真結(jié)果

表2 不同算法對算例1仿真結(jié)果

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

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

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

1—主閥左腔；2—主閥右腔彈簧；3—主閥右腔；4—過濾器；5—先導(dǎo)閥彈簧。圖3 溢流閥結(jié)構(gòu)簡圖

4 工程實例

4.1 工程實例模型

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

(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為泵至主閥左腔管道容腔中用來補償油液壓縮量及管道受壓變形量的液壓油總?cè)莘e，單位為cm3；V2為主閥右腔中用來補償油液壓縮量的液壓油總?cè)莘e，單位為cm3；P1和P2分別為主閥芯和先導(dǎo)閥閥芯的動量，單位為N·s/cm2；Y1和Y2分別為主閥芯和先導(dǎo)閥閥芯的位移，單位為cm。其余參數(shù)在表3中進行說明。

表3 溢流閥參量值

4.2 工程實例仿真結(jié)果分析

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

圖4 閥芯位移的仿真結(jié)果

圖5 入口壓力的仿真結(jié)果

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

表4 不同算法對溢流閥仿真結(jié)果

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

5 結(jié)束語

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