螢火蟲算法求解多體系統(tǒng)動力學(xué)微分-代數(shù)方程*

張笑笑 丁潔玉 ，2?

（1.青島大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，青島266071）（2.青島大學(xué)計算力學(xué)與工程仿真研究中心，青島266071）

引言

微分-代數(shù)方程是描述多體系統(tǒng)動力學(xué)的數(shù)學(xué)模型，它由微分方程和代數(shù)約束方程組成，其數(shù)值求解方法的研究是多體系統(tǒng)動力學(xué)研究的重要內(nèi)容.保持約束穩(wěn)定從而保證微分-代數(shù)方程求解的穩(wěn)定是數(shù)值方法研究的重點(diǎn)［1］.近年來逐漸發(fā)展起來的一些方法，如：坐標(biāo)縮并方法［2］、廣義-α方法［3］、辛算法［4］、能量方法［5，6］、辛保持的變分?jǐn)?shù)值積分方法［7］、Lie 群方法［8］等，均取得了較好的求解結(jié)果.

螢火蟲算法（firefly algorithm，F(xiàn)A）是XS.Yang在2008年提出的一種新穎的智能優(yōu)化算法［9］.該方法的思想是模擬螢火蟲閃爍行為，通過螢火蟲間的相互吸引來實(shí)現(xiàn)種群的迭代進(jìn)化.近年來，研究人員對標(biāo)準(zhǔn)FA算法進(jìn)行諸多改進(jìn)，較好地解決算法收斂速度慢、求解精度低、易陷入局部最優(yōu)等缺陷.例如：將變量從混沌空間變換到解空間然后再進(jìn)行搜索的混沌螢火蟲算法［10，11］；通過將FA與其他算法結(jié)合提出混合螢火蟲算法［12-14］；基于離散空間的離散螢火蟲算法［15］等.目前該類算法已經(jīng)在諸多優(yōu)化問題中得到了較好的效果，并且在機(jī)器學(xué)習(xí)、生產(chǎn)調(diào)度、機(jī)器人智能控制和圖像處理等領(lǐng)域得到了廣泛的應(yīng)用.

本文主要通過Lagrange插值，結(jié)合Gauss數(shù)值積分方法將微分-代數(shù)方程求解問題轉(zhuǎn)換成優(yōu)化問題，然后對該優(yōu)化問題進(jìn)行螢火蟲算法設(shè)計.最后對平面雙連桿系統(tǒng)進(jìn)行仿真實(shí)驗(yàn).

1 多體系統(tǒng)動力學(xué)微分-代數(shù)方程

多體系統(tǒng)動力學(xué)方程通常為指標(biāo)3的微分-代數(shù)方程組（DAEs）：

其中，q?Rn是廣義坐標(biāo)，q??Rn是廣義速度，λ?Rd是Lagrange乘子，Φ為廣義坐標(biāo)q的約束方程，Φq為約束方程的Jacobi矩陣.

將約束方程Φ(q，t)=0求兩階導(dǎo)，方程（1）可以由指標(biāo)3降為指標(biāo)1，方程形式如下：

將仿真時間[0，T]平均劃分為若干小區(qū)間[ti，ti+1]，i=0，1，…，N.在時間 [ti，ti+1]中對廣義坐標(biāo)q(t)進(jìn)行Lagrange插值：

則

為保證位移約束、速度約束和加速度約束，插值函數(shù)q(t)、q?(t)、q?(t)需要滿足如下 3d個約束方程：

通過方程組（7）將方程（2）中的變量進(jìn)行縮并，3d個 變 量 (q1，...qd；，...；，...)可 以 由3(n-d)個 變 量 (qd+1，...qn；q?d+1，...q?n；q?d+1，...q?n)來表示 .然后將q(t)、q?(t)、q?(t)代入微分-代數(shù)方程（2）得到：

若對于?t∈ [ti，ti+1]都有f(t)=0，則插值函數(shù)q(t)為式（1）在[ti，ti+1]上的精確解.顯然非線性函數(shù)（8）一般不存在精確解，只能求解相對于插值函數(shù)q(t)盡可能滿足（8）式接近于0的相對最優(yōu)解.則可以轉(zhuǎn)換成如下優(yōu)化問題：

對上式進(jìn)行Gauss數(shù)值積分，可以化成如下形式：

其中，As為Gauss積分系數(shù)，ts為Gauss積分節(jié)點(diǎn).顯然，式（10）是一個非凸優(yōu)化問題［16］，已有的研究表明，傳統(tǒng)的數(shù)學(xué)優(yōu)化方法難以對此類問題進(jìn)行有效求解，智能優(yōu)化算法是求解此類問題的有效算法.螢火蟲算法是模仿種群運(yùn)動的一種智能算法，提供了求解DAE方程的新方法.

2 螢火蟲算法設(shè)計

螢火蟲算法包含一個有M個個體的螢火蟲種群系統(tǒng) X={X1，X2，X3，…，XM}，Xi∈RD.第i個個體在時間t處的位置：Xi(t)=(xi1(t)，xi2(t)，…，xiD(t)).螢火蟲對彼此吸引的原因取決于兩個因素，即自身亮度Li和吸引度β.其中，螢火蟲的亮度取決于自身所在位置的目標(biāo)值Li=φ(Xi).吸引度與亮度相關(guān)，愈亮的螢火蟲吸引亮度比其弱的螢火蟲往這個方向移動.吸引度與距離Ri，j成反比，距離越遠(yuǎn)的螢火蟲吸引度越低.系統(tǒng)X的整體最優(yōu)位置為Xg.迭代過程中，螢火蟲的位置更新公式如下［17］：

