1、螢火蟲算法求解多體系統(tǒng)動(dòng)力學(xué)微分-代數(shù)方程摘要針對(duì)多體系統(tǒng)動(dòng)力學(xué)微分-代數(shù)方程求解問(wèn)題，研究基于螢火蟲算法的求解方法.首先將廣義坐標(biāo)和 廣義速度進(jìn)行Lagrange插值，結(jié)合Gauss數(shù)值積分方法，將微分-代數(shù)方程求解問(wèn)題轉(zhuǎn)化成求解最優(yōu)化問(wèn)題. 然后用螢火蟲算法對(duì)問(wèn)題進(jìn)行優(yōu)化求解.最后，通過(guò)對(duì)平面雙連桿機(jī)械臂的多體系統(tǒng)仿真實(shí)驗(yàn)，驗(yàn)證了螢火 蟲算法在求解動(dòng)力學(xué)方程中既保持了約束又較好地保證了能量精度.結(jié)果表明智能優(yōu)化算法在求解多體動(dòng) 力學(xué)問(wèn)題上具有較好的應(yīng)用前景.關(guān)鍵詞多體系統(tǒng)動(dòng)力學(xué)，微分代數(shù)方程，螢火蟲算法引言積分方法將微分-代數(shù)方程求解問(wèn)題轉(zhuǎn)換成優(yōu)化 問(wèn)題，然后對(duì)該優(yōu)化問(wèn)題進(jìn)行螢火蟲算法設(shè)計(jì)

2、.最微分-代數(shù)方程是描述多體系統(tǒng)動(dòng)力學(xué)的數(shù) 學(xué)模型，它由微分方程和代數(shù)約束方程組成，其數(shù) 值求解方法的研究是多體系統(tǒng)動(dòng)力學(xué)研究的重要 內(nèi)容.保持約束穩(wěn)定從而保證微分-代數(shù)方程求解 的穩(wěn)定是數(shù)值方法研究的重點(diǎn)近年來(lái)逐漸發(fā)展 起來(lái)的一些方法，如：坐標(biāo)縮并方法、廣義-a方 法、辛算法、能量方法56、辛保持的變分?jǐn)?shù)值積 分方法、Lie群方法時(shí)等，均取得了較好的求解 結(jié)果.螢火蟲算法(firefly algorithm，F(xiàn)A)是 XS.Yang 在2008年提出的一種新穎的智能優(yōu)化算法該方 法的思想是模擬螢火蟲閃爍行為，通過(guò)螢火蟲間的 相互吸引來(lái)實(shí)現(xiàn)種群的迭代進(jìn)化.近年來(lái)，研究人 員對(duì)標(biāo)準(zhǔn)FA算法進(jìn)行諸

3、多改進(jìn)，較好地解決算法 收斂速度慢、求解精度低、易陷入局部最優(yōu)等缺陷. 例如:將變量從混沌空間變換到解空間然后再進(jìn)行 搜索的混沌螢火蟲算法10,11;通過(guò)將FA與其他算 法結(jié)合提出混合螢火蟲算法12-141 ;基于離散空間的 離散螢火蟲算法15等.目前該類算法已經(jīng)在諸多優(yōu) 化問(wèn)題中得到了較好的效果，并且在機(jī)器學(xué)習(xí)、生 產(chǎn)調(diào)度、機(jī)器人智能控制和圖像處理等領(lǐng)域得到了 廣泛的應(yīng)用.本文主要通過(guò)Lagrange插值,結(jié)合Gauss數(shù)值后對(duì)平面雙連桿系統(tǒng)進(jìn)行仿真實(shí)驗(yàn).1多體系統(tǒng)動(dòng)力學(xué)微分-代數(shù)方程多體系統(tǒng)動(dòng)力學(xué)方程通常為指標(biāo)3的微分-代 數(shù)方程組(DAEs)：(q，t )q + 中泓=F (q,q,t)

