高超聲速飛行器多約束再入軌跡快速優(yōu)化

梅映雪，馮 玥，王容順，吳了泥，孫洪飛

(廈門大學航空航天學院，廈門 361102)

0 引 言

高超聲速飛行器在戰(zhàn)術(shù)意義上以其極高的飛行速度和大跨度的飛行空域，占據(jù)著得天獨厚的優(yōu)勢，目前已經(jīng)日益受到世界各國的普遍重視[1]。黃長強等[2]提出了高超聲速飛行器軌跡優(yōu)化的重點研究方向，根據(jù)兵力投送和長航程的實際作戰(zhàn)需求，勢必要考慮多預(yù)設(shè)區(qū)到達和禁飛區(qū)規(guī)避，以及如何提高軌跡優(yōu)化的計算效率等問題。因此有必要根據(jù)實際需求尋求一條適合再入滑翔飛行且高超聲速到達的最優(yōu)軌跡。

高超聲速飛行器的軌跡優(yōu)化涉及眾多非線性約束條件，是一類復(fù)雜的非線性多約束最優(yōu)控制問題[3]，其求解具有相當?shù)奶魬?zhàn)。自二十世紀后期,出現(xiàn)了多種飛行器軌跡優(yōu)化方法[4]，Gauss偽譜法作為一種基于全局插值多項式的直接配點法，它相對于一般直接配點法的優(yōu)勢是可以用較少的節(jié)點獲得較高的精度[5-7]，且計算效率較高，因此受到研究者的青睞。Reddien[8]使用Gauss偽譜法求解最優(yōu)控制問題，利用全局正交多項式在特定節(jié)點同時離散狀態(tài)量和控制量。李鐵鵬等[9]基于Gauss偽譜法以飛行時間為優(yōu)化指標對滑翔彈道進行了優(yōu)化設(shè)計。和爭春等[10]在過載限制條件下對升力體高超聲速飛行器再入彈道進行了優(yōu)化設(shè)計。上述文獻多將重點放在如何優(yōu)化得到高超聲速飛行器在各種常規(guī)約束條件下到達指定點的軌跡上，卻很少考慮飛行器在實際任務(wù)中必須面對的經(jīng)過特定任務(wù)區(qū)域并順利避開敵方攔截或勿入?yún)^(qū)域的場景。對于繞飛區(qū)這一導(dǎo)彈防御系統(tǒng)的挑戰(zhàn)，大都采用改變傾側(cè)角符號的方法單獨規(guī)劃橫向軌跡。張科南等[11]基于一定的傾側(cè)角變化規(guī)律優(yōu)化得到了高超聲速滑翔飛行器規(guī)避攔截區(qū)的再入軌跡，但計算效率和魯棒性均不高。Shen等[12]所述方法能夠快速計算出一條滿足各種過程約束和狀態(tài)約束的三維再入軌跡，但其橫向軌跡設(shè)計部分僅通過一次改變傾側(cè)角的符號進行控制，不太適用于大范圍橫向機動情況。也有學者對Gauss偽譜方法進行改進，以適用于此類考慮測控區(qū)和繞飛區(qū)的復(fù)雜多約束問題。Darby等[13]提出了一種多重區(qū)間配置法對連續(xù)時間非線性優(yōu)化問題進行研究，但未針對具體問題。Wang等[14]采用多段Gauss偽譜法通過將航路點和繞飛區(qū)分別設(shè)置為固定LG點和虛擬LG點的方法對此類問題進行研究，但設(shè)定LG點的方式無法使軌跡得到全局優(yōu)化。Jorris等[15]提出了一種實時自主優(yōu)化方法使得軌跡在滿足航路點和繞飛區(qū)等復(fù)雜多約束的同時再入時間最短，但未考慮測控區(qū)以及軌跡平滑性這樣的實際問題。因此，在如此多約束條件下進行軌跡的優(yōu)化設(shè)計是一個比較復(fù)雜的問題，需針對軌跡優(yōu)化問題的特點，采用合適的優(yōu)化策略和方法。

