基于粒子群優(yōu)化的再入飛行器在線軌跡規(guī)劃

陳上上，何英姿、2，劉賀龍

（1.北京控制工程研究所，北京 100190；2.空間智能控制技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100190）

0 引言

航天飛機(jī)再入制導(dǎo)取得了巨大成功，但其采用的參考剖面離線規(guī)劃與在線跟蹤的制導(dǎo)策略，不僅設(shè)計(jì)過程費(fèi)時(shí)費(fèi)力，而且難以自主應(yīng)對異常事件的發(fā)生，常需要人的干預(yù)［1］。為提高第二代可重復(fù)使用飛行器的自主飛行能力，從而提高安全性、可靠性，同時(shí)降低成本，美國航空航天局與美國空軍分別發(fā)起了先進(jìn)制導(dǎo)與控制計(jì)劃、綜合自適應(yīng)制導(dǎo)與控制計(jì)劃［2－3］。在線軌跡規(guī)劃是提高自主飛行能力的一個(gè)關(guān)鍵，受到了廣泛關(guān)注。軌跡規(guī)劃的基本問題是根據(jù)飛行任務(wù)的飛行條件和技術(shù)要求，尋找一條某種性能指標(biāo)最優(yōu)而又不違背各種約束的飛行軌跡?？紤]實(shí)時(shí)性要求，在線軌跡規(guī)劃時(shí)一般對再入運(yùn)動模型進(jìn)行簡化，結(jié)合解析計(jì)算與簡單數(shù)值計(jì)算得到近似最優(yōu)解。MEASE等將航天飛機(jī)的二維軌跡規(guī)劃擴(kuò)展到三維，軌跡規(guī)劃問題分為軌跡長度問題與軌跡曲率問題：解決軌跡長度問題時(shí)，先基于簡化的運(yùn)動模型得到軌跡長度的解析表達(dá)形式，再通過補(bǔ)償縱程偏差修正總航程，最后用牛頓法迭代至滿足終端約束；處理軌跡曲率問題時(shí)，用解決長度問題獲得的阻力加速度剖面，計(jì)算相關(guān)側(cè)向加速度幅值，再積分航跡方位角與橫程的微分方程，通過搜索得到最小終端橫程誤差對應(yīng)的傾側(cè)角反轉(zhuǎn)點(diǎn)［4－6］。陸平等用偽平衡滑翔條件，將強(qiáng)約束、非線性規(guī)劃問題轉(zhuǎn)為兩個(gè)一維參數(shù)搜索問題，用割線法解決了在線優(yōu)化問題［7］。SCHIERMAN等用能量多項(xiàng)式函數(shù)擬合了高度剖面，將運(yùn)動方程降階，顯著減少了數(shù)值優(yōu)化問題的搜索維數(shù)，最終用非線性規(guī)劃解決了數(shù)值優(yōu)化問題［8］。

上述數(shù)值優(yōu)化過程多采用傳統(tǒng)優(yōu)化方法，問題復(fù)雜度增加時(shí)易陷入局部極值；優(yōu)化指標(biāo)為單一的航程偏差，優(yōu)化結(jié)果對模型參數(shù)攝動、環(huán)境不確定性的適應(yīng)性較差。PSO方法實(shí)現(xiàn)簡單，全局搜索能力強(qiáng)，在再入軌跡規(guī)劃應(yīng)用中受到關(guān)注［9－11］。本文對一種基于PSO的再入飛行器在線軌跡規(guī)劃方法進(jìn)行了研究。

1 問題描述

忽略地球扁率影響，RLV再入3自由度運(yùn)動方程為

式中：r為地心距；θ為經(jīng)度；φ為緯度；v為飛行器相對地球的速度；γ為航跡傾角；ψ為航跡方位角（北偏西為正）；μ為地球引力常數(shù)；σ為傾側(cè)角；ωe為地球自轉(zhuǎn)角速度；D，L分別為阻力加速度和升力加速度，且

此處：ρ為大氣密度；M為飛行器的質(zhì)量；S為參考面積；CD為阻力系數(shù)；CL為升力系數(shù)；α為攻角。各狀態(tài)量幾何關(guān)系如圖1所示。與航天飛機(jī)類似，本文研究的飛行器氣動系數(shù)CD，CL只與α有關(guān)。

圖1 再入幾何關(guān)系Fig.1 Geometry of entry flight

采用指數(shù)大氣密度模型

式中：ρ0為海平面大氣密度；h為飛行高度；hs為大氣歸一化高度。

