孟雅哲＊

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

孟雅哲1，2，＊

1．中國科學(xué)院空間應(yīng)用工程與技術(shù)中心 太空應(yīng)用重點實驗室，北京 100094 2．中國科學(xué)院大學(xué)，北京 100049

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

改進(jìn)春分點根數(shù)；燃耗最優(yōu)軌道；直接／間接混合法；延拓；開關(guān)序列預(yù)設(shè)；步長自適應(yīng)

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

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

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

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

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

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

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

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

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

1 常值推力加速度模型

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

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

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

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

設(shè)航天器轉(zhuǎn)移軌道的初始和末端時刻分別為t0和tf，對應(yīng)軌道狀態(tài)滿足：

性能指標(biāo)（等價于燃耗最優(yōu)）為

至此，可給出始末狀態(tài)和轉(zhuǎn)移時間固定的轉(zhuǎn)移軌道優(yōu)化問題如下：

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

應(yīng)用Pontryagin極小值原理（后簡稱為極小值原理）將轉(zhuǎn)移軌道優(yōu)化問題轉(zhuǎn)化為TPBVP，過程如下：

首先構(gòu)造如式（5）所示的Hamilton函數(shù)：

式中：λ為x對應(yīng)的協(xié)態(tài)變量，需滿足式（6）所示的微分方程。

當(dāng)F＝0N時，協(xié)態(tài)變量λ可應(yīng)用附錄A中的公式進(jìn)行解析運算。

根據(jù)極小值原理，最優(yōu)解F＊（t）和α＊（t）（t∈［t0，tf］）應(yīng)當(dāng)使H最小，所需滿足的必要條件如式（7）和式（8）所示。

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

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

協(xié)態(tài)變量的等比縮放不改變?nèi)己淖顑?yōu)軌道：以x0和λ（t0）為t0時刻初值，α（t）和F （t）按式（7）和式（8）求取，積分方程如式（2）和式（6），得到任意時刻的x（t）和λ（t）。若以x0和Kλλ（t0）（Kλ為包括1在內(nèi)的任意正實數(shù)）為初值，α（t）和F （t） 按式 （7）和式 （9）求 取，積 分 如 式 （2）和式（6），則得到任意時刻的x1（t）和λ1（t），兩次積分結(jié)果之間滿足式（10）。

為描述方便，定義式（11）所示的S為開關(guān)函數(shù)。

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

文獻(xiàn)［1］以慣性直角坐標(biāo)系下位置和速度表達(dá)的二體模型為例，給出了燃耗最優(yōu)脈沖轉(zhuǎn)移軌道主矢量必要條件，其中主矢量為速度對應(yīng)的協(xié)態(tài)變量。根據(jù)極小值原理，以位置速度表達(dá)的相對運動模型和限制性三體模型等其他航天器模型，也有以速度對應(yīng)協(xié)態(tài)變量為主矢量的脈沖最優(yōu)轉(zhuǎn)移軌道必要條件。但改進(jìn)春分點根數(shù)模型下的主矢量形式與文獻(xiàn)［1］中不同，本節(jié)以文獻(xiàn)［47］中對應(yīng)不同狀態(tài)量的協(xié)態(tài)變量之間的轉(zhuǎn)換關(guān)系出發(fā)給出該改主矢量表達(dá)式。

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

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

根 據(jù) x＝Rxrv［rTvT］T，M和Mrv有 如 下關(guān)系［47］：

式中：ΦRSW→RIRF為從RSW坐標(biāo)系到RIRF坐標(biāo)系的轉(zhuǎn)換矩陣。

由式（12）和式（13）有

式中：ΦRIRF→RSW為從RIRF坐標(biāo)系到RSW坐標(biāo)系的方向余弦矩陣，為正交矩陣；λv表示速度v對應(yīng)的協(xié)態(tài)變量，是文獻(xiàn)［1］中定義的主矢量?？梢?，改進(jìn)春分點根數(shù)運動方程下與λv等價的主矢量為ΦRSW→RIRFMTλ。相應(yīng)的燃耗最優(yōu)脈沖轉(zhuǎn)移軌道主矢量必要條件如下：

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

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

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

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

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

2 推力幅值延拓方案

2．1 延拓方案與方法

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

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

2．2 雙脈沖到全程推進(jìn)

2．2．1 雙脈沖軌道春分點根數(shù)協(xié)態(tài)變量

