周敬，胡軍，*

1. 北京控制工程研究所，北京 100190 2. 空間智能控制技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100094

近年來，隨著航天基礎(chǔ)理論的不斷深入研究和工程實(shí)踐技術(shù)的快速發(fā)展，深空探測越來越受到世界各航天大國與組織的重視。根據(jù)中國航天事業(yè)發(fā)展的指南，未來航天的三個重點(diǎn)方向之一就是深空探測[1]。相比以經(jīng)典攝動二體運(yùn)動問題為基礎(chǔ)的近地空間航天活動，深空探測所涉及的運(yùn)動模型更多的是三體問題。動力學(xué)本質(zhì)上的區(qū)別使得三體問題相對二體問題更加復(fù)雜。在三體問題的研究中，平動點(diǎn)(尤其是共線平動點(diǎn))及其附近的周期/擬周期軌道所具有的獨(dú)特位置優(yōu)勢和豐富的動力學(xué)特性[2]，使其在深空探測任務(wù)設(shè)計(jì)中占據(jù)著重要的地位，已經(jīng)成為宇宙觀測、天文研究的理想場所和星際高速公路(interplanetary superhighway，IPS)的門戶站。

通常情況下，在平動點(diǎn)軌道進(jìn)入工程應(yīng)用之前，首先需要解決的便是與其有關(guān)的軌道轉(zhuǎn)移問題。按照推進(jìn)方式的不同，平動點(diǎn)軌道之間的轉(zhuǎn)移可以分為脈沖轉(zhuǎn)移和小推力轉(zhuǎn)移兩種方式，且小推力轉(zhuǎn)移已逐漸成為深空探測領(lǐng)域的研究熱點(diǎn)。

在小推力軌道轉(zhuǎn)移方面，小推力推進(jìn)系統(tǒng)具有比沖高、推力小、工作時間長、推進(jìn)劑消耗量少等特點(diǎn)，非常契合深空探測的任務(wù)特點(diǎn)，在深空探測研究中具有廣闊的應(yīng)用前景。然而，小推力的引入不可避免地使原本就復(fù)雜的三體問題的復(fù)雜性進(jìn)一步增加，導(dǎo)致相關(guān)的軌道設(shè)計(jì)與優(yōu)化面臨新的問題與挑戰(zhàn)。目前，小推力軌道優(yōu)化方法主要分為直接法、間接法和混合法[3-4]。直接法雖然收斂性較好，但求解精度較低，計(jì)算量大，對于小推力軌道優(yōu)化這類非線性、多局部極值問題，解的最優(yōu)性一般難以保證。間接法雖然可以保證解的最優(yōu)性，但協(xié)態(tài)變量無明確物理意義，難以給出合理的初始猜測值?；旌戏骖欓g接法和直接法的優(yōu)點(diǎn)，精度更高、收斂域更寬、數(shù)值穩(wěn)定性更好，但同時復(fù)雜性也更大。

在小推力研究方面，文獻(xiàn)[5]借助平動點(diǎn)軌道的相空間結(jié)構(gòu)揭示了小推力轉(zhuǎn)移的機(jī)理，對小推力方式的地月低能轉(zhuǎn)移問題進(jìn)行了研究。文獻(xiàn)[6]采用遺傳/逐次二次規(guī)劃混合優(yōu)化算法，研究了近地小推力轉(zhuǎn)移軌道的制導(dǎo)問題。文獻(xiàn)[7]基于軌道逆推思想，采用人工免疫算法對地月轉(zhuǎn)移軌道設(shè)計(jì)中的小推力捕獲軌道進(jìn)行求解。文獻(xiàn)[8]提出了基于N次逆多項(xiàng)式逼近的半解析Lambert算法，并基于該算法發(fā)展了一種轉(zhuǎn)移軌道初始設(shè)計(jì)方法，實(shí)現(xiàn)了行星際小推力軌道設(shè)計(jì)。文獻(xiàn)[9]結(jié)合小推力技術(shù)和金星借力技術(shù)，研究了火星轉(zhuǎn)移軌道設(shè)計(jì)問題。文獻(xiàn)[10]結(jié)合小推力技術(shù)和bang-bang控制技術(shù)，研究了高精度的皮納衛(wèi)星編隊(duì)構(gòu)型保持問題。

