李志華，賀英良，吳晨佳

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

0 引言

在機(jī)械、航空、航天、兵器、機(jī)器人等領(lǐng)域中，多體系統(tǒng)動力學(xué)為大量機(jī)械系統(tǒng)進(jìn)行動態(tài)性能評估和優(yōu)化提供了強(qiáng)有力的理論工具與技術(shù)支撐，是當(dāng)今力學(xué)領(lǐng)域的研究熱點和難點之一[1]。多體系統(tǒng)在運動過程中，由于不同構(gòu)件之間的特性參數(shù)存在較大差異，或者柔性體大范圍運動與構(gòu)件本身較小彈性變形之間存在耦合，使得動力學(xué)方程呈現(xiàn)剛性，對這類剛性方程進(jìn)行求解是多體系統(tǒng)動力學(xué)控制中的難點問題之一。

目前，多體系統(tǒng)動力學(xué)方程常見的數(shù)值解法有Newmark法、Wilson-θ法、Houbot法、Runge-kuta法和Gear法等[1-4]，這些方法均基于時間離散的數(shù)值積分方法。當(dāng)動力學(xué)方程為剛性，采用這些傳統(tǒng)的數(shù)值積分方法求解時，為了保證求解過程穩(wěn)定，需要強(qiáng)制使用隱式算法，因為所有顯式方法都必須顯著減少積分步長來確保數(shù)值的穩(wěn)定性，而隱式算法需要在每一個積分步中調(diào)用迭代算法來解一個2n維的非線性代數(shù)方程組，其過程繁瑣而復(fù)雜，計算成本隨系統(tǒng)規(guī)模的增大顯著增加[5]。

量化狀態(tài)系統(tǒng)(Quantized State System, QSS)方法是一種新的數(shù)值積分方法，其與傳統(tǒng)基于時間離散的積分方法明顯不同，QSS對狀態(tài)變量進(jìn)行離散，即系統(tǒng)的狀態(tài)變量以量子為單位躍遷，依次計算每次狀態(tài)變量躍遷所需要的時間，從而推進(jìn)積分。這一思想最早由Zeigler等[6]提出，并由Kofman等[7]首次實現(xiàn)。在求解一般動力學(xué)方程時，QSS不但具有穩(wěn)定性強(qiáng)、精確度高等優(yōu)點，而且計算過程完全不需要迭代，可以大大提升仿真效率[7]。

通過線性估計狀態(tài)變量躍遷所需時間的方法稱為一階QSS方法。Kofman[8]后續(xù)又給出高階的QSS2和QSS3算法，這些算法采用高次曲線估計躍遷時間來提高算法精度。然而作為一種顯式算法，QSS在求解剛性系統(tǒng)時會出現(xiàn)仿真數(shù)值振蕩現(xiàn)象(即不穩(wěn)定現(xiàn)象)[9]，因此不能完全應(yīng)用于剛性系統(tǒng)。為解決這一問題，Migoni等[10]提出一階后向QSS(Backward QSS, BQSS)算法，然而該算法精度不高(只有一階精度)，在非線性系統(tǒng)中，誤差范圍大可能會出現(xiàn)假性平衡點，使仿真無法正常進(jìn)行。因此，近年來國外學(xué)者對QSS方法的改進(jìn)從未間斷，并將其應(yīng)用于多個領(lǐng)域，例如文獻(xiàn)[11]將其應(yīng)用在建筑通風(fēng)與空調(diào)系統(tǒng)仿真計算中，得到了較好的效果。

目前，國內(nèi)只有北京航空航天大學(xué)的朱雨童等[12-13]和清華大學(xué)電機(jī)系的檀添、李帛洋、楊祎、秦建[9,14-16]等將QSS方法應(yīng)用于具有間斷和剛性的電力電子系統(tǒng)仿真求解中，并取得了一定的效果。

本文針對具有非線性和剛性特性的多體系統(tǒng)動力學(xué)，提出一種基于量化狀態(tài)系統(tǒng)的多點校正(Multi-point Correction QSS, MCQSS)顯式求解算法，該方法在保留QSS顯式計算特性的同時，引入多點校正思想對狀態(tài)變量導(dǎo)數(shù)進(jìn)行修正，有效提高了算法的精度和穩(wěn)定性。通過對多體系統(tǒng)的經(jīng)典算例雙擺進(jìn)行仿真求解，驗證了MCQSS算法的有效性，同時與傳統(tǒng)數(shù)值積分方法、QSS方法和BQSS算法進(jìn)行了性能對比。

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