一般來說，再入過程中需要實(shí)時(shí)滿足的約束主要有熱流率約束、法向氣動過載約束、動壓約束、零傾側(cè)角時(shí)的平衡滑翔約束，分別為

式中：為熱流率；cq為熱流傳遞系數(shù)；n，m為常數(shù)，一般取n＝0.5，m＝3；nz為法向氣動過載；g0為地面的重力加速度；q為動壓。

另外，再入過程還需考慮總吸熱量約束

以及航程偏差約束

式中：Rpre為解析預(yù)測航程；Rtogo為待飛航程。

本文研究的再入軌跡在線規(guī)劃問題描述為：給定初始條件、過程約束式（10）～（13）、吸熱量約束式（14）、終端位置約束式（15），考慮模型參數(shù)攝動、環(huán)境不確定性和異常事件發(fā)生等情況，在線尋找參考軌跡，使如下條件成立：

a）參考軌跡對應(yīng)的狀態(tài)變量r，θ，φ，v，γ，ψ與控制量σ，α滿足運(yùn)動方程式（1）～（6）；

b）過程約束、吸熱量約束、終端位置約束得到滿足；

c）控制量σ，α及其一階導(dǎo)數(shù)與二階導(dǎo)數(shù)不能超出姿態(tài)控制系統(tǒng)要求范圍，整個(gè)飛行過程中都能實(shí)現(xiàn)姿態(tài)配平。

2 參考剖面設(shè)計(jì)

本文軌跡規(guī)劃的方案為：

a）離線設(shè)計(jì)攻角參考剖面（α－v）；

b）在線規(guī)劃阻力加速度（D－v）參考剖面；

c）根據(jù)D－v剖面，實(shí)時(shí)得到參考升阻比，進(jìn)而得到控制量σ的大小；

d）σ的符號根據(jù)瞄準(zhǔn)誤差或通過一維搜索最優(yōu)反轉(zhuǎn)點(diǎn)確定［1、7］。

其中，D－v剖面在線規(guī)劃是軌跡規(guī)劃的核心，是本文的研究重點(diǎn)。方案中的c），d）已十分成熟，本文未做討論。

2.1 攻角剖面設(shè)計(jì)

再入α－v剖面設(shè)計(jì)主要考慮減輕防熱系統(tǒng)重量、末端能量管理段的制導(dǎo)需求，以及姿態(tài)配平要求。本文用離線優(yōu)化設(shè)計(jì)α－v剖面，再入過程中α－v剖面不再變化，根據(jù)地速直接得到參考攻角。本文對所研究的再入飛行器離線設(shè)計(jì)的α－v剖面如圖2所示。

圖2 參考攻角剖面Fig.2 Angle of attack profile

α－v剖面的函數(shù)形式為

式中：α0，vα，kα為設(shè)計(jì)參數(shù)。

α－v剖面確定后，在D－v平面內(nèi)，可直觀地用再入走廊邊界替代過程約束式（10）～（13），航程與總吸熱量也可解析表示為D與v的關(guān)系式，而D與v能通過測量系統(tǒng)或在線辨識實(shí)時(shí)獲得，因此基于阻力加速度的制導(dǎo)方法在實(shí)際工程中得到了廣泛應(yīng)用［1、4－6］。

2.2 過程約束轉(zhuǎn)化

由式（7）、（10）～（13）易得D－v平面內(nèi)各過程約束的表達(dá)式為

由上述約束形成了D－v平面內(nèi)的再入走廊，其中熱流、法向氣動過載和動壓約束共同構(gòu)成了再入走廊的上邊界（Du－v），平衡滑翔條件則獨(dú)自構(gòu)成了再入走廊的下邊界（Dl－v），如圖3所示。

圖3 再入走廊與參考阻力加速度剖面Fig.3 Reentry corridor and reference drag acceleration profile

2.3 吸熱量計(jì)算

再入的初始下降段大氣稀薄，氣動力作用很小，一般采用固定傾側(cè)角指令策略，不單獨(dú)設(shè)計(jì)參考剖面。該段結(jié)束后（高度約80km），才開始設(shè)計(jì)參考剖面。一般D－v剖面由溫控段、偽平衡滑翔段、常值阻力段和過渡段4段組成，其中溫控段包括兩個(gè)相接的二次型段。文獻(xiàn)［1］給出了各段航程的解析表達(dá)式，本文給出各段吸熱量的表達(dá)式。

吸熱量的一般表達(dá)式為

式中：t0，tF分別為初始與終端時(shí)間。再入過程中航跡傾角很小，故

