基于改進粒子群算法的月地轉移軌道優化

李遠平，蔡遠利，

(西安交通大學電子與信息學部，陜西，西安 710049)

近年來，國際上再次掀起月球探測高潮，美國、俄羅斯、中國、歐空局和日本都制定了各自的載人登月和月球基地計劃[1]。我國原計劃于2019年底實施“嫦娥五號”月球采樣返回任務，軌道器與返回器組合體從月球停泊軌道加速后，進入月地轉移軌道，在地面高度約5000 km處軌道器分離，返回器采用跳躍式再入，著陸場選擇為中緯度的內蒙古四子王旗地區。月地轉移軌道優化設計作為其中的技術難題之一，需要滿足月球停泊軌道、月地轉移和再入地球大氣層等方面的復雜約束條件，盡量減小返回窗口內的逃逸速度增量，降低返回燃料消耗，同時獲得最優的再入點參數，為跳躍式再入控制創造良好的起始條件。

月球返回軌道的設計方案可分為兩種：一是月球自由返回軌道設計[2]，采用多段拼接方法[3]，適用于短期的月球探測返回任務；二是月地快速轉移軌道設計，從月球停泊軌道加速后進入月地轉移軌道[4]，適用于各類登月返回任務。文獻[5]在圓形限制性三體模型下研究了地月空間對稱自由返回軌道的設計方法及其軌道特性。文獻[6]分析了地月自由返回時近地距的各種影響因素，包括出發時刻、轉移時間、傾角組合等，提出了一種可使返回時近地距逐漸向設定的目標值靠攏的單向搜索方法。文獻[7]分析了探測器低能量返回軌道的存在性與軌道動力學特性，提出了一種混合自適應遺傳算法，用于尋找月球探測器返回地球所需的最低能量以及對應的轉移軌道。文獻[8―9]對從月球返回地球的軌道進行了分析，給出了快速月地轉移軌道設計方法。文獻[10]采用了自起步解析算法，構造出多次點火的軌道機動序列，可實現從任意近月停泊軌道到地球再入狀態的轉移，但在某些情況下需要消耗較多的燃料。

對于復雜約束條件下的優化設計問題，采用群智能優化算法具有明顯的優勢：算法結構簡單、對初始值要求不高、不需要繁瑣的公式推導。常用的群智能優化算法包括遺傳算法[11]、模擬退火算法[12]、粒子群算法[13]、蟻群算法[14]、蜂群算法[15]、螢火蟲算法[16]和離子運動算法[17]。

本文基于月地快速轉移軌道的思想，建立了多段拼接的解析法軌道計算模型，考慮了月地轉移軌道和再入彈道參數等諸多約束條件，以月球停泊軌道的逃逸速度增量最小為優化目標，提出了一種改進粒子群算法，求解月心真赤道坐標系內的逃逸點軌道參數，搜索最優的逃逸點時刻、升交點經度和緯度幅角。

1 月地返回窗口選擇

月球停泊軌道通常采用圓軌道，軌道高度hMP和軌道傾角iMP(相對月球真赤道面)由月球軌道交會對接任務確定，月地返回窗口選擇主要是確定逃逸點時刻tMP的取值范圍，而逃逸速度增量ΔvMP、升交點經度?MP和緯度幅角uMP由下一步的月地轉移軌道優化過程來確定。

利用噴氣推進實驗室JPL的DE/LE405歷表來準確計算月球軌道位置及天平動參數，瞄準著陸場中心位置，考慮月地轉移時間、地心段軌道傾角、返回當地時間和返回航程等約束條件，在一年或者更長的時間段內采用幾何解析法進行覆蓋式搜索。

2 月地轉移軌道優化方法

2.1 軌道計算模型

在全局優化設計時，為了提高計算效率，采用雙二體引力模型來描述月地轉移軌道的運動規律。

1) 月球停泊軌道逃逸計算

先利用逃逸點時刻的月球停泊軌道參數，計算出探測器在月心真赤道坐標系內的逃逸點位置矢量rMP和速度矢量vMP，將其轉換到J2000.0月心天球慣性坐標系，再考慮逃逸速度增量ΔvMP后，得到月地轉移軌道的起始點位置矢量rM0和速度矢量vM0。

2) 月心段轉移軌道計算

探測器從月球停泊軌道逃逸后，進入一條月心雙曲線軌道。在月心天球慣性坐標系內計算探測器飛經月球引力影響球出口點的歷元時刻tM1、位置矢量rM1和速度矢量vM1，將其轉換到J2000.0地心天球慣性坐標系，得到相應的位置矢量rE1和速度矢量vE1。

3) 地心段轉移軌道計算

探測器飛出月球引力影響球后，進入一條地心橢圓軌道。在地心天球慣性坐標系內計算開普勒軌道參數并推算至再入點，得到月地轉移時間TET、地心軌道傾角iE、升交點赤經?E、近地點幅角ωE和近地點高度hEP，以及再入點歷元時刻tR、當天時間τR、位置矢量rR和速度矢量vR。

