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

孟雅哲＊

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

孟雅哲1，2，＊

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

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

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

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

協(xié)態(tài)變量初值估計(jì)是間接法和混合法的難點(diǎn)之一。Dixon和 Bartholomew－Biggs［13］給出的應(yīng)用具有物理意義的控制量和協(xié)態(tài)變量相互轉(zhuǎn)換（Adjoint－Control Transformations，ACT）對(duì)協(xié)態(tài)變量初值進(jìn)行估計(jì)的方法應(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ī)動(dòng)任務(wù)簡(jiǎn)化模型并對(duì)之解析分析以獲取協(xié)態(tài)變量初值［23］；以及應(yīng)用最優(yōu)性必要條件中協(xié)態(tài)變量相關(guān)性質(zhì)給出協(xié)態(tài)變量初值［19，24］等協(xié)態(tài)變量初值估計(jì)方法。

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

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

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

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

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

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

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

1 常值推力加速度模型

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

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

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

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

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

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

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

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

應(yīng)用Pontryagin極小值原理（后簡(jiǎn)稱為極小值原理）將轉(zhuǎn)移軌道優(yōu)化問(wèn)題轉(zhuǎn)化為T(mén)PBVP，過(guò)程如下：

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

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

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

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

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

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

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

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

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

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

文獻(xiàn)［1］中的二體運(yùn)動(dòng)方程可寫(xiě)為式（2）所示的形式，其中：狀態(tài)量對(duì)應(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個(gè)分量也為其在RIRF坐標(biāo)系3個(gè)坐標(biāo)軸上的投影。用λrv表示該模型對(duì)應(yīng)的協(xié)態(tài)變量，則改進(jìn)春分點(diǎ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)春分 點(diǎn)根 數(shù) 和 位 置速度之間的Jacobi矩陣，Rrvx的計(jì)算可參考文獻(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對(duì)應(yīng)的協(xié)態(tài)變量，是文獻(xiàn)［1］中定義的主矢量?？梢?jiàn)，改進(jìn)春分點(diǎn)根數(shù)運(yùn)動(dòng)方程下與λv等價(jià)的主矢量為Φ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在且僅在脈沖施加時(shí)刻成立。

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

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

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

2 推力幅值延拓方案

2．1 延拓方案與方法

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

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

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

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

用 Lambert 問(wèn) 題［50－53］的 求 解 等 成 熟 算法［54－55］求得雙脈沖軌道后，即可計(jì)算出改進(jìn)春分點(diǎ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）求得始末速度脈沖對(duì)應(yīng)的α0和αf。由1．2節(jié)所示的主矢量必要條件，令始末脈沖作用點(diǎn)處 MTλ ＝1，由式（7）可得

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

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

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

全程推進(jìn)軌道通過(guò)表1所示的NLP求取，表中：t1∈ ［t0，tf］為一中間時(shí)刻；x （t1）和λ（t1）為對(duì)式（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)化變量中對(duì)應(yīng)協(xié)態(tài)變量的初值，對(duì)應(yīng)Fmax取值為

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

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

t1的初值滿足

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

判斷第l－1個(gè)和第l個(gè)NLP3的開(kāi)關(guān)序列是否相同，“是”則開(kāi)關(guān)切換點(diǎn)時(shí)刻初值為第l個(gè)NLP3優(yōu)化解的對(duì)應(yīng)值，F(xiàn)max用式（24）計(jì)算。“否”則轉(zhuǎn)4）。

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

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

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

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

2．4 脈沖軌道求取

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

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

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

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

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

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

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

3．1 常值推力模型

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

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

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

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

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

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

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

協(xié)態(tài)變量等比例縮放不改變?nèi)己淖顑?yōu)軌道：以x0、m0、λ0和λm0為t0時(shí)刻初值，α（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)行同樣的積分過(guò)程，可得到x1（）t、λ1（）t、S1m（）t，有

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

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

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

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

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

