藺君，何英姿，黃盤興

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

高超聲速飛行器再入可達(dá)域又稱作再入覆蓋域，是指飛行器以給定的初始條件返回地球，到達(dá)著陸或者交班狀態(tài)，飛行器所能允許的區(qū)域范圍[1]。再入可達(dá)域是評(píng)估航程覆蓋能力的指標(biāo)，體現(xiàn)了飛行器的潛在飛行范圍，可為飛行器任務(wù)規(guī)劃和軌跡規(guī)劃提供依據(jù)。

Vinh等在忽略了軌跡約束和哥氏力等理想情況下，計(jì)算最優(yōu)落點(diǎn)方法，在此基礎(chǔ)上，Doman和Ngo提出了次優(yōu)落點(diǎn)區(qū)域的方法[2]。雍恩米等將再入可達(dá)域問題轉(zhuǎn)化為極限約束條件下的軌跡優(yōu)化問題，并采用Gauss偽譜法進(jìn)行求解，將得到的外邊界組成的多邊形作為近似可達(dá)域[3]。文獻(xiàn)[4]將再入可達(dá)域問題利用直接打靶法轉(zhuǎn)化為四類情形的彈道優(yōu)化問題，采用序列二次規(guī)劃求解滿足約束的再入可達(dá)域，文獻(xiàn)[5]利用梯度修復(fù)算法求解軌跡優(yōu)化問題，得到再入可達(dá)域。文獻(xiàn)[6]在擬平衡滑翔條件下，將再入可達(dá)域計(jì)算問題轉(zhuǎn)化為求解次優(yōu)傾側(cè)角控制律問題，進(jìn)一步將其等價(jià)為虛擬目標(biāo)逼近問題，并將過程約束轉(zhuǎn)化為傾側(cè)角約束，得到再入可達(dá)域。根據(jù)文獻(xiàn)[7]在分析攻角優(yōu)化問題的基礎(chǔ)上，基于擬平衡滑翔條件結(jié)合飛行器再入多種約束，得到再入攻角的設(shè)計(jì)空間，給出再入攻角剖面設(shè)計(jì)方法，并應(yīng)用于再入可達(dá)域優(yōu)化。文獻(xiàn)[8]利用擬平衡滑翔條件，通過零值傾側(cè)角和最大橫向航程生成虛擬目標(biāo)集，將再入強(qiáng)約束條件轉(zhuǎn)化為攻角下邊界約束條件，得到飛行器再入可達(dá)域。文獻(xiàn)[9]將可重復(fù)使用飛行器再入可達(dá)域轉(zhuǎn)化為參數(shù)優(yōu)化問題，針對(duì)該參數(shù)求解問題，結(jié)合全局和局部優(yōu)化算法的特點(diǎn)，設(shè)計(jì)了遺傳算法和模式搜索相結(jié)合的優(yōu)化算法，獲得滿足過程約束和控制約束的再入可達(dá)域。文獻(xiàn)[10]在加速度再入制導(dǎo)基礎(chǔ)上，利用阻力加速度插值和傾側(cè)角反轉(zhuǎn)，快速估算再入可達(dá)域。文獻(xiàn)[11]在能量-阻力空間里，將飛行器再入約束條件描述為阻力加速度邊界曲線，在得到實(shí)際可用的上下邊界后，利用插值計(jì)算得到阻力方案，從而生成再入可達(dá)區(qū)域。文獻(xiàn)[12]在r-V空間通過設(shè)計(jì)“初始下降段、外邊界跟蹤段和終端調(diào)整段”三段設(shè)計(jì)的軌跡優(yōu)化方法，沿走廊上沿獲得最大縱向航程，而后通過平移邊界跟蹤段獲得固定縱向航程下的最大橫向航程，得到外邊界，并沿走廊下邊界獲得內(nèi)邊界。文獻(xiàn)[13]利用粒子群優(yōu)化和傾側(cè)角反轉(zhuǎn)相結(jié)合的方法，設(shè)計(jì)參數(shù)化傾側(cè)角剖面，實(shí)現(xiàn)可達(dá)域快速計(jì)算。文獻(xiàn)[14]利用粒子群優(yōu)化算法，結(jié)合傾側(cè)角插值模型，實(shí)現(xiàn)再入可達(dá)域快速計(jì)算。