Gauss偽譜法是近年來新興的一種高效的最優(yōu)控制問題求解方法，本質(zhì)上屬于直接法，其具有三個顯著優(yōu)點(diǎn)：1)動力學(xué)約束只與當(dāng)前節(jié)點(diǎn)狀態(tài)有關(guān)，尋優(yōu)參數(shù)規(guī)模?。?)KKT(Karush-Kuhn-Tucher)條件與極大值原理中的一階最優(yōu)性條件等價，求解精度高；3)具有偽譜法共有的指數(shù)收斂特性，收斂性好。目前有關(guān)Gauss法的研究主要有：楊博等利用Gauss偽譜法研究了地球同步軌道衛(wèi)星軌道轉(zhuǎn)移問題[11]；曹喜濱等基于Gauss偽譜法將最優(yōu)控制問題離散化為NLP問題，并采用基于逆多項(xiàng)式的形狀法給出了NLP初值的計(jì)算方法[12]；尚海濱等基于Gauss偽譜法的配點(diǎn)特性，推導(dǎo)出了性能指標(biāo)和約束方程的解析雅可比矩陣，進(jìn)而解決了星際小推力轉(zhuǎn)移軌道優(yōu)化問題[13]。

由于Gauss偽譜法得到的NLP模型包含較多的等式約束，且小推力軌道的非線性較強(qiáng)，控制變量初值設(shè)置的不合理往往會導(dǎo)致NLP的求解迭代時間長、容易陷入局部極小甚至不收斂等問題，因此對Gauss偽譜法的控制變量初值確定方法進(jìn)行研究具有重要意義。目前，對初值確定方法研究最多的是形狀法，主要有：文獻(xiàn)[14-16]采用逆多項(xiàng)式函數(shù)作為小推力轉(zhuǎn)移軌道的形狀函數(shù)，對小推力軌道轉(zhuǎn)移問題進(jìn)行研究；Ellutini等利用正弦指數(shù)函數(shù)作為形狀函數(shù)，來近似小推力轉(zhuǎn)移軌道[17]；Novak等提出了一種基于偽春分軌道根數(shù)的三維形狀方法[18]，解決了地球-火星轉(zhuǎn)移軌道設(shè)計(jì)；Gondelach等提出了速度速率圖的概念，將小推力軌道轉(zhuǎn)移的速度矢量表示為時間或者極角的形狀函數(shù)，然后對此函數(shù)進(jìn)行優(yōu)化，實(shí)現(xiàn)了地球-火星、小行星、彗星和水星之間的轉(zhuǎn)移[19]；文獻(xiàn)[20-22]利用Fourier級數(shù)作為小推力轉(zhuǎn)移軌道的形狀函數(shù)，對小推力軌道轉(zhuǎn)移問題進(jìn)行了研究。

上述研究均獲得了不錯的研究成果，但基本都是針對傳統(tǒng)常見的軌道轉(zhuǎn)移問題。對于平動點(diǎn)周期軌道間的小推力轉(zhuǎn)移這一新興問題，國內(nèi)外研究很少涉及，且三體問題相對二體問題的復(fù)雜性可能使得上述研究方法無法適用于平動點(diǎn)軌道之間的小推力轉(zhuǎn)移問題。鑒于小推力技術(shù)在深空探測中具有的廣泛應(yīng)用前景，本文在充分研究三體問題特性的基礎(chǔ)上，構(gòu)造了一種新的形狀函數(shù)，并在此基礎(chǔ)上基于Gauss偽譜法對共線平動點(diǎn)周期軌道間的小推力轉(zhuǎn)移進(jìn)行了研究，以期為未來中國深空探測事業(yè)以及平動點(diǎn)軌道工程應(yīng)用奠定一定的理論基礎(chǔ)。

1 軌道轉(zhuǎn)移優(yōu)化模型

1.1 動力學(xué)模型

在三體問題中，圓型限制性三體問題(circular restricted three-body problem, CRTBP)是最簡單的三體運(yùn)動模型，常用來研究最基本的運(yùn)動規(guī)律，同時考慮到地月L1點(diǎn)作為地球通向宇宙空間的門戶站以及在此部署空間站的可能性，具有重要的工程應(yīng)用價值。因此本文以地月系統(tǒng)CRTBP下的L1點(diǎn)附近的周期軌道作為本文的研究對象。