4、()麗 q，t) = 0其中，q*是廣義坐標(biāo)，q*是廣義速度，XwR“是 Lagrange乘子，中為廣義坐標(biāo)q的約束方程，伐為 約束方程的Jacob i矩陣.將約束方程中(q，t) = 0求兩階導(dǎo)，方程(1)可 以由指標(biāo)3降為指標(biāo)1,方程形式如下：_ m e ij f, o J_a_ _-(e的)2- 2e“ -j- e(2 ) 將仿真時(shí)間0, T 平均劃分為若干小區(qū)間 t t, + ,，z = 0, 1，,N.在時(shí)間也，t, + ,中對(duì)廣義 坐標(biāo)q (t)進(jìn)行Lagrange插值：q (t) =(t) q,(3)1 = 0m (t t )l()= n t-；k = 0,k 豐 j ( t L

5、k )貝u+ 1，. qn ； + 1，. qn ； qd + 1，. qn ； qd +1，q (t) _-mtF (1(t),q (t)，t)_ A _一 q0 _-(中 qq) ,4 - 2中旗-中+,(10)(12)q( t)=項(xiàng)(t) q,.! = 0q(t)=El(t) q,.! = 0為保證位移約束、速度約束和加速度約束,插值函 數(shù)q (t)、q (t)、q (t)需要滿足如下3d個(gè)約束方程：，中(g, t) = 0 血 q,q,t) = 0曲 q,q,q,t)= 0通過(guò)方程組(7)將方程(2)中的變量進(jìn)行縮并， 3d個(gè)變量(q， qd ； q 1，qd ； q 1，q d) 可

6、以由 3(n - d)個(gè)變量(qd + 1， q ； qd + 1，qn ； qd+1，.qn) 來(lái)表示.然后將q (t)、q (t)、q (t)代入微分-代數(shù)方程 (2)得到：Vt wEt若對(duì)于Vt e t，t +1都有f (t) = 0,則插值函數(shù)q(t) 為式(1)在tt, + 1上的精確解.顯然非線性函數(shù) (8)一般不存在精確解，只能求解相對(duì)于插值函數(shù) q (t)盡可能滿足(8)式接近于0的相對(duì)最優(yōu)解.則可 以轉(zhuǎn)換成如下優(yōu)化問(wèn)題：min f f fdxq t對(duì)上式進(jìn)行Gauss數(shù)值積分，可以化成如下形式： m，in尸(ts) f (ts) dx1 s = 1 其中,為Gauss積分系數(shù)

