航天器燃耗最優軌道直接／間接混合法延拓求解

孟雅哲＊

航天器燃耗最優軌道直接／間接混合法延拓求解

孟雅哲1，2，＊

1．中國科學院空間應用工程與技術中心 太空應用重點實驗室，北京 100094 2．中國科學院大學，北京 100049

針對轉移時間和始末狀態固定的航天器燃耗最優軌道的求解，給出了一種延拓方法：以雙脈沖軌道為初值，首先求解全程推進軌道，然后逐步增加推力幅值，應用直接／間接混合法依次求解所有推力幅值下的、滿足包括開關函數在內的所有必要條件的轉移軌道，包括連續和脈沖推力軌道。通過基于開關函數曲線變化趨勢的開關序列預設方法，以及基于已有優化結果的延拓步長自適應方案，實現了延拓方法的自動運行。為實現該延拓方法，給出了適用于改進春分點根數模型的脈沖最優轉移軌道主矢量必要條件，推導了無推力軌道段改進春分點根數協態變量狀態轉移矩陣。通過3個算例對延拓求解會遇到的不同情況進行了具體說明。延拓方法可以看作現有直接／間接混合法的進一步完善與拓展，延拓過程和結果有助于對燃耗最優軌道與系統參數之間的關聯獲得更為深刻的認識。

改進春分點根數；燃耗最優軌道；直接／間接混合法；延拓；開關序列預設；步長自適應

給定轉移時間與始末狀態，求解燃耗最優轉移軌道是航天器軌道設計的基本問題，也是求解更加復雜的軌道優化問題的基礎。該方面的標志性研究可以追溯到1963年，Lawden［1］應用變分法和 Weierstrass條件（二者等價于Pontryagin極小值原理）推導了脈沖推力燃耗最優轉移軌道需滿足的最優性必要條件，建立了軌道優化的主矢量理論。基于主矢量理論，Handelsman和Lion［2］、Jezewski和 Rozendaal［3］給出了根據主矢量曲線對最優性必要條件的滿足情況逐步調整脈沖作用并最終得到脈沖最優轉移軌道的方法。主矢量理論也可用于求解連續推力轉移軌道，脈沖推進序列對應為連續推力的bang－bang控制形式，即包含推力連續作用的推進段（稱“開”段）和無推力作用的滑行段（稱“關”段）的軌道形式。基于最優性必要條件求解轉移軌道的方法一般被稱為間接法，即求解滿足必要條件的協態變量從而給出推力作用而非直接求解使性能指標最優的推力作用的方法。間接法中協態變量的求解通常采用基于梯度的算法實現，對所提供的初值十分敏感。另一大類軌道優化方法稱為直接法，其發展得益于非線性規劃（NonLinear Programming，NLP）方法和相應軟件包的發展。直接法通過連續系統離散化將軌道優化問題轉換為參數優化問題，應用NLP方法求解該問題，即對性能指標尋優直接得到推力作用，如直接轉錄法［4］、配點法［5］、偽譜法［6］等。直接法不再應用協態變量方程與主矢量，計算量較小但難以判讀解的最優性。綜合間接法與直接法，同時利用NLP與協態變量方程也可以有效地求解軌道優化問題，該方法被稱為“直接／間接混合法”（后簡稱為“混合法”）［7－8］。混合法最基本的特性是在構造NLP問題時直接應用協態變量方程等部分最優性必要條件，并通過NLP尋優盡可能彌補剩余的最優性必要條件（與充分條件）。相對于直接法，混合法所構造NLP問題的優化變量更少，求得的解可滿足更多的最優性必要條件；相對間接法而言，應用NLP尋優可在一定程度上降低協態變量敏感性，但所得結果可能并非最優。更多關于直接法、間接法和混合法的知識可參看綜述性論文［9－12］。

協態變量初值估計是間接法和混合法的難點之一。Dixon和 Bartholomew－Biggs［13］給出的應用具有物理意義的控制量和協態變量相互轉換（Adjoint－Control Transformations，ACT）對協態變量初值進行估計的方法應用較廣［14－16］。此外還有：應用已有解（如直接法所得連續推進軌道解、最優脈沖解等）給出協態變量初值［17－20］；根據已有優化結果總結協態變量初值的取值規律以給出初值［21－22］；依具體軌道機動任務簡化模型并對之解析分析以獲取協態變量初值［23］；以及應用最優性必要條件中協態變量相關性質給出協態變量初值［19，24］等協態變量初值估計方法。

當上述方法所得協態變量初值仍無法保證某復雜問題的順利求解時，可找出與該復雜問題有統一表達的、容易計算的簡單問題，該統一表達中含有一個變量，變量的兩個不同取值對應所選的簡單問題和需求解的復雜問題，此時可從簡單問題出發，逐步改變該變量的值，依次求解對應的問題直至得到所需復雜問題的解，求解過程中應用前一個問題的解給出后一個問題的初值——這就是延拓方法的基本思想。Minter和Fuller－Rowell［25］、劉滔等［26］應用延拓方法求解了最短時間轉移軌道；Shen等［27－28］應用延拓方法求解了地月限制性三體模型下雙脈沖轉移軌道，并結合延拓方法和主矢量理論求解了二體模型下多脈沖轉移軌道，還指出該方法有潛力尋求全局最優解。

本文主要處理始末狀態和轉移時間固定的bang－bang控制轉移軌道的優化問題，此問題中協態變量取值可決定其開關序列，開關序列切換造成的強非線性會在很大程度上增強軌道求解對協態變量初值的敏感性［29］。克服開關序列造成的敏感性的方法可分為兩類：① 不含延拓的預先設定開關序列；② 通過延拓求解開關序列。

不含延拓的預設開關序列的方法可分為兩類：一是事先設定開關序列進行求解，然后驗證解是否滿足必要條件，若不滿足則對開關序列進行調整，再次求解和驗證，直至得到滿足必要條件的解，Fowler和 O’Neill［30］、Enright和 Conway［31］均應用此法求得bang－bang控制轉移軌道，此類方法適用于開關切換次數較少的情況；二是經過對軌道設計要求的分析設定開關序列進行求解，但不對所得優化軌道進行開關函數等必要條件的驗證，Gao［32］、Zuiani和 Vasile［33］應用類似方法實現了多圈bang－bang控制轉移軌道求解，此類方法可以快速求得有較多開關切換的轉移軌道，但無法保證其解滿足所有最優性必要條件。

延拓方法可以求解含較多開關序列的bangbang控制轉移軌道，且可保證所得解滿足包括開關函數在內的所有最優性必要條件。較常見的求解bang－bang控制軌道的延拓方法是應用包含延拓參數的性能指標，實現從能量最優到燃耗最優的延拓求解（文獻描述中若無特定說明均為此類情況），延拓參數某一取值對應的問題用間接法求解：Bertrand和 Epenoy［34］給出了包括從能量最優到燃耗最優在內的3種性能指標延拓形式，能量最優等簡單問題的初值隨機給出。Haberkorn等［29］、Gergaud和 Haberkorn［35］分別基于改進春分點根數模型和位置速度模型，從能量最優轉移軌道出發延拓求取地球引力場內多圈燃耗最優轉移軌道，能量最優轉移軌道通過延拓初始條件求取。文獻［29］中還對延拓方法進行了簡要介紹，并分析了不同推力幅值下轉移時間和末端真經度對燃耗的影響。Jiang等［18，36］實現了位置速度模型下能量最優到燃耗最優轉移軌道的延拓，文獻［18］和文獻［36］分別應用偽譜法和粒子群算法給出了能量最優轉移軌道初值，其中文獻［36］求解的是帶有引力輔助的連續推力轉移軌道。Zhang等［37］應用延拓方法進行了限制性三體模型平動點任務的時間最短、能量最優和燃耗最優軌道的延拓求解，應用ACT作為最短時間問題的初值，降低推力幅值以逐步增加轉移時間，用最短時間軌道的解設定能量最優軌道的轉移時間，再從能量最優軌道延拓求解燃耗最優軌道。Petukhov［38－39］首先應用能量最優軌道對應的兩點邊值問題（Two Point Boundary Value Problem，TPBVP），將協態變量初值為零作為延拓的初始解，構造包含延拓參數的邊值條件實現燃耗最優軌道的求解，所給求解過程可保證轉移軌道圈數的固定，而后進行從能量最優轉移軌道到燃耗最優轉移軌道的延拓求解，文中給出了包含延拓參數的狀態和協態變量微分方程。Li和Xi［40］、Mitani和 Yamakawa［41］應用延拓方法實現了相對運動模型下的燃耗或總速度脈沖最優的bang－bang控制軌道求解：文獻［40］延續了應用從能量最優到燃耗最優的延拓，鑒于模型相對簡單，能量最優軌道的初值可解析計算；文獻［41］則是求解了包括推力方向約束在內的燃耗最優軌道，推將力方向約束作為懲罰函數寫入性能指標，延拓逐步添加此約束并求得燃耗最優解。