考慮到滑翔飛行方式可以規(guī)避攔截和進行防御，而飛行器再入初期不滿足滑翔條件，通常將再入軌跡分為初始下降段和準平衡滑翔段分別規(guī)劃。本文在此基礎(chǔ)上，考慮到實際任務(wù)中所必須面對的經(jīng)過特定任務(wù)區(qū)域并順利避開敵方攔截或勿入?yún)^(qū)域的場景，對高超聲速滑翔飛行器軌跡優(yōu)化的主要約束條件進行了分析并給出相應(yīng)的數(shù)學模型，然后對Gauss偽譜方法在求解復(fù)雜多約束條件下軌跡優(yōu)化問題時存在的主要難點進行了分析，提出了具體的軌跡分段優(yōu)化策略和方法，較好地解決了這一問題。為了便于控制系統(tǒng)操控，保證優(yōu)化效率，同時為了避免約束條件太多而減小可行域，在性能指標中引入彈道傾角和航向角相關(guān)指標的加權(quán)和代替控制角速率等相關(guān)過程約束。最后通過對比，分析了本文優(yōu)化方法的快速性和工程實用性。

1 問題描述和預(yù)備知識

1.1 問題描述

高超聲速飛行器再入段由再入點到滑翔段終端點，軌跡規(guī)劃過程需要考慮動力學方程約束，過程約束，終端約束等常規(guī)約束，還要考慮更符合實際任務(wù)的測控區(qū)約束和繞飛區(qū)約束。

1) 模型約束

(1)動力學方程約束

高超聲速飛行器在高速大氣層內(nèi)無動力返回，由于高升阻比飛行器氣動力和地球重力分別約為科氏慣性力和地球自轉(zhuǎn)引起的慣性離心力的100倍和1000倍[16]，因此可以假設(shè)地球為圓球不旋轉(zhuǎn)模型，高超聲速飛行器無量綱三自由度運動方程[17]為：

(1a)

(1b)

(1c)

(1d)

(1e)

(1f)

無量綱升力L、阻力D分別為：

(2)

式中：Sref為飛行器氣動參考面積；CL,CD分別為升力系數(shù)和阻力系數(shù)，由飛行器的攻角α和Ma決定；ρ為大氣密度。

(2)過程約束

為保證再入飛行器在結(jié)構(gòu)和熱防護上的可靠性，再入過程要求嚴格滿足駐點熱流密度、動壓和法向過載約束[18]，即:

(3a)

(3b)

(3c)

(3)控制約束

以飛行攻角α和傾側(cè)角υ為控制變量，為了滿足飛行器控制能力要求，有如下約束條件[18]：

(4)

式中：μ1,μ2為攻角取值范圍的上、下限；κ為傾側(cè)角取值范圍的上限。

(4)終端約束

為保證在一定的落角(終端彈道傾角)和速度下將飛行器引導(dǎo)至指定終端目標位置，通常要加上如下終端約束[19]：

(5)

式中：tf為終端時刻。

以上是高超聲速飛行器軌跡優(yōu)化常見的約束條件，本文考慮到工程上的實際需要，在常規(guī)約束基礎(chǔ)上還要考慮測控區(qū)約束和繞飛區(qū)約束。測控區(qū)和繞飛區(qū)約束描述如下：

Θ(λ,φ;λc,φc,Rc)d(λc,φc),Γ)-Rc<0

(6)

Ξ(λ,φ;λr,φr,Rr)d(λr,φr),Γ)-Rr≥0

(7)

2) 性能指標

在高超聲速飛行器再入軌跡優(yōu)化的研究中，多數(shù)文獻僅僅考慮某個單一的性能指標。例如文獻[20]，為了實現(xiàn)快速打擊，要求再入飛行時間最短，所以將再入過程的時間作為優(yōu)化目標：

J1=tf-t0

(8)

式中：t0為再入初始時刻。

為了兼顧到再入軌跡的平滑性，常常將彈道傾角速率的平方的積分作為性能指標[21]：

(9)

再者，為了便于控制系統(tǒng)操控，使飛行器平穩(wěn)飛行，也有學者采用航向角速率的平方的積分作為性能指標[22]：

(10)

本文將綜合考慮所設(shè)計軌跡的高效性和工程實用性，選擇式(8)，式(9)，式(10)所示性能指標的加權(quán)和作為優(yōu)化目標，通過調(diào)整權(quán)系數(shù)，能夠在時間最優(yōu)、控制量最優(yōu)和軌跡平滑性最優(yōu)之間進行權(quán)衡，以規(guī)劃出易于實現(xiàn)的軌跡。最終的性能指標式為：

