潘 迅,泮斌峰
(西北工業(yè)大學航天學院,西安 710072)
在地月系三體模型中,L2點在地月連線上,且始終位于月球背面,不僅可以作為月球探測和深空探測的中轉(zhuǎn)站,還可以在L2點建立月球觀測站,實時監(jiān)測月球背面情況,為月球背面著陸提供相關(guān)信息[1]。我國的嫦娥五號再入返回飛行試驗器的服務(wù)艙曾于2014年11月28日到達地月系L2點,并在2015年1月4日離開L2轉(zhuǎn)移至環(huán)月軌道[2-3]。
嫦娥五號再入返回試驗器利用化學發(fā)動機經(jīng)過多脈沖機動完成軌道轉(zhuǎn)移,轉(zhuǎn)移過程較為簡單,所需時間較短。然而,隨著推進技術(shù)的發(fā)展,電推進、等離子推進等持續(xù)小推力推進方式由于具有比沖大、控制精度高等優(yōu)點,能極大地提高探測器的有效載荷的比重[4],因此將在未來的航天任務(wù),特別是深空探測任務(wù)中,扮演越來越重要的角色。
由于小推力發(fā)動機的推力幅值小,持續(xù)時間長,完成轉(zhuǎn)移所需的軌道圈數(shù)多,在轉(zhuǎn)移過程中初始狀態(tài)量對轉(zhuǎn)移軌道有較大影響,即對初始值很敏感,因此優(yōu)化過程存在較大困難。此外,相比于傳統(tǒng)的二體模型,限制性三體問題的非線性更強,導致該小推力轉(zhuǎn)移軌道優(yōu)化問題的求解難度增加。小推力轉(zhuǎn)移軌道優(yōu)化方法主要可以分為直接法、間接法和混合法[4]。直接法通過將狀態(tài)量和控制量進行離散,將軌道優(yōu)化問題轉(zhuǎn)化為非線性規(guī)劃問題,直接對離散參數(shù)進行尋解[5-6];間接法基于最優(yōu)控制理論的經(jīng)典變分法推導得到最優(yōu)控制的一階必要條件,將其轉(zhuǎn)換為兩點邊值問題,通過數(shù)值方法對該兩點邊值問題進行求解,從而得到轉(zhuǎn)移軌道[7-8];混合法只利用間接法推導得到的最優(yōu)控制律,但不包含其他最優(yōu)性條件,通過參數(shù)優(yōu)化方法進行求解[9]。其中,間接法由于利用了最優(yōu)性條件,可以保證所得解的局部最優(yōu)性,且計算精度高。然而間接法將軌道優(yōu)化問題轉(zhuǎn)換成兩點邊值問題后,雖然理論上可以利用打靶法對該問題直接求解,但由于作為優(yōu)化變量的協(xié)態(tài)變量沒有物理意義,且算法的收斂域較小,因此,對于推力值小、轉(zhuǎn)移圈數(shù)多的問題,求解過程存在很大困難。
為了增加優(yōu)化算法的收斂域,降低求解難度,同倫方法被廣泛應(yīng)用于軌道優(yōu)化。同倫方法,是指從較簡單問題出發(fā),通過改變同倫參數(shù),求解對應(yīng)的一系列子問題,最終得到原問題的解。其在小推力在轉(zhuǎn)移軌道中的應(yīng)用主要有兩個方向:一是針對燃料最優(yōu)問題,通過對性能指標進行同倫,先求解連續(xù)光滑控制問題,然后通過同倫迭代,最終得到不連續(xù)的Bang-Bang控制的燃料最優(yōu)轉(zhuǎn)移軌道[10-11];二是針對推力很小的軌道優(yōu)化問題,先求解較簡單的大推力值對應(yīng)的軌道優(yōu)化問題,然后通過對推力幅值進行同倫,逐漸減小推力幅值,最終得到所需的很小推力幅值的轉(zhuǎn)移軌道[12-13]。
目前對限制性三體問題中的小推力轉(zhuǎn)移軌道優(yōu)化問題研究較少,Caillau等[12]對從航天器從GEO軌道到地月系L1點和月球軌道的轉(zhuǎn)移軌道進行相關(guān)研究,本文將其擴展到了航天器從地球GEO軌道到地月系L2點的時間最優(yōu)轉(zhuǎn)移軌道優(yōu)化問題。首先基于最優(yōu)控制理論,將地月系三體問題模型下的小推力轉(zhuǎn)移軌道問題轉(zhuǎn)換為兩點邊值問題,然后結(jié)合同倫方法對該兩點邊值問題進行求解。
考慮航天器同時受到地球引力和月球引力,忽略地球非球形攝動、太陽引力、太陽光壓等其他攝動,并假設(shè)月球圍繞地球的運行軌道為圓軌道,則可用圓限制性三體問題來描述航天器的運動。僅考慮航天器在月球軌道面內(nèi)的平面運動,并且忽略轉(zhuǎn)移過程中航天器的質(zhì)量變化,則航天器在旋轉(zhuǎn)坐標系下的動力學模型為式(1)~(2):
(1)
(2)
式中:r=[x,y]T和v=[vx,vy]T分別為航天器的位置矢量和速度矢量,m0為航天器初始質(zhì)量,Tmax為航天器可以提供的最大推力,u為轉(zhuǎn)移過程中的實際推力大小,α=[αx,αy]T為單位推力方向。則可定義航天器有效勢能函數(shù)U如式(3),g(r)和h(v)分別表示為式(4)、(5):
(3)
(4)
h(v)=[2vy,-2vx]T
(5)
需要注意的是,在動力學模型式(1)~(2)中,各變量均為單位化后的值,其中長度歸一化單位為1 DU=3.8440×105km,時間歸一化單位為1 TU=3.75676967×105s,質(zhì)量歸一化單位為1 MU=M0kg(M0為航天器實際質(zhì)量)。歸一化后,地球在旋轉(zhuǎn)坐標系的位置為(-μ,0),月球的位置為(1-μ,0)。
本文針對航天器從地球GEO軌道到L2平動點的時間最優(yōu)小推力轉(zhuǎn)移軌道進行研究。在初始時刻t0,航天器在GEO軌道上的位置自由,其狀態(tài)量可以表示為式(6):
x0=[r0cosθ0-μ,r0sinθ0,-v0sinθ0,v0cosθ0]T
(6)
時間最優(yōu)問題的性能指標可以表示為式(7):
(7)
時間最優(yōu)小推力轉(zhuǎn)移軌道優(yōu)化問題的本質(zhì),就是在滿足初始狀態(tài)約束式(6)和終端狀態(tài)約束xf的條件下,通過優(yōu)化變量u和α,使性能指標最小。根據(jù)龐德里亞金極小值原理,引入拉格朗日乘子,構(gòu)造哈密爾頓函數(shù)為式(8):
(8)
式中:λr=[λx,λy]T和λv=[λvx,λvy]T也稱為協(xié)態(tài)變量,其微分方程為式(9)~(10):
(9)
(10)
為使哈密爾頓函數(shù)為極小值,從式(8)可知,推力大小和推力方向的最優(yōu)控制律應(yīng)為式(11)~(12):
u*=1
(11)
(12)
對于航天器初始時刻的狀態(tài)量x0,由于存在優(yōu)化變量θ0,需要對橫截條件進行推導。根據(jù)最優(yōu)控制理論,可以得到對應(yīng)的橫截條件為式(13):
(x0+μ)λy0-y0λx0+vx0λvy0-vy0λvx0=0
(13)
對于時間最優(yōu)問題,轉(zhuǎn)移時間自由,對應(yīng)的橫截條件為式(14):
H(tf)=0
(14)
至此,將最優(yōu)控制問題轉(zhuǎn)換為兩點邊值問題,其初始時刻的狀態(tài)值由式(6)確定,終端時刻的約束包括狀態(tài)約束xf和橫截條件式(13)~(14)。該兩點邊值問題可利用打靶法進行求解,打靶函數(shù)表示為式(15):
F(y)=O
(15)
式中y=[θ0,λx0,λy0,λvx0,λvy0,tf]T為優(yōu)化變量。
將時間最優(yōu)小推力轉(zhuǎn)移軌道優(yōu)化問題轉(zhuǎn)換為兩點邊值問題后,雖然理論上可以利用打靶法對式(15)進行求解得到優(yōu)化變量y的值,再對初值進行積分得到轉(zhuǎn)移軌道。然而推力越小,轉(zhuǎn)移時間越長,轉(zhuǎn)移圈數(shù)越多,打靶過程對初始猜測值就越敏感,打靶的收斂域就越小。此外,相比于傳統(tǒng)的只考慮中心引力場的二體模型,限制性三體模型的非線性更強,從而增加了求解難度。因此,為了克服小推力轉(zhuǎn)移軌道優(yōu)化問題的初始值敏感問題,本文利用同倫方法進行求解。
基于所需求解的原問題,構(gòu)造容易求解的問題作為同倫初始問題,而同倫過程中具有全局收斂性,因此同倫方法可以有效增大打靶過程收斂域,減小求解難度。同倫函數(shù)的一般構(gòu)造形式可以表示為式(16):
H(y,κ)=κF(y)+(1-κ)G(y)
(16)
式中:y∈n為優(yōu)化變量,κ為同倫參數(shù),F(xiàn):n→n和G:n→n均為光滑的映射函數(shù),分別表示所需求解的原問題和構(gòu)造的同倫初始問題,H:n×→n為構(gòu)造得到的同倫函數(shù)。根據(jù)初始問題G的不同,一般可以將同倫方法劃分為牛頓同倫(Newton Homotopy),定點同倫(Fixed-point Homotopy)和尺度不變仿射同倫(Scale-invariant Affine Homotopy)[14]。針對本文的時間最優(yōu)小推力轉(zhuǎn)移軌道優(yōu)化問題,在推力幅值中引入同倫參數(shù)κ,構(gòu)造的同倫函數(shù)的動力學模型為式(17)~(18):
(17)
(18)
式中:TL?Tmax為足夠大的推力值,在該推力值下,轉(zhuǎn)移軌道優(yōu)化問題較容易直接求解,Tmax為所需求解問題的實際推力幅值。對應(yīng)于式(16)中的一般同倫形式,F(xiàn)為所需求解的推力值為Tmax的原問題,G為容易求解的推力值為TL的初始同倫問題。在時間最優(yōu)的性能指標下,得到的哈密爾頓函數(shù)為式(19):
(19)
協(xié)態(tài)變量的微分方程保持式(9)~(10)的形式不變,同時最優(yōu)推力大小u*和最優(yōu)推力方向α*和式(11)~(12)保持一致。邊界約束條件和橫截條件也與原問題保持不變。至此,同倫問題轉(zhuǎn)換成一系列隨著κ值變化的兩點邊值問題,通過求解這一系列的同倫子問題,最終可得到原問題的解。
同倫方法除了構(gòu)造合理的同倫函數(shù),還需要選取合適的同倫曲線跟蹤方法。同倫曲線的跟蹤方法主要可以分為離散同倫和連續(xù)同倫。離散同倫可以描述為,將同倫參數(shù)從0到1劃分為m個離散節(jié)點,即0=κ1<κ2<…<κm-1<κm=1,先求解κ1對應(yīng)的初始問題的解,然后將得到的解y1作為κ2對應(yīng)的同倫子問題的初始猜測值,求解得到該問題的解y2,然后再將y2作為κ3對應(yīng)的子問題的初始猜測值,依次計算,最終得到κm對應(yīng)的原問題的解ym[7]。離散同倫簡單直觀,但若相鄰兩個κ節(jié)點之間距離較遠,則以前一個節(jié)點的解作為當前節(jié)點子問題的猜測值時收斂性不能保證;更重要的是,當同倫曲線存在拐點時,離散同倫將在拐點處失效。
本文利用連續(xù)同倫中的偽弧長法跟蹤同倫路徑[15]。不同于離散同倫根據(jù)同倫參數(shù)κ來選取同倫節(jié)點,偽弧長法沿著同倫曲線的切線方程,根據(jù)偽弧長的步長Δs進行迭代計算。偽弧長法是一種預測-校正方法。首先在當前節(jié)點計算表征梯度信息的雅克比矩陣,將步長Δs分配到各個優(yōu)化變量上,得到下一個節(jié)點的信息,然后再利用牛頓法等校正方法對下一個節(jié)點的值進行修正,從而得到該節(jié)點所對應(yīng)的解[15]。如圖1所示,假設(shè)在第i個節(jié)點,(κi,yi)為非線性方程組H(κ,y)=O的已知解,定義該點處的單位切向量為式(20):
(20)
圖1 偽弧長法示意圖Fig.1 Geometric interpretation of pseudo arclength continuation method
(21)
H(κ,y)=O
(22)
由于κ值的變化方向由同倫曲線的切線方法確定,在同倫過程中κ可以增加,也可以減小,因此可以用于跟蹤存在拐點的同倫路徑。
本文的仿真算例為地月系平面限制性三體問題下航天器從GEO軌道到L2的時間最優(yōu)小推力轉(zhuǎn)移軌道,航天器的初始質(zhì)量為1500 kg,發(fā)動機最大推力值為Tmax=1 N,對應(yīng)的初始推重比為6.8027×10-5。航天器初始位置由式(6)確定,終端狀態(tài)約束為航天器終端時刻狀態(tài)量與L2點重合,GEO軌道半徑為r0=0.109689855932071 DU,航天器在GEO軌道上的初始速度為v0=3.000969693845573 DU/TU。
首先選取較大的推力值為TL=10 N,此時利用打靶法,可得初始問題的解為y0=[-7.4762,-3.7986,6.9355,-0.3092,-0.0844,1.4413]T。對應(yīng)的轉(zhuǎn)移軌道如圖2所示,其中藍色實心圓點表示地球位置,紅色實心圓點表示月球位置,航天器到達L2點所需的轉(zhuǎn)移軌道圈數(shù)為1圈。
圖2 TL=10 N時從GEO到L2點的轉(zhuǎn)移軌道Fig.2 Transfer trajectory from GEO to L2 libration point with TL=10 N
以初始解y0為同倫初始點,利用偽弧長法跟蹤同倫軌跡,可以得到從κ=0到κ=1的同倫路徑,其中優(yōu)化變量中的轉(zhuǎn)移時間tf的同倫路徑如圖3所示,初始時刻的協(xié)態(tài)變量λ0=[λx0,λy0,λvx0,λvy0]T如圖4所示,表示初始時刻航天器在GEO軌道上的位置的變量θ0的同倫曲線如圖5所示。從圖中可以看出,該同倫曲線并不是隨著κ單調(diào)遞增的,而是存在多個拐點。在同倫初始階段,κ值沿著同倫路徑增加,但同倫曲線在κ=0.8588時遇到第一個拐點,改變了同倫路徑的方向,使κ值沿著同倫路徑減小,直到在κ=0.8408時遇到第二個拐點,使κ值重新沿著同倫路徑增加,如此往復,直到同倫路徑到達κ=1,最終得到Tmax=1 N的轉(zhuǎn)移軌道。
圖3 轉(zhuǎn)移時間tf的同倫曲線Fig.3 Homotopy curve of the transfer time tf(κ)
圖4 初始時刻協(xié)態(tài)變量λ0的同倫曲線Fig.4 Homotopy curves of costate λ0(κ)
圖5 表示航天器在GEO上的位置的θ0的同倫曲線Fig.5 Homotopy curve of θ0(κ)
值得注意的是,同倫路徑上的每一個點都對應(yīng)該點κ值下的一個解,即在對應(yīng)推力值下,滿足打靶方程的一條時間最優(yōu)轉(zhuǎn)移軌道。這是因為間接法的最優(yōu)控制律和橫截條件是基于一階最優(yōu)性條件推到得到的,可以保證滿足這些約束條件的解至少是極值點,即局部最優(yōu)解,但不能保證是全局最優(yōu)解。而在軌道轉(zhuǎn)移優(yōu)化問題中,特別是轉(zhuǎn)移圈數(shù)較多的問題,一般都存在多個局部最優(yōu)解。從同倫曲線中可以看出,在第二個拐點κ=0.8408之后,每個κ值對應(yīng)的推力下都存在多個解。以κ=0.90為例,其對應(yīng)的推力值為1.9 N,在這條同倫路徑中共存在13個解,對應(yīng)的轉(zhuǎn)移時間范圍為5.0480~ 6.1190 TU,通過對性能指標即轉(zhuǎn)移時間進行對比,即可在這些局部最優(yōu)解中選取得到最優(yōu)解。對于原問題中的Tmax=1 N的轉(zhuǎn)移軌道,雖然在這條同倫路徑中只存在1個局部解,但可以預見的是,如果繼續(xù)跟蹤同倫路勁,隨著同倫過程的進行,將會得到更多的局部解。根據(jù)tf的同倫路徑的變化規(guī)律,當前得到的這個解,對應(yīng)的轉(zhuǎn)移時間將始終小于其他解,即可認為當前得到的解是最優(yōu)解。
從圖5中可知,與圖2對應(yīng)的10 N時的轉(zhuǎn)移軌道相比,Tmax=1 N的轉(zhuǎn)移圈數(shù)增加了9圈。該解對應(yīng)的轉(zhuǎn)移軌道如圖6所示,所需轉(zhuǎn)移時間為8.7113 DU,即37.8777天。轉(zhuǎn)移過程中航天器發(fā)動機的推力幅值始終保持最大值1 N,推力方向α*的變化情況如圖7所示。
圖6 Tmax=1 N時的從GEO到L2點的轉(zhuǎn)移軌道Fig.6 Transfer trajectory from GEO to L2 libration point with Tmax=1 N
圖7 Tmax=1 N的轉(zhuǎn)移軌道推力方向α*隨時間的變化曲線Fig.7 Variation of thrust direction α* along the transfer trajectory with Tmax=1 N
1) 在轉(zhuǎn)移軌道優(yōu)化問題中,針對推力大小進行同倫時,同倫曲線存在拐點,離散同倫不能得到完整的同倫曲線,因此需要選取連續(xù)同倫軌跡跟蹤方法,利用偽弧長法可以跟蹤得到連續(xù)同倫路徑,得到包括拐點在內(nèi)的完整同倫信息,有利于對同倫過程及同倫結(jié)果進行分析;
2) 對于多圈的小推力轉(zhuǎn)移軌道優(yōu)化問題,可能存在多個滿足最優(yōu)性一階必要條件的局部最優(yōu)解,需要對性能指標進行對比才能確定最優(yōu)解。若直接對小推力轉(zhuǎn)移軌道進行求解,并得到其中一個解,不能從該解本身判斷其是否為全局最優(yōu)解,但如果利用連續(xù)同倫法跟蹤完整同倫路徑得到小推力轉(zhuǎn)移軌道,則可以依據(jù)該解在同倫路徑的位置及同倫路徑變化規(guī)律來判斷該解是否為全局最優(yōu)解。