應用混合法實現單步求解的延拓方法研究較少，其中混合法所用的NLP需基于所預設開關序列進行構造，而開關序列的預設則基于對已有優化結果開關函數曲線的分析。與基于間接法的延拓不同，此法可獲得延拓參數取延拓區間任意值的、滿足包括開關函數在內所有必要條件的bang－bang控制解。Chuang等［42］進行了較短轉移時間到較長轉移時間的延拓，由開關函數曲線預測開關序列，給出了添加開關段時開關切換時刻的初值估計方法。朱小龍等［43］則針對 Hill－Clohessy－Wiltshire方程，首先構造“開－關－開”序列優化模型，逐步降低推力幅值上限以實現雙脈沖轉移軌道到全程推進轉移軌道的延拓，其中雙脈沖軌道及其協態變量可解析求解，然后構造包含所有最優性必要條件的總速度脈沖最小的直接／間接混合法模型，逐步增加推力幅值上限以實現從全程推進軌道到多脈沖軌道的延拓，文中應用主矢量曲線變化趨勢對開關序列進行預設，但對主矢量曲線變化情況說明尚可完善，且未給出程序可行的開關序列預設方法。

總結前述延拓求解航天器軌道的文獻：通過延拓參數的選取，延拓方法可以實現不同性能指標［18，29，36，38］、不 同 始 末 狀 態［29，39］、不 同 轉 移 時間［37，42］、不 同 推 進 器 參 數［43］、不 同 重 力 加 速度［25－26］、不同攝動大小、轉移軌道約束的不同約束程度［41］、bang－bang 求 解 中 不 同 開 關 序 列［42－43］所對應解的連續轉換。延拓過程在得到所需目標軌道之前，需進行多次軌道求解，無推力軌道段的解析計算有利于降低數值優化過程的計算量。Glandorf［44］給出了直角坐標系位置與速度所對應的協態變量轉移矩陣的解析形式；Fernandes［45］給出了二維極坐標系下位置速度對應協態變量的解析解；到目前為止，本文作者未見已有文獻給出改進春分點根數協態變量解析解。Xu［46］、Jamison和 Coverstone［47］對無推力軌道段結束時刻的確定進行了研究：文獻［46］分析了無推力軌道段開關函數的變化頻率，并依該頻率檢測求取開關函數零點（即無推力軌道段結束時刻）；文獻［47］進一步研究了文獻［46］的方法并對其進行了改進。此外，延拓步長越小，單步延拓可給定的協態變量初值越接近優化解，開關序列預設越準確，但延拓步長太小會增加中間問題的個數，延長延拓求解的時間，然而，延拓求解航天器軌道的文獻中幾乎沒有對如何取延拓步長進行過討論。

本文采用改進春分點根數模型，用基于混合法［7］的推力幅值延拓方法求解燃耗最優轉移軌道。與文獻［43］一致，先用雙脈沖軌道求解全程推進轉移軌道，再從全程推進轉移軌道出發，用基于主矢量曲線的變化趨勢預測開關序列，延拓求得更大推力幅值下的連續推進和脈沖轉移軌道。相對文獻［43］而言，本文對推力幅值延拓中開關序列變化情況進行了更完善的描述，給出了實現該延拓的自動運行算法，算法重點是對開關序列變化的自動判斷和簡單的延拓步長自適應規則。此外，本文采用逐個將短時大推力作用化為脈沖作用的方法求得脈沖推進轉移軌道；給出了以常值推力加速度（不考慮航天器質量變化）轉移軌道為初值延拓求解常值推力轉移軌道（考慮航天器質量變化）的延拓過程。文章給出3個算例對延拓方法進行說明，驗證了方法有效性。此外，本文基于文獻［1，47］給出了適用于改進春分點根數模型的的脈沖最優軌道主矢量必要條件，以文獻［47－48］為基礎給出了脈沖作用前后處改進春分點根數協態變量的轉化關系，推導了無推力軌道段改進春分點根數協態變量轉移矩陣的解析形式。

1 常值推力加速度模型

1．1 極小值原理必要條件

改進春分點軌道根數［48］是一類在軌道傾角和偏心率接近零時無奇點的軌道根數，各根數定義為

式中：a、e、i、ω、Ω和θ為經典軌道根數，且θ為真近點角。用x＝ ［p f g h kL ］T表示改進春分點根數向量，則航天器軌道運動微分方程為

式中：F∈ ［0 ，Fmax］為推力幅值；m≡m0為航天器質量；M和D分別為6×3和6×1的改進春分點根數的矩陣函數，文獻［7］附錄中給出了M和D及其對改進春分點根數偏導的表達式；α為推力方向矢量，3個分量為其在RSW坐標系3個坐標軸上的投影（即α為其在RSW坐標系下的表達）。RSW坐標系3個坐標軸方向定義為：R軸沿航天器矢徑方向，軌道面內垂直矢徑沿速度方向為S軸正向，W軸沿軌道角動量方向。

設航天器轉移軌道的初始和末端時刻分別為t0和tf，對應軌道狀態滿足：

性能指標（等價于燃耗最優）為

至此，可給出始末狀態和轉移時間固定的轉移軌道優化問題如下：

確定F （t）和α （t）（t∈ ［t0， tf］），使 得 式（2）表述的系統滿足邊界條件式（3），且使式（4）所示J的取值盡可能小。

應用Pontryagin極小值原理（后簡稱為極小值原理）將轉移軌道優化問題轉化為TPBVP，過程如下：

首先構造如式（5）所示的Hamilton函數：

式中：λ為x對應的協態變量，需滿足式（6）所示的微分方程。

當F＝0N時，協態變量λ可應用附錄A中的公式進行解析運算。

根據極小值原理，最優解F＊（t）和α＊（t）（t∈［t0，tf］）應當使H最小，所需滿足的必要條件如式（7）和式（8）所示。

式中：1－ MTλ 與H／F＝ （1－ MTλ ）／m0同號。不考慮奇異情況（MTλ ≡1在某一連續時段恒成立），有TPBVP如下：

確定 Fmax、λ（t0）（或 λ（tf）），使 得 式 （2）、式（6）～式（8）表述的系統滿足邊界條件式（3）。

協態變量的等比縮放不改變燃耗最優軌道：以x0和λ（t0）為t0時刻初值，α（t）和F （t）按式（7）和式（8）求取，積分方程如式（2）和式（6），得到任意時刻的x（t）和λ（t）。若以x0和Kλλ（t0）（Kλ為包括1在內的任意正實數）為初值，α（t）和F （t） 按式 （7）和式 （9）求 取，積 分 如 式 （2）和式（6），則得到任意時刻的x1（t）和λ1（t），兩次積分結果之間滿足式（10）。

為描述方便，定義式（11）所示的S為開關函數。

1．2 脈沖軌道主矢量必要條件

文獻［1］以慣性直角坐標系下位置和速度表達的二體模型為例，給出了燃耗最優脈沖轉移軌道主矢量必要條件，其中主矢量為速度對應的協態變量。根據極小值原理，以位置速度表達的相對運動模型和限制性三體模型等其他航天器模型，也有以速度對應協態變量為主矢量的脈沖最優轉移軌道必要條件。但改進春分點根數模型下的主矢量形式與文獻［1］中不同，本節以文獻［47］中對應不同狀態量的協態變量之間的轉換關系出發給出該改主矢量表達式。

