大范圍收斂的攝動Lambert問題新型解法：擬線性化-局部變分迭代法

馮浩陽,岳曉奎,汪雪川,*

1.西北工業(yè)大學 航天飛行動力學技術國家級重點實驗室,西安 710072

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

Lambert問題是航天動力學領域的經(jīng)典問題,很多學者對其進行研究并提出了不同的解決方法。在工程領域,由于攝動Lambert問題不具備解析解,因此主要通過數(shù)值方法進行求解。常用的攝動Lambert問題解法主要包括兩類,一類是將其轉(zhuǎn)換為初值問題進行迭代求解,另一類是通過離散化方法將其轉(zhuǎn)換為非線性代數(shù)方程組進行求解。

將邊值問題轉(zhuǎn)換為初值問題的方法中,最具代表性的方法即打靶法。打靶法[1-2](Shooting Method)把兩點邊值問題轉(zhuǎn)化為初值問題,并使用牛頓法(Newton’s Method)對初速度進行迭代修正,該方法雖然具有收斂速度快的優(yōu)點,但其收斂域很小且對初值非常敏感,實際應用缺陷較大。多重打靶法[3]融合了有限差分法和打靶法,它將整個求解區(qū)間劃分為多個子區(qū)間,并在各個子區(qū)間上求解初值問題,相對于簡單打靶法,該方法具有更大的收斂域,但仍具有計算不穩(wěn)定的缺陷,且在各個子區(qū)間的邊界上存在速度不連續(xù)的問題。隱式打靶法[4]在間接法思想的基礎上,借助尺度變換法和變分法,在求解包含隱式終端約束的軌道優(yōu)化問題中得到應用。直接變換法[5-6]能夠?qū)牲c邊值問題轉(zhuǎn)化為一對初值問題,計算簡便,但該方法對微分方程的形式和邊值條件有特定要求,對于一般的非線性兩點邊值問題并不適用。其他方法如Gooding給出了迭代變量的最優(yōu)起始公式[7],得到了在所有情形下都可以迅速收斂的高精度Lambert問題算法,Battin通過迭代高斯方程[8]得到Lambert問題的解,但是這些經(jīng)典方法都不適用于一般的攝動Lambert問題。

在邊值問題的離散化方法中,有限差分法[9]的應用較早,通過差分估計將轉(zhuǎn)移軌跡離散為大量的差分節(jié)點,進而對其進行代數(shù)求解。然而有限差分法需要建立大規(guī)模非線性代數(shù)矩陣方程,并且需要存儲大量節(jié)點數(shù)據(jù)。對于大范圍軌道轉(zhuǎn)移的攝動Lambert問題,有限差分法不僅求解困難,其數(shù)值計算精度也較低。配點法[10]也是求解初值問題和兩點邊值問題的重要方法,該方法使用基函數(shù)對時間域內(nèi)的配點插值,從而得到半解析解,且無需存儲大量的離散數(shù)據(jù)[11]。NK/C偽譜法[12]使用Chebyshev正交多項式作為插值基函數(shù),在求解Lambert問題、相對軌道轉(zhuǎn)移等航天動力學問題中得到應用。徑向基函數(shù)(RBFs)方法[13]則使用徑向基函數(shù)對配點插值,在求解二體問題中得到應用。配點法具有精確高效的優(yōu)點,但在計算過程中需要對雅克比矩陣求逆,這通常比較繁瑣。修正Chebyshev-Picard迭代法(MCPI)[14]利用Picard迭代構(gòu)造出修正公式,通過正交多項式插值逼近真實解,避免了矩陣求逆。反饋加速Picard迭代法[15]對修正變分迭代法[16]進行改進,并融合了配點法的思想,使用第一類Chebyshev正交多項式對配點插值,得到了一種可用于攝動軌道遞推和攝動Lambert問題的新型算法。受算法特點所限,該方法的遞推公式需要在整個時間區(qū)間上迭代,為了達到較高計算精度需要設置較多的配點和較多的基函數(shù),使計算量增加。