連續(xù)系統(tǒng)數(shù)值仿真需要離散化，例如Euler，Runge-kutta等傳統(tǒng)方法是對時間進(jìn)行離散，QSS方法是對狀態(tài)變量進(jìn)行離散，即將狀態(tài)變量離散成一個個量子，計算狀態(tài)變量由一個量化狀態(tài)躍遷到另一個量化狀態(tài)所用的時間，以此推進(jìn)積分。需要說明的是，量子是事先給定的，類似于傳統(tǒng)數(shù)值積分算法中的時間步長，量子的大小將直接影響求解的精度和速度。

連續(xù)系統(tǒng)的數(shù)學(xué)模型可用常微分方程組表示，將常微分方程組化為如下形式的狀態(tài)方程：

(1)

式中：x為狀態(tài)向量；u為輸入量。

對該系統(tǒng)進(jìn)行量化處理，可得

(2)

(3)

式中：Δqj為狀態(tài)變量的量子；εj為遲滯寬度，為了在不增加誤差的條件下減少模型震蕩，遲滯寬度εj與量子Δqj相同。

狀態(tài)變量x各分量躍遷時需要變化的值為：

(4)

式中：qj為量化向量q的第j個分量；xj為狀態(tài)向量x的第j個分量。

(5)

tj+1=tj+Δtj；

(6)

(7)

2 基于量化狀態(tài)系統(tǒng)的多點校正顯式算法

針對QSS方法無法求解剛性系統(tǒng)和BQSS算法精度不高、在非線性系統(tǒng)中仿真穩(wěn)定性差的問題，提出一種MCQSS顯式算法。該算法在QSS和BQSS的基礎(chǔ)上，引入多點校正思想對狀態(tài)變量導(dǎo)數(shù)進(jìn)行修正，大幅提高了仿真中每步時間節(jié)點的精度，從而控制了誤差范圍，避免求解非線性剛性系統(tǒng)時出現(xiàn)擾動項，有效提高了算法精度和穩(wěn)定性，同時保留了QSS方法顯式計算無需迭代的特點，提高了算法的仿真效率。

假設(shè)動力學(xué)方程最終可化為常微分方程形式

(8)

相應(yīng)經(jīng)量化處理后可表示為

(9)

qj(t)=

(10)

(11)

(12)

(13)

(14)

由3個修正值確定該段的擬合函數(shù)斜率最終值fj，

(15)

利用修正后的最終斜率確定狀態(tài)變量此步的最終時間步長：

(16)

