張 勃，泮斌峰，唐 碩

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

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

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

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

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

1 模型建立

1.1 基本假設(shè)

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

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

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

1.2 動(dòng)力學(xué)模型

在日心黃道慣性系中，探測(cè)器的軌道動(dòng)力學(xué)方程為

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

1.3 發(fā)動(dòng)機(jī)模型

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

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

其中η為發(fā)動(dòng)機(jī)的工作效率．

1.4 行星借力軌道模型

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

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

圖1 行星借力B平面示意

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

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

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

其中

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

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

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

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

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

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

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

探測(cè)器飛越借力星體時(shí)，探測(cè)器的位置必須滿足借力星體的星歷約束，即

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

當(dāng)探測(cè)器與目標(biāo)交會(huì)時(shí)，探測(cè)器的位置與速度必須滿足

其中rM，vM分別為目標(biāo)星體的位置和速度矢量．本文中的星歷均采用JPL實(shí)驗(yàn)室提供的DE405行星星歷［15］.

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

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

目標(biāo)函數(shù)(9)記為

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

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

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

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

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

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

到達(dá)目標(biāo)星體時(shí)必須滿足交會(huì)要求(12)，記為

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

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

2.3 俯仰角和偏航角的計(jì)算

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

3 高斯偽譜法基本原理及NLP問(wèn)題求解策略

3.1 高斯偽譜法基本原理

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

3.2 NLP問(wèn)題求解策略

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

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

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

3.3 自動(dòng)微分

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

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

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

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

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

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

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

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

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

4 仿真計(jì)算與分析

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

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

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

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

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

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

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

表4 自動(dòng)微分與差分方法計(jì)算導(dǎo)數(shù)對(duì)優(yōu)化過(guò)程的影響對(duì)比

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

圖3 探測(cè)器三維軌跡

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

圖5為探測(cè)器日心距和速度的變化歷程，其中探測(cè)器的速度在進(jìn)行借力后存在突變，這是由于B平面模型沒(méi)有考慮借力的過(guò)程，認(rèn)為借力是在某一點(diǎn)瞬時(shí)完成的．由于B平面模型并不影響探測(cè)器的位置，因此探測(cè)器的日心距不存在突變．

圖5 探測(cè)器日心距與速度時(shí)間歷程

圖6給出了控制角的變化歷程.同樣也是由于B平面模型的原因，俯仰角與偏航角都存在突變.借力前，俯仰角單調(diào)增加，但始終為負(fù)值，這是因?yàn)榻鹦擒壍腊霃叫∮诘厍蜍壍腊霃?，探測(cè)器必須不斷減小日心距.而借力后由于火星軌道半徑大于金星軌道半徑，探測(cè)器首先需要增加日心距，當(dāng)達(dá)到火星軌道附近時(shí)為了滿足與火星交會(huì)的要求，探測(cè)器需再對(duì)控制角進(jìn)行調(diào)整.偏航角的變化幅值較小，這是因?yàn)榻鹦?、地球和火星軌道的軌道傾角相差不大，軌道面之間的調(diào)整所需控制力較小.

圖6 控制角時(shí)間歷程

5 結(jié)語(yǔ)

本文基于自動(dòng)微分和高斯偽譜法研究了小推力借力軌道的優(yōu)化設(shè)計(jì)問(wèn)題，將發(fā)射時(shí)間的搜索與軌道的設(shè)計(jì)包含在了同一計(jì)算框架內(nèi)，避免傳統(tǒng)混合方法中隨機(jī)搜索算法收斂速度慢、精度低的缺點(diǎn).為了增強(qiáng)算法的魯棒性、加快收斂速度，本文提出了串聯(lián)優(yōu)化策略和彈性約束策略.針對(duì)有限差分法求解導(dǎo)數(shù)信息精度有限而解析導(dǎo)數(shù)又難以獲得的問(wèn)題，本文采用自動(dòng)微分的方法求解導(dǎo)數(shù)信息，保證了導(dǎo)數(shù)計(jì)算的高精度和快速收斂.通過(guò)地球-金星-火星小推力借力軌道的仿真驗(yàn)證，證明了本文算法的正確性和有效性，為未來(lái)深空探測(cè)軌道優(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］張旭輝，劉竹生.火星探測(cè)無(wú)動(dòng)力借力飛行軌道研究［J］.宇航學(xué)報(bào)，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］.計(jì)算機(jī)仿真，2013，30(1):100-103.

［9］趙遵輝，尚海濱，崔平遠(yuǎn).基于混合優(yōu)化方法的小推力借力轉(zhuǎn)移軌道設(shè)計(jì)與優(yōu)化［C］//第30屆控制年會(huì).煙臺(tái):［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］尚海濱，崔平遠(yuǎn)，徐瑞，等.結(jié)合行星借力飛行技術(shù)的小推力轉(zhuǎn)移軌道初始設(shè)計(jì)［J］.宇航學(xué)報(bào)，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ùn)載火箭上升段軌跡快速優(yōu)化研究［J］.宇航學(xué)報(bào)，2011，32(1):15-21.

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

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

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