綜上所述，現(xiàn)有的再入可達(dá)域計(jì)算問題本質(zhì)上是復(fù)雜軌跡優(yōu)化問題，并利用參數(shù)搜索方法求解。本文在此基礎(chǔ)上，通過分析傾側(cè)角剖面對(duì)航程的影響，將傾側(cè)角剖面設(shè)計(jì)為3段或5段線性分段函數(shù)，利用差分進(jìn)化算法，求解滿足航程的傾側(cè)角剖面。在得到對(duì)應(yīng)最大橫向航程傾側(cè)角剖面和最大縱向航程傾側(cè)角的基礎(chǔ)上，設(shè)計(jì)傾側(cè)角插值模型，生成傾側(cè)角指令集，近似計(jì)算可達(dá)域。

1 問題描述

1.1 再入運(yùn)動(dòng)學(xué)方程

考慮地球自轉(zhuǎn)的影響，建立高超聲速飛行器的再入無量綱運(yùn)動(dòng)學(xué)方程[15]為：

(1)

式中：r為無量綱地心距；V為飛行器無量綱飛行速度；θ和φ分別為經(jīng)度和緯度；γ和ψ分別為飛行航跡角和航向角；ωe為地球自轉(zhuǎn)角速度；σ為飛行器傾側(cè)角；L和D分別為無量綱升力加速度和阻力加速度:

L=ρR0V2SrefCL/(2m)

D=ρR0V2SrefCD/(2m)

式中:R0為地球半徑；Sref為飛行器參考面積；m為飛行器質(zhì)量；CL和CD分別為升力系數(shù)和阻力系數(shù)；ρ為大氣密度。

1.2 再入過程約束

Q=g0R0ρV2/2≤Qmax

(2)

(3)

(4)

1.3 再入終端約束

為成功保證飛行任務(wù)，飛行器接近任務(wù)終端時(shí)，需要滿足終端約束條件，通常包括[17]：

終端高度約束：

|hf-h|≤Δh

(5)

終端速度約束：

|Vf-V|≤ΔV

(6)

落點(diǎn)位置約束：

|θf-θ|≤Δθ，|φf-φ|≤Δφ

(7)

再入可達(dá)域主要采用外邊界來表征，通常為扇形區(qū)域形式，如圖1所示，E點(diǎn)是飛行器再入初始位置對(duì)應(yīng)的星下點(diǎn)，EC為最大縱程，EA和EB為最大橫程，外邊界AC和BC由再入終端約束構(gòu)成。

圖1 再入可達(dá)區(qū)域示意Fig.1 Illustrations of reentry landing footprint

2 軌跡優(yōu)化設(shè)計(jì)

高超聲速飛行器再入軌跡一般由攻角和傾側(cè)角確定。在大多數(shù)情況下，飛行器的攻角剖面根據(jù)飛行器過程的強(qiáng)約束條件，給定為速度的分段二次函數(shù)形式，如式(8)所示，飛行器再入航跡通過控制傾側(cè)角實(shí)現(xiàn)。

(8)

式中：α0為初始再入攻角；VT為確定攻角曲線開始減小的臨界速度；kα為常值系數(shù)，確定攻角下降速率。

2.1 傾側(cè)角參數(shù)化再入軌跡設(shè)計(jì)

傾側(cè)角剖面的幅值直接影響再入縱向航程，符號(hào)影響再入側(cè)向航程。將傾側(cè)角參數(shù)化為：

(9)

根據(jù)式(9)中σ1、σ2以及臨界速度V2、V3取值不同，再入傾側(cè)角剖面會(huì)具有不同的形式。當(dāng)σ1和σ2同號(hào)時(shí)，傾側(cè)角保持相同極性；當(dāng)σ1和σ2符號(hào)相反時(shí)，傾側(cè)角在[V4,V3]完成極性過渡。特別地，當(dāng)V2=V3=(V1+V4)/2，σ1= (σ0+σ2)/2時(shí)，傾側(cè)角5段分段函數(shù)則退化為3段，即：

在初始下降段，飛行器氣動(dòng)特性較弱，選取常值傾側(cè)角σ0作為控制變量。接近終端區(qū)域能量管理時(shí)，為了保證飛行器平穩(wěn)飛行，選取常值傾側(cè)角σ2作為控制變量。在度過初始下降段后，傾側(cè)角保持線性遞增或遞減規(guī)律飛行一段時(shí)間，而后保持常值傾側(cè)角σ1，調(diào)整飛行器的再入橫縱向飛行偏差，而后保持傾側(cè)角以線性遞增或遞減規(guī)律飛行。采取這樣的傾側(cè)角剖面，可以實(shí)現(xiàn)航向角的精準(zhǔn)調(diào)節(jié)。參數(shù)化傾側(cè)角剖面如圖2所示。

通過對(duì)傾側(cè)角剖面的參數(shù)化設(shè)計(jì)，再入軌跡優(yōu)化問題轉(zhuǎn)化為參數(shù)優(yōu)化問題，從而可用含有多約束的非線性優(yōu)化方法直接求解，且待求解的未知參數(shù)僅包括(σ0,σ1,σ2,V1,V2,V3,V4)，從而減小參數(shù)優(yōu)化搜索空間，提高計(jì)算效率。