J=w1J1+w2J2+w3J3

(11)

式中：w1,w2,w3為權(quán)重系數(shù)，用于調(diào)節(jié)優(yōu)化指標中飛行器飛行時間以及彈道和控制量平滑性最優(yōu)之間的權(quán)重。

注1. 對于性能指標的分析：從物理機理上，采用Gauss偽譜法將再入時間(tf-t0)作為目標函數(shù)，即要求再入時間最短，同時在性能指標中引入彈道傾角和航向角相關(guān)指標代替相關(guān)過程約束使縱向軌跡和橫向軌跡平滑，由于減少了約束條件數(shù)量，計算效率進一步提高；從控制機理上，w1,w2,w3分別為性能指標J第1、2、3項的權(quán)，當權(quán)越大，懲罰越大，對應(yīng)項的函數(shù)值越小，即對應(yīng)項越優(yōu)。對于線性系統(tǒng)，權(quán)的選取可借鑒LQR方法中加權(quán)矩陣的選取[23-25]；而對于非線性系統(tǒng)，目前更多依賴于工程實際。為了方便后續(xù)制導(dǎo)工作，相對于軌跡平滑性，再入時間是次要的，因此w2、w3都要較w1取大一些。

令x,u分別表示系統(tǒng)(1)的狀態(tài)和輸入，即：

x=[x1,x2,x3,x4,x5,x6]T=[r,λ,φ,V,θ,ψ]T；u=[u1,u2]T=[α,υ]T；再令f(x,u)為狀態(tài)方程(1a)～(1f)的右端函數(shù)。

則本文所要研究的多約束問題可描述為：

在滿足動力學方程約束(1a)～(1f)，過程約束(3a)～(3c)，控制約束(4)，終端約束(5)以及測控區(qū)約束(6)和繞飛區(qū)約束(7)的情況下，尋找合適的狀態(tài)軌跡和控制輸入，以及初末時刻(若未給定)，使得性能指標(11)最小，即：

(12)

式中：W=6,M=2；Ωζ(x,u)≤0為由過程約束(3a)～(3c)、控制約束(4)確定的q個不等式約束；Λ(x(tf))≤0為由終端狀態(tài)約束(5)確定的等式及不等式約束。

上述多約束軌跡優(yōu)化問題本質(zhì)上是最優(yōu)控制問題，工程上普遍使用GPM來求解最優(yōu)控制問題的數(shù)值解[26]。但是，由于測控區(qū)約束為局部性約束，而GPM僅對具有過程約束的最優(yōu)控制問題有效。因此，測控區(qū)甚至是多個測控區(qū)約束的存在加大了求解軌跡優(yōu)化問題的難度，本文將著重解決該問題。

1.2 預(yù)備知識

Gauss偽譜法以Legendre多項式的根為離散點，將連續(xù)最優(yōu)控制問題的狀態(tài)變量和控制變量離散化，并以離散點為節(jié)點采用全區(qū)間Lagrange插值多項式來近似狀態(tài)變量和控制變量，從而將軌跡優(yōu)化的最優(yōu)控制問題轉(zhuǎn)換為非線性規(guī)劃問題進行求解[27]。

在不考慮測控區(qū)和繞飛區(qū)的情況下，問題可以直接采用GPM求解。對于GPM的一般步驟，本文作如下歸納：

1) 對時間變量t作變換，將時間區(qū)間轉(zhuǎn)換到[-1,1]。

(13)

2)確定離散時刻點。

取N階Legendre多項式的N個根τk∈(-1,1),k=1,2,…,N，再令τ0=-1,τN+1=1，則τ0,τk,τN+1構(gòu)成彈道的N+2個離散點。

3)約束條件及性能指標的離散化。

根據(jù)上述時間變換和配點選取，構(gòu)造系統(tǒng)狀態(tài)x和輸入u關(guān)于其估計值X,U的近似表達式：

(14)