其中，γ是光吸引系數(shù)；β0是最大吸引度；α∈ (0，1)是步長因子；α為擾動項；υ為動態(tài)視覺權(quán)重，越大尋優(yōu)視野越大，越容易獲得遠(yuǎn)距離信息；Ri，j為螢火蟲i和j之間的距離.

迭代初期，希望對螢火蟲種群個體進(jìn)行較大擾動，以增強(qiáng)全局探索能力有利于跳出局部極值點(diǎn)；而在搜索后期，要求對個體減小擾動，以提高算法的搜索精度.因此，設(shè)計擾動項如下：

運(yùn)用螢火蟲算法求解區(qū)間[ti，ti+1]的插值函數(shù)q(t)，已知初始點(diǎn)q0=q(ti) ∈Rn，插值節(jié)點(diǎn)qj，j=1，2…m是待優(yōu)化的節(jié)點(diǎn)，并且qj滿足約束Φ(q，t)=0。.一般有e個約束方程，則維度d=(n-e)m.定義目標(biāo)函數(shù)如下：

如果φ(qj)=0，則qj，j=1，2…m即為所求的最優(yōu)值.由于目標(biāo)函數(shù)φ為非線性函數(shù)，一般不存在解析解，所以φ(qj)只能盡可能地接近精確解，設(shè)置收斂精度φmin和最大迭代步數(shù)gmax，若達(dá)到收斂精度或達(dá)到最大迭代步數(shù)，則輸出尋優(yōu)結(jié)果.

螢火蟲算法的求解步驟如下：

1）設(shè)置螢火蟲算法的初始參數(shù)：種群規(guī)模M，螢火蟲位置維度d，螢火蟲位置的變化范圍，最大迭代次數(shù)gmax，收斂精度φmin，光吸引系數(shù)γ；最大吸引度β0；步長因子α，初始擾動項=0，并隨機(jī)生成初始種群，定義目標(biāo)函數(shù)φ；

2）根據(jù)設(shè)定的目標(biāo)函數(shù)，分別計算每個個體的適應(yīng)度值，比較種群的適應(yīng)度，確定種群的亮度最強(qiáng)個體Xg；

4）比較移動前后的螢火蟲適應(yīng)度值，若優(yōu)于之前位置則完成位置更新，否則保持原位置不移動；

5）如果φ(Xg)<φmin或達(dá)到最大迭代步數(shù)gmax，則輸出最優(yōu)解gbest，否則重復(fù)步驟2）-步驟5），直到滿足預(yù)設(shè)的終止條件.

3 數(shù)值算例

圖1 平面雙連桿機(jī)械臂Fig.1 Two-link planar manipulator

兩點(diǎn)Gauss數(shù)值積分方法求解目標(biāo)函數(shù).

圖2-圖5為時間步長0.005時系統(tǒng)仿真結(jié)果.從連桿末端運(yùn)動軌跡和系統(tǒng)能量變化可以看出，連桿運(yùn)動過程穩(wěn)定，沒有出現(xiàn)明顯偏移.位移約束、速度級約束和加速度級約束沒有出現(xiàn)違約.

圖2 連桿末端運(yùn)動軌跡Fig.2 The trajectory of the end of rod

圖3 系統(tǒng)總能量Fig.3 Total System Energy

圖4 系統(tǒng)動能、勢能Fig.4 System Kinetic Energy and Potential Energy

圖5 雙連桿機(jī)械臂約束Fig.5 Constraints of the two-link manipulator

采用相同時間步長h=0.01對平面雙連桿系統(tǒng)進(jìn)行仿真實(shí)驗(yàn)，結(jié)果比較見表2，其中FA-DAE表示本文的求解算法，RK4表示四階Runge-Kutta方法，仿真時間為20s.

從表1可以看出，系統(tǒng)運(yùn)動過程中精確保持位移約束、速度級約束和加速度級約束，系統(tǒng)總能量最大相對誤差較小.隨著時間步長減小，系統(tǒng)總能量最大相對誤差和平均相對誤差在減小，仿真結(jié)果也更加的準(zhǔn)確.從表2可以看出，本文算法運(yùn)行時間長，但在約束和能量方面保持較好.

表1 不同時間步長的誤差結(jié)果比較Table 1 Comparison of errors for different time steps

表2 相同時間步長時各方法結(jié)果比較(h=0.01，t=20s)Table 2 Comparison of different methods with the same time step(h=0.01，t=20s)

4 小結(jié)

本文運(yùn)用螢火蟲優(yōu)化算法求解多體系統(tǒng)動力學(xué)微分-代數(shù)方程，并對平面雙連桿系統(tǒng)進(jìn)行仿真實(shí)驗(yàn).實(shí)驗(yàn)結(jié)果表明，在滿足指標(biāo)1、指標(biāo)2和指標(biāo)3約束的情況下，系統(tǒng)總能量的誤差較小.螢火蟲算法的優(yōu)化精度高、算法設(shè)計簡單且不需要目標(biāo)函數(shù)的導(dǎo)數(shù)信息，對于求解多體系統(tǒng)微分-代數(shù)方程取得了較好效果，說明智能優(yōu)化算法應(yīng)用到多體動力學(xué)仿真中的可行性.然而，該方法存在計算量大和算法運(yùn)行時間長的缺點(diǎn).今后在有效地降低運(yùn)行時間和提高算法精度方面是一個主要研究方向.