為方便描述航天器運(yùn)動和簡化計(jì)算，在CRTBP中通常需要對相關(guān)物理量進(jìn)行無量綱化處理，并以時間作為獨(dú)立變量。相應(yīng)的質(zhì)量[M]、長度[L]和時間[T]的歸一化單位取為：

(1)

式中:m1和m2為兩主天體質(zhì)量；L12為兩主天體之間的距離；G為萬有引力常數(shù)，其他相關(guān)變量的歸一化過程均可依據(jù)式(1)推導(dǎo)獲得。

為了更直觀、更清晰地描述航天器在CRTBP中的運(yùn)動，通常選擇在會合坐標(biāo)系(也稱為質(zhì)心旋轉(zhuǎn)坐標(biāo)系)下進(jìn)行研究。如圖1所示，會合坐標(biāo)系O-XYZ的原點(diǎn)位于兩主天體的公共質(zhì)心，x軸由較大主天體指向較小主天體，z軸與主天體系統(tǒng)角動量方向平行，y軸滿足右手坐標(biāo)系定則。

圖1 地月系統(tǒng)CRTBP下的會合坐標(biāo)系O-XYZ和 L1會合坐標(biāo)系L1-xyzFig.1 Synodic coordinate system O-XYZand L1-centered synodic coordinate system L1-xyz for CRTBP of Earth-Moon system

式中：F=[Fx,Fy,Fz]T為歸一化后的推力加速度矢量及三軸分量；Ω為偽勢能函數(shù)。

式中：X1為會合坐標(biāo)系下L1點(diǎn)距離主天體公共質(zhì)心的距離；γ1為L1點(diǎn)與最近主天體之間的距離。

因此，L1會合坐標(biāo)系下航天器運(yùn)動的動力學(xué)方程為：

(2)

(3)

此外，考慮到Halo軌道是共線平動點(diǎn)附近一類特殊的三維周期軌道，且常作為深空探測航天器的工作軌道，具有重要的學(xué)術(shù)研究和工程應(yīng)用價值，Richardson基于Lindstedt-Poincaré方法推導(dǎo)了Halo軌道的三階近似解析解，如下式所示，相關(guān)參數(shù)的含義及計(jì)算過程可以參考文獻(xiàn)[23]。本文以Halo 軌道間的轉(zhuǎn)移為例對小推力軌道轉(zhuǎn)移進(jìn)行研究。

(4)

1.2 優(yōu)化性能指標(biāo)和約束

1)參考實(shí)際工程應(yīng)用要求，優(yōu)化性能指標(biāo)按照燃料最優(yōu)原則，即性能指標(biāo)函數(shù)J取為:

式中：t0為軌道轉(zhuǎn)移初始時刻，通常設(shè)置為零；tf為軌道轉(zhuǎn)移結(jié)束時刻，通常根據(jù)任務(wù)需求給出取值范圍；mfuel為軌道轉(zhuǎn)移過程中的燃料消耗量；m0為航天器總質(zhì)量；Tall=‖u‖2為飛行器總的推力加速度大??；g0為海平面重力加速度；Isp為發(fā)動機(jī)比沖。

2)考慮到發(fā)動機(jī)的最大推力限制，需要對推力加速度大小進(jìn)行限制，歸一化后的約束如下:

‖u‖2≤2DU

式中：DU為無量綱單位(dimensionless unit, DU)。

2 基于Gauss偽譜法的最優(yōu)軌道 轉(zhuǎn)移

2.1 軌道轉(zhuǎn)移近似解析解

對于Gauss偽譜法下的小推力軌道轉(zhuǎn)移優(yōu)化設(shè)計(jì)，控制變量的初值選取對軌道優(yōu)化設(shè)計(jì)的迭代過程具有非常重要的影響。為獲得較快的收斂速度，合理的初始猜測值十分必要。本文在充分結(jié)合平動點(diǎn)附近軌道運(yùn)動特性的基礎(chǔ)上，構(gòu)造了一種新的、專門適用于平動點(diǎn)周期軌道間轉(zhuǎn)移的形狀函數(shù)，下面進(jìn)行詳細(xì)介紹。