文獻［1］中的二體運動方程可寫為式（2）所示的形式，其中：狀態量對應為 ［rTvT］T，r和v分別為航天器的位置矢量和速度，各分量為其在慣性直角參考（RIRF）坐標系中各坐標軸上的投影；系數矩陣Mrv＝ ［0 I3×3］T，Drv＝ ［0－μrT／r3］T（其中μ為所用RIRF坐標系下的引力常數，r＝r）；推力方向矢量α的3個分量也為其在RIRF坐標系3個坐標軸上的投影。用λrv表示該模型對應的協態變量，則改進春分點根數協態變量λ和λrv存在如式（12）所示的轉換關系［47］，該轉換可用于脈沖作用前后λ的相互變換。

式 中：Rxrv＝ ［x／rTx／vT］和 Rrvx＝［rT／x vT／x ］T為改 進春分 點根 數 和 位 置速度之間的Jacobi矩陣，Rrvx的計算可參考文獻［49］，Rxrv為其逆矩陣。

根 據 x＝Rxrv［rTvT］T，M和Mrv有 如 下關系［47］：

式中：ΦRSW→RIRF為從RSW坐標系到RIRF坐標系的轉換矩陣。

由式（12）和式（13）有

式中：ΦRIRF→RSW為從RIRF坐標系到RSW坐標系的方向余弦矩陣，為正交矩陣；λv表示速度v對應的協態變量，是文獻［1］中定義的主矢量。可見，改進春分點根數運動方程下與λv等價的主矢量為ΦRSW→RIRFMTλ。相應的燃耗最優脈沖轉移軌道主矢量必要條件如下：

1）ΦRSW→RIRFMTλ 和d（ ΦRSW→RIRFMTλ）／dt在［t0， tf］上連續。

2） MTλ ≤1（t∈ ［t0，tf］ ），其中 MTλ ＝1在且僅在脈沖施加時刻成立。

3）脈沖作用點處MTλ與α在RSW坐標系中的表達平行且方向相反（即式（7））。

4）d MTλ／dt＝0在除始末時刻以外的所有脈沖作用點處成立。

實際上，當式（8）中Fmax＝∞時，推力段近似為一個速度脈沖，即可得上述條件2）和4）。

2 推力幅值延拓方案

2．1 延拓方案與方法

延拓是求解參數連續變化的非線性問題的有效方法，本文延拓求解Fmax取區間 ［Fmin，∞］上不同取值時的燃耗最優轉移軌道：Fmax＝∞對應脈沖解；認為Fmax存在下限Fmin且Fmax＝Fmin對應全程推進軌道；Fmax＜Fmin時，在給定時間內無法實現所要求的軌道轉移。本文假定滿足最優性必要條件式（6）、式（7）和式（9）的全程推進軌道為最小推力幅值的軌道，整體延拓方案如圖1所示：延拓從簡單易求解的雙脈沖軌道出發，多數情況下，雙脈沖軌道不滿足主矢量必要條件，需先求解滿足式（6）的全程推進軌道，然后在滿足式（6）、式（7）和式（9）的前提下逐步添加滑行段，直至獲得燃耗最優脈沖軌道（這是由于任意滿足式（6）的全程推進軌道，總可以通過協態變量等比縮放化為滿足式（6）、式（7）和式（9）的軌道解。）。本文主要闡述此類情況。少數情況下，雙脈沖軌道滿足主矢量必要條件，此時可在滿足式（6）、式（7）和式（9）的前提下，首先求解短時大推力“開－關－開”軌道，然后逐步減小推力幅值、增加推進段時長，直至獲得全程推進軌道。此類情況只在4．1節用算例進行說明。

本文采用結合開關序列預設和步長自適應的DCM （Discrete Continuous Method）［29］延 拓 方法：Fmax順序取 ［Fmin，∞］內的有限個離散點，依次求解各點對應的燃耗最優軌道，燃耗最優軌道的求解采用混合法，即構造NLP問題并進行求解，所得解滿足式（6）、式（7）和式（9），構造 NLP時需要預設開關序列，前一個NLP解中各項可作為下一個NLP對應項的初值。延拓方法的核心是基于開關函數曲線變化趨勢的開關序列預設，Fmax基于開關序列預設情況進行取值，詳見2．3節。圖1中各NLP問題的定義見表1、表2、表4和表5。

2．2 雙脈沖到全程推進

2．2．1 雙脈沖軌道春分點根數協態變量

用 Lambert 問 題［50－53］的 求 解 等 成 熟 算法［54－55］求得雙脈沖軌道后，即可計算出改進春分點根數的協態變量。由于所求脈沖矢量分量表達在RIRF坐標系下，應首先給出RIRF坐標系到RSW坐標系的坐標轉換矩陣ΦRIRF→RSW為

RSW坐標系推力方向矢量α為

式中：ΔvRSW＝ΦRIRF→RSWΔv，Δv和 ΔvRSW分別為脈沖矢量在RIRF坐標系和RSW坐標系下的表達。應用式（16）求得始末速度脈沖對應的α0和αf。由1．2節所示的主矢量必要條件，令始末脈沖作用點處 MTλ ＝1，由式（7）可得

由附錄A中的式（A3）有λ（t）＝Mλλ（t0），從而有

應用式（18）求解λ（t0），其系數矩陣的秩若小于6，可求其最小二乘解；應用式（A3）求得λ（tf）。上述過程求得的λ（t0）和λ（tf）即為無推力軌道段始末的協態變量。

2．2．2 NLP1和NLP2構造說明

全程推進軌道通過表1所示的NLP求取，表中：t1∈ ［t0，tf］為一中間時刻；x （t1）和λ（t1）為對式（2）和式（6）從t0至t1數值積分所得結果；x′（t1）和λ′（t1）為從tf反向積分到t1所得結果。表1中NLP1的約束條件只保證全程推進軌道連續，NLP2則保證狀態和協態變量均連續。

2．2．3 NLP1和NLP2優化變量初值給定

NLP1的初值設定如下：雙脈沖軌道兩端協態變量初值為優化變量中對應協態變量的初值，對應Fmax取值為

表1 求解全程推進軌道所用NLP模型（NLP1、NLP2）Table 1 NLP models for full－propel trajectory optimization（NLP1，NLP2）

式中：Δv0和Δvf為雙脈沖解始末脈沖的大小；KF可從1開始以0．1為步長逐步增加，取值為使x （t1）－x′（t1） 為局部極小。

t1的初值滿足

根據實際計算經驗，依上述方法給出初值求解NLP1，并以其優化解為初值求解NLP2即可得所需全程推進軌道。若遇難以求解的問題，可應用文獻［43］中逐步降低Fmax、求解推進時長逐漸增加的“開－關－開”轉移軌道給出NLP1的初值，應用文獻［16，39］中含有延拓參數的約束條件求解NLP2。換句話說，可先嘗試以最大延拓步長對問題進行求解，若求解不成功，再降低延拓步長徐徐圖之。

2．3 全程推進到bang－bang控制

以全程推進軌道為初值求取不同Fmax取值下的bang－bang控制軌道可通過構造和求解多個NLP3實現。2．3．1節給出NLP3的結構，其優化變量包括始末協態變量和所有開關時刻點，其中協態變量初值依DCM設定為前一個NLP3優化解的對應項，開關時刻點的初值設定，即開關序列的預設，則需依前面若干個NLP3優化解對應開關函數曲線的變化進行。延拓參數Fmax的變化與開關序列的預設相關。2．3．2節敘述預設開關序列和Fmax的思路，其中將給出造成開關序列改變的各種開關函數曲線變化情況，2．3．3節將給出完整的、包括開關時刻點預設和延拓步長自適應的開關時刻點和Fmax設定算法。

2．3．1 NLP3構造說明

全程推進到滿足極小值原理必要條件的bang－bang控制軌道的延拓過程通過求解一系列如表2所示的NLP問題實現，表中各變量意義如下：t1，1，t1，2，…，t1，M1為各開 關點時 刻值，M1為開分別表示前向和反向積分的結果；約束函數除保證狀態變量和協態變量連續外，還用S （t1，j）＝0（j＝1，2，…，M1－1）保證開關點處的開關函數約束，用式（21）中dS／dt的正負保證開關函數曲線的形狀。