用 Lambert 問 題［50－53］的 求 解 等 成 熟 算法［54－55］求得雙脈沖軌道后，即可計算出改進(jìn)春分點根數(shù)的協(xié)態(tài)變量。由于所求脈沖矢量分量表達(dá)在RIRF坐標(biāo)系下，應(yīng)首先給出RIRF坐標(biāo)系到RSW坐標(biāo)系的坐標(biāo)轉(zhuǎn)換矩陣ΦRIRF→RSW為

RSW坐標(biāo)系推力方向矢量α為

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

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

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

2．2．2 NLP1和NLP2構(gòu)造說明

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

2．2．3 NLP1和NLP2優(yōu)化變量初值給定

NLP1的初值設(shè)定如下：雙脈沖軌道兩端協(xié)態(tài)變量初值為優(yōu)化變量中對應(yīng)協(xié)態(tài)變量的初值，對應(yīng)Fmax取值為

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

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

t1的初值滿足

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

2．3 全程推進(jìn)到bang－bang控制

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

2．3．1 NLP3構(gòu)造說明

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

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

式中：M 的表達(dá)式參見附錄B。

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

2．3．2 開關(guān)序列及Fmax預(yù)設(shè)

假設(shè)已成功求解l個NLP3，本節(jié)討論第l＋1個NLP3的開關(guān)序列預(yù)設(shè)和Fmax設(shè)定思路。開關(guān)序列預(yù)設(shè)通過跟蹤開關(guān)函數(shù)曲線的變化進(jìn)行，包括統(tǒng)一開關(guān)函數(shù)曲線幅值、判斷開關(guān)序列是否發(fā)生變化、預(yù)設(shè)開關(guān)點時刻3個步驟。Fmax依開關(guān)序列預(yù)設(shè)情況進(jìn)行給定。

將開關(guān)函數(shù)曲線幅值統(tǒng)一設(shè)定為10（人為設(shè)定參數(shù)，可更改）：設(shè)第n（n＝1，2，…，l）個滿足必要條件的NLP優(yōu)化解對應(yīng)開關(guān)函數(shù)曲線為Sn，用式（22）所示的Kλ對優(yōu)化解中始末協(xié)態(tài)變量進(jìn)行式（23）所示的縮放。

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

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

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

判斷是否出現(xiàn)圖2～圖4中情況的方法如下：查找當(dāng)前開關(guān)函數(shù)曲線為0的點的個數(shù)是否與開關(guān)點數(shù)目相同，不同則認(rèn)為出現(xiàn)了圖2對應(yīng)情況，相同則認(rèn)為未出現(xiàn)。

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

開關(guān)序列不變時，第l＋1個NLP3對應(yīng)Fmax依式（24）取值，所有優(yōu)化變量初值均為第l個NLP3優(yōu)化解的對應(yīng)部分。

式中：0＜δ＜1為給定常數(shù)。

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

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

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

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

2．3．3 NLP3開關(guān)切換時刻初值給定

2．3．2節(jié)中變量ε和δ的取值決定了推力幅值延拓中Fmax的變化幅度。ε和δ越小，F(xiàn)max的變化越小，開關(guān)序列和NLP初值設(shè)定越準(zhǔn)確，但延拓整體需花費時間越長。延拓算法中將根據(jù)NLP3求解情況自適應(yīng)ε和δ的值，規(guī)則可簡單描述如下：若當(dāng)前NLP求解失敗，則減小ε和δ；若NLP連續(xù)快速收斂，則適當(dāng)增大ε和δ。

設(shè)已有l(wèi)－1個轉(zhuǎn)移軌道優(yōu)化解，求解第l個NLP3后，第l＋1個NLP3的初值設(shè)定過程如下：

1）調(diào)整延拓步長參數(shù)δ和ε。

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

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

2）判斷是否出現(xiàn)圖2所示情況并進(jìn)行處理。

設(shè)S＊＝0，找到Sl上所有取值為S＊的點，若點的個數(shù)與第l個優(yōu)化解的開關(guān)切換次數(shù)不同，認(rèn)為出現(xiàn)了圖2所示情況，Sl上所有S＊＝0的點即為第l＋1個NLP3的開關(guān)序列的初值，F(xiàn)max用式（27）求解。若次數(shù)相同，則轉(zhuǎn)3）。

3）判斷第l－1次和第l次優(yōu)化解開關(guān)序列是否相同，并進(jìn)行處理。

判斷第l－1個和第l個NLP3的開關(guān)序列是否相同，“是”則開關(guān)切換點時刻初值為第l個NLP3優(yōu)化解的對應(yīng)值，F(xiàn)max用式（24）計算?！胺瘛眲t轉(zhuǎn)4）。

