潘 迅，泮斌峰

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

1 引言

在地月系三體模型中，L2點(diǎn)在地月連線上，且始終位于月球背面，不僅可以作為月球探測(cè)和深空探測(cè)的中轉(zhuǎn)站，還可以在L2點(diǎn)建立月球觀測(cè)站，實(shí)時(shí)監(jiān)測(cè)月球背面情況，為月球背面著陸提供相關(guān)信息[1]。我國(guó)的嫦娥五號(hào)再入返回飛行試驗(yàn)器的服務(wù)艙曾于2014年11月28日到達(dá)地月系L2點(diǎn)，并在2015年1月4日離開L2轉(zhuǎn)移至環(huán)月軌道[2-3]。

嫦娥五號(hào)再入返回試驗(yàn)器利用化學(xué)發(fā)動(dòng)機(jī)經(jīng)過(guò)多脈沖機(jī)動(dòng)完成軌道轉(zhuǎn)移，轉(zhuǎn)移過(guò)程較為簡(jiǎn)單，所需時(shí)間較短。然而，隨著推進(jìn)技術(shù)的發(fā)展，電推進(jìn)、等離子推進(jìn)等持續(xù)小推力推進(jìn)方式由于具有比沖大、控制精度高等優(yōu)點(diǎn)，能極大地提高探測(cè)器的有效載荷的比重[4]，因此將在未來(lái)的航天任務(wù)，特別是深空探測(cè)任務(wù)中，扮演越來(lái)越重要的角色。

由于小推力發(fā)動(dòng)機(jī)的推力幅值小，持續(xù)時(shí)間長(zhǎng)，完成轉(zhuǎn)移所需的軌道圈數(shù)多，在轉(zhuǎn)移過(guò)程中初始狀態(tài)量對(duì)轉(zhuǎn)移軌道有較大影響，即對(duì)初始值很敏感，因此優(yōu)化過(guò)程存在較大困難。此外，相比于傳統(tǒng)的二體模型，限制性三體問(wèn)題的非線性更強(qiáng)，導(dǎo)致該小推力轉(zhuǎn)移軌道優(yōu)化問(wèn)題的求解難度增加。小推力轉(zhuǎn)移軌道優(yōu)化方法主要可以分為直接法、間接法和混合法[4]。直接法通過(guò)將狀態(tài)量和控制量進(jìn)行離散，將軌道優(yōu)化問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題，直接對(duì)離散參數(shù)進(jìn)行尋解[5-6]；間接法基于最優(yōu)控制理論的經(jīng)典變分法推導(dǎo)得到最優(yōu)控制的一階必要條件，將其轉(zhuǎn)換為兩點(diǎn)邊值問(wèn)題，通過(guò)數(shù)值方法對(duì)該兩點(diǎn)邊值問(wèn)題進(jìn)行求解，從而得到轉(zhuǎn)移軌道[7-8]；混合法只利用間接法推導(dǎo)得到的最優(yōu)控制律，但不包含其他最優(yōu)性條件，通過(guò)參數(shù)優(yōu)化方法進(jìn)行求解[9]。其中，間接法由于利用了最優(yōu)性條件，可以保證所得解的局部最優(yōu)性，且計(jì)算精度高。然而間接法將軌道優(yōu)化問(wèn)題轉(zhuǎn)換成兩點(diǎn)邊值問(wèn)題后，雖然理論上可以利用打靶法對(duì)該問(wèn)題直接求解，但由于作為優(yōu)化變量的協(xié)態(tài)變量沒(méi)有物理意義，且算法的收斂域較小，因此，對(duì)于推力值小、轉(zhuǎn)移圈數(shù)多的問(wèn)題，求解過(guò)程存在很大困難。

為了增加優(yōu)化算法的收斂域，降低求解難度，同倫方法被廣泛應(yīng)用于軌道優(yōu)化。同倫方法，是指從較簡(jiǎn)單問(wèn)題出發(fā)，通過(guò)改變同倫參數(shù)，求解對(duì)應(yīng)的一系列子問(wèn)題，最終得到原問(wèn)題的解。其在小推力在轉(zhuǎn)移軌道中的應(yīng)用主要有兩個(gè)方向：一是針對(duì)燃料最優(yōu)問(wèn)題，通過(guò)對(duì)性能指標(biāo)進(jìn)行同倫，先求解連續(xù)光滑控制問(wèn)題，然后通過(guò)同倫迭代，最終得到不連續(xù)的Bang-Bang控制的燃料最優(yōu)轉(zhuǎn)移軌道[10-11]；二是針對(duì)推力很小的軌道優(yōu)化問(wèn)題，先求解較簡(jiǎn)單的大推力值對(duì)應(yīng)的軌道優(yōu)化問(wèn)題，然后通過(guò)對(duì)推力幅值進(jìn)行同倫，逐漸減小推力幅值，最終得到所需的很小推力幅值的轉(zhuǎn)移軌道[12-13]。

