一種高超聲速飛行器再入軌跡優化方法

余 躍，王宏倫

(1. 北京航空航天大學自動化科學與電氣工程學院, 北京 100191; 2. 北京航天自動控制研究所, 北京 100854; 3. 北京航空航天大學飛行器控制一體化技術重點實驗室 , 北京 100191)

0 引 言

高超聲速飛行器因具有飛行速度快、突防能力高、生存能力強、全球精確打擊和毀傷威力大等獨特優勢，已成為當今世界各國爭相研制的武器[1]。再入軌跡優化是高超聲速飛行器研究中的主要關鍵技術，是指根據飛行任務的飛行條件和技術指標，基于飛行器動力學模型求解控制量與運動狀態的最優時間序列，從而尋找一條某種性能指標最優而又滿足各種約束的飛行軌跡[2-3]。再入環境的復雜性和不確定性以及再入過程中受熱流密度、動壓、過載的嚴格約束，再加上對終端狀態的嚴格要求，給再入軌跡優化技術的研究提出了挑戰[4]。

按照是否直接對性能指標進行尋優，可將軌跡優化方法分為間接法和直接法。間接法是基于經典的變分法或者Pontryagin極小值原理，將最優控制問題轉化為哈密爾頓兩點邊值問題進行求解。與間接法不同，直接法是將最優控制問題中的變量離散并參數化，從而將最優控制問題轉化為非線性規劃(NLP)問題，再結合數值優化方法進行求解。直接法中的偽譜法采用全局正交多項式逼近狀態變量和控制變量，具有精度高且收斂性好的優點，近年來迅速成為飛行器軌跡優化領域應用最廣泛的研究方法[5]。

雍恩米等[6]針對再入軌跡優化問題，提出了基于高斯偽譜法(GPM)的串行分段優化策略，具有較高的精度和計算效率。宗群等[7]針對臨近空間飛行器上升段的軌跡優化問題，提出了將GPM和序列二次規劃相結合的求解策略。明超等[8]針對超聲速導彈爬升段的軌跡優化設計問題，利用hp自適應偽譜法構造了一種較為通用的優化求解策略。張志國等[9]將GPM運用在液體運載火箭拋罩結束到入軌飛行段的制導律設計中。盡管偽譜法具有精度高、收斂快和應用靈活的優點，但仍存在以下缺陷：1)需要目標函數的梯度值；2)相對差的全局搜索能力；3)需要合適的初始猜測值[10-11]。當面對復雜的軌跡優化問題時，對初始猜測值的敏感可能會導致沒法得到全局最優解。

近年來，各種啟發式優化算法不斷被提出，如粒子群算法(PSO)[12]、鯨魚算法(WOA)[13]、灰狼算法(GWO)[14]等，這些算法不依賴目標函數的梯度信息、對初始猜測值不敏感和具備強大的全局搜索能力，被廣泛應用于各類優化問題。樽海鞘群算法(SSA)[15]是最新提出的一種啟發式算法，通過模仿樽海鞘的捕食行為來求得全局最優解，從優化的統計數據看，SSA算法與當前頂級的啟發式算法相比具有優勢。

基于上述分析，為了充分利用啟發式優化算法和偽譜法的優勢，本文首先針對SSA算法中探索和利用的平衡過程加以改進，提出了一種新穎的動態自適應樽海鞘群算法(Dynamic Adaptation SSA, DASSA)，將DASSA算法和GPM算法相結合，提出了一種新穎的DASSAGPM算法，并運用于高超聲速飛行器再入軌跡優化問題中。DASSAGPM算法既能利用DASSA算法強大的全局搜索能力，又能利用GPM算法收斂速度快的優點，具有滿意的優化性能。

1 再入軌跡優化問題

1.1 三自由度運動方程

不考慮地球自轉的影響，建立高超聲速飛行器三自由度無量綱運動方程如下[6]：

(1)

(2)