4）判斷是否有開關(guān)函數(shù)曲線特征點過0。

判斷距0最近的特征點開關(guān)函數(shù)值是否小于臨界值ε，“否”則認(rèn)為開關(guān)序列不發(fā)生變化，初值設(shè)定同3）；“是”則轉(zhuǎn)5）。

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

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

2．4 脈沖軌道求取

2．4．1 NLP4構(gòu)造說明

bang－bang控制軌道推進(jìn)段時長足夠小時，采用表4所示的NLP問題逐個將短時大推力段轉(zhuǎn)為脈沖軌道。表中：性能指標(biāo)Δm 為應(yīng)用式（28）和式（30）估計的總?cè)剂舷模籺2，1，t2，2，…，t2，M2為所有的開關(guān)點和脈沖點時刻；v1，v2，…，vN對應(yīng)速度脈沖大小，其中N≤M2；r、v、λr和λv為正向計算的結(jié)果，r′、v′、λ′r和λ′v為反向計算的結(jié)果，可見前幾個約束是t＊2（即t2，ceil（M2／2））處位置速度及對應(yīng)協(xié)態(tài)變量連續(xù)，而系統(tǒng)方程的處理及開關(guān)函數(shù)還是按改進(jìn)春分點根數(shù)進(jìn)行計算的；剩余約束的含義同NLP3。當(dāng)t＊2為某個脈沖作用點時，約束條件是脈沖作用后的狀態(tài)和協(xié)態(tài)變量相等。

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

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

2．4．2 NLP4優(yōu)化變量初值給定

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

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

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

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

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

3．1 常值推力模型

設(shè)m 滿足微分 方 程 式 （30），則 式 （2）和式（30）組成系統(tǒng)方程。

燃耗最優(yōu)對應(yīng)目標(biāo)函數(shù)為

根據(jù)極小值原理，構(gòu)造Hamilton函數(shù)為

式中：λm為m 對應(yīng)的協(xié)態(tài)變量，其微分方程為式（33），滿足末端條件式（34）。

使得式（32）所示 Hamilton函數(shù)達(dá)極小的F＊滿足式（35）和式（36），Sm為對應(yīng)開關(guān)函數(shù)。

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

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

協(xié)態(tài)變量等比例縮放不改變?nèi)己淖顑?yōu)軌道：以x0、m0、λ0和λm0為t0時刻初值，α（t）和F（t）依式（7）、式（35）和式（36）取值，分別按式（2）、式（30）、式（6）和式（33）進(jìn)行積分，得到x（）t、λ（）t和Sm（）t；若以x0、m0、Kλmλ0和Kλmλm0為初值進(jìn)行同樣的積分過程，可得到x1（）t、λ1（）t、S1m（）t，有

3．2 NLP5和NLP5－1構(gòu)造及初值給定

表5為求解常值推力轉(zhuǎn)移軌道的NLP問題。表中：λm（t0）和λm（tf）為始末時刻質(zhì)量的協(xié)態(tài)變量；其余變量意義同NLP3。

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

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

1）構(gòu)造NLP5進(jìn)行求解，若順利求解，對優(yōu)化解進(jìn)行式（39）所示的協(xié)態(tài)變量縮放得到最優(yōu)解；若NLP5無法順利求解，則用2）。

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

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

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

首個NLP5或NLP5－1求解時，λm（t0）和λm（tf）初值設(shè)定為－1，其余優(yōu)化變量初值設(shè)定為對應(yīng)NLP5優(yōu)化解相應(yīng)量的取值。常值推力加速度到常值推力求解過程中未考慮開關(guān)序列的變化，若存在開關(guān)序列變化，可依2．3節(jié)內(nèi)容進(jìn)行開關(guān)序列預(yù)設(shè)。

4 仿真驗證

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

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

4．1 算例1

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

應(yīng)用共面共拱線橢圓軌道的雙脈沖最優(yōu)轉(zhuǎn)移［54］求得速度脈沖如式（42）所示，始末速度脈沖均為RIRF坐標(biāo)系下的表達(dá)，下文不再贅述。

應(yīng)用式（18）和式（A3）求得雙脈沖轉(zhuǎn)移軌道協(xié)態(tài)變量初始值如式（43），協(xié)態(tài)變量為無量綱單位計算結(jié)果，后面算例也是如此。

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

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

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

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