目前對(duì)限制性三體問(wèn)題中的小推力轉(zhuǎn)移軌道優(yōu)化問(wèn)題研究較少，Caillau等[12]對(duì)從航天器從GEO軌道到地月系L1點(diǎn)和月球軌道的轉(zhuǎn)移軌道進(jìn)行相關(guān)研究，本文將其擴(kuò)展到了航天器從地球GEO軌道到地月系L2點(diǎn)的時(shí)間最優(yōu)轉(zhuǎn)移軌道優(yōu)化問(wèn)題。首先基于最優(yōu)控制理論，將地月系三體問(wèn)題模型下的小推力轉(zhuǎn)移軌道問(wèn)題轉(zhuǎn)換為兩點(diǎn)邊值問(wèn)題，然后結(jié)合同倫方法對(duì)該兩點(diǎn)邊值問(wèn)題進(jìn)行求解。

2 動(dòng)力學(xué)模型

考慮航天器同時(shí)受到地球引力和月球引力，忽略地球非球形攝動(dòng)、太陽(yáng)引力、太陽(yáng)光壓等其他攝動(dòng)，并假設(shè)月球圍繞地球的運(yùn)行軌道為圓軌道，則可用圓限制性三體問(wèn)題來(lái)描述航天器的運(yùn)動(dòng)。僅考慮航天器在月球軌道面內(nèi)的平面運(yùn)動(dòng)，并且忽略轉(zhuǎn)移過(guò)程中航天器的質(zhì)量變化，則航天器在旋轉(zhuǎn)坐標(biāo)系下的動(dòng)力學(xué)模型為式(1)～(2)：

(1)

(2)

式中：r=[x,y]T和v=[vx,vy]T分別為航天器的位置矢量和速度矢量，m0為航天器初始質(zhì)量，Tmax為航天器可以提供的最大推力，u為轉(zhuǎn)移過(guò)程中的實(shí)際推力大小，α=[αx,αy]T為單位推力方向。則可定義航天器有效勢(shì)能函數(shù)U如式(3)，g(r)和h(v)分別表示為式(4)、(5)：

(3)

(4)

h(v)=[2vy,-2vx]T

(5)

需要注意的是，在動(dòng)力學(xué)模型式(1)～(2)中，各變量均為單位化后的值，其中長(zhǎng)度歸一化單位為1 DU=3.8440×105km，時(shí)間歸一化單位為1 TU=3.75676967×105s，質(zhì)量歸一化單位為1 MU=M0kg(M0為航天器實(shí)際質(zhì)量)。歸一化后，地球在旋轉(zhuǎn)坐標(biāo)系的位置為(-μ,0)，月球的位置為(1-μ,0)。

本文針對(duì)航天器從地球GEO軌道到L2平動(dòng)點(diǎn)的時(shí)間最優(yōu)小推力轉(zhuǎn)移軌道進(jìn)行研究。在初始時(shí)刻t0，航天器在GEO軌道上的位置自由，其狀態(tài)量可以表示為式(6)：

x0=[r0cosθ0-μ,r0sinθ0,-v0sinθ0,v0cosθ0]T

(6)

時(shí)間最優(yōu)問(wèn)題的性能指標(biāo)可以表示為式(7)：

(7)

時(shí)間最優(yōu)小推力轉(zhuǎn)移軌道優(yōu)化問(wèn)題的本質(zhì)，就是在滿足初始狀態(tài)約束式(6)和終端狀態(tài)約束xf的條件下，通過(guò)優(yōu)化變量u和α，使性能指標(biāo)最小。根據(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)

對(duì)于航天器初始時(shí)刻的狀態(tài)量x0，由于存在優(yōu)化變量θ0，需要對(duì)橫截條件進(jìn)行推導(dǎo)。根據(jù)最優(yōu)控制理論，可以得到對(duì)應(yīng)的橫截條件為式(13)：

(x0+μ)λy0-y0λx0+vx0λvy0-vy0λvx0=0

(13)

對(duì)于時(shí)間最優(yōu)問(wèn)題，轉(zhuǎn)移時(shí)間自由，對(duì)應(yīng)的橫截條件為式(14)：

H(tf)=0

(14)