令Xk=X(τk),Uk=U(τk)；J=w1J1+w2J2+w3J3φ(t0,tf)+η(x(t),u(t),t)dt，其中，則微分方程約束，等式/不等式約束以及目標函數(shù)經(jīng)離散化，可轉(zhuǎn)換為下述代數(shù)約束[28-29]：

(15)

Ωζ(Xk,Uk)≤0

(16)

Λ(X(tf))≤0

(17)

(18)

經(jīng)過對時間、非線性函數(shù)、微分方程等離散化后，最優(yōu)控制問題(12)(不考慮測控區(qū)和繞飛區(qū)約束)就轉(zhuǎn)化為如下離散非線性規(guī)劃問題：

求離散節(jié)點上的狀態(tài)X(τi)(i=0,…,N)和控制變量U(τi)(i=1,…,N),使得性能指標(18)最小，并滿足動力學方程約束(15)，過程約束等不等式約束(16)以及終端狀態(tài)約束(17)。

2 軌跡優(yōu)化策略

目前許多研究采用Gauss偽譜法解決高超聲速飛行器軌跡的優(yōu)化設(shè)計問題，但在面對本文提出的多約束任務(wù)場景時，考慮到測控區(qū)約束的局部性特點，無法直接求解。

2.1 多約束分段優(yōu)化方法

假設(shè)1. 飛行任務(wù)中有l(wèi)個繞飛區(qū)，n個測控區(qū)。飛行器自東向西跨經(jīng)度飛行且測控區(qū)中心軸處在不同的經(jīng)度上。同時假設(shè)所有測控區(qū)均在準平衡滑翔段。假設(shè)1具有一般性，對測控區(qū)位置的假定是方便對測控區(qū)自東向西進行編號。若部分測控區(qū)的中心軸處于同一經(jīng)度，則以該經(jīng)度對軌跡進行分段，然后按照緯度對處于相同經(jīng)度的測控區(qū)編號。不妨假設(shè)，經(jīng)過編號后的測控區(qū)按照自東向西的方向依次為測控區(qū)1、測控區(qū)2, ……，測控區(qū)n。

對于繞飛區(qū)約束和測控區(qū)約束，將問題簡化，只考慮其橫向剖面，化三維約束為繞飛圓約束和測控圓約束進行求解。則規(guī)避繞飛區(qū)，進入測控區(qū)兩個約束條件可分別嚴格描述為：

(19)

(20)

式中：下標ri和cj分別表示第i個繞飛區(qū)和第j個測控區(qū)。

注2. 本文分段點個數(shù)取決于測控區(qū)個數(shù)，自然地，測控區(qū)越多，分段數(shù)越多，對數(shù)值計算要求越高，但計算時間可接受，并不會帶來明顯的增長。

由式(19)可知，繞飛區(qū)約束是一個全局性約束條件，即：對于軌跡上所有點(λ,φ),都滿足Ξ(λ,φ;λri,φri,Rri)≥0。因此，可以直接將其設(shè)置為過程約束。則規(guī)避第i個繞飛區(qū)的約束為Ξi(λ,φ;λri,φri,Rri)≥0。

由此可見，繞飛區(qū)約束作為全局性約束條件不會給求解帶來本質(zhì)上的困難。由式(20)可知，測控區(qū)約束僅僅是對軌跡的局部性約束，此時不能將其作為過程約束來看待，這給軌跡優(yōu)化問題帶來了本質(zhì)困難。針對該棘手問題，本文將改進傳統(tǒng)Gauss偽譜法，以適用于再入飛行任務(wù)含測控區(qū)的軌跡優(yōu)化問題。具體來說，就是依據(jù)n個測控區(qū)將準平衡滑翔段分成n+1段，這樣連同初始下降段，相當于將整個軌跡分割成n+2段(如圖1所示)，然后每一段軌跡通過Gauss偽譜法來求解。為方便，自東向西依次稱每一小段軌跡為第1段、第2段，……，第n+2段。

將飛行軌跡按上述方式進行分段后，最主要的問題是分段點如何選取。分段點選取首先要考慮飛行器在每一段上的飛行要求和約束，還要兼顧到Gauss偽譜法的適用性。下面給出分段點的選取依據(jù)。

分段點1為初始下降段與準平衡滑翔段的交界點。該點可以通過平衡滑翔條件以及最大熱流密度約束(再入初期動壓和過載一般較小，可以滿足約束，因此不考慮)來決定，即滿足條件[29]：

