張 哲 代洪華 馮浩陽 汪雪川 岳曉奎

(西北工業(yè)大學(xué)航天學(xué)院,西安 710072)

(航天飛行動(dòng)力學(xué)技術(shù)國家級(jí)重點(diǎn)實(shí)驗(yàn)室,西安 710072)

引言

軌道快速計(jì)算是航天工程領(lǐng)域的基礎(chǔ)問題,廣泛存在于地球繞行軌道設(shè)計(jì)及深空探測等各種航天工程任務(wù)中,具有重要意義.軌道動(dòng)力學(xué)問題在形式上一般描述為常微分方程初值問題或邊值問題,其中最具代表性的有軌道遞推問題[1-2]、攝動(dòng)Lambert問題[3-4]和三體模型下的轉(zhuǎn)移軌道設(shè)計(jì)問題[5-6].軌道問題在早期一般通過微積分進(jìn)行解析求解,牛頓與伯努利計(jì)算了萬有引力作用下的二體軌道解析解[7],歐拉與拉格朗日求出了限制性三體問題的5 個(gè)拉格朗日點(diǎn)[8].但當(dāng)考慮更為一般的三體系統(tǒng)以至N 體系統(tǒng)時(shí),解析求解將不再適用[9].龐加萊發(fā)展了漸進(jìn)法并將其用于求解限制性三體問題,但由于三體問題具有混沌效應(yīng),其長期軌道仍然難以通過漸進(jìn)法精確預(yù)測.此外,早期對軌道動(dòng)力學(xué)問題的求解都是基于理想的均勻圓球或橢球模型,而在實(shí)際航天任務(wù)中則需要考慮地球非球形等擾動(dòng)項(xiàng),這些干擾力使得航天器軌道難以求得精確解析解.由于存在上述困難,在實(shí)際工程任務(wù)中,航天器軌道均采用數(shù)值方法進(jìn)行求解.

軌道動(dòng)力學(xué)問題最為經(jīng)典的數(shù)值解法是有限差分類方法,其中代表性的包括Euler 法、Runge-Kutta法等[10].此類算法將待求解問題劃分為多個(gè)差分節(jié)點(diǎn),逐點(diǎn)計(jì)算問題近似解.由于物理含義直觀且遞推方程較為簡單,有限差分法受到了廣泛關(guān)注,并在仿真軟件與科學(xué)計(jì)算函數(shù)中得到大量應(yīng)用[11-12].然而,有限差分法的計(jì)算精度嚴(yán)重依賴于小步長,方法性能受到限制,在諸如航天器追逃博弈以及失效航天器快速在軌維修等任務(wù)中,將難以滿足任務(wù)實(shí)時(shí)性以及長期解算過程高效高精度的要求[13-14].

為克服有限差分法的原理性缺陷,學(xué)者們發(fā)展了一系列能實(shí)現(xiàn)大步長計(jì)算的迭代算法[15-21].1963年Clenshaw 與Norton[15]首次將Chebyshev 多項(xiàng)式與Picard 迭代法相結(jié)合,用于求解常微分方程的半解析解,該算法能夠避免有限差分法逐步計(jì)算的缺點(diǎn),但其系數(shù)計(jì)算方案較為復(fù)雜,且僅適用于求解一維方程.其后,日本Fukushima[16]對Chebyshev 多項(xiàng)式與Picard 迭代法的結(jié)合方式予以改進(jìn),簡化了系數(shù)計(jì)算方案,并將求解對象拓展到高維常微分方程,但在這一改進(jìn)方案中待解函數(shù)及其導(dǎo)數(shù)被分別估計(jì),未能有效利用Chebyshev 多項(xiàng)式性質(zhì),引入的待定系數(shù)較多,計(jì)算量較大.此后,美國德州農(nóng)工大學(xué)Bai[17]提出修正Chebyshev-Picard 迭代法(modified Chebyshev-Picard iteration method,MCPI),并將其用于求解軌道動(dòng)力學(xué)問題,該算法利用Chebyshev多項(xiàng)式性質(zhì)減少了待定系數(shù),但由于在迭代公式推導(dǎo)中利用了反演公式,得到的迭代公式較為復(fù)雜,降低了算法效率.為進(jìn)一步提高M(jìn)CPI 算法計(jì)算效率,Woollands 等[18]對該算法進(jìn)行改進(jìn),將導(dǎo)函數(shù)估計(jì)值與真實(shí)值之差作為迭代修正項(xiàng),提高了算法收斂階,但引入了雅可比矩陣,增加了算法計(jì)算量.Wang 等[19]將變分迭代法與時(shí)域配點(diǎn)法相結(jié)合,提出了局部變分迭代法(local variational iteration method,LVIM),該算法將計(jì)算殘差用于誤差修正,結(jié)合時(shí)域配點(diǎn)法簡化了迭代公式推導(dǎo)過程,得到的計(jì)算公式較為簡潔,在軌道計(jì)算問題中收斂速度較快,但該算法由于對廣義拉格朗日乘子進(jìn)行Taylor展開,同樣引入了雅可比矩陣,在求解復(fù)雜高維問題時(shí)產(chǎn)生的計(jì)算量過大.2020 年,Wang 等[20]將牛頓法與Picard 迭代法相結(jié)合,提出了多步Newton-Picard法(multistep Newton-Picard method,MNP),該算法將導(dǎo)數(shù)方程在待定函數(shù)真解處進(jìn)行Taylor 展開,并依據(jù)展開式建立了一系列衍生算法,但Taylor 展開引入的偏微分算子使得算法迭代公式較為復(fù)雜,且算法收斂階較低.馮浩陽等[21]將局部變分迭代法、擬線性化法及疊加法相結(jié)合,提出了能用于求解邊值問題的擬線性化-局部變分迭代法,拓展了局部變分迭代法的適用范圍,但仍未能解決算法需要反復(fù)計(jì)算雅可比矩陣的弊端.

