基于量化狀態(tài)系統(tǒng)的多體系統(tǒng)動(dòng)力學(xué)求解方法

李志華，賀英良，吳晨佳

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

0 引言

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

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

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

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

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

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

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

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

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

(1)

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

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

(2)

(3)

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

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

(4)

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

(5)

tj+1=tj+Δtj；

(6)

(7)

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

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

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

(8)

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

(9)

qj(t)=

(10)

(11)

(12)

(13)

(14)

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

(15)

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

(16)

while(tj

end if

K1=f(tj,xj)

Δtj=Δxj/fj//本步推進(jìn)時(shí)長(zhǎng)

end while

3 算例

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

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

(2)各算法均在MATLAB R2014b平臺(tái)上實(shí)現(xiàn)。

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

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

(17)

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

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

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

(18)

(19)

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

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

(20)

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

(21)

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

(22)

H=T+V。

(23)

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

(24)

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

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

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

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

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

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

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

4 結(jié)束語(yǔ)

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

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

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

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

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

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