需要說明的是，本算例為共面轉(zhuǎn)移，可忽略h和k對應(yīng)的協(xié)態(tài)量及其微分方程，并按照h＝k＝0和α（）3＝0化簡式（2）和式（6）中的剩余矩陣元素，按照文中給出的延拓過程進(jìn)行求解，求得所需燃耗最優(yōu)軌道后，恢復(fù)所有積分?jǐn)?shù)據(jù)中h和k的值，即可得三維模型下的燃耗最優(yōu)軌道。

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

4．2 算例2

軌道轉(zhuǎn)移時間和始末改進(jìn)春分點軌道根數(shù)如表7所示，m （t0）、ge和Isp同算例1。

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

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

雙脈沖軌道協(xié)態(tài)變量初值為

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

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

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

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

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

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

4．3 算例3

軌道轉(zhuǎn)移時間和始末改進(jìn)春分點軌道根數(shù)如表9所示，m （t0）、ge和Isp同算例1。

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

設(shè)轉(zhuǎn)移軌道有20個完整軌道周期，Lambert解［55］對應(yīng)始末脈沖解如式（47）所示。

雙脈沖軌道始末協(xié)態(tài)變量如式（48）所示。

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

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

本算例中，常值推力燃耗最優(yōu)軌道需將Isp作為延拓參數(shù)進(jìn)行多個NLP5－1的延拓求解，以推力為586．23N為例，設(shè)定Isp＝603 000s－1求解NLP5－1，每個 NLP5－1求解后，令I(lǐng)sp＝ηIsp求解下一個NLP5－1，若順利求解則繼續(xù)降低Isp，若不能順利求解，則令η＝ （1 ＋ η）／2和Isp＝Ispη構(gòu)造NLP以求解，η初值為0．5，延拓結(jié)果如圖15所示，隨著Isp的減小，質(zhì)量消耗逐漸增大，總推進(jìn)時長逐漸減小，且隨著Isp的減小，η越來越大，即延拓步長越來越小。最后一個NLP5－1優(yōu)化解對應(yīng)總推進(jìn)時長為21 313．8s，總?cè)己臑?24．71kg，以其為初值對NLP5進(jìn)行求解，優(yōu)化解總推進(jìn)時長為21 312．5s。

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

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

5 結(jié) 論

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

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

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

［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］ 高揚．電火箭星際航行：技術(shù)進(jìn)展、軌道設(shè)計與綜合優(yōu)化［J］．力學(xué)學(xué)報，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］ 李俊峰，蔣方華．連續(xù)小推力航天器的深空探測軌道優(yōu)化方法綜述［J］．力學(xué)與實踐，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］ 何勝茂．多段拼接小推力轉(zhuǎn)移軌道優(yōu)化設(shè)計方法［D］．北京：中國科學(xué)院大學(xué)，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，F(xiàn)ULLER－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］ 劉滔，何兆偉，趙育善．持續(xù)推力時間最優(yōu)軌道機(jī)動問題的改進(jìn)魯棒算法［J］．宇航學(xué)報，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］ 朱小龍，劉迎春，高揚．航天器最優(yōu)受控繞飛軌跡推力幅值延拓設(shè)計方法［J］．力學(xué)學(xué)報，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］ 佘明生．軌道轉(zhuǎn)移［M］／楊嘉墀，范劍峰．航天器軌道動力學(xué)與控制．北京：中國宇航出版社，2009：333－335．SHE M S．Orbital tranfer［M］／YANG J X，F(xiàn)AN 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：無推力軌道段改進(jìn)春分點根數(shù)協(xié)態(tài)變量解析解

設(shè)［t0，t］為無推力軌道段，則該軌道段內(nèi)改進(jìn)春分點軌道根數(shù)p、f、g、h和k保持不變，真經(jīng)度從L0變?yōu)長，t0和t時刻春分點根數(shù)各協(xié)態(tài)變量如式（A1）和式（A2）所示。

當(dāng)開普勒軌道為圓時，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，上標(biāo)·表示各物理量對時間的導(dǎo)數(shù)。

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；網(wǎng)絡(luò)出版時間：2016－06－27 13：54

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

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

＊通訊作者 ．E－mail：myz＠csu．a(chǎn)c．cn

孟雅哲．航天器燃耗最優(yōu)軌道直接／間接混合法延拓求解［J］．航空學(xué)報，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．

（責(zé)任編輯：張玉）

Sciences

＊Corresponding author．E－mail：myz＠csu．a(chǎn)c．cn