將式（7）、（19）整理后代入式（18）得

式中：v0，vF分別為初始與終端速度。

因CD只與α有關(guān)，由圖2、3可知：溫控段、偽平衡滑翔段與常值阻力段阻力系數(shù)為常值，僅過渡段的阻力系數(shù)隨地速而變。

為便于表述，對D－v參考剖面的各節(jié)點(diǎn)依次編號：0為再入點(diǎn)；1為第一個(gè)二次型段起點(diǎn)；2為第一個(gè)二次型段終點(diǎn)，第二個(gè)二次型段起點(diǎn)；3為第二個(gè)二次型段終點(diǎn)，偽平衡滑翔段起點(diǎn)；4為偽平衡滑翔段終點(diǎn)，常值阻力段起點(diǎn)；5為常值阻力段終點(diǎn)，過渡段起點(diǎn)；6為過渡段終點(diǎn)。設(shè)vi，Di為節(jié)點(diǎn)i（i＝0，1，…，6）處的地速與阻力加速度。

a）溫控段

第一個(gè)二次型D－v剖面函數(shù)形式為

式中：C11，C12，C13為設(shè)計(jì)參數(shù)。

將式（21）代入式（20），可得第一個(gè)二次型段吸熱量表達(dá)式為

（a）當(dāng)C3＞0時(shí)，

（b）當(dāng)C3＝0，C2≠0時(shí)，

（c）當(dāng)C3＜0時(shí)，

第二個(gè)二次型D－v剖面函數(shù)形式為

式中：C21，C22，C23為設(shè)計(jì)參數(shù)，該段吸熱量Q2的表達(dá)式與式（22）～（24）類似。

b）偽平衡滑翔段

該段D－v剖面的函數(shù)形式為

c）常值阻力段

該段D－v剖面的函數(shù)形式為

式中：C4為設(shè)計(jì)參數(shù)，把式（28）代入式（20），可得該段吸熱量表達(dá)式為

d）過渡段

該段D－v剖面的函數(shù)形式為

式中：C5為設(shè)計(jì)參數(shù)；E為飛行器的能量，且E＝gh＋0.5v2；E6為過渡段末端能量。

求解該段吸熱量時(shí)，先擬合得h－v關(guān)系式

式中：m1，m2，m3為擬合系數(shù)。由式（30）、（31）可得

再擬合得CD－α關(guān)系式

式中：l1，l2，l3為擬合系數(shù)。

根據(jù)式（16）、（20）、（32）、（33）可得該段吸熱量表達(dá)式為

式（34）不便于進(jìn)一步解析計(jì)算，D－v剖面優(yōu)化時(shí)可通過簡單的數(shù)值積分得到該段吸熱量。

在再入走廊內(nèi)設(shè)計(jì)的D－v剖面能自動滿足各種過程約束；用本文的吸熱量表達(dá)式與文獻(xiàn)［1］給出的航程表達(dá)式，可預(yù)測吸熱量與航程誤差；若吸熱量與航程誤差不滿足要求，則需要重新設(shè)計(jì)參考剖面，重復(fù)上述過程至滿足實(shí)際需求。通常，上述迭代過程可用優(yōu)化方法完成。

3 優(yōu)化問題形成

參考剖面參數(shù)有9個(gè)，其中溫控段6個(gè)，其它3段各1個(gè)，若直接進(jìn)行優(yōu)化會十分復(fù)雜。如以各段端點(diǎn)坐標(biāo)為待優(yōu)化參數(shù)，根據(jù)任務(wù)要求及參考剖面連續(xù)性與光滑性要求，可減少優(yōu)化參數(shù)數(shù)，提高優(yōu)化效率。

3.1 參考剖面簡化設(shè)計(jì)

初始下降段采用固定傾側(cè)角飛行，通常不再調(diào)整v1，D1；為保證過渡段的航程調(diào)整能力，v5通過離線優(yōu)化設(shè)計(jì)，不做在線調(diào)整；另外，根據(jù)任務(wù)要求，一般v6，D6為固定值。

除上述5個(gè)節(jié)點(diǎn)坐標(biāo)固定外，其余參考剖面節(jié)點(diǎn)坐標(biāo)均可在線調(diào)整，調(diào)整基本原則為

a）節(jié)點(diǎn)及由節(jié)點(diǎn)連接的參考剖面應(yīng)完全在再入走廊內(nèi)；

b）再入航程與總吸熱量應(yīng)滿足任務(wù)需求；