本文提出一種基于擬線性化和局部變分迭代法的攝動Lambert問題求解方法,該方法能夠一般性的應用于地球衛(wèi)星軌道、相對運動軌道、深空軌道等不同軌道類型的轉(zhuǎn)移問題中,為地球衛(wèi)星的單圈或多圈軌道轉(zhuǎn)移、航天器的交會對接、多航天器編隊飛行、深空探測等任務提供穩(wěn)定、精確、實時的軌道轉(zhuǎn)移算法。在攝動Lambert問題中,使用擬線性化[17](Quasi Linearization)法將攝動Lambert問題轉(zhuǎn)化為迭代形式的線性邊值問題,通過疊加法將線性邊值問題分解為兩個初值問題,并得到初速度和轉(zhuǎn)移軌道的迭代公式。與傳統(tǒng)打靶法相比,其迭代格式更加簡單,計算更加穩(wěn)定,且收斂域顯著增大。在初值問題的解算方面,本文利用了局部變分迭代法[18-19](Local Variational Iteration Method)精確高效的優(yōu)點,在計算過程中不涉及非線性代數(shù)方程組的矩陣求逆操作,因此計算簡便高效。在此基礎上提出的擬線性化-局部變分迭代法(QL-LVIM),經(jīng)過較少的幾次迭代,計算結(jié)果即可達到很高精度。在初值選取方面,避免了常用的通過二體模型確定初始估計的步驟,可以通過更簡潔的方式進行確定。該方法不僅具有與牛頓打靶法相當?shù)目焖偈諗啃?還克服了傳統(tǒng)打靶法初值敏感、收斂域小的缺點,并且在同等計算精度下,能夠顯著提高計算效率。該方法的精確性、穩(wěn)定性和實時性在低軌轉(zhuǎn)移和高軌轉(zhuǎn)移情形下得到了驗證,結(jié)果表明,相對于對比方法,本方法在計算效率和收斂域方面具有明顯優(yōu)勢。通過地-月系三體轉(zhuǎn)移問題對方法的適用性進行進一步驗證。此外,文獻[17]中擬線性化方法僅用于求解一維兩點邊值問題,本文探究了該方法在多維兩點邊值問題求解中的應用。

1 問題描述

1.1 非線性邊值問題的擬線性化

考慮兩點邊值問題,對于如下非線性二階微分方程:

y″=f(t,y,y′)

(1)

滿足如下邊界條件：

y(t0)=y0,y(tf)=yf

(2)

式中：y′=dy/dt,y″=d2y/dt2。

式(1)可以改寫為

φ(t,y,y′,y″)=y″-f(t,y,y′)=0

(3)

為得到迭代方程,令yn和yn+1分別表示第n次和第n+1次迭代結(jié)果,并且均滿足φ=0,對第n次迭代,有

(4)

將第n+1次迭代結(jié)果展開有

φ(t,yn+1,y′n+1,y″n+1)=φ(t,yn,y′n,y″n)+

(5)

略去二階及以上偏導數(shù),式(5)可簡化為

(y″n+1-y″n)+…=0

(6)

將式(4)代入式(6)可得

(7)

對應的邊界條件為

yn+1(t0)=y0,yn+1(tf)=yf

(8)

則式(1)和式(2)構(gòu)成的非線性兩點邊值問題被轉(zhuǎn)化為式(7)和式(8)構(gòu)成的迭代形式的線性兩點邊值問題,通過對其多次求解和迭代,即可不斷逼近原非線性兩點邊值問題的解。在擬線性化處理中,雖然忽略了二階及以上偏導數(shù),但通過多次迭代,仍可以使計算結(jié)果達到較高精度。

1.2 攝動Lambert問題的擬線性化

航天器在繞地球運行中,會受到地球非球形攝動、大氣阻力攝動、日月攝動、太陽光壓攝動等干擾。在近地軌道,主要攝動因素為地球非球形攝動和大氣阻力攝動,在中高軌道,主要攝動因素為地球非球形攝動[20-21],這里考慮地球非球形攝動和大氣阻力攝動,忽略其他高階小攝動力的影響。在赤道慣性坐標系下,航天器的軌道動力學方程為

(9)

(10)

記軌道轉(zhuǎn)移時間為tf,邊值條件為

r(t0)=r0

(11)

r(tf)=rf

(12)

(13)

rn+1(t0)=r0

(14)