表2 求解bang－bang控制軌道所用NLP模型（NLP3）Table 2 NLP model for bang－bang control trajectory optimization（NLP3）

式中：M 的表達式參見附錄B。

其中推進段通過對式（2）和式（6）數值積分求得，推力的大小和方向滿足式（7）和式（9），滑行段改進春分點根數改變遵循開普勒軌道的規律，協態變量用式（A3）計算，延拓過程中均如此計算。

2．3．2 開關序列及Fmax預設

假設已成功求解l個NLP3，本節討論第l＋1個NLP3的開關序列預設和Fmax設定思路。開關序列預設通過跟蹤開關函數曲線的變化進行，包括統一開關函數曲線幅值、判斷開關序列是否發生變化、預設開關點時刻3個步驟。Fmax依開關序列預設情況進行給定。

將開關函數曲線幅值統一設定為10（人為設定參數，可更改）：設第n（n＝1，2，…，l）個滿足必要條件的NLP優化解對應開關函數曲線為Sn，用式（22）所示的Kλ對優化解中始末協態變量進行式（23）所示的縮放。

幅值統一有益于判斷開關函數曲線的變化趨勢，曲線的變化趨勢可由其特征點取值量化反映，特征點指曲線的始末端點、波峰和波谷。開關序列變化和特征點變化的關系見表3。

表3 特征點變化和開關序列變化對照表Table 3 Relationship between changes of feature points and switching sequences

圖2～圖4為造成開關函數變化的特征點所在的開關函數段示意圖。其中，圖2為波峰波谷出現和消失時的開關函數曲線，圖中所示的特征點出現和消失的情況在距零較遠時也可能發生，但不影響開關序列；圖3和圖4分別為特征點從0＋到0－和從0－到0＋變化時的開關函數曲線。圖中：S＊為一給定數值。

判斷是否出現圖2～圖4中情況的方法如下：查找當前開關函數曲線為0的點的個數是否與開關點數目相同，不同則認為出現了圖2對應情況，相同則認為未出現。

若未出現圖2所示情況時，查找開關函數值最接近0的特征點，判斷其是否足夠接近0，即其絕對值是否小于給定的ε（ε為一個足夠小的正實數），“否”則判定開關序列不變，“是”則查找最近連續2～3次NLP優化解中該特征點開關函數值的變化規律，若其絕對值單調遞減，則認為該特征值將過0。Sl中該特征點開關函數值的正負標識其過0方向，為負認為出現0－→0＋情況，為正則對應0＋→0－變化的情況。其中0＋和0－分別代表符號為正和負的絕對值足夠小的數。

開關序列不變時，第l＋1個NLP3對應Fmax依式（24）取值，所有優化變量初值均為第l個NLP3優化解的對應部分。

式中：0＜δ＜1為給定常數。

若判定開關序列將發生變化，則通過截取開關函數曲線預設所有開關點時刻：給定數值S＊，將當前NLP優化解對應開關函數曲線上所有取值為S＊的點預設為新的NLP問題的開關點。圖2對應情況設定S＊＝0；特征點開關函數值從0＋→0－變化時，用式（25）計算S＊（如圖3所示）；從0－→0＋變化時，則用式（26）計算S＊（如圖4所示）。

以圖4（a）所示情況為例說明單個特征點過0在整個開關函數曲線處理時的作用：如圖5所示，Sl曲線第1個波峰將出現0－→0＋，用式（26）計算常數S＊，Sl曲線上所有值為S＊的點即預設為第l＋1個問題的開關時刻點，由于Sl（t0）＜S＊，t0時刻為“開”段，開關序列預設完成，第1個波峰處開關序列變化如表3首行所示。由Sl＋1曲線可以看出，第l＋1個NLP優化解的開關切換點時刻與初值有小幅變化。開關點時刻設定成功后，應用FmaxtF不變（即性能指標式（4）不變）設定新NLP對應Fmax的值，如式（27）所示。

式中：tF為所有推進段總時長；Flmax和tFl表示第l個 NLP問題的Fmax和tF。

圖5也給出了滿足式（9）的全程推進軌道和脈沖軌道開關函數曲線示意圖。由圖5可以看出，全程推進軌道到脈沖軌道變化過程中應有如下整體趨勢：Fmax逐漸增大，滑行段總時長逐漸增加。所以當式（27）所得Fmax＜時（極少數情況下），應提高Fmax值，如設定Fmax＝以保證整體延拓趨勢。

2．3．3 NLP3開關切換時刻初值給定

2．3．2節中變量ε和δ的取值決定了推力幅值延拓中Fmax的變化幅度。ε和δ越小，Fmax的變化越小，開關序列和NLP初值設定越準確，但延拓整體需花費時間越長。延拓算法中將根據NLP3求解情況自適應ε和δ的值，規則可簡單描述如下：若當前NLP求解失敗，則減小ε和δ；若NLP連續快速收斂，則適當增大ε和δ。

設已有l－1個轉移軌道優化解，求解第l個NLP3后，第l＋1個NLP3的初值設定過程如下：

1）調整延拓步長參數δ和ε。

第l個NLP3是否成功求解，“否”則降低δ，令δ＝δ／3，若δ＜0．000 1（人為設定值）則降低ε：令ε＝ε／2，令l＝l－1轉2）設定初值。若第lnNLP3（nNLP3∈N＋為人為設定參數，建議不大于10）至第l個NLP3均成功求解且對應δ相同，則設定δ＝ （1＋δ）n′NLP3－1（其中n′NLP3≤nNLP3也為正整數）。

為保證延拓速度，可為δ設定下限值δm，如0．000 1，若連續 N1（N1為人為設定的常數，常取5～10中任意整數）個NLP對應δ均小于δm，則設定δ＝δm。

2）判斷是否出現圖2所示情況并進行處理。

設S＊＝0，找到Sl上所有取值為S＊的點，若點的個數與第l個優化解的開關切換次數不同，認為出現了圖2所示情況，Sl上所有S＊＝0的點即為第l＋1個NLP3的開關序列的初值，Fmax用式（27）求解。若次數相同，則轉3）。

3）判斷第l－1次和第l次優化解開關序列是否相同，并進行處理。

判斷第l－1個和第l個NLP3的開關序列是否相同，“是”則開關切換點時刻初值為第l個NLP3優化解的對應值，Fmax用式（24）計算。“否”則轉4）。

4）判斷是否有開關函數曲線特征點過0。

判斷距0最近的特征點開關函數值是否小于臨界值ε，“否”則認為開關序列不發生變化，初值設定同3）；“是”則轉5）。

5）若有特征點過0，判斷具體情況并進行處理。

查看之前N2（N2為人為設定常數，通常取2～5中任意整數）個NLP3的優化解的歸一化開關函數曲線特征值，判斷對應特征值是否越來越接近0，“否”則認為開關序列不變，初值設定同3）；“是”則判斷該特征點開關函數值的正負，“正”則認為出現圖3所示情況，用式（25）計算S＊，“負”則認為出現圖4所示情況，用式（26）計算S＊，Sl上所有取值為S＊的點為開關切換時刻初值，用式（27）計算Fmax。

2．4 脈沖軌道求取

2．4．1 NLP4構造說明

bang－bang控制軌道推進段時長足夠小時，采用表4所示的NLP問題逐個將短時大推力段轉為脈沖軌道。表中：性能指標Δm 為應用式（28）和式（30）估計的總燃料消耗；t2，1，t2，2，…，t2，M2為所有的開關點和脈沖點時刻；v1，v2，…，vN對應速度脈沖大小，其中N≤M2；r、v、λr和λv為正向計算的結果，r′、v′、λ′r和λ′v為反向計算的結果，可見前幾個約束是t＊2（即t2，ceil（M2／2））處位置速度及對應協態變量連續，而系統方程的處理及開關函數還是按改進春分點根數進行計算的；剩余約束的含義同NLP3。當t＊2為某個脈沖作用點時，約束條件是脈沖作用后的狀態和協態變量相等。

