基于粒子群算法的多脈沖轉(zhuǎn)移軌跡優(yōu)化*

陳 全，楊 震，羅亞中

(國防科技大學(xué)航天科學(xué)與工程學(xué)院，長(zhǎng)沙410073)

0 引言

空間脈沖軌道機(jī)動(dòng)問題是自20世紀(jì)50年代以來就得到廣泛研究的航天器軌跡優(yōu)化問題［1-3］.早期的研究主要基于開普勒軌道解析分析算法，通過函數(shù)極值條件確定脈沖轉(zhuǎn)移軌道的最優(yōu)性條件.此后隨著研究問題的復(fù)雜程度增加，單一的解析分析方法難以奏效，間接法和直接法得到了發(fā)展.前者主要是基于主矢量理論構(gòu)造優(yōu)化模型，后者則是構(gòu)造一般性非線性規(guī)劃模型.在20世紀(jì)90年代之前主要采用梯度優(yōu)化算法，此后遺傳算法等進(jìn)化算法得到了廣泛應(yīng)用.最近十多年，通過構(gòu)造一般性優(yōu)化模型，采用進(jìn)化算法來完成脈沖軌道機(jī)動(dòng)問題是廣泛采用的研究思路［3-6］.

軌道機(jī)動(dòng)問題通常包括3類問題:攔截、轉(zhuǎn)移和交會(huì)［2］.目前對(duì)于多脈沖最優(yōu)交會(huì)問題，有相對(duì)成熟的基于Lambert算法構(gòu)造的一般性多脈沖最優(yōu)交會(huì)優(yōu)化模型，并得到了廣泛應(yīng)用［4-6］;而對(duì)多脈沖最優(yōu)轉(zhuǎn)移一般性優(yōu)化模型的研究則較少，只有Abdelkhalik 和 Mortari［7］作了初步研究.

本文在文獻(xiàn)［7］和［8］的研究基礎(chǔ)上，對(duì)基于Lambert算法的多脈沖最優(yōu)轉(zhuǎn)移一般性優(yōu)化模型進(jìn)行了探討.文獻(xiàn)［7］中僅以脈沖變軌處航天器真近點(diǎn)角作為設(shè)計(jì)變量，本文同時(shí)探討了脈沖時(shí)刻作為設(shè)計(jì)變量的優(yōu)化模型，并對(duì)比分析了優(yōu)化效果;此外本文采用粒子群優(yōu)化(PSO，particle swarm optimization)算法［9］求解優(yōu)化問題，該算法新近在航天器小推力軌道優(yōu)化［10］、多脈沖交會(huì)軌道優(yōu)化［11］等中得到了成功應(yīng)用，本文通過仿真研究表明，其性能優(yōu)于文獻(xiàn)［7］中使用的遺傳算法.

1 多脈沖轉(zhuǎn)移軌跡優(yōu)化問題描述

1.1 Lambert問題

Lambert問題是軌道動(dòng)力學(xué)中的兩點(diǎn)邊值問題.在給定航天器初始位置r1、終端位置r2和轉(zhuǎn)移時(shí)間Δt的情況下，通過求解轉(zhuǎn)移時(shí)間方程獲得航天器在r1、r2處的速度v1、v2.Lagrange形式的轉(zhuǎn)移時(shí)間方程［12］可表示為

式中，μ為地球引力常數(shù)，α、β為L(zhǎng)agrange參數(shù)，a為轉(zhuǎn)移軌道半長(zhǎng)軸，s=(r1+r2+d)/2，N為轉(zhuǎn)移圈數(shù)，d=|r1-r2|，r1，r2分別為位置矢量 r1，r2的模.本文采用Vallado［12］求解 Lambert問題的普適變量算法，轉(zhuǎn)移速度公式為

若初始和終端時(shí)刻航天器速度分別為v-1，v+2，則Lambert轉(zhuǎn)移所需要的總速度增量為

式中“-”、“+”分別表示變軌前和變軌后通過式(2)求解.

