基于自動微分與偽譜法的小推力借力軌道設(shè)計

張 勃，泮斌峰，唐 碩

(西北工業(yè)大學 航天學院，710072西安)

小推力發(fā)動機具有比沖高、質(zhì)量輕的特點，在深空探測任務中，采用小推力發(fā)動機能夠有效提高探測器有效載荷的比重.在可以預見的未來，小推力發(fā)動機將成為深空探測任務中不可或缺的動力系統(tǒng).但是，由于小推力發(fā)動機產(chǎn)生的是很小的連續(xù)推力，并且工作時間很長，使得小推力轉(zhuǎn)移軌道會呈現(xiàn)出非常強的非線性［1］，導致傳統(tǒng)基于大推力脈沖變軌的軌道設(shè)計方法已不能滿足需求.而且，由于深空探測器轉(zhuǎn)移軌道必須滿足發(fā)射時間、飛行時間以及星歷等的約束，形成了一個多約束、高度非線性的最優(yōu)控制問題，使得小推力問題的求解異常復雜.

借力飛行技術(shù)也叫引力甩擺或引力輔助軌道轉(zhuǎn)移技術(shù)［2］.探測器在進入借力天體的影響球和飛出借力天體的影響球時，其相對借力天體的速度大小不變，僅改變相對速度的方向.但是探測器飛入和飛出行星影響球時相對太陽的速度大小和方向均發(fā)生了變化，從而導致探測器脫離借力天體影響球后日心軌道發(fā)生改變，達到變軌的目的.借力飛行技術(shù)不僅可以節(jié)省燃料，而且可以實現(xiàn)對借力天體的觀測，提高探測器的利用率.

近年來，由于人類深空探測活動的需要，小推力結(jié)合借力飛行的軌道設(shè)計方法受到了廣泛關(guān)注.文獻［3］用間接優(yōu)化的方法對小推力借力軌道的優(yōu)化設(shè)計進行了研究，首先用脈沖軌道轉(zhuǎn)移的方法得到借力序列和切換結(jié)構(gòu)，然后從后向前逐段優(yōu)化，再將各段拼接起來.為了簡化計算，Casalino去掉了借力高度的約束，這樣得到的結(jié)果不一定能夠滿足實際任務約束.由于基于龐特里亞金極大值原理的間接優(yōu)化算法求解過于復雜，限制了它的應用.文獻［4］提出的shape based方法采用正弦指數(shù)函數(shù)來擬合軌道的形狀，通過不斷調(diào)整相關(guān)參數(shù)，進行初值搜索，這種方法提供了一種有效的初值搜索手段.文獻［5-7］用shape based方法對小推力借力軌道的優(yōu)化設(shè)計進行了研究.但是該方法在調(diào)整函數(shù)中的參數(shù)時很困難，需要清晰理解相關(guān)參數(shù)的物理意義，而且計算規(guī)模大，計算成本高［8-9］.文獻［8］根據(jù)形狀法和偽譜法，提出一種混合優(yōu)化策略對小推力借力軌道進行了研究.文獻［9］采用微分進化、擬退火算法與序列二次規(guī)劃相結(jié)合的方法研究了小推力借力軌道的優(yōu)化問題.文獻［10-12］用隨機搜索算法與確定性算法相結(jié)合的混合方法對小推力借力軌道進行了研究，這些方法中隨機搜索算法的收斂速度慢、精度低，影響了算法的性能.文獻［13］提出了一種基于形狀逼近策略的小推力借力軌道的初始設(shè)計方法，但并沒有進行更為詳盡的優(yōu)化設(shè)計.

本文采用B平面借力模型，建立了小推力借力轉(zhuǎn)移軌道的數(shù)學模型.采用高斯偽譜法對整個軌道進行離散，把發(fā)射時間和交會條件作為端點約束，借力高度與B平面角作為借力的控制量，并約束在一定范圍內(nèi)，建立起完整的整個轉(zhuǎn)移軌道的離散模型.用SQP方法對離散后的非線性規(guī)劃問題進行求解.為了提高導數(shù)的計算精度，從而加快收斂速度，求解過程中的導數(shù)信息通過自動微分獲得.由于小推力借力軌道的空間和時間跨度很大，離散后的NLP問題規(guī)模非常大，為了使計算快速收斂，本文提出了串行優(yōu)化和彈性約束的策略.本文方法將發(fā)射時間的搜索和飛行軌跡的優(yōu)化包含在同一計算框架內(nèi)，避免了傳統(tǒng)混合算法中進化算法收斂速度慢的缺點.