假設(shè)在L1會合坐標(biāo)系下，小推力軌道轉(zhuǎn)移的起點(diǎn)狀態(tài)為x0，終點(diǎn)狀態(tài)為xf，即

根據(jù)文獻(xiàn)[23]，共線平動點(diǎn)附近運(yùn)動的近似解析解為：

根據(jù)上式，要形成周期軌道，必然使得指數(shù)運(yùn)動項(xiàng)A1eλ1t和A2e-λ1t為零，僅保留周期運(yùn)動項(xiàng)，即A1=0,A2=0,Ax≠0,Az≠0，由此便得到了共線平動點(diǎn)附近周期軌道的一階近似解析解，如式(5)所示，這也是本文提出的新的形狀函數(shù)的原始形式。

(5)

進(jìn)一步，考慮到通常情況下小推力轉(zhuǎn)移軌道的形狀為螺旋狀，如圖2所示，對于這種軌道轉(zhuǎn)移規(guī)律，Xie等和Ellutini等分別提出了經(jīng)典的正弦指數(shù)函數(shù)[16]和逆多項(xiàng)式函數(shù)來近似小推力轉(zhuǎn)移軌道[17]。對于CRTBP，參考共線平動點(diǎn)附近周期軌道的動力學(xué)特性，即式(5)，同時考慮小推力轉(zhuǎn)移軌道的螺旋特性，在周期軌道近似解析解(5)的基礎(chǔ)上，對振幅和相位進(jìn)行時變處理，類似于軌道研究中常用的常數(shù)變易法，以此來構(gòu)造一類新的形狀函數(shù)來擬合小推力轉(zhuǎn)移軌道。借鑒逆多項(xiàng)式形狀函數(shù)思想，本文假設(shè)小推力轉(zhuǎn)移軌道近似解析解的振幅和相位均按多項(xiàng)式變化，當(dāng)然也可以假設(shè)其他形式的變化規(guī)律，如指數(shù)變化等，只要使得小推力轉(zhuǎn)移軌道近似解析解滿足螺旋特性即可。

圖2 小推力轉(zhuǎn)移軌道示意Fig. 2 Sketch map of low-thrust orbit transfer

然后考慮軌道轉(zhuǎn)移中的過程約束，如推力加速度大小約束等，便可獲得其他約束方程，也可以通過試湊法配合事后驗(yàn)證的方法對多項(xiàng)式系數(shù)進(jìn)行確定，最后獲得滿足約束條件的小推力轉(zhuǎn)移軌道的近似解析解。

(6)

速度矢量表達(dá)式為：

(7)

加速度矢量表達(dá)式為：

(8)

(9)

將初始軌道和目標(biāo)軌道的x軸振幅Ax1,Ax2和z軸振幅Az1,Az2，以及轉(zhuǎn)移起點(diǎn)和終點(diǎn)在xy平面運(yùn)動的相位φ1,φ2和z平面運(yùn)動的相位Ψ1,Ψ2代入式(9)，便可以獲得小推力轉(zhuǎn)移軌道的近似解析解。這便是本文提出的一種新的、能夠深度結(jié)合平動點(diǎn)周期軌道運(yùn)動特性的形狀函數(shù)法。

此外，由于式(5)作為平動點(diǎn)附近周期軌道的通用一階近似解析解，不僅可以表示Halo軌道，還可以表示Lissajous軌道以及Lyapunov軌道等。因此，本文提出的形狀法不僅適用于Halo軌道之間的轉(zhuǎn)移，也適用于Lissajous軌道和Lyapunov軌道以及上述三種軌道相互間的轉(zhuǎn)移。

2.2 Gauss偽譜法的離散化方法

Gauss偽譜法的主要思想是采用Lagrange插值和Gauss求積公式將最優(yōu)控制問題在一系列Legendre-Gauss(LG)點(diǎn)進(jìn)行離散化處理。參照文獻(xiàn)[12]，由于LG點(diǎn)的取值區(qū)間為(-1,1)，因此需要對最優(yōu)控制問題的時間區(qū)間進(jìn)行線性變換，新的時間變量記為τ，變換方法如下：

