基于同倫方法的地月系L2點小推力轉(zhuǎn)移軌道優(yōu)化

潘 迅，泮斌峰

(西北工業(yè)大學航天學院，西安 710072)

1 引言

在地月系三體模型中，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é)合同倫方法對該兩點邊值問題進行求解。

2 動力學模型

考慮航天器同時受到地球引力和月球引力，忽略地球非球形攝動、太陽引力、太陽光壓等其他攝動，并假設(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)化變量。

3 推力幅值的同倫方法

將時間最優(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)化問題的初始值敏感問題，本文利用同倫方法進行求解。

3.1 構(gòu)造同倫函數(shù)

基于所需求解的原問題，構(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)換成一系列隨著κ值變化的兩點邊值問題，通過求解這一系列的同倫子問題，最終可得到原問題的解。

3.2 同倫曲線跟蹤方法

同倫方法除了構(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)

由于κ值的變化方向由同倫曲線的切線方法確定，在同倫過程中κ可以增加，也可以減小，因此可以用于跟蹤存在拐點的同倫路徑。

4 數(shù)值仿真驗證

本文的仿真算例為地月系平面限制性三體問題下航天器從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

5 結(jié)論

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)解。