rn+1(tf)=rf

(15)

1.3 通過疊加法求解擬線性化形式的Lambert問題

為求解式(13),使用疊加法[17]將rn+1寫為

rn+1=V+W·s

(16)

式中:V和W為分運動的位置矢量;s為飛行器在轉(zhuǎn)移軌道起點處的速度,即

(17)

在t0時刻,初值條件式(14)滿足：

rn+1(t0)=V(t0)+W(t0)·s=r0

(18)

(19)

為應用疊加法,將約束(18)和(19)做如下分解:

式中:03×3為3階零矩陣;03×1為3×1零向量,E3×3為3階單位陣。則式(13)～式(15) 構(gòu)成的兩點邊值問題被分解為如下兩個初值問題。

1) 初值問題Ⅰ

(20)

(21)

2) 初值問題Ⅱ

(22)

(23)

利用邊值條件(15)有

rn+1(tf)=V(tf)+W(tf)·s=rf

(24)

s=W(tf)-1·[rf-V(tf)]

(25)

不同于對初值敏感的牛頓法,擬線性化方法具有較大的收斂域,在確定邊值問題的初始估計解時可以較為隨意,這里以連接2個邊界點的線段作為初始估計解。另外,由于03×3和E3×3均為3×3矩陣,直接求解初值問題Ⅱ比較繁瑣,此時可以將其分解為x、y、z方向上的3個子問題,并分別求解,由3個子問題的求解結(jié)果合并得到W,注意圖1中rn,n=0表示初始估計解而非轉(zhuǎn)移軌道的初始位置。

圖1 擬線性化法求解Lambert問題的流程

1.4 局部變分迭代法對初值問題的求解

經(jīng)典的求解初值問題的方法有歐拉法、龍格庫塔方法、線性多步法等,為得到較高精度的解,這些方法需要設置較小的步長,造成計算耗時較長,并且需要存儲大量的節(jié)點數(shù)據(jù)。局部變分迭代法也是求解初值問題的方法之一,并具有精確高效的優(yōu)點,計算中選取較大的步長就可以達到較高的計算精度,且無需存儲大量節(jié)點數(shù)據(jù),這里使用該方法對以上2個初值問題進行求解。對于如下一階非線性微分方程：

(26)

x(t)=[x1(t),x2(t),…,xd(t),…,xD(t)]T為隨時間變化的D維矢量。首先給出方程解的初始估計,記為x0(τ),隨后在求解區(qū)間內(nèi)對其進行多次修正,從而不斷逼近方程的真實解,迭代公式為

τ]}dτ

(27)

式中:λ(τ)為待定的拉格朗日乘子矩陣。式(27)表明,近似解xn+1包括2部分,一部分為xn,另一部分為xn在t時刻之前累計偏差的加權(quán)平均。若使式(27)右端關于xn不變,則右端的變分為零,由此可以得到如下約束：

(28)

式中：t0≤τ≤t;I為單位陣;J(τ)=?g(xn,τ)/?xn。將λ(τ)在τ=t處做一階泰勒展開,可以得到

λ(τ)≈-I+J(t)(τ-t)

(29)

(30)

通過在時間域配點可以將式(30)離散化,記t1,t2,…,tM為配點時刻,則由式(30)可得

(31)

(32)

式中：Φd=[φd,1(t),φd,2,…,φd,N(t)]為基函數(shù)系;Ad=[αd,1(t),αd,2(t),…,αd,N(t)]T為基函數(shù)的系數(shù)矢量。當x(t)的各分量用相同的正交基函數(shù)表示時,編程會大大簡化,因此將Φd(t)統(tǒng)一寫為Φ(t)。根據(jù)式(32)有

(33)

(34)

對于第n次迭代,

(35)

式中：G(t)=[G1(t),G2(t),…,Gd(t),…,GD(t)]T為D維矢量。通過正交基函數(shù)插值,將Gd表示為

(36)

式中：Φ(t)=[φ1(t),φ2(t),…,φN(t)],Bd=[βd,1(t),βd,2(t),…,βd,N(t)]T,對xd(t)和Gd(t)的插值使用相同的基函數(shù)Φ(t)。在t1,t2,…,tM時刻,由式(36)可得