7、，ts為Gauss積分節(jié)點(diǎn).顯 然，式(10)是一個(gè)非凸優(yōu)化問(wèn)題商，已有的研究表 明，傳統(tǒng)的數(shù)學(xué)優(yōu)化方法難以對(duì)此類問(wèn)題進(jìn)行有效 求解,智能優(yōu)化算法是求解此類問(wèn)題的有效算法. 螢火蟲算法是模仿種群運(yùn)動(dòng)的一種智能算法，提供 了求解DAE方程的新方法.2螢火蟲算法設(shè)計(jì)螢火蟲算法包含一個(gè)有M個(gè)個(gè)體的螢火蟲種 群系統(tǒng) X = X,X2,X3,-,Xm,X, w RD.第 i 個(gè)個(gè)體 在時(shí)間 t 處的位置:X, (t) = ( x.1 (t),x.2 (t), ,*iD (t). 螢火蟲對(duì)彼此吸引的原因取決于兩個(gè)因素，即自身 亮度L和吸引度伙其中，螢火蟲的亮度取決于自身 所在位置的目標(biāo)值L =中(X,).

8、吸引度與亮度相關(guān)， 愈亮的螢火蟲吸引亮度比其弱的螢火蟲往這個(gè)方 向移動(dòng).吸引度與距離Rhl成反比，距離越遠(yuǎn)的螢火 蟲吸引度越低.系統(tǒng)X的整體最優(yōu)位置為X，.迭代 過(guò)程中，螢火蟲的位置更新公式如下17：Xt+1 = XI + 6(X，- X；)+ an(11)B = 60 e-R其中,Y是光吸引系數(shù);60是最大吸引度;a w (0, 1) 是步長(zhǎng)因子；ani為擾動(dòng)項(xiàng);u為動(dòng)態(tài)視覺(jué)權(quán)重，越 大尋優(yōu)視野越大，越容易獲得遠(yuǎn)距離信息;R“,為螢 火蟲i 和j之間的距離.迭代初期,希望對(duì)螢火蟲種群個(gè)體進(jìn)行較大擾 動(dòng)，以增強(qiáng)全局探索能力有利于跳出局部極值點(diǎn)； 而在搜索后期，要求對(duì)個(gè)體減小擾動(dòng)，以提高算法 的

9、搜索精度.因此，設(shè)計(jì)擾動(dòng)項(xiàng)如下：n+1 = ant + 6 (X，- x;)= 0運(yùn)用螢火蟲算法求解區(qū)間ti, ti+1的插值函數(shù) q (t),已知初始點(diǎn)q 0 = q (t,) w Rn，插值節(jié)點(diǎn)q,= 1, 2 m.是待優(yōu)化的節(jié)點(diǎn)，并且q,滿足約束0(q, t)= 0。. 一般有e個(gè)約束方程,則維度d = (n-e)m.定 義目標(biāo)函數(shù)如下：中(q,) = AfT( ts )f (ts)(13)如果中(q,) =0,則q, j = 1,2-m.即為所求的最優(yōu) 值.由于目標(biāo)函數(shù)中為非線性函數(shù)，一般不存在解 析解，所以中(q,)只能盡可能地接近精確解,設(shè)置收 斂精度J和最大迭代步數(shù)M，若達(dá)到收斂

10、精度 或達(dá)到最大迭代步數(shù)，則輸出尋優(yōu)結(jié)果.螢火蟲算法的求解步驟如下：設(shè)置螢火蟲算法的初始參數(shù):種群規(guī)模M, 螢火蟲位置維度d,螢火蟲位置的變化范圍，最大 迭代次數(shù)，收斂精度 J，光吸引系數(shù)Y ；最大吸 引度60；步長(zhǎng)因子a,初始擾動(dòng)項(xiàng)n1 = 0,并隨機(jī)生 成初始種群,定義目標(biāo)函數(shù)中；根據(jù)設(shè)定的目標(biāo)函數(shù)，分別計(jì)算每個(gè)個(gè)體的 適應(yīng)度值，比較種群的適應(yīng)度，確定種群的亮度最 強(qiáng)個(gè)體X，；計(jì)算螢火蟲Xi到亮度最強(qiáng)個(gè)體X，的距離 R、，,確定X，對(duì)所有個(gè)體的吸引度6,根據(jù)公式(10) 對(duì)所有螢火蟲進(jìn)行移動(dòng)操作，更新擾動(dòng)項(xiàng)n+,；比較移動(dòng)前后的螢火蟲適應(yīng)度值，若優(yōu)于之 前位置則完成位置更新,否則保持原位置

11、不移動(dòng)；如果中(X，)(Pmin或達(dá)到最大迭代步數(shù)gma,，則輸出最優(yōu)解gbest，否則重復(fù)步驟2)-步驟5),直 到滿足預(yù)設(shè)的終止條件.3數(shù)值算例下面以平面雙連桿機(jī)械臂系統(tǒng)為例，對(duì)上述方 法進(jìn)行驗(yàn)證和分析.如圖1所示，設(shè)桿1質(zhì)量m- 1 = g，桿長(zhǎng)L1 = 1m，質(zhì)心坐標(biāo)為3，).桿2的質(zhì)量 m2 = 2kg，桿長(zhǎng)L2 = a/3 m，質(zhì)心坐標(biāo)為(*2, y2).兩 個(gè)剛性連桿可以繞端點(diǎn)自由轉(zhuǎn)動(dòng)，取X軸為零勢(shì)能 面，假設(shè)系統(tǒng)只受重力作用，重力加速度g = 9.81. 系統(tǒng)狀態(tài)變量q = (xl，yl，0l，X2，y2, 02)，則系統(tǒng)動(dòng)力 學(xué)微分-代數(shù)方程為：L1 TOC o 1-5 h