4) 再入彈道參數計算

計算再入相對速度VR、再入角θR、再入點地理經度λR、再入點地理緯度φR、再入彈道面傾角iR和再入段航程LR。

2.2 優化目標函數

通常多目標優化問題可描述如下：

式中：x為決策向量，x∈Rm；J為目標函數向量，J∈Rn；fj(x)為目標函數；gj(x)為約束條件。

對于月地轉移軌道優化問題，目標函數可簡化為：

式中，Δv1、Δv2分別為月球逃逸速度增量的下界和上界。

約束條件如下：

式中：hEP0為參考近地點高度；TET0為標稱月地轉移時間；ΔTET為月地轉移時間偏差允許值；τR0為標稱再入點當天時間；ΔRτ為再入點當天時間偏差允許值；iE0為標稱地心軌道傾角；ΔiE為地心軌道傾角偏差允許值；ΔiR為再入彈道面傾角偏差允許值；ΩE0為標稱地心軌道升交點赤經；ΔEΩ為地心軌道升交點赤經偏差允許值；ωE0為標稱地心軌道近地點幅角；ΔEω為地心軌道近地點幅角偏差允許值；θR0為標稱再入角；ΔRθ為再入角偏差允許值；LR0為標稱再入航程；ΔLR為再入航程偏差允許值。

為了求解上述多約束條件下的優化問題，定義如下罰函數：

上述問題轉化為：

式中：χ為收斂域條件標識，0為不滿足，1為滿足。

2.3 優化求解方法

由于搜索空間大，約束條件多，目標函數復雜，本文提出了一種改進粒子群算法。

2.3.1 標準粒子群算法

粒子群算法(Particle Swarm Optimization, PSO)的基本概念源于對于鳥群捕食行為的簡化社會模型的模擬，在解空間追隨最優的粒子進行搜索[18]。

在PSO算法中，每個粒子在飛行過程所經歷過的最好位置，就是粒子本身找到的最優解，稱之為個體極值(pi)。整個粒子群所經歷過的最好位置，就是整個群體找到的最優解，稱之為全局極值(pg)。實際操作中，通過適應度值(fitness value)來評價粒子的“好壞”程度。

設粒子的群體規模為N，第j次迭代時第i個粒子的位置表示為，所經歷過的“最好”位置記為，其速度表示為，群體中所有粒子所經歷過的“最好”位置表示為。在標準PSO算法中，第i個粒子的速度和位置更新規則如下：

式中：j為迭代次數；c1、c2均為學習因子，取值為常數；r1、r2均為[0，1]區間上的隨機數；w為慣性權重(Inertia Weight)。

設適應度函數為F(X)，第i個粒子個體極值的更新公式為：

所有粒子的全局極值更新公式為：

標準PSO算法在應用時，存在收斂速度慢和早熟的缺陷[19]。

2.3.2 改進粒子群算法

為了克服標準PSO算法的缺點，針對上述月地轉移軌道優化問題，提出一種改進粒子群算法(Modified Particle Swarm Optimization, MPSO)，利用適應度函數F(X)來更新慣性權重w，并對粒子的速度加以約束。

對于月地轉移軌道優化問題，可將適應度函數F(X)的最大化變換為目標函數J的最小化。

對不同的粒子位置，采用不同的慣性權重

式中，JL為標函數的上限值，取值為10。

隨著目標函數J的最小化，慣性權重→0，加快了粒子位置的收斂速度。

定義粒子的位置向量X為：

其搜索空間如下：

式中：[t1,t2]為月球停泊軌道逃逸點時間窗口；[Δv1,Δv2]為月球停泊軌道逃逸速度增量區間，[850,860]m/s；[?1,?2]為月球停泊軌道升交點經度區間，[0, 2π]；[u1,u2]為月球停泊軌道緯度幅角區間，[0,2π]。

粒子位置的初值為各維參數區間內的隨機數組合，即：

式中，r為[0，1]區間上的隨機數。

粒子速度的初值為：

粒子速度的最大約束值為：

為進一步提高MPSO算法的收斂能力，利用月球停泊軌道逃逸點時刻tMP與再入點當天時間τR的相關性，當約束條件g2(x)≤0滿足時，對第一維位置參數引入下面的隨機反饋控制：

2.3.3 算法收斂性分析

利用式(24)進行位置更新的粒子，最終會收斂到一個穩定的位置[20]：

將式(35)代入式(24)中，可得：

式中，c=c1r1+c2r2。

設yj=pc-Xj，算法的速度和位置更新公式可變換為：

式(37)中，系數矩陣的特征方程為：

求解式(38)可得：

算法的收斂條件為max{|λ1|, |λ2|}＜1，代入式(39)可得：