(21)

分段點2～(n+1)根據(jù)測控區(qū)約束選取，將第j個測控區(qū)約束設(shè)置為第j+1段軌跡的終端約束。根據(jù)測控區(qū)位置及大小，動態(tài)設(shè)置軌跡上點Aj為分段點。設(shè)Aj點所在經(jīng)緯度坐標為(λAj,φAj)，若Aj點位于測控區(qū)內(nèi)部，即：進入第j個測控區(qū)的約束為Θj(λAj,φAj;λcj,φcj,Rcj)<0，則軌跡成功經(jīng)過測控區(qū)域，符合任務(wù)要求。

令*(p)表示第p段彈道對應(yīng)參數(shù)，p∈{1,2,…,n+2}。

在每段軌跡上選擇合適數(shù)量的離散點將狀態(tài)變量和控制變量進行離散化，將整個軌跡分段離散化后，將所有的離散點狀態(tài)變量、控制變量作為優(yōu)化參數(shù)同時進行優(yōu)化設(shè)計，最終通過插值獲得對應(yīng)的最優(yōu)控制及最優(yōu)軌跡，第p段彈道優(yōu)化問題描述如下:

當然，為確保狀態(tài)量和控制量的連續(xù)性，在軌跡優(yōu)化過程中各段軌跡在分段點處需要滿足如下連接條件：

(22)

2.2 優(yōu)化流程

本文采用GPM作為每階段優(yōu)化的基本方法。根據(jù)上述優(yōu)化策略，為充分利用再入軌跡不同階段特性,將再入軌跡分為初始下降段和滑翔段分別優(yōu)化，同時根據(jù)測控區(qū)約束條件將滑翔段分段，將一般最優(yōu)控制問題轉(zhuǎn)換為多段最優(yōu)控制問題，以求得滿足多約束條件且具有一定精度的最優(yōu)軌跡。圖2描述了本文提出的分段優(yōu)化策略的基本流程。

3 仿真結(jié)果及分析

在MATLAB中進行仿真，仿真的目的是檢驗任意給定測控區(qū)和繞飛區(qū)的情況下，所設(shè)計的方法是否能夠快速得到優(yōu)化軌跡使飛行器從再入點快速平滑地過渡到終端點，且成功進入測控區(qū)并規(guī)避所設(shè)繞飛區(qū)；同時就優(yōu)化效果與已有研究進行對比，檢驗本文所設(shè)計的改進方法的高效性和實用性。

3.1 算例

以通用航空器CAV為對象，再入飛行器氣動特性參數(shù)采用通用航空器CAV的相關(guān)參數(shù)[30]，其最大升阻比為3.4，氣動參考面積為0.8 m2，質(zhì)量為900 kg，最大升阻比攻角為12°。

根據(jù)實際要求，仿真參數(shù)設(shè)置如下：

假設(shè)飛行器要進入1個測控區(qū)并繞過3個繞飛區(qū)，即，n=1,p=3。

給定測控區(qū)1中心位于(136.11°E，44.74°N),半徑為100 km；繞飛區(qū)1中心位于(140.03°E，35.53°N)，半徑為600 km；繞飛區(qū)2中心位于(123.91°E，53.98°N)，半徑為500 km；繞飛區(qū)3中心位于(116.37°E，39.690°N)，半徑為1000 km。

性能指標的權(quán)重系數(shù)取為w1=0.2,w2=0.4,w3=0.4。

圖3～6為采用本文提出的改進Gauss偽譜法對算例進行的有效性驗證曲線。由圖3可知，再入過程中熱流密度、法向過載以及動壓都嚴格低于最大峰值，滿足常規(guī)過程約束。圖4為高度、速度關(guān)于時間變化的曲線，從圖中可以看出，曲線平滑且滿足終端狀態(tài)約束等條件。圖5表明控制變量在要求范圍內(nèi)變化，且變化較為平緩。但根據(jù)實際工程需要，由于考慮攻角幅度限制，導(dǎo)致出現(xiàn)了兩次較長時間的攻角飽和。