式中：mbefore和mafter分別為脈沖作用前和作用后的航天器質量；v為脈沖大小；ge和Isp分別為地球表面重力加速度和推進器比沖，本文中二者均為常數。

表4 求解含有脈沖的軌道所用NLP模型（NLP4）Table 4 NLP model for impulse trajectory optimization（NLP4）

2．4．2 NLP4優化變量初值給定

每個NLP4優化后，判斷最短推進段時長t2，k＋1－t2，k是否小于設定參數timp，“否”則開關序列和脈沖作用情況不變，用式（24）增加Fmax的值進行下一個NLP4的求解；“是”則將此推進段化為脈沖，脈沖時刻初值為 （t2，k＋t2，k＋1）／2，脈沖點個數加1，新添加速度脈沖的大小vnew如式（29）所示。

式中：縮放系數Kv∈［0．5，1］，取值應使得 NLP6中約束條件對應項－r′和v－盡可能小，為當前所添加脈沖的時刻點。通常從1逐步減小進行搜索，Kv＝1時vnew對應脈沖燃耗與所代推力段相同。

timp為人為設定的參數，取值范圍為40～90s。若總推進時長小于時NLP3優化解對應各推進段時長相差不大，可直接將所有推進段均化為脈沖可一次求解。常將設為0．1～0．2tu（tu為時間對應的歸一化單位，將在第4節中介紹）。

3 考慮航天器質量變化的bang－bang控制軌道

本節用常值推力加速度模型某一Fmax對應優化解作為初值，構造NLP問題求解考慮航天器質量變化（常值推力）的Fmax燃耗最優軌道。實際上，常值推力模型可直接用于推力幅值延拓，但實際運行中常值加速度模型對應的NLP問題較常值推力模型對應的NLP問題更容易求解（模型相對簡單），且對常值加速度和常值推力模型的分析比較，有利于對主矢量和開關函數關系的深刻理解。

3．1 常值推力模型

設m 滿足微分 方 程 式 （30），則 式 （2）和式（30）組成系統方程。

燃耗最優對應目標函數為

根據極小值原理，構造Hamilton函數為

式中：λm為m 對應的協態變量，其微分方程為式（33），滿足末端條件式（34）。

使得式（32）所示 Hamilton函數達極小的F＊滿足式（35）和式（36），Sm為對應開關函數。

由式（7）、式（30）、式（33）和式（36）可知：

與式（11）進行對比可知Sm與S，即m≠0與m＝0的開關函數曲線增減性一致，但過0時刻會有差別，差別的大小取決于質量變化的多少。

協態變量等比例縮放不改變燃耗最優軌道：以x0、m0、λ0和λm0為t0時刻初值，α（t）和F（t）依式（7）、式（35）和式（36）取值，分別按式（2）、式（30）、式（6）和式（33）進行積分，得到x（）t、λ（）t和Sm（）t；若以x0、m0、Kλmλ0和Kλmλm0為初值進行同樣的積分過程，可得到x1（）t、λ1（）t、S1m（）t，有

3．2 NLP5和NLP5－1構造及初值給定

表5為求解常值推力轉移軌道的NLP問題。表中：λm（t0）和λm（tf）為始末時刻質量的協態變量；其余變量意義同NLP3。

NLP5中忽略約束條件λm（tf）＝－1，需對其優化結果進行如式（39）所示變換才能得到滿足所有必要條件的解。

以常值加速度優化軌道為初值，求常值推力優化軌道的過程如下：

1）構造NLP5進行求解，若順利求解，對優化解進行式（39）所示的協態變量縮放得到最優解；若NLP5無法順利求解，則用2）。

2）構造 NLP5－1，將 NLP5－1的解作為初值按1）處理；若NLP5－1無法順利求解，則用3）。

表5 考慮航天器質量變化的連續推力軌道優化所用NLP模型（NLP5、NLP5－1）Table 5 NLP models for continuous thrust trajectory optimization with variable spacecraft mass （NLP5，NLP5－1）

3）Isp取從大到小的一系列離散值，各值對應NLP5－1延拓求解；Isp為模型設定值的 NLP5－1的解作為初值按1）處理。

首個NLP5或NLP5－1求解時，λm（t0）和λm（tf）初值設定為－1，其余優化變量初值設定為對應NLP5優化解相應量的取值。常值推力加速度到常值推力求解過程中未考慮開關序列的變化，若存在開關序列變化，可依2．3節內容進行開關序列預設。

4 仿真驗證

算例1（4．1節）屬于共面共拱線軌道的雙切轉移，用于說明雙脈沖軌道為脈沖最優軌道的情況；算例2（4．2節）和算例3（4．3節）為非共面轉移，其中算例2中的開關函數曲線出現如圖2所示變化情況，算例3則是一個包含12個軌道周期的算例。文中將簡要敘述各算例延拓過程和延拓結果，并對各軌道根數的變化趨勢進行簡要分析。算例1和算例3為地心系轉移軌道，算例2為日心系轉移軌道。3個算例均包含了從常值加速度到常值推力模型的延拓過程。

實際計算中長度、速度和時間采用無量綱單位au、vu和tu，無量綱單位對應引力常數為1。地心系和日心系的引力常數及所用無量綱單位分別如式（40）和式（41）所示

4．1 算例1

軌道轉移時間和始末改進春分點軌道根數如表6所示，設t0時刻的航天器質量為m t（）0＝2 000kg，地球重力加速度為 ge＝9．806 7×10－3km／s2，比沖Isp＝3 000s為常量。

應用共面共拱線橢圓軌道的雙脈沖最優轉移［54］求得速度脈沖如式（42）所示，始末速度脈沖均為RIRF坐標系下的表達，下文不再贅述。

應用式（18）和式（A3）求得雙脈沖轉移軌道協態變量初始值如式（43），協態變量為無量綱單位計算結果，后面算例也是如此。

表6 算例1中的始末改進春分點軌道根數Table 6 Initial and final modified equinoctial orbit elements for Example 1

繪制雙脈沖軌道主矢量曲線可以看出，雙脈沖轉移軌道滿足所有開關函數必要條件，故以雙脈沖軌道為初值，首先用短時大推力段取代始末脈沖，構造NLP3求解“開－關－開”轉移軌道，而后逐步減小Fmax的值獲得一系列“關”段時長逐漸減小的“開－關－開”轉移軌道，當“關”段時長足夠小時，應用NLP2求得全程推進軌道。首個NLP3初值依2．2．3節內容給出，其中：Fmax＝58．422N（對應KF＝30，此處KF取值較大以保證首個NLP3的順利求解），以燃耗不變為原則給出始、末推進段時長初值分別為47．340s和47．331s。

推力幅值延拓過程中，Fmax的變化規律為解，η初值設定為0．5。延拓過程所得全部優化解的推力作用分布和燃耗隨Fmax的變化如圖6所示。圖6中的燃耗由）／（geIsp）求得，算例2和算例3的LP3延拓中的質量變化也是相同估計，此處不再一一說明。圖7為延拓結果中若干開關函數曲線變化示意圖，曲線的不同顏色標識不同推力大小，虛線表示滑行段，實線表示推進段，可見開關函數曲線形狀基本不變，推進段時長隨推力減小逐漸增大。

算例1中，轉移軌道質量消耗在0．02kg左右，遠小于m （t0），直接構造NLP5可得常值推力燃耗最優軌道，Fmax＝3．948 8N時對應常值加速度解的開關點時刻分別為1 340．914s和1 503．528s，常值推力解的開關點時刻分別為1 340．902s和1 503．543s，可見 m 變化只造成開關點時刻極小幅度的偏差，“開”段總時長增長。常值加速度和常值推力全程推進解Fmax分別為3．933 86N 和3．933 88N。

需要說明的是，本算例為共面轉移，可忽略h和k對應的協態量及其微分方程，并按照h＝k＝0和α（）3＝0化簡式（2）和式（6）中的剩余矩陣元素，按照文中給出的延拓過程進行求解，求得所需燃耗最優軌道后，恢復所有積分數據中h和k的值，即可得三維模型下的燃耗最優軌道。