1 模型建立

1.1 基本假設(shè)

探測器在太陽系內(nèi)運動時，受到太陽和各個行星引力的作用，是一個多體問題.在對探測器軌道進行精確設(shè)計之前，一般忽略次要因素，以二體問題為基礎(chǔ)設(shè)計軌道的初步方案.同樣，探測器進行借力的過程也是一個復雜的過程，在初步方案設(shè)計階段予以簡化.本文基于如下簡化假設(shè).

1)探測器飛離地球影響球，飛向目標星體的星際航行段僅考慮太陽作為中心天體的引力作用，此段軌道通常為橢圓軌道.探測器進入借力天體影響球內(nèi)時僅考慮借力天體為中心引力體，其飛行軌跡相對于借力天體為雙曲線軌道.

2)由于探測器日心轉(zhuǎn)移飛行時間較長，借力時間可以忽略不計(以金星為例，高度10 000 km的圓軌道周期大約為6.24 h，而整個日心轉(zhuǎn)移時間至少需要上百天)，即將借力飛行看作是瞬時脈沖(無需消耗工質(zhì))，借力飛行前后的探測器日心位置沒有變化.

1.2 動力學模型

在日心黃道慣性系中，探測器的軌道動力學方程為

其中:r為探測器的位置矢量;v為速度矢量;μ為太陽引力常數(shù);T為發(fā)動機的推力大小;m為探測器質(zhì)量;α為發(fā)動機推力的單位方向矢量;g0為地球重力加速度;Ⅰsp為發(fā)動機比沖．為提高計算精度和收斂速度，需要用正則天文單位對動力學方程進行無量綱化，無量綱化后的動力學方程與方程(1)的形式完全相同．

1.3 發(fā)動機模型

本文中的探測器采用太陽能電離子發(fā)動機.假定太陽能帆板受照面積為常值，則發(fā)動機的輸入功率P隨著探測器日心距離的增大而減小，它與探測器日心距的平方成反比，即

其中P0為探測器日心距離為1 AU時的輸入功率，r采用天文單位，發(fā)動機推力大小T為

其中η為發(fā)動機的工作效率．

1.4 行星借力軌道模型

探測器的行星借力飛行使得整個動力學過程的非線性更強，解空間變得更為復雜.選擇合適的借力軌道模型不僅可以簡化計算，而且可以獲得較高的精度.本文采用B平面模型［14］對借力軌道進行分析，引入了B平面角和借力高度作為控制量，降低了軌道對控制變量的敏感度，提高了軌跡優(yōu)化設(shè)計的魯棒性.

行星借力B平面示意如圖1所示．B平面垂直于探測器飛入行星引力影響球時的雙曲線超速矢量v-∞，并且經(jīng)過飛越行星的引力中心．定義B矢量為B平面與探測器軌道面的交線，由行星引力中心指向探測器方向，記為B;T矢量定義為B平面與日心黃道坐標系的交線，由行星引力中心指向探測器方向，記為T．定義B平面角γ為矢量B與矢量T的夾角．

圖1 行星借力B平面示意

根據(jù)借力飛行的原理，借力前后探測器相對于飛越行星的相對速度大小不變，只改變相對速度矢量的指向，相對速度矢量的方向變化角度由式(6)給出:

其中rp為飛越半徑，μp為飛越行星的引力常數(shù)．

根據(jù)B平面的幾何關(guān)系，可以得到探測器經(jīng)行星借力后飛出行星影響球的雙曲線超速矢量為

其中