當(dāng)式(1)中的轉(zhuǎn)移時(shí)間Δt較短時(shí)，Lambert轉(zhuǎn)移軌道僅有一條，此時(shí)Lambert轉(zhuǎn)移軌道就是燃料最省軌道.當(dāng)飛行時(shí)間較長(zhǎng)時(shí)，Lambert轉(zhuǎn)移軌道有2N+1條［12］，直接采用式(2)算法獲得的轉(zhuǎn)移軌道是其中一個(gè)解，該轉(zhuǎn)移軌道呈現(xiàn)較大的偏心率，即航天器要沿這條軌道飛行很長(zhǎng)的路程才能完成轉(zhuǎn)移，這樣的燃料消耗巨大，為實(shí)際工程所不允許，此時(shí)需要求解多圈Lambert問題，從其2N+1個(gè)解［12］中尋找最省時(shí)間或最省燃料解.

1.2 多脈沖軌道轉(zhuǎn)移問題

工程應(yīng)用中航天器軌道轉(zhuǎn)移的脈沖數(shù)目一般大于2，多脈沖軌道轉(zhuǎn)移示意如圖1所示.二體假設(shè)下，航天器軌道動(dòng)力學(xué)方程為

方程(4)的解可表示為 r(t+Δt)=p(r(t)，v(t)，t，t+Δt)和 v(t+Δt)=q(r(t)，v(t)，t，t+Δt).

圖1 N脈沖轉(zhuǎn)移軌道示意圖Fig.1 N-impulse orbit transfer

脈沖假設(shè)下，航天器變軌前后狀態(tài)可表示為

假設(shè)航天器在初始t0時(shí)刻狀態(tài)為r0，v0，終端0+Δt時(shí)刻狀態(tài)為 rf，vf，則該問題即尋找ti，Δvi，(i=1，2，…，n)在滿足如下約束的條件下:

使得總速度增量最小

對(duì)轉(zhuǎn)移時(shí)間較短的兩脈沖變軌，式(6)存在唯一解，可通過求解Lambert問題獲得最省燃料轉(zhuǎn)移軌道;對(duì)轉(zhuǎn)移時(shí)間較長(zhǎng)的兩脈沖變軌或n＞2的多脈沖變軌，式(6)的解不唯一，需要通過求解多圈Lambert問題或優(yōu)化多脈沖變軌方案以尋找滿足約束的最省燃料轉(zhuǎn)移軌道.式(6)～(7)所示優(yōu)化問題的參數(shù)化處理方法可以分為可行解迭代算法和非可行解迭代算法.前者迭代過程中無需約束處理，所產(chǎn)生的每一個(gè)解均為可行解;而后者則需要算法的最終收斂才能獲得可行解［3，5］.試驗(yàn)表明，可行解迭代算法的收斂效率比非可行解好，但其優(yōu)化變量需要進(jìn)行精心設(shè)計(jì).

1.3 基于Lambert算法的多脈沖轉(zhuǎn)移軌跡優(yōu)化模型

方程(6)所示的優(yōu)化模型適用于一般的多脈沖共面或異面軌道轉(zhuǎn)移問題.二體假設(shè)下，航天器所在的初始軌道可由5個(gè)軌道根數(shù)確定:半長(zhǎng)軸aI，偏心率eI，軌道傾角iI，升交點(diǎn)赤經(jīng)ΩI和近地點(diǎn)幅角ωI確定.同理，終端軌道由aF，eF，iF，ΩF和 ωF確定.本文考慮的軌道轉(zhuǎn)移問題不限制航天器到達(dá)終端軌道的時(shí)間及位置，即終端真近點(diǎn)角θF、轉(zhuǎn)移時(shí)間Δt為自由設(shè)計(jì)變量.