式中：K=0.5R0S/m；ρ為大氣密度，CL和CD分別為升力系數和阻力系數，S為飛行器參考面積，m為飛行器質量。

1.2 約束條件

1.2.1過程約束

(3)

1.2.2控制邊界

為使飛行器穩定飛行，控制量攻角α和傾側角σ需滿足如下條件：

(4)

1.2.3終端約束

為了滿足與末制導段的交接班要求，再入終端狀態需要滿足終端約束條件，一般為地心距和經緯度組成的位置約束，某些飛行器對速度、航跡傾角和航跡偏角也有終端約束要求。用公式表示為：

(5)

式中:下標“f”表示終端狀態。

1.3 目標函數

通常情況下，目標函數根據特定任務選定。對高超聲速飛行器來說，可選為最小化熱載、最大化航程、最大化橫程或者最小化到達時間等。本文選擇最小化熱載為優化目標，同時考慮彈道平滑的因素，取綜合性能指標為：

(6)

式中：c為加權系數，本文優化的主要指標為熱載，故取c=0.8。

1.4 問題描述

基于上述分析，再入軌跡優化問題可描述為：在滿足動力學方程(1)以及多種約束(3)～(5)的條件下，尋找最優控制量攻角α和傾側角σ使得目標函數(6)最小。更為一般化的最優控制問題可描述為：尋找控制變量u(t)，使得具有一般性的Bolza型性能指標最小

J=Φ(x(τ0),t0,x(τf),tf)+

(7)

且滿足動力學微分方程約束

(8)

邊界條件約束

φ(x(τ0),t0,x(τf),tf)=0

(9)

以及不等式約束

C(x(τ),u(τ),τ;t0,tf)≤0

(10)

2 樽海鞘群算法和高斯偽譜法

2.1 樽海鞘群算法

從海洋生物樽海鞘的航行和捕食行為中受到啟發，Mirjalili提出了SSA算法，并用于解決優化問題。具體來說，將整個樽海鞘種群分為領導者和跟隨者，領導者負責帶領整個種群，跟隨者彼此跟隨。通過這種方式，整個種群逐步向食物位置(即全局最優點)移動。領導者的位置通過以下公式更新：

(11)

(12)

2.2 高斯偽譜法

首先簡單介紹一下高斯偽譜法的基本原理。在一系列Gauss點上離散運動學方程中的狀態變量和控制變量，并以這些離散點為節點構造Lagrange插值多項式來近似狀態變量和控制變量。狀態變量對時間的導數通過對全局插值多項式求導來近似，從而將微分方程約束轉化為一組代數約束。終端狀態可以由初始狀態和高斯積分聯合表示，目標函數中的積分部分也可以由高斯積分表示。通過變化，可將最優控制問題轉化為具有一系列代數約束的NLP問題。

1)時域變換

采用GPM求解軌跡優化問題時，由于所涉及的正交多項式的正交區間是τ∈[-1,1]，因此需要將時間區間[t0,tf]轉換到[-1,1]，可對時間變量t做如下變換：

(13)

2)狀態和控制變量近似

GPM的配點τk(k=1,2,…,N)為N個LG(Legendre-Gauss)點，即N階Legendre多項式的根。配點分布區間為(-1,1)，加上邊界點τ0=-1，作為區間[-1,1)的N+1個節點。因此，狀態變量x(τ)可由N+1個Lagrange插值多項式的組合來近似：

(14)

類似地，控制變量u(τ)可由N個Lagrange插值多項式的組合來近似：

(15)

3)微分方程約束

對式(14)求微分，可得

(16)

式中：微分矩陣Dki可以離線確定

(17)

式中：PN(τ)為N階Legendre多項式。從而將動力學微分方程約束轉化為代數約束：

(k=1,2,…,N)

(18)

4)終端狀態約束

由于式(14)中狀態變量的近似忽略了終端時刻τf，因此需要補充終端狀態約束：

(19)