至此，將最優(yōu)控制問(wèn)題轉(zhuǎn)換為兩點(diǎn)邊值問(wèn)題，其初始時(shí)刻的狀態(tài)值由式(6)確定，終端時(shí)刻的約束包括狀態(tài)約束xf和橫截條件式(13)～(14)。該兩點(diǎn)邊值問(wèn)題可利用打靶法進(jìn)行求解，打靶函數(shù)表示為式(15)：

F(y)=O

(15)

式中y=[θ0,λx0,λy0,λvx0,λvy0,tf]T為優(yōu)化變量。

3 推力幅值的同倫方法

將時(shí)間最優(yōu)小推力轉(zhuǎn)移軌道優(yōu)化問(wèn)題轉(zhuǎn)換為兩點(diǎn)邊值問(wèn)題后，雖然理論上可以利用打靶法對(duì)式(15)進(jìn)行求解得到優(yōu)化變量y的值，再對(duì)初值進(jìn)行積分得到轉(zhuǎn)移軌道。然而推力越小，轉(zhuǎn)移時(shí)間越長(zhǎng)，轉(zhuǎn)移圈數(shù)越多，打靶過(guò)程對(duì)初始猜測(cè)值就越敏感，打靶的收斂域就越小。此外，相比于傳統(tǒng)的只考慮中心引力場(chǎng)的二體模型，限制性三體模型的非線性更強(qiáng)，從而增加了求解難度。因此，為了克服小推力轉(zhuǎn)移軌道優(yōu)化問(wèn)題的初始值敏感問(wèn)題，本文利用同倫方法進(jìn)行求解。

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

基于所需求解的原問(wèn)題，構(gòu)造容易求解的問(wèn)題作為同倫初始問(wèn)題，而同倫過(guò)程中具有全局收斂性，因此同倫方法可以有效增大打靶過(guò)程收斂域，減小求解難度。同倫函數(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ù)，分別表示所需求解的原問(wèn)題和構(gòu)造的同倫初始問(wèn)題，H:n×→n為構(gòu)造得到的同倫函數(shù)。根據(jù)初始問(wèn)題G的不同，一般可以將同倫方法劃分為牛頓同倫(Newton Homotopy)，定點(diǎn)同倫(Fixed-point Homotopy)和尺度不變仿射同倫(Scale-invariant Affine Homotopy)[14]。針對(duì)本文的時(shí)間最優(yōu)小推力轉(zhuǎn)移軌道優(yōu)化問(wèn)題，在推力幅值中引入同倫參數(shù)κ，構(gòu)造的同倫函數(shù)的動(dòng)力學(xué)模型為式(17)～(18)：

(17)

(18)

式中：TL?Tmax為足夠大的推力值，在該推力值下，轉(zhuǎn)移軌道優(yōu)化問(wèn)題較容易直接求解，Tmax為所需求解問(wèn)題的實(shí)際推力幅值。對(duì)應(yīng)于式(16)中的一般同倫形式，F(xiàn)為所需求解的推力值為Tmax的原問(wèn)題，G為容易求解的推力值為TL的初始同倫問(wèn)題。在時(shí)間最優(yōu)的性能指標(biāo)下，得到的哈密爾頓函數(shù)為式(19)：

(19)

協(xié)態(tài)變量的微分方程保持式(9)～(10)的形式不變，同時(shí)最優(yōu)推力大小u*和最優(yōu)推力方向α*和式(11)～(12)保持一致。邊界約束條件和橫截條件也與原問(wèn)題保持不變。至此，同倫問(wèn)題轉(zhuǎn)換成一系列隨著κ值變化的兩點(diǎn)邊值問(wèn)題，通過(guò)求解這一系列的同倫子問(wèn)題，最終可得到原問(wèn)題的解。

3.2 同倫曲線跟蹤方法

同倫方法除了構(gòu)造合理的同倫函數(shù)，還需要選取合適的同倫曲線跟蹤方法。同倫曲線的跟蹤方法主要可以分為離散同倫和連續(xù)同倫。離散同倫可以描述為，將同倫參數(shù)從0到1劃分為m個(gè)離散節(jié)點(diǎn)，即0=κ1<κ2<…<κm-1<κm=1，先求解κ1對(duì)應(yīng)的初始問(wèn)題的解，然后將得到的解y1作為κ2對(duì)應(yīng)的同倫子問(wèn)題的初始猜測(cè)值，求解得到該問(wèn)題的解y2，然后再將y2作為κ3對(duì)應(yīng)的子問(wèn)題的初始猜測(cè)值，依次計(jì)算，最終得到κm對(duì)應(yīng)的原問(wèn)題的解ym[7]。離散同倫簡(jiǎn)單直觀，但若相鄰兩個(gè)κ節(jié)點(diǎn)之間距離較遠(yuǎn)，則以前一個(gè)節(jié)點(diǎn)的解作為當(dāng)前節(jié)點(diǎn)子問(wèn)題的猜測(cè)值時(shí)收斂性不能保證；更重要的是，當(dāng)同倫曲線存在拐點(diǎn)時(shí)，離散同倫將在拐點(diǎn)處失效。