為了瞄準(zhǔn)終端位置速度，至少需要施加2次機(jī)動(dòng).當(dāng)機(jī)動(dòng)脈沖數(shù)目n≥2時(shí)，設(shè)計(jì)變量一共4n個(gè)，包括脈沖施加時(shí)刻(n個(gè))和脈沖分量(3n個(gè)).為了建立可行解迭代優(yōu)化模型，本文首先選取最后2次脈沖(共6個(gè)脈沖分量)用于滿足終端約束條件，通過求解Lambert問題獲得;其余4n-6個(gè)設(shè)計(jì)變量通過粒子群優(yōu)化算法給出.當(dāng)不限制航天器到達(dá)終端軌道的時(shí)間及位置時(shí)，優(yōu)化變量一共有4n-5.

對(duì)二體軌道，可分別采用時(shí)間或真近點(diǎn)角來表示航天器在軌位置，則軌道轉(zhuǎn)移的脈沖施加位置也可分別采用變軌時(shí)刻和變軌點(diǎn)真近點(diǎn)角給出:

1)以變軌點(diǎn)真近點(diǎn)角θi表示脈沖施加位置.此時(shí)優(yōu)化變量可表示為

式中Δt為最后兩次脈沖的間隔時(shí)間，即Lambert變軌的轉(zhuǎn)移時(shí)間.對(duì)軌道轉(zhuǎn)移問題，不限制其到達(dá)終端軌道的時(shí)間，因此將Δt也作為優(yōu)化變量，這與時(shí)間固定的交會(huì)問題有所不同.

2)以變軌點(diǎn)時(shí)刻ti表示脈沖施加位置.此時(shí)優(yōu)化變量可表示為

式中Δθ為最后兩次脈沖的間隔角距，即Lambert變軌的轉(zhuǎn)移角度.對(duì)軌道轉(zhuǎn)移問題，不限制其到達(dá)終端軌道的位置，因此將Δθ也作為優(yōu)化變量.

綜上，基于Lambert算法的多脈沖轉(zhuǎn)移軌跡優(yōu)化模型可描述為

對(duì)n脈沖軌道轉(zhuǎn)移，共有n-1段轉(zhuǎn)移軌道.實(shí)際任務(wù)中，航天器在某一段轉(zhuǎn)移軌道上的停泊可能會(huì)超過一圈.此時(shí)需要通過設(shè)計(jì)變軌點(diǎn)真近點(diǎn)角θi或變軌時(shí)刻ti的范圍來控制停泊圈數(shù)，對(duì)最后一段轉(zhuǎn)移軌道，需要通過求解多圈Lambert問題［3］確定.

2 粒子群算法

2.1 算法基本模型

粒子群算法基于群體智能，采用全局搜索策略尋優(yōu).每個(gè)粒子根據(jù)全局最優(yōu)位置和自身最優(yōu)位置更新速度.粒子群算法相對(duì)遺傳算法具有更加高效的信息共享機(jī)制，即啟發(fā)性更強(qiáng).因此，理論上粒子群算法收斂速度更快、全局搜索能力更強(qiáng)［9］.

設(shè)目標(biāo)搜索空間為M維，粒子總數(shù)為I，總迭代次數(shù)為K.第i個(gè)粒子在M維空間的位置矢量為Xi=(xi1，xi2，…，xiM)，飛行速度矢量為 Vi=(vi1，vi2，…，viM)，第i個(gè)粒子的歷史最優(yōu)位置為 Pi=(pi1，pi2，…，piM)，全局歷史最優(yōu)位置為 Pg=(pg1，pg2，…，pgM).根據(jù)如下公式更新粒子速度和位置:

式中，w為慣性因子，c1，c2為加速因子(常數(shù))，R1，R2，R3為［0，1］上的隨機(jī)數(shù).迭代中若xim，vim超出邊界，則取邊界值.為保證算法有較好的收斂性和搜索能力，可令w為隨迭代次數(shù)變化的線性減小函數(shù).對(duì)于該多脈沖轉(zhuǎn)移優(yōu)化問題，適應(yīng)度函數(shù)F即目標(biāo)函數(shù)Δv.

2.2 優(yōu)化流程