c）各段相接的節(jié)點(diǎn)處盡量保證連續(xù)性與光滑性。

根據(jù)上述原則，為保證初始下降段與第一個(gè)二次型段的平滑過渡，設(shè)計(jì)D－v剖面時(shí)應(yīng)保證兩段交點(diǎn)處斜率相等，設(shè)均為k1。因初始下降段不做調(diào)整，故k1為已知。

同理，為保證其它各段的連續(xù)性與光滑性，則有

式中：

由式（35）～（37）可發(fā)現(xiàn)，參考剖面的9個(gè)參數(shù)能由X唯一確定：給定k1，v1，D1，v5，v6，D6，X后，由式（35）可得C11，C12，C13，之后由式（36）可得C21，C22，C23，K3，D3，由式（37）可得v4。

3.2 參考剖面優(yōu)化

當(dāng)α－v剖面確定后，相應(yīng)的再入走廊也隨之確定。但由于模型參數(shù)攝動及環(huán)境不確定性的存在，設(shè)計(jì)再入走廊時(shí)應(yīng)保守一些，另外為保證一定的側(cè)向機(jī)動能力，參考剖面不應(yīng)過于接近再入走廊下邊界。

由上述分析已知：4個(gè)節(jié)點(diǎn)參數(shù)可唯一確定D－v剖面，因此也決定了參考剖面與標(biāo)稱再入走廊邊界的最小距離。本文在優(yōu)化指標(biāo)中考慮此最小距離、總吸熱量和航程偏差。

由式（17）可得再入走廊上下邊界的D－v函數(shù)分別為

由前述已知，參考剖面參數(shù)由X唯一確定，參考阻力加速度可表示為X，v的函數(shù)，即

取設(shè)計(jì)優(yōu)化目標(biāo)為

式中：Qall為總吸熱量，且Qall＝Q1＋Q2＋Q3＋Q4＋Q5；ΔDu為參考剖面與標(biāo)稱再入走廊上邊界的距離，且 ΔDu＝min ［fu（v）－f（X，v）］；ΔDl為與下邊界的距離，且ΔDl＝min ［f（X，v）－fl（v） ］；kf1，kf2，kf3，kf4，kf5為設(shè)計(jì)參數(shù)。最后兩項(xiàng)設(shè)計(jì)為指數(shù)形式，其值隨參考剖面與走廊邊界距離的減小而迅速增大，而參考剖面與走廊邊界的距離增大到一定程度后，這兩項(xiàng)不再起作用。實(shí)際應(yīng)用時(shí)，根據(jù)任務(wù)情況與飛行器自身特性可合理選取式（42）中的某幾項(xiàng)組成優(yōu)化目標(biāo)。考慮再入走廊邊界約束以及節(jié)點(diǎn)坐標(biāo)之間的大小關(guān)系，最終的優(yōu)化問題為

4 基于改進(jìn)PSO的在線剖面優(yōu)化

圖4 二維PSO示意圖Fig.4 Two－dimensional PSO

PSO算法中的m個(gè)粒子一直并行搜索。在每次迭代過程中，記錄每個(gè)粒子的迄今最優(yōu)位置pi（k），通過比較獲得整個(gè)粒子群迄今最優(yōu)位置pg（k）。通過速度更新與位置更新得到下一步的搜索方向與位置分別為

本文定義使優(yōu)化目標(biāo)式（42）數(shù)值最小的位置為最優(yōu)位置。參數(shù)搜索的更新速度應(yīng)與該參數(shù)范圍相匹配，否則，更新速度太快易引發(fā)搜索振蕩，更新速度太慢搜索時(shí)間太長，不利于在線應(yīng)用。為此，本文將式（44）改進(jìn)為

式中：λj，ηj分別為控制個(gè)體認(rèn)知分量和群體社會分量相對貢獻(xiàn)的第j維學(xué)習(xí)率，可根據(jù)優(yōu)化量的搜索范圍設(shè)計(jì)，改進(jìn)后每維待優(yōu)化量都有與之匹配的學(xué)習(xí)率。

算法流程如圖5所示。

圖5 基于PSO的在線剖面優(yōu)化算法流程Fig.5 Flowchart of onboard profile planning algorithm based on PSO

5 仿真

對本文使用的RLV及相關(guān)飛行任務(wù)，過渡段參考剖面易接近再入走廊下邊界，從而導(dǎo)致末端側(cè)向機(jī)動能力不足，因此需在優(yōu)化指標(biāo)中考慮過渡段參考剖面與走廊下邊界的距離。取優(yōu)化指標(biāo)