本文將誤差反饋機(jī)制與Picard 迭代法相結(jié)合,提出一種新的反饋迭代算法,避免對原函數(shù)進(jìn)行泰勒展開,消除迭代公式中的雅可比矩陣,以保證計(jì)算效率;利用配點(diǎn)法與Chebyshev 多項(xiàng)式性質(zhì),將迭代公式推導(dǎo)中涉及的復(fù)雜符號(hào)運(yùn)算轉(zhuǎn)化為線性代數(shù)運(yùn)算,簡化推導(dǎo)過程,降低迭代公式的復(fù)雜度;結(jié)合擬線性化法與疊加法,拓展算法適用范圍,使之能夠求解兩點(diǎn)邊值問題;將變參數(shù)計(jì)算方案引入局部配點(diǎn)反饋迭代法,使算法能夠在運(yùn)行過程中自適應(yīng)選擇計(jì)算參數(shù),進(jìn)一步提高算法計(jì)算效率與計(jì)算精度,從而為在軌航天器的長期軌道快速預(yù)測提供一種高性能數(shù)值計(jì)算方法.

1 局部配點(diǎn)反饋迭代法

1.1 修正Chebyshev-Picard 迭代法

修正Chebyshev-Picard 迭代法是用于迭代求解常微分方程初值問題的數(shù)值算法[17],該算法求解常微分方程的具體過程簡述如下.

對于具有如下形式的一階常微分方程

假設(shè)未知函數(shù)及其導(dǎo)數(shù)可以用正交多項(xiàng)式的線性組合來近似,即

式中,Ti(t)為一組正交多項(xiàng)式組的第i 項(xiàng),βi及Fi為Ti(t)的系數(shù).

Picard 迭代法計(jì)算公式為

將式(2)與式(3)代入式(4),并令式(4)兩邊正交多項(xiàng)式對應(yīng)項(xiàng)系數(shù)在一系列給定時(shí)刻相等,進(jìn)一步整理可得未知函數(shù)y(t)的迭代計(jì)算公式為

1.2 局部配點(diǎn)反饋迭代法

在式(4)中,迭代公式通過對導(dǎo)數(shù)向量積分得到函數(shù)值,利用新的函數(shù)值更新導(dǎo)數(shù)向量并再次進(jìn)行積分迭代,直到達(dá)到目標(biāo)精度.本節(jié)將誤差反饋機(jī)制引入迭代過程,結(jié)合時(shí)域配點(diǎn)法與局部離散化思想,建立一種新的局部配點(diǎn)反饋迭代法.

Picard 迭代公式(4)與下式等價(jià)

對式(6)進(jìn)行整理,以導(dǎo)數(shù)估計(jì)值 y˙n(t) 與利用狀態(tài)估計(jì)值計(jì)算得到的導(dǎo)數(shù)值 g(yn(t),t) 之差作為修正項(xiàng),將式(6)改寫為如下迭代公式

稱式(7)對應(yīng)的迭代方案為反饋迭代法.

對于[t0,tf]范圍內(nèi)的任意時(shí)刻tm,可以使用下式將其轉(zhuǎn)化到區(qū)間[-1,1]范圍內(nèi)

通過式(8) 將求解區(qū)間進(jìn)行轉(zhuǎn)化后,使用Chebyshev 正交多項(xiàng)式的線性組合作為未知函數(shù)y(t)的估計(jì)解,即近似認(rèn)為

式中,cos(karccost) 為第k 項(xiàng)Chebyshev 多項(xiàng)式在t 時(shí)刻的值,lk為第k 項(xiàng)Chebyshev 多項(xiàng)式的系數(shù).

為求解式(9)中N 個(gè)未知多項(xiàng)式的系數(shù)lk(k=1,2,···,N),對式(9)采用時(shí)域配點(diǎn)法[22],即令正交多項(xiàng)式線性組合所構(gòu)成的函數(shù)估計(jì)解在N 個(gè)給定時(shí)刻t1,t2,···,tN與函數(shù)真值相等,從而構(gòu)造由N 個(gè)線性方程構(gòu)成的線性代數(shù)方程組

將式(10)記為如下矩陣形式

式中,y 為N 個(gè)時(shí)刻 t1,t2,···,tN對應(yīng)的函數(shù)真值構(gòu)成的真值向量,C 為Chebyshev 多項(xiàng)式矩陣,其中第i 行j 列元素Cij=cos(jarccosti),L 為系數(shù)向量L=[l1,l2,···,lN]T.

在給定時(shí)刻,未知函數(shù)的真值與Chebyshev 多項(xiàng)式矩陣C 可以預(yù)先求得,因此可以利用下式求解系數(shù)向量L

將式(12)得到的系數(shù)向量代入式(11),再將式(11)代入式(7),得到配點(diǎn)反饋迭代法的迭代公式

由于Chebyshev 多項(xiàng)式形式已知,對其進(jìn)行微積分運(yùn)算可得

由式(14)及式(15)可以得到Chebyshev 多項(xiàng)式矩陣的微分形式及積分形式.據(jù)此可將式(5)推導(dǎo)中涉及的符號(hào)運(yùn)算全部轉(zhuǎn)化為代數(shù)運(yùn)算,最終得到配點(diǎn)反饋迭代法的迭代計(jì)算式

