基于人工免疫算法的地球-火星小推力轉移軌道優化

彭坤 徐世杰 果琳麗 王平

（1中國空間技術研究院載人航天總體部，北京100094）（2北京航空航天大學宇航學院，北京100191）

1 引言

目前，小推力軌道的優化方法可分為間接法和直接法兩大類。對于間接法，文獻［1］中應用主矢量理論求解最優推力方向；文獻［2］中分別用多點打靶法和最小邊界條件法求解了近地小推力軌道優化問題。文獻［3］中用多重打靶軟件BNDSCO求解了多種燃料最省的小推力地火轉移軌道。但間接法對初值敏感，難以收斂。直接法的初值均有物理意義，容易猜測，同時其通用性強。直接法的不足是尋優結果的最優性不能嚴格保證；尋優精度越高，計算量越大。在直接法研究中，文獻［4］通過直接轉錄法求解最省燃料的地月小推力轉移軌道。文獻［5］中用序列二次規劃（Sequential Quadric Programming，SQP）方法求解了地球到火星的燃料最省小推力轉移軌道。但是SQP會遇到病態梯度、初始點敏感和局部收斂問題。文獻［6］中利用退火遺傳算法求解了定常推力地球-木星小推力軌道轉移問題，有效避免了病態梯度、初始點敏感問題，但遺傳算法本身容易陷入早熟收斂。

本文依據直接法的思想，將小推力軌道優化問題轉化為非線性規劃（Nonlinear Programming，NLP）問題，采用一種新興的優化算法——人工免疫算法（Artificial Immune Algorithm，AIA）求解NLP。該算法具有個體多樣性好的特點，能有效避免早熟收斂。同時，對基本型AIA進行改進，在算法中加入引導因子以提高算法的局部尋優能力，得到引導型人工免疫算法（Guiding Artificial Immune Algorithm，GAIA），將其用于求解地球到火星的燃料最省小推力轉移軌道。最后將GAIA尋優結果與間接法最優解進行對比，驗證其最優性。

2 地球-火星小推力軌道轉移問題

2.1 系統模型

地球到火星的小推力軌道轉移過程可分3個階段：地球逃逸段，星際巡航段和火星捕獲段。在地球逃逸段和火星捕獲段，最優推力方向可以用切向推力來近似，因此研究重點可放在星際巡航段的燃料消耗上。在星際巡航過程中，地球引力主導段主要位于巡航初期，約占整個過程的0.32%；火星引力主導段主要位于巡航末期，約占星際巡航過程的0.17%。故飛行器在星際巡航段大部分時間主要受太陽引力作用，可將地球和火星的引力攝動忽略。同時，地球軌道和火星軌道的偏心率均小于0.1，兩軌道平面夾角小于2°，可將地球軌道和火星軌道簡化為同平面的圓軌道。因此，可建立二維極坐標系（如圖1所示）描述飛行器的運動，選擇太陽為極點O，定義由太陽指向地球初始位置的射線為極軸Ox。飛行器的運動可由二維極坐標系下的位置速度攝動方程來描述，其形式如下：

圖1 地球-火星軌道轉移極坐標系Fig.1 Polar coordinate system for Earth-Mars orbit transfer

式中r、θ、vr、vθ和m分別是飛行器的極半徑 （日心距）、極角、徑向速度、橫向速度和質量；μs為太陽引力常數；Ft為發動機推力，其幅值假設為常數；u為推力方向角；w為排氣速度。由于忽略地球和火星的引力，飛行器只受太陽引力Fs和小推力發動機推力Ft的作用。

巡航段中地球到火星的軌道轉移可設為從地球軌道到火星軌道。初始條件為飛行器在地球軌道上運動，終端約束為飛行器到達火星軌道。令變軌初始時刻為t0＝0，終端時刻tf自由，相應的邊界條件為

式中r0為地球軌道的日心距；θ0為飛行器在地球軌道上的逃逸點極角；m0為飛行器的初始質量；rf為火星軌道的日心距。

2.2 模型歸一化

在小推力軌道優化過程中，由于各個狀態變量的量級相差較大，尋優過程可能會丟失有效位數而難以收斂。為此，可將系統模型進行歸一化處理。