while(tj

end if

K1=f(tj,xj)

Δtj=Δxj/fj//本步推進(jìn)時長

end while

3 算例

本章以典型多體系統(tǒng)雙擺為例，將MCQSS算法與傳統(tǒng)數(shù)值積分方法(ode45，ode23s，ode15s)、QSS方法和BQSS算法從仿真精度與仿真效率兩方面進(jìn)行性能對比。仿真條件如下：

(1)計算機(jī)硬件配置為Intel i5-3427U處理器，CPU運行速度為1.80 GHz，操作系統(tǒng)為Windows 8.1。

(2)各算法均在MATLAB R2014b平臺上實現(xiàn)。

(3)對所有的算法設(shè)置相同的誤差容限，即Rel.Tol=Abs.Tol=10-3，同時將MCQSS算法中的量子大小分別設(shè)為Δq=10-2和Δq=10-3。

(4)總能量H的相對誤差計算公式為

(17)

式中：u[k]為各算法求得總能量H的值；udassl[k]為微分代數(shù)系統(tǒng)求解器(Differential Algebraic System Solver, DASSL)算法在很小誤差設(shè)定(10-9)下求的總能量H的值，在此作為基準(zhǔn)值[17]。DASSL算法是一種使用后向差分公式的多步隱式求解方法，是求解剛性方程常用的時間積分方法。因為該方法需要使用牛頓法對差分方程進(jìn)行迭代求解，雖然需要占用較大的資源且時間較長，但是精度較高[18]，所以大多數(shù)仿真工具使用DASSL或其變種作為剛性方程的默認(rèn)求解器[11]。

如圖1所示的雙擺系統(tǒng)，若m1和m2的質(zhì)量相差懸殊，則可以使雙擺系統(tǒng)動力學(xué)方程呈現(xiàn)剛性[19]。這里取m1=1 kg，m2=500 kg，兩根擺桿長度l1=l2=1 m，不計重量。

采用拉格朗日方法對其建模，取系統(tǒng)的廣義坐標(biāo)q=[θ1θ2]T,則廣義質(zhì)量矩陣與廣義力矩陣分別為：

(18)

(19)

式中θ1和θ2為兩根擺桿與豎直方向所成的角度。

雙擺系統(tǒng)的動力學(xué)方程為

(20)

系統(tǒng)的動能T、勢能V和總能量H分別為[20]：

(21)

V=-m1gl1cosθ1-m2g(l1cosθ1+l2cosθ2)；

(22)

H=T+V。

(23)

將該雙擺系統(tǒng)的動力學(xué)方程式(20)化為如式(8)所示的常微分方程形式：

(24)

對于本文的理想雙擺系統(tǒng)(即不計空氣阻力、鉸鏈約束無摩擦)，根據(jù)能量守恒定律，由給定的初始條件可知系統(tǒng)總能量始終保持為零。采用不同算法仿真雙擺系統(tǒng)的總能量H，得到如圖3所示的軌跡，其中QSS無解。由圖3可見，DASSL的求解結(jié)果最好，MCQSS算法的結(jié)果較理想，其他算法精度較差。隨著仿真時間的增加，它們的累積誤差越來越大，特別是ode45，不但精度差，而且還會出現(xiàn)不穩(wěn)定現(xiàn)象；MCQSS算法能將誤差控制在一個較小的范圍內(nèi),說明其精度和穩(wěn)定性較好。

圖3中，DASSL的求解結(jié)果最好，MCQSS算法的結(jié)果較理想，它們之間的異同之處如下：

(1)DASSL算法是基于時間離散的數(shù)值積分方法，MCQSS算法是基于狀態(tài)變量離散的積分方法。

(2)DASSL算法為隱式方法，在很小誤差設(shè)定(10-9)下，得到圖3中最好的結(jié)果(即精度最高)，但計算時間長；本文方法為顯式方法，在誤差設(shè)定為10-3時，得到圖3中較理想的結(jié)果(即精度較高)，且計算時間短。

各算法對雙擺系統(tǒng)總能量的求解結(jié)果如表1所示，可見QSS無解，說明其不適合求解剛性系統(tǒng)；相比ode45，MCQSS算法精度提高了20倍，效率提高了近3倍；相比ode15s和ode23s，MCQSS算法精度提高了近7倍，效率提高了1倍多；相比BQSS算法，在CPU運行時間相差不大的情況下，MCQSS算法精度提高了3倍多。

表1 各算法結(jié)果比較

為了研究量子大小對MCQSS算法性能的影響，分別設(shè)置了兩個不同的量子值來仿真求解。由表1可見，當(dāng)量子值減小時，CPU運行時間相差并不大，但是仿真精度卻提升較多。

4 結(jié)束語

針對多體系統(tǒng)動力學(xué)剛性方程求解問題，本文嘗試采用基于狀態(tài)變量離散的積分方法，提出一種MCQSS顯式算法，通過求解雙擺系統(tǒng)動力學(xué)剛性方程得到如下結(jié)論：

(1)ode45是顯式算法，穩(wěn)定性較差，很難兼顧仿真的精度和效率，且很難滿足剛性系統(tǒng)在實際求解中的需要。

(2)ode15s和ode23s是傳統(tǒng)數(shù)值積分方法中求解剛性系統(tǒng)較好的算法，本文MCQSS算法在仿真精度與仿真效率兩方面比ode15s和ode23s更具優(yōu)勢。

(3)QSS方法不適合求解剛性系統(tǒng)，BQSS算法盡管能求解剛性系統(tǒng)，但其精度不如MCQSS算法，可見MCQSS算法中的導(dǎo)數(shù)修正值對提升精度是有效的。

(4)適當(dāng)減小量子值能在幾乎不增加仿真運行時間的同時有效提高仿真精度。

下一步將開展基于量化狀態(tài)的混合系統(tǒng)求解方法研究，以擴(kuò)展QSS方法的應(yīng)用范圍。