式中,P為積分形式的Chebyshev 多項(xiàng)式矩陣,計(jì)算公式為P=

在計(jì)算過程中對式(16)反復(fù)迭代,直至計(jì)算結(jié)果精度達(dá)到目標(biāo)要求.此時(shí),得到的系數(shù)向量L 中每一項(xiàng)即為使用式(9)作為未知函數(shù)估計(jì)解時(shí),滿足精度要求的Chebyshev 多項(xiàng)式對應(yīng)項(xiàng)系數(shù).

Cuthrell 和Biegler[23-24]的相關(guān)研究結(jié)果指出,將整個(gè)函數(shù)待解區(qū)域作為配點(diǎn)區(qū)間,直接得到全局估計(jì)解的方法僅適用于求解光滑問題,而對于非光滑或難以直接全局近似估計(jì)的問題,應(yīng)當(dāng)對全局區(qū)域進(jìn)行離散,并在離散后得到的每個(gè)較小子區(qū)間內(nèi)使用正交多項(xiàng)式進(jìn)行分段估計(jì).考慮到將區(qū)間離散化后在每個(gè)較小子區(qū)間內(nèi)對未知函數(shù)進(jìn)行估計(jì)能夠取得更好的估計(jì)效果,在應(yīng)用配點(diǎn)反饋迭代法時(shí),首先將整個(gè)待解區(qū)間離散化為多個(gè)子區(qū)間,從第一個(gè)子區(qū)間開始對未知函數(shù)進(jìn)行迭代計(jì)算,之后將第一個(gè)子區(qū)間終點(diǎn)函數(shù)值作為下一子區(qū)間初值,再次使用配點(diǎn)反饋迭代法估計(jì)該區(qū)間內(nèi)的未知函數(shù)數(shù)值解,重復(fù)上述過程,直至求得未知函數(shù)在全部待解區(qū)間內(nèi)的數(shù)值解.這一結(jié)合局部離散思想的配點(diǎn)反饋迭代法稱為局部配點(diǎn)反饋迭代法(local collocation feedback iteration method,LCFI),算法計(jì)算過程示意圖見圖1.

圖1 局部配點(diǎn)反饋迭代法示意圖Fig.1 Illustration of LCFI

本節(jié)提出的局部配點(diǎn)反饋迭代法與Picard 迭代法在數(shù)學(xué)上雖然等價(jià),但式(7)通過結(jié)合誤差反饋,將導(dǎo)數(shù)殘差作為修正項(xiàng)引入迭代公式,實(shí)際計(jì)算效率具有明顯優(yōu)勢.此外,在Bai[17]提出的修正Chebyshev-Picard 迭代算法中,由于其利用多項(xiàng)式正交性計(jì)算系數(shù),推導(dǎo)過程十分繁瑣,且迭代公式較為復(fù)雜,增加了編程工作量,降低了算法的計(jì)算效率.本文通過結(jié)合配點(diǎn)法,直接計(jì)算正交多項(xiàng)式系數(shù)向量,在實(shí)際計(jì)算中效率更高.汪雪川[25]通過相關(guān)數(shù)值實(shí)驗(yàn)證明了局部變分迭代法的計(jì)算效率高于修正Chebyshev-Picard 迭代算法,在本文下一節(jié)的數(shù)值仿真中可以看到,本文提出的局部配點(diǎn)反饋迭代法計(jì)算效率遠(yuǎn)高于局部變分迭代算法.上述事實(shí)證明相比于修正Chebyshev-Picard 迭代算法及變分迭代法,本文算法在計(jì)算效率上具有顯著的優(yōu)越性.同時(shí),通過比較本文算法迭代公式與文獻(xiàn)[25]中局部變分迭代算法迭代公式不難發(fā)現(xiàn),本文算法消除了雅可比矩陣,這一特點(diǎn)使得算法能夠更容易的被應(yīng)用于求解復(fù)雜高維問題.

2 擬線性化-局部配點(diǎn)反饋迭代法

局部配點(diǎn)反饋迭代法能夠計(jì)算常微分方程初值問題,但不能直接求解Lambert 問題等邊值問題,本節(jié)介紹一種將擬線性化法及疊加法與局部配點(diǎn)反饋迭代法相結(jié)合的方案,使算法能夠處理兩點(diǎn)邊值約束下的軌道計(jì)算問題.

2.1 擬線性化法

擬線性化法是一種將非線性微分方程近似線性化的數(shù)學(xué)方法[21].

對具有如下形式的二階非線性微分方程邊值問題

式(17)可以改寫為

令y0為未知函數(shù)y 的初始估計(jì)解,yn和yn+1分別表示第n 次和第n+1 次迭代結(jié)果,將第n+1 次迭代結(jié)果在處展開,略去二階及以上偏導(dǎo)數(shù),可以得到將非線性二階微分方程(17)擬線性化后的迭代計(jì)算公式

式(19)對應(yīng)的邊界條件為 yn+1(a)=ya,yn+1(b)=yb.

利用式(19)迭代求解未知函數(shù)y,求解過程中,對于第n 次迭代得到的yn,在代入式(19)進(jìn)行第n+1次迭代時(shí),若將yn視為已知量,則式(19)可視為僅關(guān)于yn+1的二階微分方程,求解該方程即可得到y(tǒng)n+1.重復(fù)上述迭代過程,當(dāng)yn+1=yn時(shí),迭代結(jié)束,此時(shí)yn+1即為微分方程(17)的解.