表1 優(yōu)化參數(shù)初值Tab.1 Initial values of optimized parameter

采用上述設(shè)計(jì)參數(shù)，遍歷考查異常時(shí)間發(fā)生時(shí)飛行任務(wù)更改為短航程與長航程時(shí)參考剖面優(yōu)化結(jié)果。短航程下的指標(biāo)優(yōu)化結(jié)果為Qall＝0.717Qmax，ΔR＝0.310km，ΔDl＝1.398m／s2；長航程下的指標(biāo)優(yōu)化結(jié)果為Qall＝0.777Qmax，ΔR＝1.211km，ΔDl＝1.166m／s2。優(yōu)化的D－v剖面分別如圖6、7所示。由圖可知：兩種情況下過渡段參考剖面對應(yīng)最小傾側(cè)角均為30°左右，側(cè)向機(jī)動能力得到了保證，設(shè)計(jì)的參考剖面與再入走廊上下邊界都保持一定距離，在異常事件發(fā)生時(shí)能應(yīng)對一定的模型參數(shù)攝動和環(huán)境不確定性。

在普通計(jì)算機(jī)Windows XP系統(tǒng)中用Matlab仿真，優(yōu)化算法平均耗時(shí)6s，能滿足在線規(guī)劃要求。另外，根據(jù)實(shí)際需求，在C環(huán)境中對程序進(jìn)行優(yōu)化可進(jìn)一步提高該方法的實(shí)時(shí)性。

圖6 短航程優(yōu)化結(jié)果Fig.6 Optimization result of short range simulation

圖7 長航程優(yōu)化結(jié)果Fig.7 Optimization result of long range simulation

6 結(jié)束語

本文對基于改進(jìn)PSO的再入飛行器的在線軌跡規(guī)劃方法進(jìn)行了研究?？紤]了再入過程約束、終端位置約束、總吸熱量約束、模型參數(shù)攝動與環(huán)境不確定性影響，保證了再入末端的側(cè)向機(jī)動能力，提高了異常事件發(fā)生時(shí)優(yōu)化軌跡對模型參數(shù)攝動以及環(huán)境不確定性的適應(yīng)能力。仿真驗(yàn)證了方法的實(shí)時(shí)性。

［1］ HARPOLD J C，GRAVES C A.Shuttle entry guidance［J］.Journal of the Astronautical Sciences，1979，27（3）：239－268.

［2］ HANSON J M.A plan for advanced guidance and control.technology for 2nd generation reusable launch vehicles［C］／／ AIAA Guidance，Navigation，and Control Conference and Exhibit. Monterey：AIAA，2002：1－9.

［3］ SCHIERMAN J D，WARD D G，HULL J R，et al.Integrated adaptive guidance and control for re－entry vehicles with flight test results［J］.Journal of Guidance，Control，and Dynamics，2004，27（6）：975－988.

［4］ MEASE K D，CHEN D T，TEUFEL P，et al.Reduced－order entry trajectory planning for acceleration guidance［J］.Journal of Guidance，Control and Dynamics，2002，25（2）：257－266.

［5］ SARAF A，LEAVITT A，CHEN D T，et al.Design and evaluation of an acceleration guidance algorithm for entry［J］.Journal of Spacecraft and Rockets，2004，41（6）：986－996.

［6］ LEAVITT A，MEASE K D.Feasible trajectory generation for atmospheric entry guidance［J］.Journal of Guidance，Control and Dynamics，2007，30（2）：473－481.

［7］ SHEN Z J，LU P.Onboard generation of three－dimensional constrained entry trajectories［J］.Journal of Guidance，Control and Dynamics，2003，26（1）：111－121.

［8］ SCHIERMAN J D，HULL J R.In－flight entry trajectory optimization for reusable launch vehicles［C］／／AIAA Guidance，Navigation，and Control Conference and Exhibit.San Francisco：AIAA，2005：1－18.

［9］ 孫增圻，鄧志東，張?jiān)倥d.智能控制理論與技術(shù)［M］.北京：清華大學(xué)出版社，2011.

［10］ RAHIMI A，KUMAR K D，ALIGHANBARI H.Particle swarm optimization applied to spacecraft reentry trajectory［J］.Journal of Guidance，Control and Dynamics，2013，36（1）：307－310.

［11］ 謝富強(qiáng)，吳 浩，唐靈靈.基于粒子群算法的飛行器再入軌跡優(yōu)化［J］.計(jì)算技術(shù)與自動化，2008，27（4）：72－75.