根據(jù)上述借力模型，對于給定的借力飛越半徑rp和B平面角γ，就可以由探測器進入行星引力影響球時的雙曲線超速矢量v-∞計算得到探測器飛出行星引力影響球的雙曲線超速矢量v+∞．在借力軌道的設(shè)計與優(yōu)化過程中可以通過調(diào)節(jié)飛越半徑rp和B平面角γ來控制v+∞．

2 軌道優(yōu)化問題描述

2.1 性能指標函數(shù)及約束條件

根據(jù)深空探測任務的性質(zhì)，燃料越省，就意味著更長的使用壽命和更高的有效載荷.因此性能指標取為

其中:t0為從地球發(fā)射的時間;tf為到達目標星體的時間;m0為探測器的初始質(zhì)量．

探測器的發(fā)射時間不定，通常被要求在某一時間范圍內(nèi)完成，而探測器到達目標星體的時間亦不定．但是探測器的發(fā)射時間和到達時間必須滿足星歷約束．對于任何發(fā)射時間t0，探測器必須滿足地球星歷約束，即

其中:r，v是發(fā)射時刻探測器在日心黃道坐標系中的位置向量和速度向量;rE，vE是此刻地球的位置向量和速度向量;v∞是地球逃逸雙曲線速度．本文中v∞=0，即探測器以雙曲線超速為零逃逸出地球．

探測器飛越借力星體時，探測器的位置必須滿足借力星體的星歷約束，即

其中tm為借力時刻，rF為借力行星的位置矢量．

當探測器與目標交會時，探測器的位置與速度必須滿足

其中rM，vM分別為目標星體的位置和速度矢量．本文中的星歷均采用JPL實驗室提供的DE405行星星歷［15］.

2.2 最優(yōu)控制問題

為了便于最優(yōu)控制問題的描述，令狀態(tài)變量x(t)=[rvm]T，x(t)∈R7，推力控制量u(t)=α，u(t)∈R3．由式(1)可以得動力學方程的形式為

目標函數(shù)(9)記為

在飛行過程中，狀態(tài)變量必須處于合理的取值區(qū)間內(nèi)，即

推力的控制變量為單位方向矢量，因此必須滿足

探測器發(fā)射時刻必須滿足地球星歷約束(10)，記為

借力控制量飛躍半徑rp和B平面角γ的取值必須滿足

借力時刻tm，探測器的速度必須滿足借力前后的關(guān)系(7)，記為

探測器的位置矢量必須滿足借力行星星歷約束(11)，記為

到達目標星體時必須滿足交會要求(12)，記為

最后，總的飛行時間必須限制在可接受的范圍內(nèi)

小推力借力軌道的優(yōu)化設(shè)計問題可以描述為:尋找推力控制量u*(t)∈R3，借力控制量和γ*，使得目標函數(shù)J最小，同時滿足動力學方程(13)以及約束(15)～(22)．可見小推力借力軌道的優(yōu)化設(shè)計問題是一個高度非線性，多邊值約束，具有內(nèi)點約束并且不連續(xù)(速度存在突變)的最優(yōu)控制問題，導致問題的求解異常復雜．

2.3 俯仰角和偏航角的計算

在前面的動力學模型中，發(fā)動機推力的方向矢量定義在日心黃道坐標系中，不便于軌道機動的實施．實際上，發(fā)動機推力的方向由探測器軌道坐標系中的兩個角度來確定，即俯仰角和偏航角．俯仰角α定義為推力矢量在探測器軌道平面內(nèi)的投影與當?shù)厮矫娴膴A角，偏航角β定義為推力矢量與探測器軌道平面之間的夾角．俯仰角和偏航角的大小分別由下式給出:

3 高斯偽譜法基本原理及NLP問題求解策略

3.1 高斯偽譜法基本原理