1）首先定義兩個參考變量r＊和m＊［2］，其表達式為r＊＝r0，m＊＝m0。

2）然后定義相關的歸一化變量為

所有變量歸一化后均為量綱為1的變量，且幾乎處于同一數量級上。極角本身是量綱為1的變量，不進行歸一化處理。

3）相應的狀態方程變為

由式（2）可知，極角在狀態方程中是解耦的，可先不考慮的狀態變化。

4）相應的邊界條件為

地球到火星的小推力轉移軌道有很多條，為盡可能增加有效載荷的質量，需要設計一條燃料消耗最省的轉移軌道，即要求以下性能指標達到最大：

3 參數化及約束處理

3.1 控制函數參數化

在小推力軌道優化問題中，控制函數的搜索空間是一個泛函空間，不能直接應用優化方法進行求解，必須先將控制函數參數化。目前主要有3種參數化方法［7］：直接離散法、多段參數插值方法和函數逼近方法。這3種方法各有優缺點，本文采用一種新的方法——插值函數逼近法，其步驟如下：

1）采用函數逼近方法逼近控制變量函數。由維爾斯特拉斯 （Weierstrass）逼近定理可知，對在閉區間內連續的任何函數，都可用多項式級數一致逼近。本文的控制變量函數為閉區間內連續的函數，故可采用多項式函數進行逼近，其形式如下：

式中p＝1，2，…，N；N為冪函數的個數。

2）采用插值擬合的方法確定多項式函數的系數。

3）確定插值的特征點位置。采用直接離散方法，將最優控制問題的積分區間等分為M-1段，將M個時間節點上的控制向量設為特征點。

這種插值函數逼近方法的待優化變量少，具有明確的物理意義，容易猜測取值范圍，同時計算速度快，函數逼近的精度高。

3.2 約束處理

小推力軌道優化問題是具有終端狀態約束的最優控制問題，如何處理約束關系到直接法求解小推力軌道優化問題的成敗。目前處理約束的方法主要有修復不可行解法和罰函數法。對于求解小推力軌道優化這類高維的非線性優化問題，罰函數法的應用效果更好。本文采用罰函數的形式定義評價函數faff（X）：

式中X為待優化量；σ1，σ2，σ3分別為3個懲罰項的權重系數。如果權重系數太小，起不到懲罰的作用，而權重系數太大又會影響全局尋優。因此，在尋優初期應放寬約束的懲罰，得到盡可能多的可行解以便尋找最優解；而在尋優后期應加大約束的懲罰，使不滿足約束的可行解被淘汰，從而得到滿足約束的最優解。為此，可吸取模擬退火算法的思想對權重系數進行自適應調整，其自適應調整規則如下：

式中λ＝1，2，3；β為溫度下降系數，決定了權重系數的變化率，β∈［0，1］；Tλ為與σλ對應的初始溫度，決定了權重系數最終值；g為當前迭代代數；gmax為最大迭代代數。為保證權重系數在尋優初期不會因數值太小而失去約束功能，將退火曲線值與常數1取最大值作為權重系數。

4 人工免疫算法

經過參數化和約束處理，小推力軌道優化問題轉化為高維非線性規劃問題。目前，非線性規劃算法應用最多的是SQP。但SQP在求解優化問題時會出現病態梯度、初始點敏感和局部收斂問題。近年來，遺傳算法（Genetic Algorithm，GA）由于不依賴梯度信息而在軌跡優化上得到廣泛應用，但遺傳算法容易陷入早熟收斂，局部搜索能力較差。本文采用AIA求解非線性規劃問題，該算法具有尋優成功率高、個體多樣性好的特點，不易陷入局部最優。

4.1 人工免疫算法原理

人工免疫算法模擬生物免疫系統，將待優化的問題對應抗原，可行解對應抗體，可行解的質量對應抗體與抗原的親和度，將尋優過程與生物免疫系統識別抗原并實現抗體進化的過程對應起來，形成一種智能優化算法。

4.2 人工免疫算法設計

按照人工免疫算法原理，設計適合小推力軌道優化問題的人工免疫算法：

1）計算親和度：由于AIA求解的是最優問題的最大值，因此可將親和度取為式（3）所示的形式。這里，X為待優化抗體，在本文的小推力軌道優化問題中X＝［u1，u2，…，uM，］T。