每個(gè)優(yōu)化變量的上下界可由軌道轉(zhuǎn)移問題的動(dòng)力學(xué)背景給出.同時(shí)采用式(10)模型獲得的每個(gè)解均為可行解，因而不涉及等式或不等式約束處理.基于粒子群算法的優(yōu)化流程如圖2所示.

圖2 粒子群算法流程圖Fig.2 Flowchart of particle swarm optimization

1)設(shè)置粒子總數(shù)I和最大迭代次數(shù)K，設(shè)置航天器初始軌道要素E0和終端軌道要素Ef，優(yōu)化變量上下界Dlow，Dup，根據(jù)式(11)第三式隨機(jī)初始化各個(gè)粒子速度和位置，迭代次數(shù)k=1.

2)求解式(10)所示優(yōu)化模型，獲得變軌總速度增量，計(jì)算各個(gè)粒子適應(yīng)度值.

3)對(duì)每個(gè)粒子，若其適應(yīng)度值比當(dāng)前個(gè)體極值更優(yōu)，則更新當(dāng)前個(gè)體極值;若其適應(yīng)度值比當(dāng)前全局極值更優(yōu)，則更新當(dāng)前全局極值.

4)根據(jù)式(11)前兩式更新各個(gè)粒子速度和位置.

5)若k＜K，k=k+1.更新慣性因子，返回步驟2;否則，輸出優(yōu)化結(jié)果Δv、D或D'.

3 算例分析

首先基于本文建立的多脈沖軌道轉(zhuǎn)移優(yōu)化模型，分別對(duì)兩脈沖、三脈沖軌道轉(zhuǎn)移問題進(jìn)行了優(yōu)化求解，通過與解析解、現(xiàn)有文獻(xiàn)研究結(jié)果對(duì)比，驗(yàn)證了優(yōu)化模型和粒子群算法的正確性和高效性;其次分別采用變軌位置和變軌時(shí)間為設(shè)計(jì)變量，分析其對(duì)優(yōu)化結(jié)果的影響.以下各算例均為單圈轉(zhuǎn)移，Lambert問題采用式(2)求解.

3.1 兩脈沖共面軌道轉(zhuǎn)移

采用與文獻(xiàn)［7-8］相同的兩脈沖共面圓軌道轉(zhuǎn)移算例.中心天體為火星 (μMars=4.283×104km3/s2)，初末軌道半徑分別為r1=8 000 km，rF=15 000 km.該算例最優(yōu)解為霍曼轉(zhuǎn)移，最優(yōu)總速度增量為 0.609 km/s［13］.

粒子群算法參數(shù)設(shè)置為:粒子總數(shù)15，最大迭代次數(shù)30，加速因子c1=2，c2=2，初始慣性因子w1=0.5，30次迭代時(shí)慣性因子w2=0.2.選用真近點(diǎn)角θ作為設(shè)計(jì)變量表示航天器變軌位置.優(yōu)化變量的取值范圍為:0≤θi≤2π(i=1，2)，0≤Δt≤30 000 s.

圖3給出了迭代過程中目標(biāo)函數(shù)的變化曲線，并將其與文獻(xiàn)［7］中遺傳算法的進(jìn)化歷史對(duì)比.由圖3及表1可知，粒子群算法在前5代中迅速下降，通過30次迭代求得最優(yōu)值為0.609 15 km/s，比遺傳算法更快收斂到最小值，并得到了更優(yōu)的結(jié)果，說明該算法具有較好的優(yōu)化效率.

表1 粒子群算法與遺傳算法優(yōu)化結(jié)果比較Tab.1 Comparison of optimal two-impulse solutions obtained by PSO and GAs respectively

圖3 粒子群算法進(jìn)化歷史Fig.3 Evolution history of particle swarm optimization

表2給出了粒子群算法優(yōu)化得到的兩脈沖軌道轉(zhuǎn)移最優(yōu)變軌位置和變軌脈沖.轉(zhuǎn)移角度為179.97°，非常接近霍曼轉(zhuǎn)移結(jié)果.圖4給出了轉(zhuǎn)移軌跡.