偽譜法是近年提出的一種求解最優(yōu)控制問題的直接方法，具有收斂半徑大、收斂速度快、初值不敏感等特性，在航天航空領(lǐng)域得到了廣泛的應用［16-18］.高斯偽譜法在勒讓德多項式的零點(即Legendre Gauss點)上對狀態(tài)變量和控制變量進行離散，以離散點上的狀態(tài)變量和控制變量構(gòu)造拉格朗日多項式對連續(xù)狀態(tài)變量和連續(xù)控制變量進行逼近，通過對全局差值多項式求導來近似狀態(tài)變量的導數(shù)，從而將微分方程約束轉(zhuǎn)化為一組代數(shù)約束.性能指標中的積分項和終端狀態(tài)由高斯積分得到.經(jīng)過上述變換，就可將最優(yōu)控制問題轉(zhuǎn)化為具有一系列代數(shù)約束的非線性規(guī)劃(NLP)問題.

3.2 NLP問題求解策略

在高精度的優(yōu)化設(shè)計中，離散點的數(shù)量通常很大，導致求解變量數(shù)量也隨之劇增，若直接對大規(guī)模NLP問題進行求解可能會導致收斂速度非常慢，甚至不能收斂.針對這一問題，本文提出串行優(yōu)化和彈性約束的求解策略(如圖2所示).串行優(yōu)化有兩層含義:1)將優(yōu)化問題分成若干子段，分段進行優(yōu)化，以前一段的結(jié)果作為后一段的初值，求出每一段各自部分最優(yōu)解，再將各段優(yōu)化結(jié)果結(jié)合起來得到次優(yōu)解，并作為原優(yōu)化問題的初始猜想，最終求得整體最優(yōu)解;2)先以較少的離散點來對軌道進行參數(shù)化，求得優(yōu)化結(jié)果;對優(yōu)化結(jié)果進行插值，將優(yōu)化結(jié)果映射到新的計算網(wǎng)格上作為初始猜想，然后逐步增加離散點，提高計算精度.彈性約束指先忽略(或者放寬)某些約束，使得NLP問題能夠快速收斂，當離散點逐步增多時增加(或收緊)被忽略(或放寬)的約束，這樣逐級進行計算，不斷提高精度，最終求得滿足約束的精確解.

圖2 兩種不同意義下的串聯(lián)優(yōu)化策略

在進行計算時，不僅兩種意義下的串行優(yōu)化策略可以互相包含，交叉使用，而且串行優(yōu)化策略可以與彈性約束策略相互耦合.靈活地運用串行優(yōu)化和彈性約束策略不僅可以增強計算的魯棒性，而且可以顯著加快收斂速度.

3.3 自動微分

采用SQP方法對得到的NLP問題進行求解時，導數(shù)的計算精度關(guān)系到最終結(jié)果的精度以及NLP問題收斂的快慢，導數(shù)計算精度不高甚至會導致NLP問題的求解過程不能收斂.常用的導數(shù)計算方法一般有解析導數(shù)、有限差分和自動微分3種.解析導數(shù)只適用于導數(shù)所依賴的自變量關(guān)系較為簡單的情況，當導數(shù)所依賴的自變量關(guān)系較為復雜時則很難得到解析導數(shù)，因此解析導數(shù)的適用范圍有限.有限差分法是普遍采用的一種數(shù)值計算導數(shù)的方法，分為前向差分、后向差分以及中心差分，有限差分實現(xiàn)簡單，但計算精度不高.自動微分是機械地運用鏈式求導法則對計算機程序形式的函數(shù)求導的一組技術(shù)［19］.自動微分的基本思想是無論描述函數(shù)的程序多么復雜，其本質(zhì)都是執(zhí)行一系列的元代數(shù)運算或元函數(shù)運算.對這些初等運算迭代地運用鏈式求導法則，就可以自動、精確地得到目標函數(shù)的任意階導數(shù)［20］.

根據(jù)鏈式求導法則累加的方式不同，自動微分分為正向模式和逆向模式.如果從中間變量到獨立變量按照程序的執(zhí)行順序進行求導稱為正向模式;反之，如果從依賴變量到中間變量按照與程序執(zhí)行順序相反的方向進行求導稱為逆向模式.由于向量函數(shù)的求導最終需要轉(zhuǎn)換為標量函數(shù)的求導，因此下面以標量函數(shù)的求導為例來說明正向模式和逆向模式的實現(xiàn)方式.