(37)

(38)

類似地可以得到

(39)

(40)

可以得到迭代公式為

(41)

2 求解驗證

2.1 計算精度和效率分析

本節(jié)使用QL-LVIM對攝動Lambert問題進行求解,并與另外4種求解兩點邊值問題的方法進行對比。打靶法是求解兩點邊值問題的經(jīng)典方法之一,利用牛頓法可以得到初速度的迭代公式,并將兩點邊值問題轉(zhuǎn)化為兩個初值問題,若通過經(jīng)典的四階龍格庫塔方法(RK4)對初值問題進行求解,就得到了牛頓-四階龍格庫塔方法(Newton-RK4),這里將其作為第1種對比方法。同時,將牛頓法與局部變分迭代法結(jié)合,構(gòu)造Newton-LVIM作為第2種對比方法。牛頓法具有初值敏感性,當初始估計與實際初速度偏差不是很大時,才可以保證算法收斂。為了觀察大范圍收斂情況下不同方法的精度和效率,構(gòu)造擬線性化-四階龍格庫塔方法(QL-RK4)作為第3種對比方法。最后,與文獻[15]提出的反饋加速Picard迭代(FAPI)方法對比,該方法具有精度高、效率高的特點。本文使用聯(lián)想筆記本R480(CPU： Intel Core i5-8250U；RAM：8.00 G)安裝的MATLAB R2017a軟件進行數(shù)值仿真,未采用GPU加速和并行計算。

表1 攝動Lambert問題的邊值條件和轉(zhuǎn)移時間

表2 計算參數(shù)

初始估計解可以通過多種方法確定,這里選取連接初始位置和末位置的勻速直線軌道,這種初始估計也被稱作“冷啟動”(詳見附錄A)。求解后得到低軌和高軌的轉(zhuǎn)移軌道如圖2所示,作為參照,調(diào)用MATLAB內(nèi)置的ODE45函數(shù)(將相對誤差和絕對誤差分別設置為1×10-13和1×10-15),根據(jù)QL-LVIM求出的初速度做軌道遞推,遞推結(jié)果也在圖2中給出,可以看到QL-LVIM的計算結(jié)果和ODE45函數(shù)的遞推結(jié)果能夠很好的擬合。

圖2 擬線性化-局部變分迭代法的求解結(jié)果

為進一步分析QL-LVIM的計算精度和計算效率,將其與Newton-RK4、Newton-LVIM、QL-RK4、FAPI的計算誤差和計算時間進行對比,相關參數(shù)和初始估計見表2。求解兩點邊值問題的關鍵在于求出準確的初速度,初速度的精度反映了算法的精度,因此采用如下方式定義誤差：首先求出轉(zhuǎn)移軌道的初速度,再根據(jù)起點位置和初速度使用ODE45函數(shù)(調(diào)至最高精度)做軌道遞推,將遞推出的終點位置與實際終點位置rf進行比較,則二者的偏差反映了初速度的精度和算法的精度,將該偏差定義為誤差。5種方法求出的初速度如表3所示,計算誤差和計算時間如表4所示,其中的計算時間為10次計算的平均值。

表3顯示,不同算法求出的初速度非常接近,5種算法的有效性得到互相驗證。在2種軌道類型下,QL-LVIM、Newton-RK4、Newton-LVIM、FAPI算法求解結(jié)果的小數(shù)點后6位、后8位完全一致,QL-RK4求解結(jié)果的小數(shù)點后2位與其他算法一致,表明QL-LVIM、Newton-RK4、Newton-LVIM、FAPI算法精度較高,QL-RK4精度較低。

表4顯示,在表2的參數(shù)條件下,可以將QL-LVIM、Newton-RK4、Newton-LVIM、FAPI方法的計算誤差都調(diào)整至1×10-6級別,這能夠使不同算法的計算時間更具有可比性。在2種軌道類型下,QL-LVIM的計算時間均最短,其計算速度不同程度的快于其他算法,表明QL-LVIM具有更高的計算效率。而采用QL-RK4方法在經(jīng)過150 s以上的運算后,精度也只能達到個位數(shù)級別,由于更長的計算時間很可能無法滿足航天工程中的實時性需要,因此沒有必要繼續(xù)延長計算時間以提高精度。這一結(jié)論也與表3的結(jié)果吻合。

