基于人工免疫算法的地球-火星小推力轉(zhuǎn)移軌道優(yōu)化

彭坤 徐世杰 果琳麗 王平

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

1 引言

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

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

2 地球-火星小推力軌道轉(zhuǎn)移問題

2.1 系統(tǒng)模型

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

圖1 地球-火星軌道轉(zhuǎn)移極坐標(biāo)系Fig.1 Polar coordinate system for Earth-Mars orbit transfer

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

巡航段中地球到火星的軌道轉(zhuǎn)移可設(shè)為從地球軌道到火星軌道。初始條件為飛行器在地球軌道上運(yùn)動(dòng)，終端約束為飛行器到達(dá)火星軌道。令變軌初始時(shí)刻為t0＝0，終端時(shí)刻tf自由，相應(yīng)的邊界條件為

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

2.2 模型歸一化

在小推力軌道優(yōu)化過程中，由于各個(gè)狀態(tài)變量的量級相差較大，尋優(yōu)過程可能會(huì)丟失有效位數(shù)而難以收斂。為此，可將系統(tǒng)模型進(jìn)行歸一化處理。

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

2）然后定義相關(guān)的歸一化變量為

所有變量歸一化后均為量綱為1的變量，且?guī)缀跆幱谕粩?shù)量級上。極角本身是量綱為1的變量，不進(jìn)行歸一化處理。

3）相應(yīng)的狀態(tài)方程變?yōu)?/p>

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

4）相應(yīng)的邊界條件為

地球到火星的小推力轉(zhuǎn)移軌道有很多條，為盡可能增加有效載荷的質(zhì)量，需要設(shè)計(jì)一條燃料消耗最省的轉(zhuǎn)移軌道，即要求以下性能指標(biāo)達(dá)到最大：

3 參數(shù)化及約束處理

3.1 控制函數(shù)參數(shù)化

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

1）采用函數(shù)逼近方法逼近控制變量函數(shù)。由維爾斯特拉斯 （Weierstrass）逼近定理可知，對在閉區(qū)間內(nèi)連續(xù)的任何函數(shù)，都可用多項(xiàng)式級數(shù)一致逼近。本文的控制變量函數(shù)為閉區(qū)間內(nèi)連續(xù)的函數(shù)，故可采用多項(xiàng)式函數(shù)進(jìn)行逼近，其形式如下：

式中p＝1，2，…，N；N為冪函數(shù)的個(gè)數(shù)。

2）采用插值擬合的方法確定多項(xiàng)式函數(shù)的系數(shù)。

3）確定插值的特征點(diǎn)位置。采用直接離散方法，將最優(yōu)控制問題的積分區(qū)間等分為M-1段，將M個(gè)時(shí)間節(jié)點(diǎn)上的控制向量設(shè)為特征點(diǎn)。

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

3.2 約束處理

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

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

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

4 人工免疫算法

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

4.1 人工免疫算法原理

人工免疫算法模擬生物免疫系統(tǒng)，將待優(yōu)化的問題對應(yīng)抗原，可行解對應(yīng)抗體，可行解的質(zhì)量對應(yīng)抗體與抗原的親和度，將尋優(yōu)過程與生物免疫系統(tǒng)識別抗原并實(shí)現(xiàn)抗體進(jìn)化的過程對應(yīng)起來，形成一種智能優(yōu)化算法。

4.2 人工免疫算法設(shè)計(jì)

按照人工免疫算法原理，設(shè)計(jì)適合小推力軌道優(yōu)化問題的人工免疫算法：

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

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

式中n為種群中抗體個(gè)數(shù)；XI為種群中的第I個(gè)抗體；fbff（XI，Xi）為抗體XI與抗體Xi的相似度，其表達(dá)式如下：

式中XI，j和Xi，j分別是抗體XI和抗體Xi的第j個(gè)變量；L為抗體個(gè)體的變量個(gè)數(shù)。

激勵(lì)度計(jì)算方法［8］為

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

3）免疫操作：首先根據(jù)激勵(lì)度進(jìn)行免疫選擇，然后進(jìn)行m次克隆、變異操作，實(shí)現(xiàn)局部搜索。變異操作采用如下規(guī)則：

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

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

4.3 引導(dǎo)型局部搜索策略

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

式中Xbj為當(dāng)前種群的最優(yōu)解Xb的第j個(gè)變量；τ為克隆數(shù)目；G為引導(dǎo)因子，取值為［0，1］，表示Xb對克隆體的引導(dǎo)強(qiáng)度。優(yōu)化過程初期應(yīng)取較小的引導(dǎo)因子，有利于全局搜索；優(yōu)化過程末期應(yīng)取較大的引導(dǎo)因子，有利于局部搜索。為此，定義引導(dǎo)因子的自適應(yīng)調(diào)整公式如下：

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

5 仿真及分析

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

由地球-火星小推力軌道轉(zhuǎn)移的優(yōu)化經(jīng)驗(yàn)可知，推力角的前5個(gè)特征值在 ［0，π］內(nèi)，后5個(gè)特征值在［π，2π］內(nèi)。轉(zhuǎn)移時(shí)間一般在200天以上，可設(shè)f取值范圍為［3.5，4］。

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

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

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

表1 GAIA和AGA的尋優(yōu)結(jié)果Tab.1 Optimization results obtained by GAIA and AGA

圖2 地球-火星小推力軌道轉(zhuǎn)移最優(yōu)軌跡Fig.2 Optimal Earth-Mars low-thrust trajectory

圖3 狀態(tài)變量變化曲線Fig.3 Time history of state variables

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

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

6 結(jié)束語

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

［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］尚海濱，崔平遠(yuǎn)，欒恩杰.地球-火星的燃料最省小推力轉(zhuǎn)移軌道的設(shè)計(jì)與優(yōu)化 ［J］.宇航學(xué)報(bào)，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］任遠(yuǎn)，崔平遠(yuǎn)，欒恩杰.基于退火遺傳算法的小推力軌道優(yōu)化問題研究 ［J］.宇航學(xué)報(bào)，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］陳剛，萬自明，徐敏，等.飛行器軌跡優(yōu)化應(yīng)用遺傳算法的參數(shù)化與約束處理方法研究 ［J］.系統(tǒng)仿真學(xué)報(bào)，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］孫寧.人工免疫優(yōu)化算法及其應(yīng)用研究 ［D］.哈爾濱：哈爾濱工業(yè)大學(xué)，2006.SUN NING.Artificial immune optimization algorithm and applications［D］.Harbin：Harbin Institute of Technology，2006.