將終端約束條件離散并用Gauss積分來近似，可得：

U(τk),τk;t0,tf)

(20)

5)目標函數近似

將目標函數(7)中的積分項用Gauss積分近似，可得：

(21)

邊界條件約束(9)和不等式約束(10)離散為：

φ(X(τ0),t0,X(τf),tf)=0

(22)

C(X(τk),U(τk),τk;t0,tf)≤0

(23)

經過上述變換，原最優控制問題轉化為一個NLP問題。

3 基于DASSAGPM的再入軌跡優化

3.1 動態自適應樽海鞘群算法

本質上來說，智能優化算法都會面臨陷入局部最優的問題，因此改進策略都是圍繞著平衡算法尋優的探索和利用過程、避免陷入局部最優來展開。

從式(11)可以看出，系數c1(l)是平衡優化過程中探索和利用的主要參數。隨著迭代次數的增加，c1(l)不斷減小，可見優化前期重探索，更加注重全局搜索能力，后期重利用，更加注重局部搜索能力。雖然這一變化趨勢是合理的，但是這種變化機制是被動的，沒法根據當前種群優化狀況進行動態精確的調整。為了使種群更加合理地進化，本節在傳統SSA算法的基礎上，借鑒經典控制理論中的“反饋”思想，引入種群改善率作為反饋量來動態自適應更新系數c1(l)，提出了一種新穎的DASSA優化算法。

定義種群改善率R=Ms/M，其中Ms為當前種群中比上一代更接近全局最優解的個體數，M為種群大小。根據1/5原理，算法參數應該動態自適應變化以使種群改善率為20%[17]。所以，當R<0.2時，表明探索能力強而利用能力弱，此時應該減小c1(l)來提高搜索精度。當R>0.2時，種群主要在進行局部尋優，此時應該增大c1(l)來提高全局搜索能力。當R=0.2時，探索和利用達到了合理的平衡，此時不需要對c1(l)進行動態自適應調整。上述動態自適應規則可以表示為：

(24)

式中：c1(l-1)為上次迭代的調整系數，0<ζ<1為學習因子。

需要說明的是，式(24)中將不需要調整c1(l)的范圍由R=0.2擴大至0.2≤R≤0.3，這是由于R=0.2出現的概率很低，會導致c1(l)頻繁變化，不利于參數自調整的穩定。

3.2 基于罰函數法的約束條件的處理

SSA優化算法最初被提出是用來解決不帶約束的靜態優化問題，而再入軌跡優化問題包含多種約束，沒法直接利用SSA算法求解，而是需要對約束進行處理。一般來說，罰函數法是最為普遍的約束處理方法，具有原理簡單且容易實現的優點[18]。由此，在利用DASSA算法優化過程中，適應度函數可定義為目標函數加上罰函數，表示如下：

F(x)=J+P(x)

(25)

式中：J為目標函數，P(x)為罰函數，可表示為：

pb(x)+pf(x)]

(26)

式中：N為軌跡優化過程中的離散點個數，ph(x)，pd(x)，pl(x)，pb(x)和pf(x)分別為超出熱流密度、動壓、過載、邊界條件和終端狀態約束的罰函數，定義如下：

(27)

(28)

(29)

(30)

(31)

(32)

(33)

3.3 DASSAGPM算法的優化程序

針對高超聲速飛行器再入軌跡優化問題，利用DASSAGPM混合優化算法求解最優控制量。具體來說，將整個優化過程分為兩個階段。在第一階段，利用DASSA優化算法對控制量進行全局尋優，這樣能充分利用智能優化算法的全局搜索能力，避免陷入局部最優。當適應度函數變化值達到預先設置的停止標準時，停止DASSA算法的尋優過程，并將求得的近似全局最優解作為GPM優化的初始猜測值。第二階段，利用GPM對控制量進行全局尋優，這樣便于利用GPM收斂速度快、精度高的優點。DASSAGPM優化算法的詳細步驟如下：