由表6可知，本算例目標為增加p，圖8給出了NLP3延拓過程中不同推力下改進春分點根數隨時間變化的情況。推力作用方向可用偏航角和俯仰角表達，其中偏航角為推力和推力在RS平面上投影的夾角，俯仰角為推力在RS平面的投影與其在S方向上投影的夾角。對于共面轉移軌道，偏航角為0°。圖9給出了本算例俯仰角隨時間變化的情況。由圖6、圖8和圖9可以看出，Fmax較大時，轉移軌道燃耗較小，推進段內p單調增加、f和g變化幅度較低，推力主要在遠近地點附近，俯仰角接近0°，即推力集中在S方向；隨著Fmax逐漸減小，轉移軌道燃耗逐漸增加，推進段內p隨時間增加率降低，全程推進時甚至出現p單調下降的時段，f和g變化幅度逐漸增大，推力作用從遠近地點逐漸擴散于整條轉移軌道上，俯仰角逐漸增加，即R（航天器矢徑）方向上的推力分量逐漸增加。

4．2 算例2

軌道轉移時間和始末改進春分點軌道根數如表7所示，m （t0）、ge和Isp同算例1。

表7 算例2中的始末改進春分點軌道根數Table 7 Initial and final modified equinoctial orbit elements for Example 2

設轉移軌道含1個完整軌道周期，采用文獻［55］的方法求解Lambert問題得始末脈沖如式（44）所示。

雙脈沖軌道協態變量初值為

NLP1初值設定如下：Fmax＝0．327 39N（KF＝1．8），t1＝391．09d，經 NLP1和 NLP2得到全程推進解。從全程推進轉移軌道延拓求解至得到多脈沖轉移軌道過程中所有優化解推力作用分布燃耗隨推力幅值變化情況如圖10所示。圖10中黑色線段為延拓過程中全程推進和所有bang－bang控制軌道解推力作用的時間分布情況，末端紅色星號標識了脈沖作用時刻，對應藍色線條為這些解對應的燃耗，推進段逐個化為脈沖時燃耗逐次減小，三脈沖轉移軌道不含連續推力作用，所對應推力大小數據無意義。圖11繪制了若干Fmax下的開關函數曲線以描述曲線的形狀及變化情況，表8則用具體數據對其進行了說明。

以第1行數據為例對表8進行說明：全程推進軌道燃耗最優解對應Fmax＝0．226 107N，對應開關函數曲線101．89d時刻處波峰出現如圖4（a）所示從0－到0＋的情況，按照式（26）給出S＊的值，用S＊截取整條開關曲線得到兩個開關切換時刻初值，按照“開－關－開”序列設定NLP3，按式（27）給出下一個 NLP3的Fmax為0．226 109N，用MATLAB中的fmincon對其進行求解得到滿足必要條件的“開－關－開”軌道，而后Fmax逐漸增加，開關序列保持不變不變，至Fmax＝0．226 12N時，圖4（a）情況再次出現，依照表8中的第2行信息構造NLP3以實現開關序列的變換。當圖2（a）所示的情況出現時，表8中給出了相應開關函數曲線的波谷波峰時刻，在Fmax至26．226N后，構造NLP4依次將3個推進段化為脈沖，最終脈沖作用時刻和脈沖大小如式（46）所示。

全程推進常值加速度到常值推力轉移軌道的延拓通過構造NLP5－1和NLP5完成，對應Fmax分別為0．226 11N 和0．217 99N。以 Fmax＝0．215 56N 對應 NLP3優 化 解為初 值，構 造NLP5求常值推力軌道，常值加速度和常值推力軌道的推進段總時長分別為6．021 2d和5．871 4d。

表8 算例2中Fmax延拓過程中的開關序列變化情況Table 8 Variation of switching sequences during Fmaxcontinuation for Example 2

圖12 給出了延拓過程中俯仰角和偏航角的變化情況，可以看出，Fmax從0．226 11N上升到0．282 03N的過程中，俯仰角和偏航角的最大值大幅縮減，對照圖10可發現此段燃耗快速降低，Fmax從0．282 03N繼續上升時，俯仰角基本維持在±5°范圍內，偏航角最大值持續降低，燃耗降低幅度較小。

本算例始末時刻6個軌道根數均不同，幾個改進春分點根數的變化規律不如算例1明確，但滑行段的添加均發生在某個春分點根數隨時間變化率接近零的時刻附近，其變化過程符合燃耗最優的目的，由于篇幅所限，本文未給出軌道根數變化圖。

4．3 算例3

軌道轉移時間和始末改進春分點軌道根數如表9所示，m （t0）、ge和Isp同算例1。

表9 算例3中的始末改進春分點軌道根數Table 9 Initial and final modified equinoctial orbit elements for Example 3

設轉移軌道有20個完整軌道周期，Lambert解［55］對應始末脈沖解如式（47）所示。

雙脈沖軌道始末協態變量如式（48）所示。

初值設定中 KF＝1．8，F1max＝294．07N，t1＝89 277s，構造NLP1和NLP2所示模型得到全程推進解。圖13給出了整個延拓過程——從全程推進解通過求解一系列NLP3得到不同Fmax下的bang－bang控制軌道，而后應用NLP4求得脈沖解——中所有優化解的推力作用分布和燃耗和推力幅值的對應圖。整體而言，每個軌道周期內（真經度L從0到2π）均有若干次開關序列的變化，步長參數初值設為ε＝0．01，δ＝0．01，δm＝0．000 1，程序共有65次開關序列調整，但開關序列最多的燃耗最優軌道含有22個推進段和21個滑行段，即為“關－開－關－…－開－關”的序列，這是由于圈數較多時會累積積分誤差，從而S曲線特征點開關函數值有小的擾動，會影響開關序列預設，但對整體求解趨勢無影響；本算例中的最終大推力段位置均位于所在軌道的遠地點或近地點，且各推進段時長相差不大，可應用一個NLP4完成脈沖軌道求解；NLP4初值設定時，式（29）中的Kv均取為1。脈沖作用時刻和脈沖大小如式（49）所示。

整體延拓過程中只出現了圖3和圖4所示的特征點變化情形，S曲線形狀不變，只是呈上移趨勢，圖14則給出了拓過程中Fmax下的若干S曲線。

本算例中，常值推力燃耗最優軌道需將Isp作為延拓參數進行多個NLP5－1的延拓求解，以推力為586．23N為例，設定Isp＝603 000s－1求解NLP5－1，每個 NLP5－1求解后，令Isp＝ηIsp求解下一個NLP5－1，若順利求解則繼續降低Isp，若不能順利求解，則令η＝ （1 ＋ η）／2和Isp＝Ispη構造NLP以求解，η初值為0．5，延拓結果如圖15所示，隨著Isp的減小，質量消耗逐漸增大，總推進時長逐漸減小，且隨著Isp的減小，η越來越大，即延拓步長越來越小。最后一個NLP5－1優化解對應總推進時長為21 313．8s，總燃耗為424．71kg，以其為初值對NLP5進行求解，優化解總推進時長為21 312．5s。

圖16為圖14中各推力下軌道的春分點根數變化情況。此算例的主要變化量為h和k，軌道變化規律為首先增加軌道半長軸和偏心率，然后改變軌道傾角和升交點赤經，而后減小軌道半長軸和偏心率，故p、f和g均有期望值以外的變化。

篇幅所限文中不再給俯仰角和偏航角的變化的示意圖，但同前兩個算例一樣，隨Fmax的增加，推力逐漸集中于對軌道根數改變效率最高的位置附近；且明顯可看出前4個軌道周期的推力角度基本一致，后4個軌道周期的推力角度基本一致，圖16中這幾個軌道周期各軌道根數變化規律也基本一致。

5 結 論

1）延拓方法是求解復雜非線性問題的一種有效數值方法，延拓方法應用的重點和難點在于延拓方案的選擇：延拓方案應從容易求解的問題，如有成熟的計算方法的雙脈沖軌道，出發；延拓過程中單個NLP問題收斂域越大，延拓過程越容易進行、整體延拓時間越短；需要特別說明的是，在實際應用中應優先選取較為簡單快速的延拓策略，若難以求解再考慮更加復雜的延拓過程，如添加質量變化模型部分的設置。為嚴謹起見用延拓方法進行一類問題的求解，可嘗試從理論上論證解的存在性，譬如文獻［29，35］的分析。對延拓結果的分析有助于對問題的進一步認知。