表2 最優(yōu)變軌方案Tab.2 The optimal maneuver plan

圖4 兩脈沖轉(zhuǎn)移軌跡Fig.4 Trajectory for two-impulse transfer

由于粒子群算法具有隨機(jī)性，因此進(jìn)行了50次仿真.將結(jié)果與文獻(xiàn)［7］中遺傳算法求得最優(yōu)解比較，如表1所示，本文采用的粒子群算法優(yōu)化結(jié)果更接近霍曼變軌.且相比遺傳算法［7-8］，在相同迭代次數(shù)下粒子群算法使用粒子總數(shù)更少，獲得最優(yōu)解所需要的進(jìn)化代數(shù)更少，計(jì)算目標(biāo)函數(shù)的次數(shù)大大減少，因此具有更高的優(yōu)化效率與優(yōu)化性能.

3.2 三脈沖軌道轉(zhuǎn)移

對(duì)三脈沖轉(zhuǎn)移軌道，采用與文獻(xiàn)［7］算例相同的初始條件和邊界條件，中心天體為地球初始軌道和終端軌道為共面圓軌道，初末軌道半徑分別為r1=8 000 km，rF=15 000 km.

當(dāng)變軌次數(shù)n=3時(shí)，由式(10)可知，設(shè)計(jì)變量數(shù)為7.粒子群算法參數(shù)設(shè)置為:粒子總數(shù)30，迭代次數(shù)100，加速因子c1=2，c2=2，初始慣性因子w1=0.4 ，100 次迭代時(shí)慣性因子w2=0.2 .優(yōu)化變量的取值范圍為:0 ≤ θi≤ 2π(i=1，2，3)，0 ≤Δt≤86 400 s，-3 km/s≤Δv1i≤3 km/s(i=x，y，z).優(yōu)化結(jié)果如表3所示.

對(duì)于此共面圓軌道轉(zhuǎn)移，其最優(yōu)解為霍曼轉(zhuǎn)移.由圖5和表3中數(shù)據(jù)可知，采用三脈沖轉(zhuǎn)移模型求解該問題，施加兩次脈沖后航天器已經(jīng)轉(zhuǎn)移至終端軌道，第三次速度增量接近0，也可認(rèn)為是兩脈沖轉(zhuǎn)移.而前兩次脈沖分別發(fā)生在第一轉(zhuǎn)移軌道的近地點(diǎn)和遠(yuǎn)地點(diǎn)，即最優(yōu)解仍為霍曼轉(zhuǎn)移.

表3 三脈沖轉(zhuǎn)移軌道優(yōu)化結(jié)果對(duì)比Tab.3 Comparison of optimal three-impulse transfer solutions

圖5 三脈沖轉(zhuǎn)移軌跡Fig.5 Trajectory for three-impulse transfer

表3中給出了粒子群算法優(yōu)化結(jié)果與文獻(xiàn)［7］中遺傳算法優(yōu)化結(jié)果對(duì)比.可以看出，在迭代次數(shù)較少的情況下，粒子群算法最優(yōu)解，更接近霍曼轉(zhuǎn)移.

3.3 不同優(yōu)化模型的影響分析

對(duì)于3.2中算例，使用式(10)中不同設(shè)計(jì)變量D，D'，各采用粒子群算法仿真100次，算法參數(shù)設(shè)置同3.2中算例.同時(shí)，統(tǒng)計(jì)不同設(shè)計(jì)變量下目標(biāo)函數(shù)值下降到3.800 km/s時(shí)的平均迭代次數(shù)，結(jié)果如表4所示(霍曼轉(zhuǎn)移最優(yōu)解為3.770 725 km/s).

由于收斂速度與算法參數(shù)中粒子速度有關(guān)，這里把解的每一維粒子速度最大值都設(shè)為該維度解的區(qū)間長(zhǎng)，以消除粒子速度的影響.由表4可知，對(duì)同一問題選用不同的設(shè)計(jì)變量，算法收斂速度和最優(yōu)解會(huì)有所不同.以變軌點(diǎn)真近點(diǎn)角為設(shè)計(jì)變量時(shí)結(jié)果更優(yōu)，收斂速度更快.