階段1：

1)將帶復雜約束的高超聲速再入軌跡優化問題公式化。

2)設定DASSAGPM優化算法的參數：最大迭代次數L，種群大小M，離散LG點數N，DASSA算法停止標準eDASSA，GPM算法停止標準eGPM。

3)將時間區間[t0,tf]用N個LG點等距離散。將控制變量u=[α,σ]T作為待優化的變量，每個種群個體的位置即代表一組控制變量的優化值。種群的完整信息可用矩陣表示為：

X=[X1,…,Xn,…,XN]

式中:離散點n的優化矩陣為：

式中:αnm表示第n個離散點上第m個種群個體尋優的攻角值，σnm表示第n個離散點上第m個種群個體尋優的傾側角值。

4)在約束范圍內隨機初始化M個種群個體的位置。

5)基于種群個體的離散位置信息通過插值法計算連續控制量u，結合控制量對動力學方程(1)積分求得系統狀態和狀態的微分，接著計算每個種群個體的適應度函數值。

6)更新所有種群個體的位置，并更新全局最優值和局部最優值。

7)如果達到最大迭代次數或者DASSA算法停止標準eDASSA，則停止DASSA算法優化過程，并保存目前獲得的最優控制量以及相關狀態量，然后，進入8)；否則返回5)。

階段2：

8)將離散點處DASSA算法優化得到的控制量值作為GPM的初始猜測值。GPM優化過程中的相關參數，如變量范圍、邊界范圍和離散LG點都保持不變。

9)通過GPM將再入軌跡優化問題轉化為NLP問題后，利用NLP求解器求解。如果達到GPM算法停止標準eGPM，則停止GPM算法優化過程。

10)利用得到的最優控制量u*對動力學方程(1)積分求得系統狀態量，由此可得再入過程中的最優控制變量和狀態變量。

11)如果計算時間和代價函數值滿足條件，結束仿真過程；否則，調整離散LG點數N并返回3)。

DASSAGPM算法的優化流程圖如圖1所示。需要指出的是，軌跡優化過程是離線完成的，DASSAGPM算法比GPM算法的步驟多，實際上并不會影響再入過程的實時性。另外，本文算法雖然是以高超聲速飛行器為背景，針對在存在復雜約束條件下傳統GPM算法在解決高超聲速飛行器再入軌跡優化問題時對初始猜測值敏感的不足而提出的。但是，該算法本質上是一種優化方法，因此不僅僅適用于高超聲速飛行器的再入軌跡優化問題，同時也適用于任何飛行軌跡優化甚至一般的優化問題。

圖1 DASSAGPM算法優化流程圖Fig.1 Optimization flow chart of DASSAGPM algorithm

4 仿真校驗

本節針對DASSA優化算法和DASSAGPM軌跡優化算法進行了大量的仿真，以驗證所提算法的可行性和優越性。

4.1 DASSA的優化性能

為了驗證DASSA算法的優化特性，利用DASSA算法對文獻[13]中的部分基準函數進行迭代尋優。作為比較，同時采用WOA、SSA、PSO和GWO對基準函數尋優。為保證對比的公平性，統一設置參數為L=500，M=30。表1給出了基準函數描述，其中，基準函數表達式為：

表1 基準函數Table 1 Benchmark functions

圖2～圖5給出了五種優化算法對基準函數優化過程中的函數值收斂曲線?？梢钥闯?，SSA和DASSA算法在收斂速度上具有明顯優勢，同時也更容易達到全局最優。進一步觀察可知，由于在SSA算法中引入了動態自適應機制，使得改進后的DASSA算法在尋優過程中探索和利用的平衡更為合理，收斂速度和精度進一步得到改善。

圖2 f1優化過程Fig.2 Optimization process of f1 function

圖3 f2優化過程Fig.3 Optimization process of f2 function

圖4 f3優化過程Fig.4 Optimization process of f3 function