2.2 疊加法

疊加法能夠?qū)⑽⒎址匠踢呏祮栴}轉(zhuǎn)化為初值問題[26].對于如下形式的微分方程邊值問題

假設(shè)解的形式為

將式(21)代入式(20),可將邊值問題(20)轉(zhuǎn)化為求解使如下兩個(gè)二階常微分方程成立的未知函數(shù)y1與y2

滿足邊界條件的y1與y2可由如下兩組微分方程初值問題分別求得

滿足邊界條件的μ 由下式計(jì)算得到

在計(jì)算時(shí),首先求解式(24)和式(25)兩個(gè)微分方程初值問題得到y(tǒng)1(t)與y2(t),之后將b 點(diǎn)的兩個(gè)函數(shù)值y1(b)與y2(b)代入式(26)計(jì)算μ,最后,將μ及y1(x)與y2(x)代入式(21),得到未知函數(shù)y(x)的解.

對于Lambert 問題等兩點(diǎn)邊值條件約束下的軌道動(dòng)力學(xué)問題,首先利用擬線性化法將其轉(zhuǎn)化為二階線性微分方程邊值問題,再利用疊加法將二階線性微分方程邊值問題轉(zhuǎn)化為一組初值問題,使用局部配點(diǎn)反饋迭代法對初值問題分別求解,再依據(jù)式(21)即可得到兩點(diǎn)邊值問題的解.

3 數(shù)值仿真結(jié)果

3.1 軌道遞推

軌道遞推問題的動(dòng)力學(xué)方程為

式中,r=[x,y,z]T,表示航天器的位置矢量,μ=3 986 004.418 × 108m3/s2,表示地球引力常數(shù),r=‖r‖,表示地心與航天器質(zhì)心間的距離,t0為軌道遞推初始時(shí)刻.aJ為攝動(dòng)項(xiàng),本算例中考慮J2項(xiàng)攝動(dòng).

軌道遞推為初值問題,在給定初值條件的情況下,可以使用局部配點(diǎn)反饋迭代法直接求解.本文通過連續(xù)100 次仿真實(shí)驗(yàn)比較算法性能.為了與同類方法的計(jì)算性能進(jìn)行對比,本算例使用文獻(xiàn)[21]中軌道遞推問題初始位置,相關(guān)計(jì)算參數(shù)見表1,100 次仿真平均單次計(jì)算時(shí)間比較結(jié)果見表2.

表1 遞推軌道計(jì)算參數(shù)Table 1 Calculation parameters of LCFI in orbit propagation

表2 軌道遞推計(jì)算效率對比Table 2 Comparation of efficiency on orbit propagation

由表2 可見,在相同配點(diǎn)條件及計(jì)算精度下,針對軌道遞推問題,局部配點(diǎn)反饋迭代算法(LCFI)的計(jì)算速度為局部變分迭代法(LVIM)的1.5 倍以上(1.78 倍),在相同精度要求下,LCFI 算法的計(jì)算效率更高.進(jìn)一步的仿真結(jié)果表明,當(dāng)配點(diǎn)個(gè)數(shù)上升到32 個(gè)時(shí),在表1 計(jì)算精度條件下,LCFI 算法計(jì)算步長可以達(dá)到12 000 s,而基于有限差分原理的ode113 算法與ode45 算法計(jì)算步長均小于1200 s(最大步長分別為889.9 s 與1 147.3 s),本文算法在上述問題中能將計(jì)算步長提高到有限差分類算法的10 倍以上.