考慮函數(shù)y=f(x)，f:Rn→R，可以由下面的程序?qū)崿F(xiàn):

其中函數(shù)fi依賴于已經(jīng)計算出的xj的值，j∈Ji，Ji?{1，2，…，i-1}，i=1，2，…，m，即f由m-n個基本初等函數(shù)fi復合而成．基本初等函數(shù)fi的導數(shù)▽fi=(?fi/?xj)j∈Ji可以精確求得．由鏈式求導法則可以得到正向自動微分的實現(xiàn)方式:

其中ei為R空間中第i個笛卡爾基向量．

為了描述自動微分的逆向模式，記與中間變量xi相關(guān)的導數(shù)為≡?xm/?xi．由定義可以得到1，?f(x)/?xi=i，i=1，2，…，n．由鏈式求導法則可以得到其中Tj≡{i≤m:j∈Ji}，進而得到逆向模式的實現(xiàn)方式:

其中(·)ni=1表示i從1增長到n．當取初始向量=0，γ=1時，計算所得向量g即為梯度▽f．

為提高計算精度和收斂速度，本文采用自動微分計算約束和目標函數(shù)的一階導數(shù)，以此提供給SQP方法進行迭代運算.假設(shè)離散后的動力學方程在第k個離散點處的狀態(tài)為Xk，約束向量為Ck，則采用前向自動微分求解約束Ck對Xk的一階導數(shù)可以表示為其中ckm為離散約束向量Ck的元素，xij為離散狀態(tài)向量Xk的元素．均為基本初等函數(shù)的求導．無論ckm的形式多么復雜，經(jīng)過有限步總能分解為初等函數(shù)的導數(shù)相乘的形式，而初等函數(shù)的導數(shù)都具有解析表達式，可以精確求解．把等的導數(shù)值相乘即可得到導數(shù)的值，進而求得的值．針對配點法(偽譜法屬于全局配點法)得到的非線性規(guī)劃問題，離散點處的約束通常只與當前點的離散狀態(tài)變量以及離散控制變量相關(guān)，因此當前約束對非當前點的離散量的導數(shù)為0．所以通過預先的判斷，可以避免不必要的微分計算，減小計算量．

目前運用算子重載方法的自動微分軟件的典型 代 表 有ADOL-C、ADF、FADBAD/TADIFF、CppAD、ADC和MAD等，基于代碼轉(zhuǎn)化方法實現(xiàn)的典型軟件包有ADIC、TAPENADE、ADIFOR和OpenAD等.本文采用MAD軟件包作為自動微分工具，它是TOMLAB軟件中的一個MATLAB工具箱，采用前向自動微分模式.

4 仿真計算與分析

本文以地球-金星-火星交會任務的燃料最省小推力借力轉(zhuǎn)移軌道進行設(shè)計.探測器在2021年1月1日到2021年12月31日之間從地球發(fā)射，假設(shè)出發(fā)時刻探測器的日心位置和速度與地球相同(即探測器的地球逃逸速度為零);飛越金星進行借力，借力高度范圍為200～10 000 km，最終與火星交會，總飛行時間約束在500 d到1 500 d之間.行星借力模型采用B平面模型，采用高斯偽譜法分別對地球到金星和金星到火星的軌道進行離散化.由于整個軌道的空間和時間跨度很大，為了得到較高的精度就必須采用足夠多的配點.整個軌道采用的配點數(shù)N=200(地球到金星、金星到火星分別100個配點)，使用基于MTALAB的TOMLAB對NLP問題進行編程求解，導數(shù)信息通過MAD用自動微分法進行計算.探測器的主要參數(shù)見表1.

表1 探測器主要參數(shù)