圖5 f4優化過程Fig.5 Optimization process of f4 function

4.2 DASSAGPM的軌跡優化性能

本節將通過大量仿真驗證DASSAGPM算法在飛行器軌跡優化上的優越性。在利用GPM進行軌跡優化時，本文采用開源軟件GPOPS求解最優軌跡，其中SNOPT軟件作為NLP求解器。為了展示DASSAGPM算法的優越性，將GPM、SSAGPM和DASSAGPM三種算法的軌跡優化性能進行對比。出于對比的公平性考慮，參數統一設置為L=500，M=40，N=31，eSSA=eDASSA=100，eGPM=1×10-5。熱流密度約束為1.5 MW/m2，動壓約束為200 kPa，過載約束為4.5g。狀態變量的初始值和終端約束如表2所示。

表2 狀態變量的初始值和終端約束Table 2 Initial values and terminal constraints of state variables

圖6～圖9給出了飛行器再入過程中的三維軌跡、速度、航跡傾角和航跡偏角曲線，可以看出，三種算法均能使高度、經緯度、速度、航跡傾角和航跡偏角滿足終端約束要求。進一步觀察可知，GPM算法優化得到的軌跡跳躍最劇烈，SSAGPM次之，DASSAGPM算法的軌跡曲線最為平滑平緩，說明DASSAGPM算法的優化解最接近全局最優解，代價函數值最小。

圖6 三維軌跡曲線Fig.6 Three dimensional trajectory

圖7 速度曲線Fig.7 Velocity profile

圖8 航跡傾角曲線Fig.8 Flight-path angle profile

圖9 航跡偏角曲線Fig.9 Heading angle profile

圖10和圖11分別給出了控制量攻角和傾側角的曲線對比?？梢钥闯觯N優化算法的控制量都在控制邊界內，DASSAGPM算法優化得到的控制量最光滑，而GPM算法的控制量最不光滑。

圖10 攻角曲線Fig.10 Angle of attack profile

圖11 傾側角曲線Fig.11 Bank angle profile

圖12～圖14分別給出了三種算法的熱流密度、動壓和過載曲線對比，可以看出均滿足約束條件，其中DASSAGPM算法的振幅最小，再入過程中的熱載最小，優化效果最好。

圖12 熱流密度曲線Fig.12 Heat flux density profile

圖13 動壓曲線Fig.13 Dynamic pressure profile

圖14 過載曲線Fig.14 Overload profile

為了進一步驗證DASSAGPM算法在軌跡優化方面的優越性，對三種算法進行50次蒙特卡羅仿真試驗，所得到的代價函數值如圖15所示。可以明顯看出，DASSAGPM算法的全局尋優能力最強，整體代價函數值最小，次好的是SSAGPM算法，相比之下，GPM算法的全局尋優能力最差。表3給出了代價函數值統計數據對比，DASSAGPM算法的代價函數值均值最小，全局尋優能力最好。

圖15 蒙特卡羅仿真代價函數值曲線Fig.15 Fitness value curves of Monte Carlo simulation

表3 代價函數值統計Table 3 Fitness value statistics

5 結 論

針對高超聲速飛行器面臨復雜的再入環境以及大的不確定性下的軌跡優化問題，提出了一種基于DASSAGPM算法的解決方案。

1)為了使SSA優化算法探索和利用之間的平衡更加地合理，提出了一種新穎的DASSA優化算法。仿真結果表明DASSA算法在基準函數優化上具有收斂速度快、更易獲取全局最優值的特點。

2)在利用DASSA算法進行軌跡優化時，采用罰函數法處理再入過程中的各種約束。

3)針對傳統GPM在軌跡優化過程中對初始猜測值敏感的不足，借助DASSA算法強大的全局搜索能力，提出了DASSAGPM再入軌跡優化方法。仿真結果表明所提出的DASSAGPM算法在求解再入軌跡優化問題時收斂速度快、精度高且更容易搜尋到全局最優解。