本文利用連續(xù)同倫中的偽弧長(zhǎng)法跟蹤同倫路徑[15]。不同于離散同倫根據(jù)同倫參數(shù)κ來(lái)選取同倫節(jié)點(diǎn)，偽弧長(zhǎng)法沿著同倫曲線的切線方程，根據(jù)偽弧長(zhǎng)的步長(zhǎng)Δs進(jìn)行迭代計(jì)算。偽弧長(zhǎng)法是一種預(yù)測(cè)-校正方法。首先在當(dāng)前節(jié)點(diǎn)計(jì)算表征梯度信息的雅克比矩陣，將步長(zhǎng)Δs分配到各個(gè)優(yōu)化變量上，得到下一個(gè)節(jié)點(diǎn)的信息，然后再利用牛頓法等校正方法對(duì)下一個(gè)節(jié)點(diǎn)的值進(jìn)行修正，從而得到該節(jié)點(diǎn)所對(duì)應(yīng)的解[15]。如圖1所示，假設(shè)在第i個(gè)節(jié)點(diǎn)，(κi,yi)為非線性方程組H(κ,y)=O的已知解，定義該點(diǎn)處的單位切向量為式(20)：

(20)

圖1 偽弧長(zhǎng)法示意圖Fig.1 Geometric interpretation of pseudo arclength continuation method

(21)

H(κ,y)=O

(22)

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

4 數(shù)值仿真驗(yàn)證

本文的仿真算例為地月系平面限制性三體問(wèn)題下航天器從GEO軌道到L2的時(shí)間最優(yōu)小推力轉(zhuǎn)移軌道，航天器的初始質(zhì)量為1500 kg，發(fā)動(dòng)機(jī)最大推力值為Tmax=1 N，對(duì)應(yīng)的初始推重比為6.8027×10-5。航天器初始位置由式(6)確定，終端狀態(tài)約束為航天器終端時(shí)刻狀態(tài)量與L2點(diǎn)重合，GEO軌道半徑為r0=0.109689855932071 DU，航天器在GEO軌道上的初始速度為v0=3.000969693845573 DU/TU。

首先選取較大的推力值為TL=10 N，此時(shí)利用打靶法，可得初始問(wèn)題的解為y0=[-7.4762,-3.7986,6.9355,-0.3092,-0.0844,1.4413]T。對(duì)應(yīng)的轉(zhuǎn)移軌道如圖2所示，其中藍(lán)色實(shí)心圓點(diǎn)表示地球位置，紅色實(shí)心圓點(diǎn)表示月球位置，航天器到達(dá)L2點(diǎn)所需的轉(zhuǎn)移軌道圈數(shù)為1圈。

圖2 TL=10 N時(shí)從GEO到L2點(diǎn)的轉(zhuǎn)移軌道Fig.2 Transfer trajectory from GEO to L2 libration point with TL=10 N

以初始解y0為同倫初始點(diǎn)，利用偽弧長(zhǎng)法跟蹤同倫軌跡，可以得到從κ=0到κ=1的同倫路徑，其中優(yōu)化變量中的轉(zhuǎn)移時(shí)間tf的同倫路徑如圖3所示，初始時(shí)刻的協(xié)態(tài)變量λ0=[λx0,λy0,λvx0,λvy0]T如圖4所示，表示初始時(shí)刻航天器在GEO軌道上的位置的變量θ0的同倫曲線如圖5所示。從圖中可以看出，該同倫曲線并不是隨著κ單調(diào)遞增的，而是存在多個(gè)拐點(diǎn)。在同倫初始階段，κ值沿著同倫路徑增加，但同倫曲線在κ=0.8588時(shí)遇到第一個(gè)拐點(diǎn)，改變了同倫路徑的方向，使κ值沿著同倫路徑減小，直到在κ=0.8408時(shí)遇到第二個(gè)拐點(diǎn)，使κ值重新沿著同倫路徑增加，如此往復(fù)，直到同倫路徑到達(dá)κ=1，最終得到Tmax=1 N的轉(zhuǎn)移軌道。

圖3 轉(zhuǎn)移時(shí)間tf的同倫曲線Fig.3 Homotopy curve of the transfer time tf(κ)