盡管高斯偽譜法具有較強的魯棒性，但是由于小推力借力軌道的優(yōu)化設(shè)計非常復雜，如果直接對大規(guī)模的NLP問題進行求解會導致收斂速度比較慢，甚至不能收斂，因此本文采用前面所述的串行優(yōu)化和彈性約束策略.首先將軌跡優(yōu)化問題分解為地球到金星和金星到火星的小推力軌道轉(zhuǎn)移問題，求出次優(yōu)解;然后以次優(yōu)解作為初始猜想，分別采用不同的高斯點進行串聯(lián)求解，即以較少離散點的優(yōu)化結(jié)果作為下一步優(yōu)化的初始猜想，逐級進行優(yōu)化.表2給出了不同串聯(lián)級數(shù)的優(yōu)化耗時，當配點數(shù)N=20、40、60時不考慮借力高度約束，N=80、100、200時增加對借力高度的約束，使NLP問題能夠很快收斂.

可以看到，隨著串聯(lián)級數(shù)的增加，總的優(yōu)化耗時逐漸減少，但是在一定串聯(lián)級數(shù)之后優(yōu)化耗時隨串聯(lián)級數(shù)的增加不再明顯減少，這是由于一定串聯(lián)級數(shù)之后優(yōu)化結(jié)果已經(jīng)接近最終優(yōu)化結(jié)果.N=200時離散變量個數(shù)達到了2 019個，約束為1 618個，對于如此大規(guī)模的NLP問題，通過串聯(lián)優(yōu)化仍然能夠快速收斂，精確得到了滿足上述約束條件的小推力借力軌道，證明了本文求解策略的有效性.主要優(yōu)化結(jié)果見表3.

表2 不同串聯(lián)優(yōu)化級數(shù)求解耗時對比

表3 小推力借力軌道優(yōu)化結(jié)果

導數(shù)信息的計算精度對依賴于導數(shù)信息的NLP問題求解方法的收斂速度有很大影響，精度越高，收斂速度越快.為了排除實現(xiàn)語言和實現(xiàn)方式對計算時間的影響，以約束函數(shù)(目標函數(shù)過于簡單，無比較意義)的計算次數(shù)和約束雅可比矩陣的計算次數(shù)作為指標對自動微分方法和差分方法對優(yōu)化過程的影響進行了比對，如表4所示，結(jié)果顯示自動微分方法能夠明顯加快優(yōu)化的收斂速度.

表4 自動微分與差分方法計算導數(shù)對優(yōu)化過程的影響對比

圖3給出了探測器的整個飛行軌跡，可以看到優(yōu)化結(jié)果很好地滿足了借力和交會條件，驗證了本文算法的正確性.圖4為飛行軌跡在x-y平面的投影，清晰地反映了軌跡與各行星之間的關(guān)系.

圖3 探測器三維軌跡

圖4 探測器軌跡在x-y平面的投影

圖5為探測器日心距和速度的變化歷程，其中探測器的速度在進行借力后存在突變，這是由于B平面模型沒有考慮借力的過程，認為借力是在某一點瞬時完成的．由于B平面模型并不影響探測器的位置，因此探測器的日心距不存在突變．

圖5 探測器日心距與速度時間歷程

圖6給出了控制角的變化歷程.同樣也是由于B平面模型的原因，俯仰角與偏航角都存在突變.借力前，俯仰角單調(diào)增加，但始終為負值，這是因為金星軌道半徑小于地球軌道半徑，探測器必須不斷減小日心距.而借力后由于火星軌道半徑大于金星軌道半徑，探測器首先需要增加日心距，當達到火星軌道附近時為了滿足與火星交會的要求，探測器需再對控制角進行調(diào)整.偏航角的變化幅值較小，這是因為金星、地球和火星軌道的軌道傾角相差不大，軌道面之間的調(diào)整所需控制力較小.

圖6 控制角時間歷程

5 結(jié)語

本文基于自動微分和高斯偽譜法研究了小推力借力軌道的優(yōu)化設(shè)計問題，將發(fā)射時間的搜索與軌道的設(shè)計包含在了同一計算框架內(nèi)，避免傳統(tǒng)混合方法中隨機搜索算法收斂速度慢、精度低的缺點.為了增強算法的魯棒性、加快收斂速度，本文提出了串聯(lián)優(yōu)化策略和彈性約束策略.針對有限差分法求解導數(shù)信息精度有限而解析導數(shù)又難以獲得的問題，本文采用自動微分的方法求解導數(shù)信息，保證了導數(shù)計算的高精度和快速收斂.通過地球-金星-火星小推力借力軌道的仿真驗證，證明了本文算法的正確性和有效性，為未來深空探測軌道優(yōu)化提供了有益的參考.