算法的收斂區域見圖1。

圖1 改進粒子群算法的收斂區域Fig.1 The convergence region of MPSO algorithm

對于MPSO算法，由于w→0，當滿足c1+c2≤2時，可保證收斂性。

3 算例及分析

利用本文所提出的MPSO算法進行月地轉移軌道優化設計，所使用的計算機配置為2.53 GHz單核CPU、4G內存，編程環境為MATLAB 7.11.0。

設定月球停泊軌道高度hMP為210 km、軌道傾角iMP為45°，著陸場中心點地理經度λL=111.43°、地心緯度ΦL=42.16°，標稱月地轉移時間TET0=3.0 d、偏差允許值ΔTET=0.5 d，標稱地心軌道傾角iE0=45.0°、偏差允許值ΔiE=3.0°，標稱升交點赤經?E0=156.6°、偏差允許值ΔEΩ=10.0°，標稱近地點幅角ωE0=10.8°、偏差允許值ΔEω=10.0°，參考近地點高度hEP0=45.0 km，標稱再入點當天時間τR0=22:00:00(UTC)、偏差允許值ΔRτ=0.5 h，標稱再入角θR0=-6.0°、偏差允許值ΔRθ=0.5°，標稱再入航程LR0=6500 km、偏差允許值ΔLR=2000 km，再入彈道面傾角偏差允許值ΔiR=0.5°。

月地返回時間窗口為2017-01-30 03:00:00～2017-01-31 03:00:00(UTC)，粒子的群體規模數N為20。采用標準PSO算法，慣性權重w=0.73，學習因子c1=c2=0.75，搜索步數設置為300，運行一次計算程序耗時為54 s，部分粒子的目標函數值一度陷入局部極值，經過228步才收斂，但最終結果仍不滿足全部約束條件，目標函數值搜索狀態見圖2。

采用MPSO改進算法，搜索步數設置為100，運行一次計算程序耗時為18 s。每次運行可獲得一個優化結果，10次運行結果見表1。

綜合MPSO改進算法的運行結果，最小收斂步數為21，最大收斂步數為68，平均收斂步數為38，月地轉移軌道優化結果可歸納為兩組，包含所有局部最優解，滿足全部約束條件，見表2。兩組優化結果的目標函數值相差僅為0.008，選取逃逸速度增量最小的第1組結果作為全局最優解。

圖2 目標函數值搜索狀態(標準PSO算法)Fig.2 The search status of objective function(PSO)

一次軌道優化過程的目標函數值搜索狀態見圖3，逃逸速度增量搜索狀態見圖4。其中，經過25步搜索后，全部20個粒子收斂一致。

為驗證上述優化結果的正確性，繪制出目標函數的全局等高線和曲面見圖5、圖6。

從圖5和圖6中，可以明顯看到存在兩組局部最優解，使目標函數達到局部極小值，與MPSO算法的優化結果一致，表明該算法對解決搜索空間大、約束條件多的月地轉移軌道優化問題是正確有效的。

表1 10次MPSO算法運行結果Table 1 Ten computational results of MPSO algorithm

表2 月地轉移軌道優化結果Table 2 The optimized result of moon-to-earth transfer trajectory

圖3 目標函數值搜索狀態(MPSO算法)Fig.3 The search status of objective function(MPSO)

圖4 逃逸速度增量搜索狀態(MPSO算法)Fig.4 The search status of escape velocity increment(MPSO)

圖5 目標函數全局等高線Fig.5 The contour map of objective function

圖6 目標函數全局曲面Fig.6 The curved surface of objective function

4 結論

本文結合我國月球探測三期工程的實際需求，針對返回中緯度著陸場，考慮月球停泊軌道、月地轉移和跳躍式再入等方面的復雜約束條件，提出了一種改進的粒子群算法，實現了返回窗口內逃逸速度增量的最小化，降低了返回燃料消耗，有效解決了多約束條件下的月地轉移軌道優化問題：

(1) 通過罰函數法將多約束條件下的月地轉移軌道優化問題轉化為無約束優化問題，在收斂域內外使用不同的目標函數，便于采用粒子群算法進行優化求解。

(2) 用標準PSO算法求解月地轉移軌道優化問題時收斂較慢且最終結果仍不滿足全部約束條件，采用MPSO改進算法后能夠較快收斂且滿足全部約束條件，表明改進算法對解決復雜優化問題的有效性。

(3) 在預定月地返回時間窗口內，月地轉移軌道優化的多次運行結果可歸納為兩組，逃逸速度增量最小的那組結果是全局最優解。

(4) 通過目標函數的全局等高線和曲面圖分析，驗證了優化結果的正確性。

本文提出的優化方法，具有求解效率高、收斂性好、適應性強等優點，能夠有效避免陷入局部極值，可適用于解空間有界的多參數優化問題。