表3 5種方法計算出的初速度對比

表4 5種方法的計算誤差與計算效率對比

QL-LVIM和FAPI都使用了變分迭代的思想,但前者計算速度更快,這是由于在應用FAPI方法時,必須在整個轉(zhuǎn)移區(qū)間上迭代,由于時間跨度較大,迭代所需的配點個數(shù)和基函數(shù)個數(shù)較多,參數(shù)選擇為M=N=32,使得計算量較大,且迭代次數(shù)較多。而本文方法克服了這一困難,在求解初值問題時可以將整個轉(zhuǎn)移區(qū)間劃分為眾多子區(qū)間,在每個子區(qū)間分別迭代求解,在M=N=10時即可達到同樣的計算精度。

QL-LVIM和QL-RK4都使用了擬線性化方法,但二者的計算精度和計算效率差異較大,這是由于RK4和LVIM具有不同的計算特點所致。具體來說,2種方法在對初速度進行迭代時,都用到了兩個初值問題在終點時刻的計算結(jié)果,在RK4方法中,終點之前所有節(jié)點的誤差都會累積至終點,這造成累積誤差嚴重,計算精度較低,而在LVIM中,前一步長的誤差不會累積到下一步長,即終點處的計算結(jié)果不受終點之前各節(jié)點計算誤差的累積影響,因此精度較高。計算中RK4是利用小步長做單點遞推,LVIM則是在大步長內(nèi)對多個配點同時迭代,因此計算效率大大提高。這說明在應用擬線性化方法之后求解初值問題時,LVIM比RK4更有優(yōu)勢。

Newton-RK4與QL-RK4的計算效率和計算精度也差異較大,這是由于在Newton-RK4中,第1個初值問題的微分方程是原非線性微分方程,第2個初值問題的微分方程是原非線性微分方程關于初速度的偏導[17],非線性性質(zhì)沒有損失,因此計算精度較高。而QL-RK4方法中,2個初值問題求解的均為原非線性微分方程擬線性化后的微分方程,非線性性質(zhì)的丟失造成計算精度的損失,為了提高初速度的精度,就要進一步縮小RK4方法的計算步長,這造成計算量和計算時間增加,但即便如此,計算結(jié)果也無法獲得令人滿意的精度。

雖然Newton-LVIM方法也有較好的計算精度和效率,但該方法僅在初速度估計值比較接近真實值的情況下才有效,這也是Newton打靶法的主要缺陷。相對于Newton法,QL方法的一大優(yōu)勢是具有大收斂域,為分析QL法和Newton法的收斂性,將Newton-LVIM與QL-LVIM進行對比,兩種方法均采用LVIM對初值問題進行求解,以排除初值問題求解方法差異對收斂域的影響,對比結(jié)果在2.2節(jié)給出。

2.2 擬線性化方法與牛頓法的收斂性分析

在工程應用中,算法的收斂性是衡量算法性能優(yōu)劣的重要標準之一,本節(jié)對QL-LVIM和Newton-LVIM的收斂域進行對比分析。由于2種方法的收斂域難以直接精確求出,這里采用蒙特卡羅模擬方法,為2種方法隨機選取大量初值,考察在這些初值條件下兩種方法的收斂情況。如前文所述,QL法通過對初始估計解(一段運行軌跡)不斷迭代得到真實解,Newton法通過對初速度估計值不斷迭代得到精確初速度。由于牛頓法不存在通用的初速度估計方法,這里假設,通過二體模型或歷史數(shù)據(jù)大致知道真實初速度的數(shù)量級,在此基礎上,在以原點為球心、真實初速度的模為半徑的球面上隨機選取Ns個點,以球心到這Ns個點的矢徑做為Ns個初速度估計值,并分別使用Newton-LVIM迭代求解,以驗證各個初速度估計值的收斂性。為增強可比性,利用這Ns個初速度和初始位置r0、運行時間tf構(gòu)造Ns種“冷啟動”(詳見附錄A),做為QL-LVIM的初始估計解,并分析收斂性,相關參數(shù)和模擬結(jié)果如表5和圖3所示。