進(jìn)行上述變換后，本文所研究的最優(yōu)控制問題將轉(zhuǎn)化為如下形式，性能指標(biāo)函數(shù)J為：

滿足約束條件：

(1)約束方程離散化

Gauss偽譜法對于狀態(tài)變量和控制變量的逼近通過全局Lagrange插值實(shí)現(xiàn)。假設(shè)離散化的LG點(diǎn)的個數(shù)為K，則采用初始狀態(tài)變量x(τ0)和K個離散節(jié)點(diǎn)處的狀態(tài)x(τi)=X(τi)，可得到K+1階Lagrange插值多項(xiàng)式，利用該多項(xiàng)式對系統(tǒng)的狀態(tài)變量進(jìn)行近似：

(10)

式中：X(τ)為Lagrange插值多項(xiàng)式；Li(τ)為Lagrange插值基函數(shù)，其計(jì)算公式如下：

式中：g(τ)為以各插值節(jié)點(diǎn)為根的多項(xiàng)式，即

顯然式(9)中滿足x(τi)=X(τi)，對式(10)進(jìn)行求導(dǎo)，得到狀態(tài)變量在LG點(diǎn)的導(dǎo)數(shù)為：

(11)

(12)

將式(12)代入式(11)，則微分形式的狀態(tài)方程轉(zhuǎn)化為以各個離散點(diǎn)表示的代數(shù)約束方程:

(13)

通過上式共得到K個等式約束，其中X(τk),U(τk)分別為離散點(diǎn)處的狀態(tài)變量和控制變量，作為后續(xù)非線性規(guī)劃的優(yōu)化參數(shù)。

(2)末端約束離散化

在以上過程中，并未考慮末端狀態(tài)約束，因此需要額外附加約束。采用Gauss求積公式對動力學(xué)方程進(jìn)行積分：

則末端狀態(tài)可形成如下代數(shù)約束：

(14)

式中：ωi為Gauss求積系數(shù)，計(jì)算方式如下：

式中：P(τ)為n階Legendre多項(xiàng)式。

(3)路徑約束離散化

根據(jù)原始問題中控制變量的約束范圍，將控

制變量的路徑約束離散化為：

(15)

式中：k=1,…,K。

(4)目標(biāo)函數(shù)離散化

采用Gauss求積公式近似式(6)中的積分項(xiàng)，則目標(biāo)函數(shù)J可以離散化為：

(16)

至此，便實(shí)現(xiàn)了小推力轉(zhuǎn)移軌道優(yōu)化問題的離散化，即將最優(yōu)控制問題轉(zhuǎn)化為了NLP問題。根據(jù)Gauss偽譜法的求解過程，取如下非線性規(guī)劃的變量作為優(yōu)化變量：

Z={X1,…,XK,U1,…,UK,t1,tf}

(17)

式中：Xk和Uk分別為離散LG點(diǎn)處的六維狀態(tài)變量和三維控制變量。

通過以上離散化方法，構(gòu)成了以式(16)為性能指標(biāo)，以式(17)中的參數(shù)為優(yōu)化變量，以式(13)～(15)為約束方程的NLP問題。采用Gauss偽譜法的Matlab優(yōu)化工具箱GPOPS對此NLP問題進(jìn)行求解，并以形狀法獲得的控制加速度近似解作為初值，便可以獲得要求的小推力轉(zhuǎn)移軌道。

3 仿真分析

本節(jié)以地月CRTBP中L1點(diǎn)附近的Halo軌道作為研究對象，利用本文所提方法對Halo軌道之間的小推力軌道轉(zhuǎn)移進(jìn)行仿真研究。初始Halo軌道的振幅記為Ax1,Az1，目標(biāo)Halo軌道的振幅記為Ax2,Az2，轉(zhuǎn)移起點(diǎn)和終點(diǎn)的狀態(tài)矢量分別記為x0,xf，以上參數(shù)的具體數(shù)值在L1會合坐標(biāo)系下的取值如表1所示。

表1 Halo軌道間小推力軌道轉(zhuǎn)移相關(guān)參數(shù)