12、z cos 0121L1 sin 012101L 2L1 cos 01 cos 02.L 2L1 sin 01 sin 0202(16)in 02(16)in 0201 01L1.cos 01 0121 1圖1平面雙連桿機(jī)械臂Fig. 1 Twolink planar manipulator01.L 2 .-L1 sins 01 01 sin 002l_L1 cos 01 01 cos02 0210_(LLcos 000sin 00LL- cos 000cos 000Lcos 000 L sins 00-Lsin000+ L cos00(17)L fin2-sin 0002 +COS 0022

13、 2cos 000in 00_ M _ M F-中q0 _A_(14)(19)(20)其中，廣義力矢量 F = (0, - m 1 g，0, 0, - m-2g，0)， 廣義質(zhì)量矩陣M = diag(m 1, m-1, J1,m2, m2, J),連桿 轉(zhuǎn)動(dòng)慣量J1 =二m 1L2, J2 =二m2L2.采用三點(diǎn)Lagrange 插值逼近 q (t),由初值 q(t0) = q0，q(t0) = q0礙： TOC o 1-5 h z / hq0+ 3q+ 92、了mq】=(4)，h= tl + 1 一 t (15)為滿足約束方程中二0,吊二0, = 0,插值函數(shù)q (t)中變量可以縮并表示為:

14、-/(18) 兩點(diǎn)Gauss數(shù)值積分方法求解目標(biāo)函數(shù). 咿(#尸(t)f (t叫、hA1 = a,=122-a/3 hti + 1 + ti11 = +1322_ V3 . hti + 1 + ti23 22取時(shí)間步長(zhǎng)h = 0.002s，仿真時(shí)間20s.連桿初始狀態(tài)：01=n ,0 _=- n，廣義速度q =0.螢火蟲算36法參數(shù)設(shè)置:種群規(guī)模M = 60,光吸引系數(shù)y = 1, 最大吸引度仇=0.7,步長(zhǎng)因子a = 0.7，收斂精度 J = 10-15，最大迭代步數(shù)g“ =1000.仿真結(jié)果如 圖所示：E.N)sJoua -Els12.74612.745圖2連桿末端運(yùn)動(dòng)軌跡Fig. 2 T

15、he trajectory of the end of rod圖4系統(tǒng)動(dòng)能、勢(shì)能Fig. 4 System Kinetic Energy and Potential Energy12.74912.7481274712.744圖3系統(tǒng)總能量Fig. 3 Total System Energy圖5雙連桿機(jī)械臂約束Fig. 5 Constraints of the two-link manipulatorE.N)sJoua -Els12.74612.745圖2連桿末端運(yùn)動(dòng)軌跡Fig. 2 The trajectory of the end of rod圖4系統(tǒng)動(dòng)能、勢(shì)能Fig. 4 System K

16、inetic Energy and Potential Energy12.74912.7481274712.744圖3系統(tǒng)總能量Fig. 3 Total System Energy圖5雙連桿機(jī)械臂約束Fig. 5 Constraints of the two-link manipulator采用不同時(shí)間步長(zhǎng)對(duì)平面雙連桿系統(tǒng)進(jìn)行仿 真實(shí)驗(yàn)，結(jié)果比較見(jiàn)表1.其中H = T + U為系統(tǒng)總 能量lH (t，)-H (10)|為系統(tǒng)總能量的 能量，8 (H ) = max為系統(tǒng)總能量的0 eH (t。)1 x | H(t,) - H(10) | 最大相對(duì)誤差，MRE(H)=1史 為QH (t0) 系統(tǒng)

17、總能量的平均相對(duì)誤差，8(中)、8傍)、8傍)分 別為系統(tǒng)位移約束、速度約束和加速度約束的最大 誤差，仿真時(shí)間為20s.采用相同時(shí)間步長(zhǎng)h = 0.01對(duì)平面雙連桿系 統(tǒng)進(jìn)行仿真實(shí)驗(yàn)，結(jié)果比較見(jiàn)表2,其中FA-DAE 表示本文的求解算法，RK4表示四階Runge-Kutta 方法，仿真時(shí)間為20s.表1不同時(shí)間步長(zhǎng)的誤差結(jié)果比較Table 1 Comparison of errors for different time stepsStep8(H)MRE(H)8(08(擊)8(市)Runtime T/sh = 0.0022.8782 x 10-43.2808 x 10-51.1102 x 10

18、-168.8817 x 10-162.8421 x 10-142.9 x 104h = 0.0051.8150 x 10-33.7283 x 10-41.1102 x 10-168.8817 x 10-162.8421 x 10-141.2 x 104h = 0.017.4496 x 10-31.3995 x 10-31.1102 x 10-168.8817 x 10-162.8421 x 10-148.8 x 103從表1可以看出，系統(tǒng)運(yùn)動(dòng)過(guò)程中精確保持位 最大相對(duì)誤差較小.隨著時(shí)間步長(zhǎng)減小，系統(tǒng)總能 移約束、速度級(jí)約束和加速度級(jí)約束，系統(tǒng)總能量 量最大相對(duì)誤差和平均相對(duì)誤差在減小，仿真結(jié)果表2相同時(shí)間步長(zhǎng)時(shí)各方法結(jié)果比較(h=0.01,t=20s)Table 2 Comparison of different methods with the same time step (h = 0.01, t=20s)Step研H)MRE(H)以0研擊)E(b)Runtime T/sFA-DAE7.4496 x 10-31.3995 x10-31.1102 x10-168.8817 x10-16RK49.0827 x 10-22.8464 x10-24.0164 x10-24.2045 x10-32.8421 x 10-14