步增加，α先上升維持飽和值30°使得V先下降以滿足測控/繞飛區(qū)約束以及終端高度等約束，后下降以滿足終端速度等約束，同時υ需配合α呈先上升后下降的趨勢以滿足測控區(qū)約束等多約束條件。

圖6為軌跡的橫側(cè)向和三維軌跡圖，目的是便于觀察軌跡的橫側(cè)向走向，圖6(a)、圖6(b)表明軌跡成功進入測控區(qū)并且有效規(guī)避了繞飛區(qū)，最后到達終端點，且在分段點處狀態(tài)量和控制量平滑銜接。

為驗證測控區(qū)/繞飛區(qū)數(shù)目對模型維數(shù)的影響，在上述算例的基礎(chǔ)上增加兩個測控區(qū)，中心分別位于(130.11°E，40.74°N)、(126.04°E，32.83°N),半徑都為100 km，分別記為測控區(qū)2和測控區(qū)3，統(tǒng)計計算時間，如表1、表2、表3。

表1 測控區(qū)個數(shù)對計算時間的影響Table 1 Influence of the number of monitoring zones on the calculation time

表2 繞飛區(qū)個數(shù)對計算時間的影響Table 2 Influence of the number of no-fly zones on the calculation time

表3 繞飛區(qū)分布/大小對計算時間的影響Table 3 Influence of the distribution of no-fly zones around the calculation time

從表3中可以看出：1)測控區(qū)個數(shù)的增加會延長計算時間t；2)繞飛區(qū)個數(shù)對計算時間有影響，但不顯著，繞飛區(qū)的分布及大小對計算時間的影響更大。

這說明由于本文分段點個數(shù)取決于測控區(qū)個數(shù)、自然地、測控區(qū)越多，分段數(shù)越多，對數(shù)值計算要求越高，但計算時間可接受，并不會帶來明顯的增長；而繞飛區(qū)只是作為過程約束，繞飛區(qū)的數(shù)目確實對模型的維數(shù)造成了影響，但計算時間更大程度上會受到繞飛區(qū)分布及其大小的影響，若分布過于緊密或者繞飛區(qū)較大，軌跡曲度可能發(fā)生較大變化，對在線規(guī)劃勢必會造成影響。

文獻[27]考慮到快速打擊和軌跡平滑性，引入再入時間和彈道傾角相關(guān)指標的加權(quán)和作為目標函數(shù)。本文在此基礎(chǔ)上考慮到軌跡快速優(yōu)化要求，力求控制量變化平緩，同時為了避免約束條件增多減小可行域，在目標函數(shù)中引入航向角相關(guān)指標代替在過程約束中增加控制角速率約束，取再入時間、彈道傾角和航向角相關(guān)指標三項指標加權(quán)和作為目標函數(shù)。圖7至圖8為采用文獻[27]中性能指標與本文改進性能指標下對原算例(n=1,p=3)進行的優(yōu)化對比曲線。圖7表明，兩種情況下常規(guī)過程約束都能被控制在峰值之內(nèi)，但改進方法使得熱流率、法向過載和動壓在再入過程中明顯低于未改進方法，從而減小了飛行器承受能力，提高了飛行器性能和使用壽命。圖8表明改進方法下的控制變量變化曲線較未改進方法平滑，波動明顯更小，從而便于控制系統(tǒng)操控，極大地提高了計算效率。表4為采用兩種優(yōu)化方法所需計算時間的對比表，改進方法計算效率較未改進方法大大提高，更加滿足再入飛行任務(wù)快速規(guī)劃的需求。

計算方法計算時間/s改進38.974302未改進453.015268

4 結(jié) 論

本文針對高超聲速飛行器再入過程中面臨多個測控區(qū)和繞飛區(qū)的再入軌跡設(shè)計問題，分析了采用傳統(tǒng)Gauss偽譜法設(shè)計時存在的主要問題，引入了分段優(yōu)化的軌跡設(shè)計思路，并提出了具體的軌跡分段策略。以再入時間、彈道傾角以及航向角相關(guān)指標的加權(quán)和最小為性能指標設(shè)計了考慮各種約束條件下的最優(yōu)軌跡。仿真結(jié)果表明本文所設(shè)計的改進方法能夠快速有效的用于解決任意給定測控區(qū)和繞飛區(qū)約束的軌跡優(yōu)化問題。