首先采用形狀法獲得小推力轉(zhuǎn)移軌道的近似解析解，在轉(zhuǎn)移起點(diǎn)和終點(diǎn)的狀態(tài)矢量x0,xf以及轉(zhuǎn)移時間tf的范圍確定之后，便可以根據(jù)式(6)～(9)獲得小推力轉(zhuǎn)移軌道的近似解析解，即獲得了任意時刻的位置矢量、速度矢量以及控制加速度矢量的近似解，然后按照2.1節(jié)所述過程進(jìn)行解算和處理，將獲得的控制加速度矢量作為Gauss偽譜法的控制變量的初始猜測值，利用GPOPS優(yōu)化工具箱對離散得到的NLP問題進(jìn)行求解，最終獲得最優(yōu)小推力轉(zhuǎn)移軌道。為滿足不同的任務(wù)背景需求，下面分別針對不同轉(zhuǎn)移時間的情形進(jìn)行了研究。

3.1 轉(zhuǎn)移時間較短

GPOPS無初始猜測值和有初始猜測值時的仿真結(jié)果如表2所示。由表2可以看出，在轉(zhuǎn)移時間較短的情況下，在不同控制變量初始猜測值情況下，Gauss偽譜法獲得了相同的最優(yōu)性能指標(biāo)0.7744和轉(zhuǎn)移時間1.3818，說明不同控制變量初始猜測值不會對Gauss偽譜法的最終結(jié)果產(chǎn)生影響。然而，當(dāng)形狀法為Gauss偽譜法提供控制加速度的初始猜測值時，迭代次數(shù)由3469降至1558，降幅為55.1%，仿真時間由166.3s降至34.8s，降幅為79%，說明本文提出的形狀法可以為Gauss偽譜法提供較好的初始猜測值，可以顯著加快Gauss偽譜法的迭代過程。

表2 轉(zhuǎn)移時間較短時不同初始猜測值下的 GPOPS仿真結(jié)果

圖 3 轉(zhuǎn)移時間較短時的小推力轉(zhuǎn)移軌道Fig.3 Low-thrust transfer trajectories with short transfer time

圖 4 轉(zhuǎn)移時間較短時的三軸控制加速度Fig.4 Three-axis controlled acceleration with short transfer time

形狀法和Gauss偽譜法下小推力轉(zhuǎn)移軌道的軌道如圖3所示，對應(yīng)的三軸控制加速度的時間序列如圖4所示。根據(jù)圖3，采用形狀法獲得的小推力轉(zhuǎn)移軌道與Gauss偽譜法優(yōu)化得到的最優(yōu)轉(zhuǎn)移軌道基本一致，說明了本文提出的形狀法的正確性。根據(jù)圖4，Gauss偽譜法優(yōu)化得到的最優(yōu)推力加速度矢量與基于形狀法獲得的近似推力加速度矢量大體一致，且均滿足推力加速度大小約束，進(jìn)一步說明形狀法具備為Gauss偽譜法提供合理初始猜測值的能力。

此外，形狀法獲得的加速度曲線比較光滑，而Gauss偽譜法獲得的加速度曲線則有幾處突變，說明形狀法相比Gauss偽譜法，所求解的轉(zhuǎn)移軌道具有更好的穩(wěn)定性。

3.2 轉(zhuǎn)移時間較長

在一定任務(wù)背景下，當(dāng)需要進(jìn)行較慢的軌道轉(zhuǎn)移時，此時軌道轉(zhuǎn)移時間較長，轉(zhuǎn)移軌道可以形成多圈(即圈數(shù)N≥1)，以單圈為例，此時可以設(shè)置轉(zhuǎn)移時間tf的取值范圍為[T,2T]。根據(jù)3.1節(jié)所述方法，獲得Gauss偽譜法下的最優(yōu)小推力轉(zhuǎn)移軌道。