［1］BETTS J T.Very low-thrust trajectory optimization using a direct SQP method［J］.Journal of Computational and Applied Mathematics，2000，120(1):27-40.

［2］張旭輝，劉竹生.火星探測無動力借力飛行軌道研究［J］.宇航學報，2008，29(6):1739-1746.

［3］CASALINO L，COLASURDO G，SENTINELLA M R.Low-thrust trajectories to mercury with multiple gravity assists［C］//43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference＆Exhibit.Cincinnati，OH:［s.n.］，2007.

［4］PETROPOULOS A E.A shape-based approach to automated，low-thrust，gravity-assist trajectory design［D］.West Lafayette:Purdue University，2001.

［5］MCCONAGHY T T，DEBBAN T J，PETROPOULOS A E，etal.Design and optimization of low-thrust trajectories with gravity assists［J］.Journal of Spacecraft and Rockets，2003，40(3):380-387.

［6］YAM C H，MCCONAGHY T T，CHEN K J，et al.Preliminary design of nuclear electric propulsion missions to the outer planets［C］//AIAA/AAS Astrodynamics Specialist Conference and Exhibit.Providence，Rhode Island:［s.n.］，2004.

［7］PETROPOULOS A E，LONGUSKI J M，F(xiàn)GLIO E P B.Trajectories to jupiter via gravity assists from Venus，Earth，and Mars［J］.Journal of Spacecraft and Rockets，2000，37(6):776-783.

［8］李小玉，鄭建華.基于形狀法和偽譜法的小推力借力優(yōu)化研究［J］.計算機仿真，2013，30(1):100-103.

［9］趙遵輝，尚海濱，崔平遠.基于混合優(yōu)化方法的小推力借力轉(zhuǎn)移軌道設(shè)計與優(yōu)化［C］//第30屆控制年會.煙臺:［s.n.］，2007.

［10］WOO B，COVERSTONE V L，CUPOLES M.Lowthrust trajectory optimization procedure for gravity-assist，outer-planet missions［J］.Journal of Spacecraft and Rockets，2006，43(1):121-129.

［11］de PASCALE P，VASILE M.Preliminary design of lowthrust multiple gravity-assist trajectories［J］.Journal of Spacecraft and Rockets，2006，43(5):1065-1076.

［12］CARNELLI I，DACHWALD B，VASILE M.Evolutionary neurocontrol a novel method for low-thrust gravity-assist trajectory optimization［J］.JournalofGuidance，Control，and Dynamics，2009，32(2):612-624.

［13］尚海濱，崔平遠，徐瑞，等.結(jié)合行星借力飛行技術(shù)的小推力轉(zhuǎn)移軌道初始設(shè)計［J］.宇航學報，2011，

32(1):29-37.

［14］MCCONAGHY T T.Design and optimization of interplanetary spacecraft trajectories［D］.West Lafayette:Purdue University，2004.

［15］CURTIS H D.Orbital mechanics for engineering students［M］.Oxford:Butterworth-Heinemann，2005.

［16］BENSON D.A gauss pseudospectral transcription for optimal control［D］.Cambridge:Massachusetts Institute of Technology，2005.

［17］楊希祥，張為華.基于Gauss偽譜法的固體運載火箭上升段軌跡快速優(yōu)化研究［J］.宇航學報，2011，32(1):15-21.

［18］宗群，田柏苓，竇立謙.基于Gauss偽譜法的臨近空間飛行器上升段軌跡優(yōu)化［J］.宇航學報，2010，31(7):1775-1781.

［19］李翔.基于自動微分算法的過程系統(tǒng)優(yōu)化［D］.杭州:浙江大學，2003.

［20］蔣占四，吳義忠，蔣慧.敏度分析的數(shù)值方法比較研究［J］.計算機與數(shù)字工程，2009，37(5):1-5.