周敬，胡軍，*

1. 北京控制工程研究所，北京 100190 2. 空間智能控制技術重點實驗室，北京 100094

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

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

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

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

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

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

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

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

1.1 動力學模型

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

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

(1)

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

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

圖1 地月系統(tǒng)CRTBP下的會合坐標系O-XYZ和 L1會合坐標系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為會合坐標系下L1點距離主天體公共質心的距離；γ1為L1點與最近主天體之間的距離。

因此，L1會合坐標系下航天器運動的動力學方程為：

(2)

(3)

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

(4)

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

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

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

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

‖u‖2≤2DU

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

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

2.1 軌道轉移近似解析解

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

假設在L1會合坐標系下，小推力軌道轉移的起點狀態(tài)為x0，終點狀態(tài)為xf，即

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

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

(5)

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

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

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

(6)

速度矢量表達式為：

(7)

加速度矢量表達式為：

(8)

(9)

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

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

2.2 Gauss偽譜法的離散化方法

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

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

滿足約束條件：

(1)約束方程離散化

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

(10)

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

式中：g(τ)為以各插值節(jié)點為根的多項式，即

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

(11)

(12)

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

(13)

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

(2)末端約束離散化

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

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

(14)

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

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

(3)路徑約束離散化

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

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

(15)

式中：k=1,…,K。

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

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

(16)

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

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

(17)

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

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

3 仿真分析

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

表1 Halo軌道間小推力軌道轉移相關參數(shù)

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

3.1 轉移時間較短

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

表2 轉移時間較短時不同初始猜測值下的 GPOPS仿真結果

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

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

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

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

3.2 轉移時間較長

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

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

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

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

表3 轉移時間較長時不同初始猜測值下的 GPOPS仿真結果

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

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

4 結束語

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

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

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

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

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