GPOPS無初始猜測值和有初始猜測值時的仿真結(jié)果如表3所示。由表3可以看出，在轉(zhuǎn)移時間較長，且控制變量初始猜測值不同的情況下，Gauss偽譜法獲得了相同的最優(yōu)性能指標(biāo)0.8859和轉(zhuǎn)移時間4.1444，同樣說明不同控制變量初始猜測值不會對Gauss偽譜法的最終結(jié)果產(chǎn)生影響。同樣地，當(dāng)形狀法為Gauss偽譜法提供控制加速度的初始猜測值時，迭代次數(shù)由3657降至1251，降幅為65.8%，仿真時間由171.5s降至30.7s，降幅為82.1%，說明本文提出的形狀法可以為Gauss偽譜法提供較好的初始猜測值，可以顯著加快Gauss偽譜法的迭代過程，并且轉(zhuǎn)移時間越長，加速效果越顯著。

轉(zhuǎn)移時間較長時，形狀法和Gauss偽譜法下小推力轉(zhuǎn)移軌道的軌跡如圖5所示，對應(yīng)的三軸控制加速度的時間序列如圖6所示。根據(jù)圖5可知，采用形狀法獲得的小推力轉(zhuǎn)移軌道與Gauss偽譜法優(yōu)化得到的最優(yōu)轉(zhuǎn)移軌道基本一致，說明了本文所提出的形狀法的正確性。根據(jù)圖6可知，Gauss偽譜法優(yōu)化得到的最優(yōu)推力加速度矢量與基于形狀法獲得的推力加速度矢量基本一致，且均滿足推力加速度大小約束，同樣進(jìn)一步說明了形狀法具備為Gauss偽譜法提供合理的初始猜測值的能力。

同樣，形狀法獲得的加速度曲線相比Gauss偽譜法依然較為光滑，反映到轉(zhuǎn)移軌跡上，形狀法轉(zhuǎn)移軌跡變化平穩(wěn)、且一直處于目標(biāo)軌跡的邊界內(nèi)，而Gauss偽譜法轉(zhuǎn)移軌跡變化略大、且有超出目標(biāo)軌跡的邊界，進(jìn)一步說明形狀法相比Gauss偽譜法，所求解的轉(zhuǎn)移軌道具有更好的穩(wěn)定性。

表3 轉(zhuǎn)移時間較長時不同初始猜測值下的 GPOPS仿真結(jié)果

圖5 轉(zhuǎn)移時間較長時的小推力轉(zhuǎn)移軌道Fig.5 Low-thrust transfer trajectories with long transfer time

圖6 轉(zhuǎn)移時間較長時的三軸控制加速度Fig.6 Three-axis controlled acceleration with long transfer time

4 結(jié)束語

本文在深度研究平動點(diǎn)附近周期軌道特性的基礎(chǔ)上構(gòu)造了一種新的形狀函數(shù)，并在此基礎(chǔ)上提出了一種基于Gauss偽譜法的小推力軌道轉(zhuǎn)移研究方法，解決了三體問題下共線平動點(diǎn)附近周期軌道間的小推力軌道轉(zhuǎn)移問題。主要結(jié)論如下：

1)根據(jù)初始軌道和目標(biāo)軌道的軌道類型，結(jié)合小推力轉(zhuǎn)移軌道的螺旋特性和提出的振幅和相位按多項(xiàng)式變化的假設(shè)，獲得的小推力轉(zhuǎn)移軌道的近似解析解具備有效性和一定的普適性，可以為Gauss偽譜法提供較為有效的控制變量初始猜測值；

2)對小推力轉(zhuǎn)移軌道的近似解析解進(jìn)行解算和處理所獲得的控制加速度估計(jì)值，作為Gauss偽譜法中控制變量的初始猜測值，可以顯著提高Gauss偽譜法的迭代速度；

3)采用Gauss偽譜法將小推力軌道轉(zhuǎn)移的最優(yōu)控制問題離散化為NLP問題，并結(jié)合形狀法提供的控制變量初始猜測值，可有效解決共線平動點(diǎn)附近周期軌道間的小推力軌道轉(zhuǎn)移問題，同時也可為中國的深空探測事業(yè)和平動點(diǎn)軌道的工程應(yīng)用奠定一定的理論基礎(chǔ)。

本文提出的近似解析解的系數(shù)通過試湊法配合后期過程約束檢驗(yàn)的方法確定，缺乏一定的理論指導(dǎo)，因此后續(xù)工作可以對過程約束指導(dǎo)下的近似解析解進(jìn)行研究。