圖2 參數(shù)化傾側(cè)角剖面Fig.2 Profile illustrations of parameterized bank angles

2.2 差分進(jìn)化算法

差分進(jìn)化算法于1997年由Storn和Price在遺傳算法等進(jìn)化思想的基礎(chǔ)上提出，本質(zhì)是一種多目標(biāo)(連續(xù)變量)優(yōu)化算法，用于求解多維空間中整體最優(yōu)解。差分進(jìn)化思想的來源是遺傳算法，模擬遺傳學(xué)中的雜交、變異、復(fù)制來設(shè)計(jì)遺傳算子。不同之處在于遺傳算法是根據(jù)適應(yīng)度值來控制父代雜交，變異后產(chǎn)生的子代被選擇的概率值，在最大化問題中適應(yīng)值大的個(gè)體被選擇的概率相應(yīng)也會(huì)大一些。而差分進(jìn)化算法變異向量是由父代差分向量生成，并與父代個(gè)體向量交叉生成新個(gè)體向量，直接與其父代個(gè)體進(jìn)行選擇。這種選擇更接近實(shí)際且逼近效果也更好。

定義{x1,x2,…,xD}為待優(yōu)化的D個(gè)未知參數(shù)，且滿足:

基于差分進(jìn)化算法，實(shí)現(xiàn)參數(shù)優(yōu)化的流程為：

1)初始化種群。隨機(jī)產(chǎn)生初始種群滿足

(10)

2)變異操作。通過差分策略實(shí)現(xiàn)種群個(gè)體變異，生成中間體，即：

vi(k+1)=xr1(k)+λ(xr2(k)-xr3(k))

(11)

式中：vi為變異操作生產(chǎn)的隨機(jī)個(gè)體；r1，r2，r3為從第k代種群中隨機(jī)選擇的三個(gè)個(gè)體，i≠r1≠r2≠r3；λ為變異因子。

3)交叉操作。對(duì)第k代種群{xi(k)}及其隨機(jī)個(gè)體{vi(k)}進(jìn)行交叉操作：

(12)

式中：uj,i為通過概率方式生產(chǎn)的隨機(jī)個(gè)體；cr為交叉概率；jrand為[1,2,…,D]的隨機(jī)整數(shù)。

4)選擇操作。通過貪婪算法選擇進(jìn)入下一代的種群個(gè)體：

xi(k+1)=

(13)

軌跡優(yōu)化問題包含復(fù)雜的終端約束和過程約束，有效處理再入過程的等式和不等式約束是保證再入優(yōu)化軌跡可行的關(guān)鍵。本文將這些約束條件以罰函數(shù)形式，添加到適應(yīng)度函數(shù)中，為：

(14)

式中：ωm≥0(m=1,2,…,p)為罰函數(shù)權(quán)重；?m(x)為與優(yōu)化參數(shù)x相關(guān)的等式約束。

對(duì)于不等式約束，包括動(dòng)壓約束式(2)、熱流密度約束式(3)或過載約束式(4)，當(dāng)這些強(qiáng)約束條件中的任一個(gè)超出給定的最大范圍時(shí)，則將該個(gè)體對(duì)應(yīng)的適應(yīng)度函數(shù)J進(jìn)行極大化，從而保證通過父本得到的新個(gè)體能夠滿足硬約束條件。

2.3 傾側(cè)角插值模型

基于傾側(cè)角插值模型的再入可達(dá)區(qū)域計(jì)算可歸結(jié)為：1)求解兩條具有最大橫向航程的再入滑翔軌跡；2)將最大橫向航程軌跡對(duì)應(yīng)的傾側(cè)角剖面指定為線性插值模型的基準(zhǔn)剖面，并選取合適的插值模型系數(shù)，計(jì)算傾側(cè)角插值剖面集合；3)將傾側(cè)角剖面集合中的控制指令代入三自由度再入運(yùn)動(dòng)學(xué)模型，生成具有不同機(jī)動(dòng)能力的再入軌跡集合，完成再入可達(dá)區(qū)域的近似計(jì)算。

設(shè)計(jì)傾側(cè)角線性插值模型為：

σ(V)=σup(V)+η(σdown(V)-σup(V))

(15)

式中:σup(V)和σdown(V)分別為給定的基準(zhǔn)傾側(cè)角剖面;η∈[0, 1]表示傾側(cè)角插值模型系數(shù)。

