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

張 勃，泮斌峰，唐 碩

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

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

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

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

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

1 模型建立

1.1 基本假設

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

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

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

1.2 動力學模型

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

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

1.3 發動機模型

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

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

其中η為發動機的工作效率．

1.4 行星借力軌道模型

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

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

圖1 行星借力B平面示意

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

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

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

其中

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

2 軌道優化問題描述

2.1 性能指標函數及約束條件

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

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

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

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

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

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

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

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

2.2 最優控制問題

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

目標函數(9)記為

在飛行過程中，狀態變量必須處于合理的取值區間內，即

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

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

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

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

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

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

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

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

2.3 俯仰角和偏航角的計算

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

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

3.1 高斯偽譜法基本原理

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

3.2 NLP問題求解策略

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

圖2 兩種不同意義下的串聯優化策略

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

3.3 自動微分

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

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

考慮函數y=f(x)，f:Rn→R，可以由下面的程序實現:

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

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

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

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

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

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

4 仿真計算與分析

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

表1 探測器主要參數

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

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

表2 不同串聯優化級數求解耗時對比

表3 小推力借力軌道優化結果

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

表4 自動微分與差分方法計算導數對優化過程的影響對比

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

圖3 探測器三維軌跡

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

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

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

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

圖6 控制角時間歷程

5 結語

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

［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，FGLIO 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］李小玉，鄭建華.基于形狀法和偽譜法的小推力借力優化研究［J］.計算機仿真，2013，30(1):100-103.

［9］趙遵輝，尚海濱，崔平遠.基于混合優化方法的小推力借力轉移軌道設計與優化［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］尚海濱，崔平遠，徐瑞，等.結合行星借力飛行技術的小推力轉移軌道初始設計［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偽譜法的固體運載火箭上升段軌跡快速優化研究［J］.宇航學報，2011，32(1):15-21.

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

［19］李翔.基于自動微分算法的過程系統優化［D］.杭州:浙江大學，2003.

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