表5中,對于QL-LVIM方法,在LEO情形下,2 000個初始估計中有1 606個收斂,經(jīng)計算這些收斂結(jié)果的方差為[5.7×10-21, 4.8×10-23, 6.7×-23],表明這些情形均收斂到了真實初速度且具有較高的精度,在HEO情形下,有1 933 個初始估計收斂到了真實初速度,這些結(jié)果的方差為[1.35×10-21, 3.36×10-25, 1.96×10-21],2種情形收斂初值的比例均達到80%以上。由于受大氣阻力影響更大,LEO情形的初值收斂比例略低。在Newton-LVIM方法下,隨機選取的初速度收斂的比例在20%左右,表明在上述隨機采樣方式下,QL法比Newton法具有更大的收斂域。其原因在于,不同于只對初速度進行迭代的Newton法,QL算法對整條估計軌道上的節(jié)點都進行迭代,這使得QL方法收斂域更大。圖3 直觀展示了2種方法下收斂初值的空間分布。

表5 擬線性化法與牛頓法的收斂性對比

圖3 擬線性化法與牛頓法收斂初值分布對比

3 QL-LVIM在三體問題求解中的應用

為驗證本文方法在其他軌道轉(zhuǎn)移問題中的適用性,本節(jié)對地-月系的圓型限制性三體問題進行求解。

從地球到月球Halo軌道的轉(zhuǎn)移通常借助Halo軌道的穩(wěn)定流形[23],以地球到月球L1點Halo軌道的轉(zhuǎn)移為例,通常包括2個階段,第1階段是航天器從地球停泊軌道到穩(wěn)定流形的轉(zhuǎn)移,第2階段是航天器沿著穩(wěn)定流形的運動,這一階段不消耗或較少消耗能量。其中,地球軌道到穩(wěn)定流形的轉(zhuǎn)移軌道要通過大量計算才能確定,傳統(tǒng)方法通常先借助蒙特卡羅模擬或經(jīng)驗方法確定一條估計軌道,再通過打靶法將其修正到精確結(jié)果。但由于邊界條件和轉(zhuǎn)移時間是任意的,這樣的軌道搜尋過程較為耗時和不便,這里使用QL-LVIM方法對該問題進行求解。假設轉(zhuǎn)移軌道的起點位于185 km高度的地球停泊軌道,其他初始條件和計算參數(shù)見表6。在地-月系統(tǒng)的旋轉(zhuǎn)坐標系下,將地球、月球、航天器分別記為P1、P2、P3,則無量綱形式的航天器運動方程為

表6 地-月Halo軌道轉(zhuǎn)移問題參數(shù)

(42)

圖4 地-月系三體轉(zhuǎn)移問題求解結(jié)果

4 結(jié) 論

設計了一種航天器軌道轉(zhuǎn)移問題的新型解法——擬線性化-局部變分迭代法(QL-LVIM)。該方法兼具擬線性化的大收斂域優(yōu)勢和局部變分迭代法的快速收斂優(yōu)勢,能夠?qū)壍擂D(zhuǎn)移問題進行快速、精確、穩(wěn)定求解,較好地克服了牛頓打靶法收斂域小的缺陷和有限差分法計算效率低的缺陷。

1) 對LEO和HEO情形下的攝動Lambert問題的求解結(jié)果表明,QL-LVIM在計算效率方面具有明顯優(yōu)勢。在同等計算精度下,該方法的計算耗時遠低于幾類對比方法。

2) 通過蒙特卡羅模擬對收斂域進行分析,結(jié)果表明本文方法的收斂域遠大于牛頓法,該結(jié)論在三維圖形中得到展示和驗證。

3) 本文方法在計算效率和收斂域方面均具有優(yōu)勢,即使在計算能力較弱的星載計算機上也能快速計算出變軌結(jié)果,較少消耗星載計算機的計算資源。在未來工作中,將使用本方法解決相對運動、深空探測等更復雜空間任務中的邊值問題,并對方法進行優(yōu)化以進一步提高計算性能。

4) 除連續(xù)和可導,本文方法對于待求解的微分方程沒有更多要求,因此對其他領域的多維兩點邊值問題也具有適用性。