通過選取一定間隔的插值模型系數(shù)η(如η=0.2, 0.4, …, 0.8)，能夠直接生成一組具有不同幅值和不同臨界速度的傾側(cè)角剖面集合,如圖4所示，當(dāng)傾側(cè)角系數(shù)分別取η= 0和η= 1時(shí)，對(duì)應(yīng)的傾側(cè)角剖面即為插值模型式(15)的上邊界σup(V)和下邊界σdown(V)。

3 仿真分析

飛行器氣動(dòng)參數(shù)[19]可近似為：

CL=-0.234 2+0.051 36α+0.294 3e-0.100 7Ma

CD=0.024 67+0.000 714 3α2+0.325 2e-0.279 0Ma

給定再入攻角剖面為再入初始攻角α0= 30°，臨界速度VT=4 500 m/s，kα=1.429 837 9×10-6，傾側(cè)角取值范圍約束為-60° ～ 60°，初始下降段選取常值傾側(cè)角為σ0= 0°。

3.1 差分進(jìn)化算法再入軌跡求解方案

給定5組不同的高超聲速飛行器再入初始條件，如表1所示，并指定期望再入終端約束為hf=24 km，Vf=1 400 m/s，θf= 45°，φf= 30°。

表1 飛行器再入初始條件

利用差分進(jìn)化算法，針對(duì)表1中給出的初始條件，可求解滿足終端約束、過程約束的傾側(cè)角剖面，即參數(shù)：(σ0,σ1,σ2,V1,V2,V3,V4)。

圖5給出了不同再入條件下飛行器的經(jīng)緯度軌跡，5組仿真結(jié)果都滿足預(yù)定的經(jīng)緯度精度要求。圖6給出了求解得到的傾側(cè)角剖面。

由圖6可見，組4和組5相較于組1～組3多進(jìn)行了一次傾側(cè)角反轉(zhuǎn)，用于調(diào)整飛行器航向角。參數(shù)化的傾側(cè)角剖面求解結(jié)果如表2所示。

圖5 再入經(jīng)緯度軌跡Fig.5 Reentry longitude and latitude trajectory

圖6 再入傾側(cè)角剖面Fig.6 Reentry bank angle profile

表2 參數(shù)化傾側(cè)角剖面求解結(jié)果

3.2 再入可達(dá)域計(jì)算

考慮最大航向航程再入軌跡，選取性能指標(biāo)為：

minJ=

并選取飛行器再入初始條件如表3所示，給定終端約束為hf= 24 km，Vf= 1 500 m/s。

表3 飛行器再入初始條件

線性插值模型參數(shù)η在區(qū)間[0, 1]上以0.1間隔依次取值，在得到飛行器最大橫向航程的基礎(chǔ)上，生成再入傾側(cè)角指令集，如圖7所示。

如圖8所示,ηr= 1和ηl= 1分別對(duì)應(yīng)了高超聲速飛行器最大橫向航程時(shí)的軌跡，且對(duì)應(yīng)于圖7中的ηr= 1和ηl= 1表示的傾側(cè)角剖面。以η= 0和ηr= 1以及η= 0和ηl= 1分別構(gòu)建出式(15)中的下邊界σdown(V)和上邊界σup(V)，根據(jù)η取值生成傾側(cè)角插值剖面，將其代入飛行器運(yùn)動(dòng)學(xué)模型，可計(jì)算得到飛行器的再入可達(dá)域。

圖7 傾側(cè)角插值剖面Fig.7 Bank angle interpolation profile

圖8 再入可達(dá)域Fig.8 Reentry landing footprint

圖9給出了本文方法與Gauss偽譜法求解再入可達(dá)域的對(duì)比。由圖9可見，本文方法和偽譜法得到的再入可達(dá)域的外邊界形狀基本保持一致，但略小于偽譜法得到的優(yōu)化結(jié)果。偽譜法需要逐條計(jì)算和優(yōu)化再入軌跡，插值法僅需要一條最大橫向航程軌跡則可近似計(jì)算。

圖9 可達(dá)域?qū)Ρ菷ig.9 Comparisons of landing footprints

4 結(jié)束語

本文針對(duì)升力式高超聲速飛行器再入可達(dá)域問題，提出帶有罰函數(shù)的差分進(jìn)化算法和傾側(cè)角插值相結(jié)合的求解方案，研究分析和仿真結(jié)果表明：

1)帶有罰函數(shù)的差分進(jìn)化算法可簡化傾側(cè)角剖面，從而減小搜索空間，提高求解效率；

2)傾側(cè)角插值模型直接生成指令集，快速計(jì)算得到再入可達(dá)域具有保守型，在需要快速估計(jì)飛行器能力時(shí)，給出的可達(dá)域具有參考性。