表4 不同優(yōu)化模型結(jié)果對(duì)比Tab.4 Comparison of different optimization models

當(dāng)選用變軌時(shí)刻為設(shè)計(jì)變量時(shí)，在轉(zhuǎn)移軌道確定前，軌道周期是未知的，且不能確定航天器在軌運(yùn)行圈數(shù)，故選取變軌點(diǎn)真近點(diǎn)角為設(shè)計(jì)變量更易實(shí)現(xiàn).但在某些特殊情況下任務(wù)對(duì)轉(zhuǎn)移時(shí)間有特定要求時(shí)，以變軌時(shí)刻為設(shè)計(jì)變量更合適.

4 結(jié)論

對(duì)于多脈沖軌道轉(zhuǎn)移問題，本文建立了基于Lambert算法的可行解迭代優(yōu)化模型，采用粒子群算法進(jìn)行求解，并對(duì)分別采用變軌點(diǎn)真近點(diǎn)角和變軌時(shí)刻作為設(shè)計(jì)變量的優(yōu)化結(jié)果進(jìn)行了對(duì)比分析.算例分析表明:1)本文基于粒子群算法的多脈沖轉(zhuǎn)移軌跡優(yōu)化模型合理可行，優(yōu)化結(jié)果更優(yōu)，效率更高;2)當(dāng)采用真近點(diǎn)角表示航天器位置時(shí)，算法更易實(shí)現(xiàn)，優(yōu)化結(jié)果更優(yōu).

［1］GOBETZ F W，DOLL J R.A survey of impulsive trajectories［J］.AIAA Journal，1969，7(5):801-834.

［2］楊嘉墀.航天器軌道動(dòng)力學(xué)與控制(下)［M］.北京:宇航出版社，2002:61-137.

［3］唐國金，羅亞中，雍恩米.航天器軌跡優(yōu)化理論、方法及應(yīng)用［M］.北京:科學(xué)出版社，2012:167-195.

［4］PRUSSING J E，CHIU J H.Optimal multiple-impulse time-fixed rendezvous between circular orbits［J］.Journal of Guidance Control and Dynamics，1986，9(1):17-22.

［5］HUGHES S P，MAILHE L M，GUZMAN J J.A comparison of trajectory optimization methods for the impulsive minimum fuel rendezvous problem ［J］.Advances in the Astronautical Sciences，2003，113:85-104.

［6］LUO Y Z，TANG G J，LEI Y J，et al.Optimization of multiple-impulse multiple-revolution rendezvous phasing maneuvers［J］.Journal of Guidance，Control and Dynamics，2007，30(4):946-952.

［7］ABDELKHALIK O，MORTARI D.N-impulse orbit transfer using genetic algorithms［J］.Journal of Spacecraft and Rockets，2007，44(2):456-460.

［8］ABDELKHALIK O.Orbit design and estimation for surveillance missions using genetic algorithms［D］.Austin:University of Texas，2005.

［9］KENNEDY J，EBERHART R.Particle swarm optimization［C］//IEEE International Conference on Neural Networks.New York:IEEE，1995:1942-1948.

［10］ZHU K J，JIANG F H，LI J F，et al.Trajectory optimization of multi-asteroids exploration with low thrust［J］.Transactions of the Japan Society for Aeronautical and Space Sciences，2009，52:47-54.

［11］PONTANI M，GHOSH P，CONWAY B A.Particle swarm optimization of multiple-burn rendezvous trajectories［J］.Journal of Guidance，Control，and Dynamics，2012，35(4):1192-1207.

［12］VALLADO D A.Fundamentals of astrodynamics and applications［M］.2nd ed.EI Segundo，California:Microcosm Press，2001:303-395.

［13］BROWN C D.Spacecraft mission design［M］.2nd ed.Washington D.C.:AIAA，1998:39-51.