2）計算濃度和激勵度：在尋優過程中，AIA優化算法對濃度過高的抗體進行抑制以保持個體的多樣性?？贵w濃度fden（XI）計算方法為

式中n為種群中抗體個數；XI為種群中的第I個抗體；fbff（XI，Xi）為抗體XI與抗體Xi的相似度，其表達式如下：

式中XI，j和Xi，j分別是抗體XI和抗體Xi的第j個變量；L為抗體個體的變量個數。

激勵度計算方法［8］為

式中faff（XI）為抗體XI的評價函數，即親和度；fden（XI）為抗體XI的濃度；a為計算參數，可根據實際情況確定。激勵度是對抗體質量的最終評價結果，它綜合考慮了抗體的親和度和濃度。個體的親和度越大，濃度越小，激勵度越大。

3）免疫操作：首先根據激勵度進行免疫選擇，然后進行m次克隆、變異操作，實現局部搜索。變異操作采用如下規則：

式中XI，j，k是抗體XI的第k個克隆體的第j個變量；δ為定義的鄰域范圍；Pr是0和1之間的隨機數；Pm為變異概率。最后進行克隆抑制，找出變異后的抗體及源抗體中親和度最高的抗體AI替代源抗體XI，使得種群中抗體個數不變。

4）種群刷新：對激勵度低的抗體，AIA將進行刪除并隨機生成新抗體Bi進行替代。

4.3 引導型局部搜索策略

基本型AIA的局部搜索機制可簡單描述為克隆→變異→克隆抑制，盲目性較大，沒有信息交換，導致AIA的局部搜索能力不強。本文將當前群體搜索到的最優解引入局部搜索中，得到了引導型人工免疫算法GAIA。相對基本型AIA，其變異算子改寫為

式中Xbj為當前種群的最優解Xb的第j個變量；τ為克隆數目；G為引導因子，取值為［0，1］，表示Xb對克隆體的引導強度。優化過程初期應取較小的引導因子，有利于全局搜索；優化過程末期應取較大的引導因子，有利于局部搜索。為此，定義引導因子的自適應調整公式如下：

式中b為調整參數；fb為當前種群最優抗體的親和度；f0是引導因子G的跳變閾值。當fb接近f0時，引導因子G從較小值跳變到較大值，參數b可調整引導因子G跳變處曲線的變化趨勢。優化前期f0取為1.01fb，優化后期f0取為0.99fb。若種群最優抗體的親和度fb不再增大，則f0維持不變。

5 仿真及分析

地球軌道的日心距r0＝1AU（AU 為天文單位，1AU＝1.496×1011m）；太陽引力常數μs＝1.327 15×1020m3／s2；火星軌道的日心距為rf＝2.279 4×1011m；設飛行器在地球軌道的逃逸點極角為θ0＝-10°，初始質量為m0＝800kg，發動機推力為Ft＝0.5N，排氣速度為w＝3.726 527×104m／s。在參數化方法中，取M＝10，N＝6。在罰函數法中，取β＝0.005，Tλ＝0.1 （λ＝1，2，3）。

由地球-火星小推力軌道轉移的優化經驗可知，推力角的前5個特征值在 ［0，π］內，后5個特征值在［π，2π］內。轉移時間一般在200天以上，可設f取值范圍為［3.5，4］。

人工免疫算法中，種群規模n＝100，最大迭代代數gmax＝100，克隆規模τ＝20，變異概率Pm＝0.8，激勵度計算參數a＝2，引導因子G調整參數b＝150。由于GAIA是一種隨機搜索算法，本文進行10次尋優仿真，并與自適應遺傳算法 （Adaptive Genetic Algorithm，AGA）的尋優結果進行對比，如表1所示。將GAIA得到的最優歸一化數據還原到真實模型，可得燃料最省的地球-火星小推力轉移軌道，如圖2和圖3所示。

圖2為GAIA求得的最優地球-火星小推力轉移軌道的飛行軌跡。為驗證GAIA尋優結果的最優性，將其與Pontryagin極大值原理 （間接法）求得的結果進行對比，如圖3所示。由圖3可知，采用這兩種算法得到的小推力軌道轉移結果基本相同。其中日心距和極角的變化曲線幾乎重合，而徑向速度和橫向速度稍有差異但變化趨勢相同，說明GAIA的優化精度高。GAIA求得的最優燃料消耗僅比間接法的結果多1.55%，其原因是多項式函數對最優控制函數的擬合存在一定誤差 （如圖4所示）。