LCFI 與LVIM 計(jì)算結(jié)果及誤差見圖2～ 圖4,對于7 日內(nèi)的軌道遞推結(jié)果,兩種算法所得目標(biāo)軌道最大絕對誤差小于0.4 m,最大相對誤差小于4.9 ×10-8,這樣的計(jì)算精度在工程實(shí)際中是可以接受的.同時(shí),從圖3 絕對誤差變化曲線與圖4 相對誤差變化曲線可以看出,算法在迭代過程中能夠?qū)φ`差進(jìn)行周期性自動(dòng)修正,從而降低誤差發(fā)散速度.

圖2 LCFI 與LVIM 所得遞推軌道對比Fig.2 Comparison between orbits propagated via LCFI and LVIM

圖3 LCFI 與LVIM 所得遞推軌道絕對誤差Fig.3 Absolute error of the orbits propagated via LCFI and LVIM

圖4 LCFI 與LVIM 所得遞推軌道相對誤差Fig.4 Relative error of the orbits propagated via LCFI and LVIM

3.2 攝動(dòng)Lambert 軌道

Lambert 軌道轉(zhuǎn)移問題是經(jīng)典的軌道力學(xué)問題,也是航天器軌道動(dòng)力學(xué)快速計(jì)算方法研究中常用的求解對象[21,25].本節(jié)對攝動(dòng)Lambert 轉(zhuǎn)移軌道及轉(zhuǎn)移初速度進(jìn)行求解,并與擬線性化-局部變分迭代算法求解結(jié)果及計(jì)算效率進(jìn)行對比,驗(yàn)證LCFI 算法的精確性及快速性.

攝動(dòng)Lambert 問題的動(dòng)力學(xué)方程與式(27)所示二體軌道動(dòng)力學(xué)系統(tǒng)的遞推問題相同,將式(27)右側(cè)記為g(r),擬線性化處理后的Lambert 問題動(dòng)力學(xué)方程及其邊界條件表達(dá)式為

依據(jù)疊加法,將 rn+1寫為如下形式

其中 r1與 r2利用如下兩個(gè)微分方程初值問題求得

式(29)中μ 由下式求得

為與同類方法的計(jì)算性能進(jìn)行對比,本算例同樣使用文獻(xiàn)[21]中Lambert 轉(zhuǎn)移問題對應(yīng)的始末位置條件,坐標(biāo)系采用地球固連坐標(biāo)系,經(jīng)過連續(xù)100 次仿真,單次平均用時(shí)見表3,初始條件及相關(guān)計(jì)算參數(shù)見表4,計(jì)算所得軌道結(jié)果及不同算法所得軌道計(jì)算誤差見圖5～ 圖7.

圖5 LCFI 與LVIM 所得Lambert 轉(zhuǎn)移軌道結(jié)果對比Fig.5 Comparison between Lambert orbits obtained via LCFI and LVIM

圖6 LCFI 與LVIM 所得Lambert 轉(zhuǎn)移軌道絕對誤差Fig.6 Absolute error of Lambert orbits obtained via LCFI and LVIM

圖7 LCFI 與LVIM 所得Lambert 轉(zhuǎn)移軌道相對誤差Fig.7 Relative error of Lambert orbits obtained via LCFI and LVIM

表3 攝動(dòng)Lambert 問題計(jì)算效率對比Table 3 Comparation of efficiency on Lambert problem

表4 攝動(dòng)Lambert 轉(zhuǎn)移軌道計(jì)算參數(shù)Table 4 Calculation parameters of LCFI in perturbed Lambert problem

表3 表明,在相同計(jì)算參數(shù)及計(jì)算精度條件下,針對攝動(dòng)Lambert 軌道轉(zhuǎn)移問題,局部配點(diǎn)反饋迭代法(LCFI)計(jì)算速度為文獻(xiàn)[21]中擬線性化-局部變分迭代算法(QL-LVIM) 的2 倍以上(2.68 倍).

在計(jì)算結(jié)果精度方面,LCFI 算法計(jì)算得到的軌道轉(zhuǎn)移初速度為[-1 687.309 900 000,-0.024 275 234,2 922.608 000 000] m/s,LVIM 算法計(jì)算得到的軌道轉(zhuǎn)移初速度為[-1 687.309 900 000,-0.024 305 875,2 922.607 900 000] m/s,兩種算法求得的轉(zhuǎn)移初速度絕對誤差二范數(shù)小于1.05 × 10-4m/s,相對誤差二范數(shù)小于3.10 × 10-8,速度計(jì)算誤差遠(yuǎn)小于實(shí)際測量與控制精度,這樣的計(jì)算誤差在實(shí)際應(yīng)用中是可以忽略的.在計(jì)算誤差變化趨勢方面,由圖6 及圖7 可見,針對本算例中的攝動(dòng)Lambert 軌道,兩種算法對位置計(jì)算結(jié)果的最大絕對誤差不超過0.83 m,最大相對誤差不超過2.1×10-8,且整個(gè)計(jì)算過程中誤差經(jīng)過初期增長后受到抑制,在后續(xù)計(jì)算過程中呈下降趨勢.

3.3 圓型限制性三體模型約束下轉(zhuǎn)移軌道

當(dāng)兩個(gè)主要天體圍繞其系統(tǒng)質(zhì)心做圓周運(yùn)動(dòng),同時(shí)存在一個(gè)質(zhì)量可忽略的第三天體在兩個(gè)主要天體運(yùn)動(dòng)的公共平面內(nèi)運(yùn)動(dòng),這樣的問題就構(gòu)成了平面圓形限制性三體問題[27].月球探測航天器在地-月轉(zhuǎn)移過程中,運(yùn)動(dòng)范圍始終位于地球影響球內(nèi),而地-月系本身圍繞系統(tǒng)質(zhì)心旋轉(zhuǎn),因此在初步設(shè)計(jì)地-月間循環(huán)軌道時(shí),可以將航天器的運(yùn)動(dòng)簡化為地-月-航天器圓形限制性三體問題模型(circular restricted three body problem,CR3BP),該模型也是描述星際探測航天器軌道轉(zhuǎn)移問題的經(jīng)典非線性模型[28-33],本節(jié)通過解算該模型驗(yàn)證LCFI 算法的有效性.

對于地球、月球、航天器構(gòu)成的三體系統(tǒng),令地-月系質(zhì)心為坐標(biāo)系原點(diǎn),x 軸正方向由地球質(zhì)心指向月球質(zhì)心,地-月系軌道平面為(x,y)平面,z 軸與地月系統(tǒng)軌道面垂直,建立空間右手直角坐標(biāo)系.設(shè)地球質(zhì)量為m1,月球質(zhì)量為m2,地球到坐標(biāo)系原點(diǎn)距離為l1,月球到坐標(biāo)系原點(diǎn)距離為l2,將地月間距離及地月系質(zhì)量均無量綱化為1,即令月球的無量綱質(zhì)量為μm,地球的無量綱質(zhì)量為1-μm,地球到原點(diǎn)的距離為μ,月球到原點(diǎn)的距離為1-μ,航天器的無量綱形式運(yùn)動(dòng)方程為[33]

式(33)表征的圓形限制性三體問題動(dòng)力學(xué)模型中,在其L1,L2以及L3拉格朗日點(diǎn)附近存在對應(yīng)于系統(tǒng)某一周期解的Halo 軌道.由于Halo 軌道附近存在著類似管道的不變流形,航天器在流向Halo 軌道的穩(wěn)定流形中運(yùn)動(dòng)時(shí)幾乎不需要消耗能量,能夠顯著節(jié)約燃料[34],因此Halo 軌道及其不變流形常常被用于設(shè)計(jì)低成本登月軌道[9,35-38].

在由近地軌道或繞月軌道向不變流形轉(zhuǎn)移的過程中,需要對航天器運(yùn)動(dòng)狀態(tài)進(jìn)行多次調(diào)整,這一過程中的轉(zhuǎn)移軌道需要精確計(jì)算[21].本節(jié)算例通過計(jì)算探月航天器從185 km 高度地球停泊軌道向流向L1點(diǎn)的穩(wěn)定流形機(jī)動(dòng)時(shí)的轉(zhuǎn)移軌道,驗(yàn)證算法在圓形限制性三體問題動(dòng)力學(xué)模型中軌道計(jì)算的的高效性.為了與同類方法的計(jì)算性能進(jìn)行對比,本算例使用文獻(xiàn)[21]中該算例計(jì)算參數(shù),具體計(jì)算參數(shù)見表5,其中轉(zhuǎn)移時(shí)間以地月系統(tǒng)運(yùn)動(dòng)周期作為系統(tǒng)的無量綱時(shí)間2π,相關(guān)轉(zhuǎn)移時(shí)間的計(jì)算以此為基準(zhǔn),連續(xù)100 次計(jì)算時(shí)間比較結(jié)果見表6.計(jì)算所得軌道結(jié)果及不同算法所得軌道計(jì)算誤差見圖8～ 圖10.

表5 圓型限制性三體模型約束下轉(zhuǎn)移軌道計(jì)算參數(shù)Table 5 Calculation parameters of LCFI in transfer orbit calculation problem with the constraint of circular restricted three-body model

表6 圓形限制性三體模型約束下轉(zhuǎn)移軌道計(jì)算效率對比Table 6 Comparation of efficiency on transfer orbit calculation problem with the constraint of three-body model

圖8 LCFI 與LVIM 軌道計(jì)算結(jié)果對比Fig.8 Comparison between orbits obtained via LCFI and LVIM

圖9 LCFI 與LVIM 對平面圓形限制性三體模型約束下轉(zhuǎn)移軌道計(jì)算絕對誤差Fig.9 Absolute error of transfer orbits under the constraint of circular restricted three body model obtained via LCFI and LVIM

圖10 LCFI 與LVIM 對平面圓形限制性三體模型約束下轉(zhuǎn)移軌道計(jì)算相對誤差Fig.10 Relative error of transfer orbits under the constraint of circular restricted three body model obtained via LCFI and LVIM

由表6 可見,在相同計(jì)算參數(shù)條件下,針對平面圓形限制性三體模型約束下的軌道轉(zhuǎn)移問題,LCFI 的計(jì)算速度為LVIM 的6 倍以上(6.92 倍),即在相同精度要求下,LCFI 算法的計(jì)算效率更高.從圖9 及圖10 誤差變化趨勢來看,誤差在計(jì)算初期增長后迅速受到抑制,最大絕對誤差小于279.5 m,最大相對誤差小于5.2 × 10-6.

4 變參數(shù)計(jì)算方案及仿真實(shí)例

本文提出的局部配點(diǎn)反饋迭代法相比于傳統(tǒng)數(shù)值積分方法具有大步長計(jì)算特性,為發(fā)揮算法大步長收斂優(yōu)勢,在滿足工程計(jì)算精度要求的基礎(chǔ)上盡可能減少計(jì)算量,本文借鑒ph 網(wǎng)格細(xì)化法(ph mesh refinement method)[39-40]計(jì)算思想及針對該方法的相關(guān)改進(jìn)策略,將變參數(shù)計(jì)算方案引入本文算法,從而使得局部配點(diǎn)反饋迭代法能夠自適應(yīng)選擇具有更高計(jì)算效率的參數(shù)組合.

4.1 變參數(shù)計(jì)算方案

參數(shù)變化算法流程圖見圖11,具體的自適應(yīng)參數(shù)調(diào)節(jié)過程分為兩部分.

圖11 變參數(shù)局部變分迭代算法流程圖Fig.11 Flow chart of adaptive local collocation feedback iteration method

(1)在一個(gè)子區(qū)間內(nèi),迭代計(jì)算精度大于迭代終止誤差限,但迭代次數(shù)超過預(yù)定的迭代計(jì)算終止次數(shù)Nmax時(shí),將當(dāng)前區(qū)間的計(jì)算步長減小為原來的a 倍(a＜1),并對當(dāng)前計(jì)算區(qū)間內(nèi)各配點(diǎn)時(shí)刻對應(yīng)的函數(shù)值重新進(jìn)行計(jì)算.

(2)在一個(gè)子區(qū)間內(nèi),迭代計(jì)算精度達(dá)到迭代終止誤差限,且迭代次數(shù)小于預(yù)設(shè)定的迭代計(jì)算終止次數(shù)Nmax時(shí),首先對子區(qū)間內(nèi)M 個(gè)節(jié)點(diǎn)狀態(tài)進(jìn)行插值,利用插值函數(shù)計(jì)算子區(qū)間配置M+1 個(gè)節(jié)點(diǎn)時(shí)各節(jié)點(diǎn)狀態(tài)估計(jì)值.其次對子區(qū)間內(nèi)M 個(gè)節(jié)點(diǎn)的導(dǎo)數(shù)估計(jì)結(jié)果進(jìn)行插值,并利用插值函數(shù)計(jì)算子區(qū)間配置M+1 個(gè)節(jié)點(diǎn)時(shí)各節(jié)點(diǎn)的導(dǎo)數(shù)估計(jì)值,對導(dǎo)數(shù)估計(jì)值積分后得到子區(qū)間內(nèi)M+1 個(gè)節(jié)點(diǎn)函數(shù)狀態(tài)估計(jì)值.根據(jù)相關(guān)研究結(jié)果,兩次估計(jì)結(jié)果差的模值與當(dāng)前對M 個(gè)節(jié)點(diǎn)的估計(jì)結(jié)果的誤差基本相等[39],因此可以使用該方法對當(dāng)前計(jì)算結(jié)果的誤差進(jìn)行估計(jì),上述相對誤差計(jì)算過程的具體計(jì)算公式如下

當(dāng)計(jì)算結(jié)果的相對誤差向量er中值最大的元素高于規(guī)定的誤差上限e1時(shí),此時(shí)計(jì)算結(jié)果誤差過大,需要增加計(jì)算精度以滿足計(jì)算精度要求.當(dāng)計(jì)算結(jié)果的相對誤差向量er中最大元素低于規(guī)定的誤差下限e2時(shí),需要適當(dāng)降低計(jì)算精度以提高計(jì)算速度.依據(jù)ph 網(wǎng)格細(xì)化法,需用節(jié)點(diǎn)個(gè)數(shù)計(jì)算公式為

式中,er＞e1時(shí)P=lg(er/e1),er＜e2時(shí)P=lg(er/e2),M2為計(jì)算所得的需用節(jié)點(diǎn)個(gè)數(shù),M1為當(dāng)前的節(jié)點(diǎn)個(gè)數(shù).

當(dāng)計(jì)算結(jié)果的相對誤差er大于規(guī)定的誤差上限e1且需用節(jié)點(diǎn)個(gè)數(shù)大于節(jié)點(diǎn)數(shù)上限Mmax時(shí),令當(dāng)前計(jì)算區(qū)間的步長降為當(dāng)前步長的b 倍(b＜1),節(jié)點(diǎn)數(shù)等于最小節(jié)點(diǎn)個(gè)數(shù),并重新計(jì)算當(dāng)前區(qū)間內(nèi)各個(gè)節(jié)點(diǎn)函數(shù)值.當(dāng)計(jì)算結(jié)果的相對誤差er大于規(guī)定的誤差上限e1且需用節(jié)點(diǎn)個(gè)數(shù)小于節(jié)點(diǎn)數(shù)下限Mmax時(shí),令節(jié)點(diǎn)數(shù)等于需用節(jié)點(diǎn)個(gè)數(shù),并重新計(jì)算當(dāng)前區(qū)間結(jié)果.當(dāng)計(jì)算結(jié)果的相對誤差er低于規(guī)定的誤差下限e2且需用節(jié)點(diǎn)個(gè)數(shù)小于節(jié)點(diǎn)數(shù)下限Mmin時(shí),令當(dāng)前計(jì)算區(qū)間的步長增加到當(dāng)前步長的c 倍(c＞1),節(jié)點(diǎn)數(shù)等于最小節(jié)點(diǎn)個(gè)數(shù),并重新計(jì)算當(dāng)前區(qū)間結(jié)果.當(dāng)計(jì)算結(jié)果的相對誤差er低于規(guī)定的誤差下限e2且需用節(jié)點(diǎn)個(gè)數(shù)大于節(jié)點(diǎn)數(shù)下限Mmin時(shí),令節(jié)點(diǎn)數(shù)等于需用節(jié)點(diǎn)個(gè)數(shù),并令當(dāng)前計(jì)算區(qū)間的步長增加到當(dāng)前步長的d 倍(d＞1),重新計(jì)算當(dāng)前區(qū)間結(jié)果.當(dāng)計(jì)算結(jié)果的相對誤差er高于規(guī)定的誤差下限e2且低于規(guī)定的誤差上限e1時(shí),保留當(dāng)前區(qū)間計(jì)算結(jié)果與計(jì)算參數(shù)設(shè)置情況,并進(jìn)行下一區(qū)間的計(jì)算.

4.2 仿真實(shí)例

近年來,微小衛(wèi)星在航天工程領(lǐng)域受到了廣泛關(guān)注.缺乏推進(jìn)系統(tǒng)的微小衛(wèi)星在軌運(yùn)行狀況發(fā)射前需要依賴軌道動(dòng)力學(xué)模型經(jīng)過數(shù)值計(jì)算得到.目前,微小衛(wèi)星大多運(yùn)行于近地軌道,在該軌道上航天器受到的最主要干擾力為地球非球形攝動(dòng),其大小比其他攝動(dòng)力高出3 個(gè)量級(jí)以上[41],本文基于40 階EGM2008 球諧重力場模型,在相同初始計(jì)算參數(shù)及迭代計(jì)算精度下,計(jì)算兩組不同初始條件下衛(wèi)星軌道經(jīng)過12 960 400 s (15 天)后的目標(biāo)位置.初始計(jì)算參數(shù)見表7 及表8,變參數(shù)LCFI 相關(guān)計(jì)算參數(shù)見表9,計(jì)算用時(shí)見表10,不同算法計(jì)算結(jié)果見圖12 及圖13,計(jì)算過程中參數(shù)變化情況見圖14～ 圖17.

表7 圓軌道遞推初始計(jì)算參數(shù)Table 7 Calculation parameters of orbit propagation

表8 橢圓軌道遞推初始計(jì)算參數(shù)Table 8 Calculation parameters of orbit propagation

表9 變參數(shù)LCFI 遞推軌道計(jì)算參數(shù)Table 9 Calculation parameters of adaptive LCFI in orbit propagation

表10 不同算法計(jì)算時(shí)間Table 10 Time cost of different numerical method

圖12 變參數(shù)LCFI 及定參數(shù)LCFI 圓軌道遞推計(jì)算結(jié)果Fig.12 Circular orbit obtained via Adaptive LCFI and LCFI

圖13 變參數(shù)LCFI 及定參數(shù)LCFI 橢圓軌道遞推計(jì)算結(jié)果Fig.13 Elliptical orbit obtained via adaptive LCFI and LCFI

圖14 變參數(shù)LCFI 圓軌道遞推計(jì)算步長變化Fig.14 Step size of ALCFI in the propagation of a circular orbit

圖15 變參數(shù)LCFI 圓軌道遞推計(jì)算配點(diǎn)數(shù)變化Fig.15 Variation of the collocation point number in the propagation of a circular orbit

圖16 橢圓軌道遞推計(jì)算步長變化Fig.16 Step size of ALCFI in the propagation of an elliptical orbit

圖17 橢圓軌道遞推計(jì)算配點(diǎn)數(shù)變化Fig.17 Variation of the collocation point number in the propagation of an elliptical orbit

根據(jù)圖14～圖17 可見,針對考慮40 階EGM2008地球球諧重力場模型的軌道遞推問題,在相同初始計(jì)算參數(shù)及迭代計(jì)算精度條件下,變參數(shù)LCFI 算法能夠自適應(yīng)選擇更優(yōu)的計(jì)算步長及基函數(shù)階次,從而降低計(jì)算量.

根據(jù)表10 可見,變參數(shù)LCFI 計(jì)算速度是定參數(shù)LCFI 的6 倍以上(6.57 倍與7.04 倍),定參數(shù)LCFI 由于初始計(jì)算參數(shù)設(shè)置的合理性不足,限制了算法自身大步長高速計(jì)算特性的發(fā)揮,該問題通過本節(jié)變參數(shù)計(jì)算方案得到了有效解決.

從圖12 及圖13 可以看出,由于變參數(shù)計(jì)算方案中引入了誤差評估機(jī)制,在航天器高精度軌道預(yù)測問題中,相同迭代計(jì)算精度下,變參數(shù)局部配點(diǎn)反饋迭代法的(ALCFI)計(jì)算結(jié)果基本落在精確參考軌道附近,而采用定參數(shù)局部配點(diǎn)反饋迭代法的(LCFI)的軌道計(jì)算結(jié)果則產(chǎn)生了更大程度的偏離(圖12 及圖15 中精確參考軌道均使用Matlab 軟件中ode45 程序計(jì)算得到,ode45 設(shè)置精確參考軌道計(jì)算相對誤差2.3 × 10-14,絕對誤差1 × 10-14m),表明變參數(shù)計(jì)算方案引入的誤差評估與調(diào)節(jié)機(jī)制能夠進(jìn)一步提高局部配點(diǎn)反饋迭代法計(jì)算結(jié)果精度.

5 結(jié)論

針對在軌航天器軌道動(dòng)力學(xué)系統(tǒng)快速解算需求,本文基于積分補(bǔ)償?shù)枷?提出了一種適用于軌道動(dòng)力學(xué)快速計(jì)算的高性能數(shù)值積分方法,有效克服了傳統(tǒng)數(shù)值差分類方法的缺陷.主要工作與結(jié)論總結(jié)如下.

(1) 結(jié)合Picard 迭代法、誤差反饋思想和時(shí)域配點(diǎn)法,提出了一種能夠高效求解微分方程初值問題的局部配點(diǎn)反饋迭代法.

(2) 將局部配點(diǎn)反饋迭代法和擬線性化法及疊加法相結(jié)合,求解了軌道轉(zhuǎn)移兩點(diǎn)邊值問題.

(3) 通過解算遞推軌道、攝動(dòng)Lambert 軌道轉(zhuǎn)移問題以及圓型限制性三體模型約束下的轉(zhuǎn)移軌道問題驗(yàn)證了本文算法的有效性.計(jì)算結(jié)果顯示在相同迭代精度及計(jì)算參數(shù)設(shè)置條件下,本文算法在上述軌道問題中計(jì)算效率比局部變分迭代算法提高50%以上.

(4) 將變參數(shù)計(jì)算方案引入局部配點(diǎn)反饋迭代法,使算法能夠自適應(yīng)調(diào)節(jié)計(jì)算參數(shù).相關(guān)計(jì)算結(jié)果表明,本文引入的變參數(shù)計(jì)算方案能夠在保證計(jì)算精度條件下發(fā)揮算法大步長快速收斂優(yōu)勢,并進(jìn)一步提高算法計(jì)算精度.

局部配點(diǎn)反饋迭代法對軌道動(dòng)力學(xué)初值問題及兩點(diǎn)邊值問題的計(jì)算效率相較于現(xiàn)有方法具有明顯優(yōu)勢,在計(jì)算能力較弱的星載計(jì)算機(jī)中[42],該優(yōu)勢對計(jì)算時(shí)間的節(jié)約程度將更為顯著,這對執(zhí)行失效衛(wèi)星抓捕等空間快速機(jī)動(dòng)任務(wù)具有重要意義.在后續(xù)工作中,將探究局部配點(diǎn)反饋迭代法的并行算法設(shè)計(jì)等改進(jìn)方案,以進(jìn)一步提高算法性能.