2）構(gòu)造 NLP5－1，將 NLP5－1的解作為初值按1）處理；若NLP5－1無(wú)法順利求解，則用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取從大到小的一系列離散值，各值對(duì)應(yīng)NLP5－1延拓求解；Isp為模型設(shè)定值的 NLP5－1的解作為初值按1）處理。

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

4 仿真驗(yàn)證

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

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

4．1 算例1

軌道轉(zhuǎn)移時(shí)間和始末改進(jìn)春分點(diǎn)軌道根數(shù)如表6所示，設(shè)t0時(shí)刻的航天器質(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)變量為無(wú)量綱單位計(jì)算結(jié)果，后面算例也是如此。

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

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

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

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

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

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

4．2 算例2

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

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

設(shè)轉(zhuǎn)移軌道含1個(gè)完整軌道周期，采用文獻(xiàn)［55］的方法求解Lambert問(wèn)題得始末脈沖如式（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)移軌道過(guò)程中所有優(yōu)化解推力作用分布燃耗隨推力幅值變化情況如圖10所示。圖10中黑色線段為延拓過(guò)程中全程推進(jìn)和所有bang－bang控制軌道解推力作用的時(shí)間分布情況，末端紅色星號(hào)標(biāo)識(shí)了脈沖作用時(shí)刻，對(duì)應(yīng)藍(lán)色線條為這些解對(duì)應(yīng)的燃耗，推進(jìn)段逐個(gè)化為脈沖時(shí)燃耗逐次減小，三脈沖轉(zhuǎn)移軌道不含連續(xù)推力作用，所對(duì)應(yīng)推力大小數(shù)據(jù)無(wú)意義。圖11繪制了若干Fmax下的開(kāi)關(guān)函數(shù)曲線以描述曲線的形狀及變化情況，表8則用具體數(shù)據(jù)對(duì)其進(jìn)行了說(shuō)明。

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

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

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

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

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

4．3 算例3

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

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

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

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

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

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

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

5 結(jié) 論

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

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

3）文中給出的開(kāi)關(guān)序列預(yù)測(cè)方法，可用于引言中總結(jié)的各種延拓參數(shù)對(duì)應(yīng)的延拓過(guò)程，有望求解更多類型更為復(fù)雜的最優(yōu)轉(zhuǎn)移軌道。本文是試圖用推力幅值延拓求解多圈問(wèn)題的第一步，就算法本身而言，下一步的工作是提升算法的速度和能力，如引入梯度信息、提升開(kāi)關(guān)序列和優(yōu)化參數(shù)的預(yù)測(cè)精度、引入有效的分叉點(diǎn)驗(yàn)證和處理策略等。參 考 文 獻(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］ 高揚(yáng)．電火箭星際航行：技術(shù)進(jìn)展、軌道設(shè)計(jì)與綜合優(yōu)化［J］．力學(xué)學(xué)報(bào)，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ù)小推力航天器的深空探測(cè)軌道優(yōu)化方法綜述［J］．力學(xué)與實(shí)踐，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è)計(jì)方法［D］．北京：中國(guó)科學(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ù)推力時(shí)間最優(yōu)軌道機(jī)動(dòng)問(wèn)題的改進(jìn)魯棒算法［J］．宇航學(xué)報(bào)，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áng)．航天器最優(yōu)受控繞飛軌跡推力幅值延拓設(shè)計(jì)方法［J］．力學(xué)學(xué)報(bào)，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］／楊嘉墀，范劍峰．航天器軌道動(dòng)力學(xué)與控制．北京：中國(guó)宇航出版社，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：無(wú)推力軌道段改進(jìn)春分點(diǎn)根數(shù)協(xié)態(tài)變量解析解

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

當(dāng)開(kāi)普勒軌道為圓時(shí)，Mλ為

開(kāi)普勒軌道不為圓時(shí)，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)·表示各物理量對(duì)時(shí)間的導(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ò)出版時(shí)間：2016－06－27 13：54

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

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

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

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