圖5為GAIA尋優過程中最大親和度隨迭代代數變化的曲線。由圖5可知，GAIA在第10代就搜索到最優值附近，第40代完全搜尋到最優值。這說明GAIA的收斂速度比較快，同時局部尋優能力比較強，能夠有效跳出局部最優點。

表1 GAIA和AGA的尋優結果Tab.1 Optimization results obtained by GAIA and AGA

圖2 地球-火星小推力軌道轉移最優軌跡Fig.2 Optimal Earth-Mars low-thrust trajectory

圖3 狀態變量變化曲線Fig.3 Time history of state variables

圖4 推力方向角變化曲線Fig.4 Time history of direction angle of the thrust

圖5 最大親和度隨迭代代數變化的曲線Fig.5 Maximum affinity value of every generation

6 結束語

本文首先通過插值函數逼近法和自適應權重罰函數，將地火小推力轉移軌道優化問題轉化為NLP；然后在基本型AIA算法中加入引導因子以提高局部收斂能力，得到GAIA，并將其應用到地球-火星小推力轉移軌道優化中。仿真結果證明，與AGA相比，GAIA能搜索到滿足終端狀態約束的最優地球-火星小推力轉移軌道，局部尋優能力強；同時，該算法求得的最優軌道與間接法求得的結果近似，尋優精度高，且其求解難度比間接法小，不需要猜測初值，易于編程實現，魯棒性好。通過仿真，證明了GAIA具有良好的全局和局部搜索能力，為其他優化問題提供了一種新的優化方法。

［1］REDDING D，BREAKWELL J V.Optimal Low-Thrust Transfers to Synchronous Orbit［J］.Journal of Guidance，Control，and Dynamics，1984，7（2）：148-155.

［2］CHUANG C H，GOODSON T D，HANSON J.Computation of Optimal Low-and Medium-Thrust Orbit Transfer［C］.AIAA Guidance，Navigation and Control Conference，Monterey，CA，Aug.9-11，1993：1391-1402.

［3］OBERLE H J，TAUBERT K.Existence and Multiple Solution of the Minimum-Fuel Orbit Transfer Problem ［J］.Journal of Optimization Theory and Applications，1997，95（2）：243-262.

［4］BETTS J T，ERB S O.Optimal Low Thrust Trajectories to the Moon ［J］.SIAM Journal of Applied Dynamical Systems，2003，2（2）：144-170.

［5］尚海濱，崔平遠，欒恩杰.地球-火星的燃料最省小推力轉移軌道的設計與優化 ［J］.宇航學報，2006，27（6）：1168-1173.SHANG HAIBIN，CUI PINGYUAN，LUAN ENJIE.Design and Optimization of Earth-Mars Optimal-Fuel Low-Thrust Trajectory ［J］.Journal of Astronautics，2006，27（6）：1168-1173.

［6］任遠，崔平遠，欒恩杰.基于退火遺傳算法的小推力軌道優化問題研究 ［J］.宇航學報，2007，28（1）：162-166.REN YUAN，CUI PINGYUAN，LUAN ENJIE.Low-Thrust Trajectory Optimization Based on Annealing-Genetic Algorithm ［J］.Journal of Astronautics，2007，28（1）：162-166.

［7］陳剛，萬自明，徐敏，等.飛行器軌跡優化應用遺傳算法的參數化與約束處理方法研究 ［J］.系統仿真學報，2005，17（11）：2737-2740.CHEN GANG，WAN ZIMING，XU MIN，et al.Parameterization Method and Constraints TransformationMethod for RLV Reentry Trajectory Optimization Using Genetic Algorithm ［J］.Journal of System Simulation，2005，17（11）：2737-2740.

［8］孫寧.人工免疫優化算法及其應用研究 ［D］.哈爾濱：哈爾濱工業大學，2006.SUN NING.Artificial immune optimization algorithm and applications［D］.Harbin：Harbin Institute of Technology，2006.