2）延拓方法的不足之處在于計算量較大，如本文所述基于混合法的延拓過程，不論轉移軌道有多少整體圈數，每一圈軌道從全程推進到脈沖延拓過程中都會經歷若干次開關序列變化，且每次開關序列變化都至少需要求解2～3個NLP問題，從而隨著轉移軌道圈數的增加，需計算的NLP問題的數目有大幅增加，而本文所給bangbang控制軌道若用直接法求解也可求出性質相差不大的軌道，文章所給脈沖軌道若用其他方式求解也可能更為快捷。

3）文中給出的開關序列預測方法，可用于引言中總結的各種延拓參數對應的延拓過程，有望求解更多類型更為復雜的最優轉移軌道。本文是試圖用推力幅值延拓求解多圈問題的第一步，就算法本身而言，下一步的工作是提升算法的速度和能力，如引入梯度信息、提升開關序列和優化參數的預測精度、引入有效的分叉點驗證和處理策略等。參 考 文 獻

［1］ LAWDEN D F．Optimal trajectories for space navigation［M］．London：Butterworths，1963：54－69．

［2］ HANDELSMAN M，LION P M．Primer vector on fixedtime impulsive trajectories［J］．AIAA Journal，1967，6（1）：127－132．

［3］ JEZEWSKI D J，ROZENDAAL H L．An efficient method for calculating optimal free－space n－impulse trajectories［J］．AIAA Journal，1968，6（11）：2160－2165．

［4］ HULL D G．Conversion of optimal control problems into parameter optimization problems［J］．Journal of Guidance，Control，and Dynamics，1997，20（1）：57－60．

［5］ HARGRAVES C R，PARIS S W．Direct trajectory optimization using nonlinear programming and collocation［J］．Journal of Guidance，Control，and Dynamics，1987，10（4）：338－342．

［6］ FAHROO F，ROSS I M．Direct trajectory optimization by a Chebyshev pseudospectral method［J］．Journal of Guidance，Control，and Dynamics，2002，25（1）：160－166．

［7］ GAO Y，KLUEVER C A．Low－thrust interplanetary orbit transfers using hybrid trajectory optimization method with multiple shooting［C］／Proceedings of AIAA／AAS Astrodynamics Specialist Conference．Reston：AIAA，2004：726－747．

［8］ KLUEVER C A，PIERSON B L．Optimal earth－moon trajectories using nuclear electric propulsion［J］．Journal of Guidance，Control，and Dynamics，1997，20（2）：239－245．

［9］ BETTS J T．Survey of numerical methods for trajectory optimization［J］．Journal of Guidance，Control，and Dynamics，1998，21（2）：193－207．

［10］ 高揚．電火箭星際航行：技術進展、軌道設計與綜合優化［J］．力學學報，2011，43（6）：991－1019．GAO Y．Interplanetary travel with electric propulsion：Technological progress，trajectory design，and comprehensive optimization［J］．Chinese Journal of Theoretical and Applied Mechanics，2011，43（6）：991－1019（in Chinese）．

［11］ 李俊峰，蔣方華．連續小推力航天器的深空探測軌道優化方法綜述［J］．力學與實踐，2011，33（3）：1－6．LI J F，JIANG F H．Survey of low－thrust trajectory optimization methods for deep space exploration［J］．Mechanics in Engineering，2011，33（3）：1－6（in Chinese）．

［12］ CONWAY B A．A survey of methods available for the numerical optimization of continuous dynamic systems［J］．Journal of Optimization Theory and Applications，2012，152（2）：271－306．

［13］ DIXON L C W，BARTHOLOMEW－BIGGS M C．Adjoint－control transformations for solving practical optimal control problems［J］．Optimal Control Applications and Methods，1981，2：365－381．

［14］ KLUEVER C A．Optimal earth－moon trajectories using combined chemical－electric propulsion［J］． Journal of Guidance，Control，and Dynamics，1997，20（2）：253－258．

［15］ RUSSEL R P．Primer vector theory applied to global lowthrust trade studies［J］．Journal of Guidance，Control，and Dynamics，2007，30（2）：460－472．

［16］ 何勝茂．多段拼接小推力轉移軌道優化設計方法［D］．北京：中國科學院大學，2012：26－31．HE S M．Design and optimization of low－thrust transfer trajectories with multiple patched orbital segments［D］．Beijing：University of Chinese Academy of Science，2012：26－31（in Chinese）．

［17］ FAHROO F，ROSS I M．Costate estimation by a legendre pseudospectral method［J］．Journal of Guidance，Control，and Dynamics，2001，24（2）：270－277．

［18］ GUO T，JIANG F，LI J．Homotopic approach and pseudospectral method applied jointly to low thrust trajectory optimization［J］．Acta Astronautica，2012，71（71）：38－50．

［19］ SEYWALD H，KUMAR R R．Finite difference scheme for automatic costate calculation［J］．Journal of Guidance，Control，and Dynamics，1996，19（1）：231－239．

［20］ PINES S．Constants of the motion for optimum thrust trajectories in a central force field［J］．AIAA Journal，1964，2（11）：2010－2014．

［21］ RANIERI C L，OCAMPO C A．Indirect optimization of spiral trajectories［J］．Journal of Guidance，Control，and Dynamics，2006，29（6）：1360－1366．

［22］ LEE D，BANG H．Efficient initial costates estimation for optimal spiral orbit transfer trajectories design［J］．Journal of Guidance，Control，and Dynamics，2009，32（6）：1943－1947．

［23］ THORNE J D，HALL C D．Approximate initial lagrange costates for continuous－thrust spacecraft［J］．Journal of Guidance，Control，and Dynamics，1996，19（2）：283－288．

［24］ YAN H，WU H．Initial adjoint－variable guess technique and its application in optimal orbital transfer［J］．Journal of Guidance，Control，and Dynamics，1999，22（3）：490－492．

［25］ MINTER C，FULLER－ROWELL T．A robust algorithm for solving unconstrained two－point boundary value problems（AAS 05－129）［J］．Advances in the Astronautical Sciences，2005，120（1）：409－444．

［26］ 劉滔，何兆偉，趙育善．持續推力時間最優軌道機動問題的改進魯棒算法［J］．宇航學報，2008，29（4）：1216－1221．LIU T，HE Z W，ZHAO Y S．Continuous－thrust orbit maneuver optimization using modified robust algorithm［J］．Journal of Astronautics，2008，29（4）：1216－1221（in Chinese）．

［27］ SHEN H X，CASALINO L，LI H Y．Adjoints estimation methods for impulsive moon－to－earth trajectories in the restricted three－body problem［J］．Optimal Control Applications and Methods，2014，36（4）：463－474．

［28］ SHEN H X，CASALINO L，LUO Y Z．Global search capabilities of indirect methods for impulsive transfers［J］．The Journal of the Astronautical Sciences，2015，62（3）：212－232．

［29］ HABERKORN T，MARTINON P，GERGAUD J．Low thrust minimum－fuel orbital transfer：An homotopic approach［J］．Journal of Guidance，Control，and Dynamics，2004，27（6）：1046－1060．

［30］ FOWLER W T，O’NEILL P M．Relationship between coast arc length and switching function value during optimization［J］．Journal of Spacecraft and Rockets，1976，13（7）：445－446．

［31］ ENRIGHT P J，CONWAY B A．Discrete approximations to optimal trajectories using direct transcription and nonlinear programming［J］．Journal of Guidance，Control，and Dynamics，1992，15（4）：994－1002．

［32］ GAO Y．Near－optimal very low－thrust earth－orbit transfers and guidance schemes［J］．Journal of Guidance，Control，and Dynamics，2007，30（2）：529－539．

［33］ ZUIANI F，VASILE M．Extended analytical formulas for the perturbed keplerian motion under a constant control acceleration［J］．Celestial Mechanics and Dynamical Astronomy，2015，121（3）：275－300．