圖4 初始時(shí)刻協(xié)態(tài)變量λ0的同倫曲線Fig.4 Homotopy curves of costate λ0(κ)

圖5 表示航天器在GEO上的位置的θ0的同倫曲線Fig.5 Homotopy curve of θ0(κ)

值得注意的是，同倫路徑上的每一個(gè)點(diǎn)都對(duì)應(yīng)該點(diǎn)κ值下的一個(gè)解，即在對(duì)應(yīng)推力值下，滿足打靶方程的一條時(shí)間最優(yōu)轉(zhuǎn)移軌道。這是因?yàn)殚g接法的最優(yōu)控制律和橫截條件是基于一階最優(yōu)性條件推到得到的，可以保證滿足這些約束條件的解至少是極值點(diǎn)，即局部最優(yōu)解，但不能保證是全局最優(yōu)解。而在軌道轉(zhuǎn)移優(yōu)化問(wèn)題中，特別是轉(zhuǎn)移圈數(shù)較多的問(wèn)題，一般都存在多個(gè)局部最優(yōu)解。從同倫曲線中可以看出，在第二個(gè)拐點(diǎn)κ=0.8408之后，每個(gè)κ值對(duì)應(yīng)的推力下都存在多個(gè)解。以κ=0.90為例，其對(duì)應(yīng)的推力值為1.9 N，在這條同倫路徑中共存在13個(gè)解，對(duì)應(yīng)的轉(zhuǎn)移時(shí)間范圍為5.0480～ 6.1190 TU，通過(guò)對(duì)性能指標(biāo)即轉(zhuǎn)移時(shí)間進(jìn)行對(duì)比，即可在這些局部最優(yōu)解中選取得到最優(yōu)解。對(duì)于原問(wèn)題中的Tmax=1 N的轉(zhuǎn)移軌道，雖然在這條同倫路徑中只存在1個(gè)局部解，但可以預(yù)見的是，如果繼續(xù)跟蹤同倫路勁，隨著同倫過(guò)程的進(jìn)行，將會(huì)得到更多的局部解。根據(jù)tf的同倫路徑的變化規(guī)律，當(dāng)前得到的這個(gè)解，對(duì)應(yīng)的轉(zhuǎn)移時(shí)間將始終小于其他解，即可認(rèn)為當(dāng)前得到的解是最優(yōu)解。

從圖5中可知，與圖2對(duì)應(yīng)的10 N時(shí)的轉(zhuǎn)移軌道相比，Tmax=1 N的轉(zhuǎn)移圈數(shù)增加了9圈。該解對(duì)應(yīng)的轉(zhuǎn)移軌道如圖6所示，所需轉(zhuǎn)移時(shí)間為8.7113 DU，即37.8777天。轉(zhuǎn)移過(guò)程中航天器發(fā)動(dòng)機(jī)的推力幅值始終保持最大值1 N，推力方向α*的變化情況如圖7所示。

圖6 Tmax=1 N時(shí)的從GEO到L2點(diǎn)的轉(zhuǎn)移軌道Fig.6 Transfer trajectory from GEO to L2 libration point with Tmax=1 N

圖7 Tmax=1 N的轉(zhuǎn)移軌道推力方向α*隨時(shí)間的變化曲線Fig.7 Variation of thrust direction α* along the transfer trajectory with Tmax=1 N

5 結(jié)論

1) 在轉(zhuǎn)移軌道優(yōu)化問(wèn)題中，針對(duì)推力大小進(jìn)行同倫時(shí)，同倫曲線存在拐點(diǎn)，離散同倫不能得到完整的同倫曲線，因此需要選取連續(xù)同倫軌跡跟蹤方法，利用偽弧長(zhǎng)法可以跟蹤得到連續(xù)同倫路徑，得到包括拐點(diǎn)在內(nèi)的完整同倫信息，有利于對(duì)同倫過(guò)程及同倫結(jié)果進(jìn)行分析；

2) 對(duì)于多圈的小推力轉(zhuǎn)移軌道優(yōu)化問(wèn)題，可能存在多個(gè)滿足最優(yōu)性一階必要條件的局部最優(yōu)解，需要對(duì)性能指標(biāo)進(jìn)行對(duì)比才能確定最優(yōu)解。若直接對(duì)小推力轉(zhuǎn)移軌道進(jìn)行求解，并得到其中一個(gè)解，不能從該解本身判斷其是否為全局最優(yōu)解，但如果利用連續(xù)同倫法跟蹤完整同倫路徑得到小推力轉(zhuǎn)移軌道，則可以依據(jù)該解在同倫路徑的位置及同倫路徑變化規(guī)律來(lái)判斷該解是否為全局最優(yōu)解。