［34］ BERTRAND R，EPENOY R．New smoothing techniques for solving bang－bang optimal control problems—Numerical results and statistical interpretation［J］．Optimal Control Applications and Methods，2002，23（4）：171－197．

［35］ GERGAUD J，HABERKORN T．Homotopy method for minimum consumption orbit transfer problem［J］．ESAIM：Control，Optimisation and Calculus of Variations，2006，12（2）：294－310．

［36］ JIANG F，BAOYIN H，LI J．Practical techniques for low－thrust trajectory optimization with homotopic approach［J］．Journal of Guidance，Control，and Dynamics，2012，35（1）：245－258．

［37］ ZHANG C，TOPPUTO F，BERNELLI－ZAZZERA F，et al．Low－thrust minimum－fuel optimization in the circular restricted three－body problem［J］．Journal of Guidance，Control，and Dynamics，2015，38（8）：1501－1510．

［38］ PETUKHOV V．Optimization of interplanetary trajectories for spacecraft with ideally regulated engines using the continuation method［J］．Cosmic Research，2008，46（3）：219－232．

［39］ PETUKHOV V．Method of continuation for optimization of interplanetary low－thrust trajectories［J］．Cosmic Research，2012，50（3）：249－261．

［40］ LI J，XI X N．Fuel－optimal low－thrust reconfiguration of formation－flying satellites via homotopic approach［J］．Journal of Guidance，Control，and Dynamics，2012，35（6）：1709－1717．

［41］ MITANI S，YAMAKAWA H．Continuous－thrust transfer with control magnitude and direction constraints using smoothing techniques［J］．Journal of Guidance，Control，and Dynamics，2012，36（1）：163－174．

［42］ CHUANG C H，TROY G，JOHN H．Fuel－optimal，lowand medium－thrust orbit transfers in large numbers of burns：AIAA－1994－3650［R］．Reston：AIAA，1994．

［43］ 朱小龍，劉迎春，高揚．航天器最優受控繞飛軌跡推力幅值延拓設計方法［J］．力學學報，2014，46（5）：756－769．ZHU X L，LIU Y C，GAO Y．Thrust－amplitude continuation design approach for solving spacecraft optimal controlled fly－around trajectory［J］．Chinese Journal of Theoretical and Applied Mechanics，2014，46（5）：756－769（in Chinese）．

［44］ GLANDORF D R．Lagrange multipliers and the state transition matrix for coasting arcs［J］．AIAA Journal，1969，7（2）：363－365．

［45］ FERNANDES S D S．Universal closed－form of lagrangian multipliers for coast－arcs of optimum space trajectories［J］．Journal of the Brazilian Society of Mechanical Sciences ＆Engineering，2003，25（4）：336－340．

［46］ XU Y．Enhancement in optimal multiple－burn trajectory computation by switching function analysis［J］．Journal of Spacecraft and Rockets，2006，44（1）：264－272．

［47］ JAMISON B R，COVERSTONE V．Analytical study of the primer vector and orbit transfer switching function［J］．Journal of Guidance，Control，and Dynamics，2010，33（1）：235－245．

［48］ WALKER M J H，IRELAND B，OWENS J．A set modified equinoctial orbit elements［J］．Celestial Mechanics，1985，36（4）：409－419．

［49］ LIU H，TONGUE B H．Indirect spacecraft trajectory optimization using modified equinoctial elements［J］．Journal of Guidance，Control，and Dynamics，2010，33（2）：619－623．

［50］ VOLK O．Johann heinrich lambert and the determination of orbits for planets and comets［J］．Celestial Mechanics，1980，21（2）：237－250．

［51］ GAUSS C F．Theoriy of the motion of the heavenly bodies moving about the sun in conic sections：A translation of carl frdr．Gauss“theoria motus”：With an appendix．By ch．H．Davis［M］．Boston：Little，Brown and Comp．，1857：161－233．

［52］ BATTIN R H．An introduction to the mathematics and methods of astrodynamics［M］．Reston：AIAA，1999：141－236．

［53］ VALLADO D A．Fundamentals of astrodynamics and applications［M］．New York：Springer Science ＆ Business Media，2007：419－495．

［54］ 佘明生．軌道轉移［M］／楊嘉墀，范劍峰．航天器軌道動力學與控制．北京：中國宇航出版社，2009：333－335．SHE M S．Orbital tranfer［M］／YANG J X，FAN J F．Spacecraft orbital dynamics and control．Beijing：China Astronautic Publishing House，2009：333－335 （in Chinese）．

［55］ IZZO D．Revisiting lambert’s problem［J］．Celestial Mechanics and Dynamical Astronomy，2015，121（1）：1－15．

附錄A：無推力軌道段改進春分點根數協態變量解析解

設［t0，t］為無推力軌道段，則該軌道段內改進春分點軌道根數p、f、g、h和k保持不變，真經度從L0變為L，t0和t時刻春分點根數各協態變量如式（A1）和式（A2）所示。

當開普勒軌道為圓時，Mλ為

開普勒軌道不為圓時，Mλ為

附錄B：

其中：s2＝1＋h2＋k2，s＝2hh＋2kk，w＝fcos L－fLsin L＋gsin L＋gLcos L，w＝1＋fcos L＋gsin L，上標·表示各物理量對時間的導數。

Minimum－fuel spacecraft transfer trajectories solved by direct／indirect hybrid continuation method

MENG Yazhe1，2，＊
1．Key Laboratory of Space Utilization，Technology and Engineering Center for Space Utilization，Chinese Academy of Sciences，Beijing 100094，China
2．University of Chinese Academy of Sciences，Beijing 100049，China

A continuation method is proposed for fixed－time minimum－fuel spacecraft trajectory optimization given initial and terminal states．The continuation method，based on the direct／indirect hybrid method，starts with a transfer solution with two impulses，followed by calculating firstly the full－propelling transfers and then the fuel－optimal transfers（including continuous and impulsive thrust ones）with continuation on thrust amplitude．All fuel－optimal solutions solved satisfy the necessary optimality conditions derived from the primer vector theory．A thrust switching presetting method based on the variational trends of switching function curves and a simple step－length adaptation rule based on the previous calculation results are proposed to enable the continuation automatic．The necessary optimality conditions for the impulsive trajectory expressed by the modified equinoctial orbit elements are given and verified，and the state transition matrix of costates of modified equinoctial orbit elements for coast arcs is derived．Three numerical examples are given to represent various transfer scenarios．Continuation can be regarded as improvement and extension of the direct／indirect hybrid trajectory optimization method．Continuation procedure and results can help to improve our understanding of the relations of the optimal control trajectories to system parameters．

modified equinoctial orbit elements；fuel－optimal trajectory；direct／indirect hybrid method；continuation；thrust switching sequence presetting；step－length adaptation

2016－02－26；Revised：2016－03－29；Accepted：2016－06－06；Published online：2016－06－27 13：54

URL：www．cnki．net／kcms／detail／11．1929．V．20160627．1354．006．html

s：National Natural Science Foundation of China（11372311）；Project of the Space Science Academy，Chinese Academy of

V412．4＋1

A

1000－6893（2017）01－320168－22

http：／hkxb．buaa．edu．cn hkxb＠buaa．edu．cn

10．7527／S1000－6893．2016．0183

2016－02－26；退修日期：2016－03－29；錄用日期：2016－06－06；網絡出版時間：2016－06－27 13：54

www．cnki．net／kcms／detail／11．1929．V．20160627．1354．006．html

國家自然科學基金 （11372311）；中科院空間科學研究院培育項目

＊通訊作者 ．E－mail：myz＠csu．ac．cn

孟雅哲．航天器燃耗最優軌道直接／間接混合法延拓求解［J］．航空學報，2017，38（1）：320168．MENG Y Z．Minimum－fuel spacecraft transfer trajectories solved by direct／indirect hybrid continuation method［J］．Acta Aeronautica et Astronautica Sinica，2017，38（1）：320168．

（責任編輯：張玉）

Sciences

＊Corresponding author．E－mail：myz＠csu．ac．cn