高超聲速伸縮式變形飛行器再入軌跡快速優化

岳彩紅, 唐勝景, 郭 杰,*, 王 肖, 張浩強

(1.北京理工大學宇航學院, 北京 100081; 2.中國運載火箭技術研究院戰術武器事業部, 北京 100076;3.中國北方工業有限公司, 北京 100053)

0 引 言

近年來,復雜的飛行環境和飛行任務對具有任務靈活性和多功能性的飛行器的需求,以及新型材料、新型結構等技術的進步,使得能夠根據飛行環境、任務需求等自適應改變外形,且飛行性能始終最優的變形飛行器重新受到了各軍事強國的重視[1-6]。高超聲速滑翔飛行器速度快、航程遠、強突防等優勢,使其成為當前高超聲速領域的研究熱點之一[7-11]。然而現有的固定外形的氣動布局,使得跨大空域、全速域飛行的滑翔飛行器無法達到整個飛行階段的性能最優。高超聲速飛行器中引入變形技術,可以有效改善現有滑翔飛行器的性能局限,使飛行器的潛力得到更好的發揮,具有以下優勢:

(1)通過變形提高飛行器的環境適應能力,實現整個飛行過程的性能最優;

(2)通過變形改變飛行器的氣動特性,進而結合攻角和傾側角控制實現協調控制,增加飛行器的軌跡控制能力;

(3)通過變形使彈道形式多變,增強飛行器的突防能力。

在變形飛行器軌跡優化方面,文獻[12]研究了飛行器變形與任務設計相耦合的優化問題。文獻[13]對襟翼可變形的空中客機的巡航性能進行了優化研究。然而,上述研究均是針對低速飛行器,對于高超聲速變形飛行器軌跡優化問題,文獻[14]采用基于分解的多目標優化算法優化了高超聲速變展長飛行器的再入軌跡,結果表明變形飛行器的射程增加且總吸熱量降低。文獻[15]采用偽譜法研究了連續變形飛行器的再入飛行能力。文獻[16]基于改進的高斯偽譜法優化求解了具有伸縮小翼的臨近空間飛行器的爬升段軌跡,仿真結果表明伸縮小翼的引入有利于節省燃料。

高超聲速變形飛行器再入軌跡優化與傳統高超聲速飛行器軌跡優化一樣,均為最優控制問題,主要有直接法和間接法兩種求解方法[17]。間接法精度較高,但推導過程復雜,對初值敏感,收斂性差。直接法通過直接離散控制量和/或狀態量將最優控制問題轉化為非線性規劃(nonlinear programming, NLP)問題。偽譜法作為一種離散控制變量和狀態變量的直接配點法,具有初值不敏感、易收斂、精度高等優點,受到了學者們的廣泛研究[18-21]。Gauss偽譜法區別于其他偽譜法的顯著優點是其離散的NLP問題的KKT(Karush-Kuhn-Tucker)條件與原最優控制問題的一階必要條件的離散形式等價,能夠保證解的精度[22]。

NLP問題需要借助參數優化方法求解,常用的參數優化方法包括序列二次規劃(sequential quadratic programming, SQP)和內點法。使用SQP算法或內點法求解NLP時,需要提供目標函數和約束對設計變量的偏導數,由于NLP問題的約束和設計變量數目很多,使用有限差分法計算時,偏導數求解非常耗時,優化效率很低。文獻[23-26]分別研究了局部配點法、Radau偽譜法偏導數矩陣的稀疏性,并將NLP梯度中非零項的求解轉化為原最優控制問題偏導數的求解,仿真結果表明可有效提高NLP的求解效率,然而目前對于Gauss偽譜法,還沒有文獻基于偏導數矩陣的稀疏性給出其梯度的高效計算方法。

本文針對一種高超聲速伸縮式連續變形飛行器,建立了變形飛行器再入軌跡優化模型,采用Gauss偽譜法將最優控制問題轉化為NLP問題,基于NLP偏導數的稀疏性,推導了目標函數梯度和約束Jacobian矩陣的高效計算方法,優化求解了變形飛行器的最大橫向航程、再入可達區、最大終端速度和最小飛行時間,仿真結果表明推導的梯度求解方法可有效提高優化效率,變形飛行器相對于固定外形飛行器的性能更加優越。

1 高超聲速伸縮式變形飛行器

本文研究對象為翼身組合伸縮式變形飛行器,如圖1所示。梯形升力彈翼兩端通過套筒結構連接了兩級伸縮彈翼[27],兩側的伸縮翼為對稱連續變化。相對于傳統的翼身組合式飛行器,該變形飛行器可通過機翼的伸縮改變機翼面積、機翼展長和展弦比,進而改變作用在飛行器上的氣動力。

圖1 變形飛行器構型

圖1中,W1和W2分別為一級伸縮翼和二級伸縮翼,面1為梯形翼與一級伸縮翼的連接面,面2為二級伸縮翼與一級伸縮翼的連接面。W1的翼梢位于面1時,一級伸縮翼為完全收縮狀態,W1的翼根位于面1時,一級伸縮翼為完全伸長狀態。類似地,W2的翼梢位于面2時,二級伸縮翼為完全收縮狀態,W2的翼根位于面2時,二級伸縮翼為完全伸長狀態。定義W1和W2均完全收縮的狀態為變形飛行器的原外形,W1完全伸長、W2完全收縮的狀態為變形飛行器的變形1,W1和W2均完全伸長的狀態為變形飛行器的變形2,其余狀態為中間狀態。設梯形翼的展長為l0,W1和W2的伸長量分別為l1和l2,變形量[15]定義為η=(l0+l1+l2)/l0,原外形、變形1、變形2這3種狀態下的變形量分別為η=1,η=2,η=3。

2 高超聲速變形飛行器再入軌跡優化問題

2.1 變形飛行器無動力再入運動模型

2.1.1 氣動力模型

變形飛行器氣動系數是攻角、馬赫數和變形量的函數。首先利用氣動工程計算軟件DATCOM計算原外形、變形1、變形2這3種外形在不同攻角和馬赫數下的升力系數和阻力系數,然后將得到的3種外形的氣動系數擬合成如下所示的函數[28]:

(1)

式中:CL為升力系數;CD為阻力系數;p0～p5和q0～q5分別為升力系數和阻力系數表達式的常值參數;α為攻角;Ma為馬赫數。

表1中給出了3種外形下的擬合結果。方程的確定系數RS均接近于1,表明了擬合結果的高精度。變形飛行器的氣動數據可通過不同外形下的氣動系數對變形量進行插值得到,例如原外形、變形1、變形2這3種外形下的p0值,對變形量進行插值可得到變形飛行器的p0值。

表1 氣動系數擬合結果

2.1.2 無量綱動力學方程

假設地球為均勻圓球,飛行器為無動力飛行的質點,不考慮地球自轉,高超聲速變形飛行器的三自由度無量綱再入運動方程[29]為

(2)

2.2 再入軌跡約束模型

2.2.1 初始條件

高超聲速變形飛行器再入的初始條件為

(3)

式中:(·)0為初始時刻的狀態量;t0為初始時刻。

2.2.2 路徑約束

為保證再入過程中飛行器的安全性,需滿足如下路徑約束:

(4)

(5)

(6)

2.2.3 終端約束

終端約束包括終端位置約束和終端速度約束,即

(7)

式中:(·)f為終端狀態量;(·)c為期望的終端狀態量。根據實際的任務需求,終端約束一般取部分約束。

2.2.4 控制量約束

變形飛行器的控制量包括攻角、傾側角和展長變形量,滿足如下約束:

(8)

2.3 目標函數

為了研究變形飛行器的再入飛行性能,本文的目標函數主要涉及最大橫向航程、最大縱向航程、最大終端速度和最小飛行時間。本文初始條件中,λ0=0,φ0=0,所以目標函數為

(9)

2.4 最優控制問題

高超聲速變形飛行器再入軌跡優化問題可表示為連續時間最優控制問題,尋找控制變量u=[α,β,η]T,最小化Bolza型目標函數:

(10)

滿足微分方程約束:

(11)

式中:f∈R6為6維向量函數,表達式如式(2)所示。

邊界約束:

Bmin≤B(x(t0),x(tf),t0,tf)≤Bmax

(12)

式中:B為邊界約束;Bmin,Bmax分別為邊界約束的下界和上界,表達式如式(3)和式(7)所示。

路徑約束:

Cmin≤C(x(t),u(t),t)≤Cmax

(13)

式中:C為路徑約束;Cmin,Cmax分別為路徑約束的下界和上界,表達式如式(4)～式(6)所示。

3 Gauss偽譜法及NLP梯度計算

3.1 Gauss偽譜法

Gauss偽譜法將狀態變量和控制變量在LG(Legendre-Gauss)點上進行離散,采用基于LG點的全局插值多項式對節點上的狀態進行近似。對LG點處的狀態求導,將微分方程約束轉化為一系列代數約束,終端狀態通過初始值和右函數的積分求得,目標函數中的積分項可由高斯積分近似,將最優控制問題轉換為NLP問題。

Gauss偽譜法離散步驟如下。

步驟 1時域轉換

LG點位于(-1,1)之間,故需要將原優化問題的時間區間轉換到(-1,1)之間:

(14)

步驟 2全局插值多項式近似狀態變量和控制變量

Gauss偽譜法的N個配點(N階Legendre多項式的根){τ1,τ2,…,τN},加上τ0=-1,構成N+1個節點,狀態變量由基于N+1個節點的Lagrange多項式近似表示：

(15)

式中:x(τ)為τ時刻實際狀態量;X(τ)為τ時刻近似狀態量;X(τi)為τi時刻的離散狀態量;Li(τ)為Lagrange插值基函數，可表示為

(16)

控制變量一般也采用Lagrange多項式近似表示:

(17)

步驟 3動力學微分方程約束轉換為代數約束

在配點處對式(15)求導可得

(18)

式中:k=1,2,…,N;D∈RN×(N+1)為微分差分矩陣,可以由下面公式離線求得:

(19)

令式(18)求得的導數等于狀態變量的右端函數,可將微分方程約束轉化為代數約束:

(20)

步驟 4終端約束

終端狀態需要滿足動力學約束,可由初始狀態加上右函數積分求得,積分部分由Gauss積分近似:

(21)

步驟 5目標函數轉化

將Bolza型目標函數中積分項用Gauss積分近似:

(22)

同時邊界約束和路徑約束離散為如下形式:

B(X0,t0,Xf,tf)=0

(23)

C(Xk,Uk,τk;t0,tf)≤0,k=1,2,…,N

(24)

經過以上轉換,Gauss偽譜法將高超聲速變形飛行器再入軌跡優化問題轉化為如下形式的NLP問題:

minf(x)

(25)

式中:x為NLP問題的設計變量,包括所有離散點處的狀態變量、控制變量和端點時間t0和tf(若時間不固定);gi(x)為NLP的不等式約束,包括式(24),hj(x)為NLP的等式約束,包括式(20)、式(21)和式(23)。

3.2 NLP梯度的高效計算方法

Gauss偽譜法將變形飛行器再入軌跡優化問題轉化為式(25)所示的NLP問題,NLP問題通過常用的參數優化方法如SQP和內點法求解時需要提供梯度信息。然而NLP問題的設計變量和約束數目繁多,例如,當N=30時,設計變量有284個,約束有276個,若用有限差分法計算NLP的梯度,計算將非常耗時,優化效率很低。考慮到目標函數和約束函數往往是部分或所有配點處設計變量函數的疊加,疊加的某一部分只與部分配點處的變量有關,并不是和所有設計變量均相關,NLP的梯度矩陣和約束的Jacobian矩陣是非常稀疏的,即矩陣中有很多零元素,計算梯度時利用其稀疏性省略零元素會提高計算效率,矩陣中非零項的求解可轉化為求解原最優控制問題的偏導數,考慮到原最優控制問題的自變量和約束數目均很少,故將NLP問題的偏導數轉化原最優控制問題的偏導數會進一步提高計算效率。基于第3.1節內容,采用Gauss偽譜法將高超聲速變形飛行器再入軌跡優化問題轉化為NLP問題,本節利用NLP梯度的稀疏性,并將梯度中的非零項轉化為求原最優控制問題的偏導數[25],推導NLP梯度的高效計算方法。

為了推導目標函數的梯度和約束的Jacobian矩陣,首先進行如下定義:

(26)

式中:X∈R(N+2)×6;U∈RN×3;ω∈RN×1;g∈RN×1;F∈RN×6;C∈RN×3;gi(i=1,2,…,N),fi(i=1,2,…,N),Ci(i=1,2,…,N)分別為配點i處的積分型目標函數、右函數值、路徑約束值。

NLP的設計變量為

y=[X:,1,…,X:,6,U:,1,…,U:,3,t0,tf]T

(27)

式中:X:,i(i=1,2,…,6)為X的第i列;U:,i(i=1,2,3)為U的第i列。

目標函數式(22)可寫成如下形式:

(28)

等式約束式(20)可寫成如下形式:

(29)

式中:X(1:N+1)表示X的前N+1行。

NLP的約束函數為

l=[Δ:,1,…,Δ:,6,B,C:,1,…,C:,3]T

(30)

3.2.1 目標函數的梯度

目標函數的梯度可分解為Mayer函數的梯度和積分型函數的梯度,Mayer函數的梯度為

(31)

式中:

(32)

記積分型目標函數為

(33)

則S的梯度為

(34)

式中:

(35)

目標函數的梯度為

(36)

3.2.2 約束的Jacobian矩陣

約束的梯度矩陣為

(37)

式中:

(38)

式中:Δ:, j(j=1,2,…,6)對X:,i(i=1,2,…,6)、U:,l(l=1,2,3)、t0和tf的偏導數的表達式為

(39)

Bm(m=1,2,…,6)對X:,i(i=1,2,…,6)、U:,l(l=1,2,3)、t0和tf的偏導數的表達式為

(40)

C:,z(z=1,2,3)對X:,i(i=1,…,6)、U:,l(l=1,2,3)、t0和tf的偏導數的表達式為

(41)

使用內點法求解NLP時,需要提供目標函數和約束對設計變量的偏導數,分別采用本文的式(31)～式(36)計算目標函數梯度,式(37)～式(41)計算約束Jacobian矩陣,可大大提高計算效率,減少優化時間。

4 仿真結果及分析

本節以第1節設計的高超聲速變形飛行器為研究對象求解其再入軌跡優化問題,氣動模型為第2.1節求得的氣動模型。NLP問題求解采用內點法工具包IPOPT。控制量邊界取值為αmin=0°,αmax=30°,σmin=-90°,σmax=90°,ηmin=1,ηmax=3過程約束為

(42)

邊界條件設置如表2所示。

4.1 梯度計算效率驗證

為了驗證本文推導的梯度求解方法的快速性,分別以最大緯度、最大速度和最小飛行時間為目標函數,對于每一個目標函數,配點數分別取N=10和N=20,采用本文推導的梯度計算方法和有限差分法進行梯度計算,兩種方法的其他仿真條件均相同,收斂到最優解的時間如表3所示。

表3 優化收斂時間

由表3可知,本文方法的優化收斂時間均小于有限差分法的收斂時間。經計算,當目標函數為最大緯度,N=10和N=20時,本文方法的優化收斂時間分別是有限差分法收斂時間的3.29%和1.15%;當目標函數為最大速度,N=10和N=20時,本文方法的優化收斂時間分別是有限差分法收斂時間的3.20%和1.33%;當目標函數為最小飛行時間,本文方法的優化收斂時間分別是有限差分法收斂時間的2.87%和1.25%。仿真算例表明,采用本文推導的梯度計算方法可大大節省優化時間,且配點數目越多,相對效率越高。

4.2 最大橫向航程軌跡和可達區

4.2.1 最大橫向航程

最大橫向航程可以表征飛行器的再入機動能力,為了研究變形飛行器相對固定外形飛行器機動性能的變化,采用Gauss偽譜法離散加IPOPT優化的求解步驟分別求解變形飛行器和固定外形飛行器的最大橫向航程,其中固定外形飛行器取為所設計的變形飛行器的原外形狀態,即η=1。梯度的計算采用本文第3.2節推導的方法,配點N=60。仿真結果如圖2～圖10所示。

圖2 最大橫向航程高度變化曲線

圖3 最大橫向航程速度變化曲線

圖4 最大橫向航程地面軌跡

圖5 熱流曲線

圖6 動壓曲線

圖7 過載曲線

圖8 最大橫向航程攻角曲線

圖9 傾斜角曲線

圖10 變形量曲線

圖2～圖10為變形飛行器和固定外形飛行器以最大橫向航程為目標函數的優化結果。圖2～圖7分別給出了高度、速度和地面軌跡等狀態變量,以及熱流、動壓和過載等過程約束的優化結果,離散最優解和沿優化控制變量積分動力學方程的數值積分軌跡高度重合,表明了優化結果的可靠性和高精度。由圖2、圖5、圖6和圖7可知,兩種飛行器初始段采用平滑軌跡以滿足熱流約束,后段通過跳躍軌跡使得橫向航程最大,再入過程嚴格滿足過程約束,再次表明優化方法的精度。圖8～圖10分別為攻角、傾側角和變形量等控制量的優化結果,變形飛行器的攻角和傾側角幅值整體小于固定外形飛行器,表明變形飛行器可通過增大變形量減小需用的攻角和傾側角,有利于緩解舵機壓力。圖10中,固定外形飛行器的變形量始終為1是因為本文取η=1為原外形,變形飛行器大部分時間處于最大變形狀態,大變形能夠增加飛行器的升阻比,進而提高飛行器的射程。飛行器再入末段,高度和速度逐漸接近相應的約束,變形量減小,一方面減小阻力使飛行器速度緩慢減小以增加橫向航程,另一方面減小升力使飛行器下降以滿足高度約束。

變形飛行器和固定外形飛行器的優化終端結果如表4所示。

表4 優化結果

由表4可知,變形飛行器相對于固定外形飛行器的最大緯度增加了52.23%,其機動能力顯著增強。變形飛行器的航程相對于固定外形飛行器增加了1 402.5 km,這是由于變形提高了飛行器的升阻比,使得航程增加,與上文理論分析結果一致。

4.2.2 可達區

可達區是指再入飛行器能夠到達的地球表面或者交班狀態的邊界曲線的范圍[30]。再入可達區的確定有利于任務規劃,可達區的外邊界能夠表征再入飛行器的機動能力。

本節以第1節設計的變形飛行器為對象,利用Gauss偽譜法優化求解其可達區,并與固定外形飛行器的可達區作對比。可達區的求解可轉化為求解一系列優化問題,求解流程如下。

步驟 1求解射向的最大縱程,目標函數為

J=min[-λ(tf)]

(43)

步驟 2求解射向的最小縱程,目標函數為

J=minλ(tf)

(44)

步驟 3取最大縱程和最小縱程之間的多個縱程,分別優化求解其最大橫向航程,目標函數為

J=maxφ(tf)

(45)

步驟 4將上述求得的末端軌跡相連可得一側的可達區,將可達區對稱可得另一側可達區(由于不考慮地球自轉,兩側可達區對稱)。

仿真結果如圖11～圖13所示。

圖11 變形飛行器左半側可達區邊界軌跡示例

圖12 固定外形飛行器左半側可達區邊界軌跡示例

圖13 兩種飛行器的可達區邊界

圖11和圖12分別給出了變形飛行器和固定外形飛行器的左側可達區的若干條邊界再入軌跡,由仿真結果可知,離散最優解與數值積分軌跡高度重合,表明優化結果的可靠性和精度。為了便于顯示,圖11和圖12中給出了再入軌跡的平面投影軌跡和可達區邊界。

按照上述可達區求解步驟,將優化得到的一系列落點相連可分別得到變形飛行器和固定外形飛行器的可達區,如圖13所示。兩種飛行器的可達區均呈橢圓形,固定外形飛行器的可達區位于變形飛行器的可達區內部,明顯小于變形飛行器的可達區。固定外形飛行器的縱向航程范圍λ=4.41°～63.97°,橫向航程范圍φ取值為-25.58°～25.58°,變形飛行器的縱向航程范圍λ取值為-0.56°～82.28°,橫向航程范圍φ取值為-38.94°～38.94°,經計算,相對于固定外形飛行器,變形飛行器的縱向航程覆蓋范圍提高了39.1%,橫向航程覆蓋范圍提高了52.23%,可達區增大,飛行器機動能力增強。

4.3 最大終端速度

滑翔段的終端速度越大越有利于末端打擊,為了研究變形飛行器相對固定外形飛行器最大終端速度的變化,在滿足終端經度、緯度約束的前提下求解最大終端速度,優化設置同第4.2節。優化結果如圖14～圖19所示。

圖14 最大終端速度高度曲線

圖15 最大終端速度速度曲線

圖16 最大終端速度地面軌跡

圖17 最大終端速度攻角曲線

圖18 最大終端速度傾側角曲線

圖19 最大終端速度變形量曲線

圖14～圖19給出了以最大終端速度為目標函數的優化結果。圖14～圖16分別為高度、速度和地面軌跡曲線,離散最優解與數值積分軌跡高度重合,表明優化結果的可靠性和精度。圖17～圖19分別為攻角、傾側角和變形量的優化結果。變形飛行器的終端速度為4 515.6 m/s,飛行時間為888.7 s,固定外形飛行器的終端速度為3 810 m/s,飛行時間為951.53 s,變形飛行器的終端速度提高了705.6 m/s,飛行時間減小了62.83 s,這說明變形飛行器可以較大的速度、較短的時間飛完相同的航程,具有更加快速的到達能力。由圖17可知,變形飛行器的攻角明顯小于固定外形飛行器,說明增加變形量這個控制量可以顯著改善攻角特性,緩解舵機控制的壓力,且有利于隔熱系統的設計。

4.4 最小飛行時間

飛行時間最小有利于飛行器的再入突防,本節以最小飛行時間為目標函數進行優化求解。考慮終端經度和終端緯度約束,優化設置同第4.2節,仿真結果如圖20～圖25所示。

圖20 最小飛行時間高度曲線

圖21 最小飛行時間速度曲線

圖22 最小飛行時間地面軌跡

圖23 最小飛行時間攻角曲線

圖24 最小飛行時間傾側角曲線

圖25 最小飛行時間變形量曲線

圖20～圖25給出了以最小飛行時間為目標函數的優化結果。圖20～圖22分別為高度、速度和地面軌跡曲線,離散最優解與數值積分軌跡高度重合,表明優化結果的可靠性和精度。圖23～圖25為攻角、傾側角和變形量的優化結果,變形飛行器的攻角始終小于固定外形飛行器,變形量幾乎處于最大變形狀態,表明變形量的引入緩解了舵機控制系統的壓力。變形飛行器的飛行時間為887 s,相對于固定外形飛行器的948 s減小了6.88%,說明變形飛行器具有更加快速的到達能力。

5 結 論

本文針對一種高超聲速伸縮式連續變形飛行器,采用Gauss偽譜法將高超聲速變形飛行器再入軌跡優化問題轉化為NLP問題,基于偏導數的稀疏性推導了NLP目標函數梯度和約束的Jacobian的高效計算方法,優化求解了變形飛行器的最大橫向航程、再入可達區、最大終端速度和最小飛行時間,仿真結果表明:

(1)基于稀疏性推導的梯度計算方法可以很大程度地節省優化時間,提高優化效率;

(2)高超聲速變形飛行器相比固定外形飛行器在飛行性能方面具有顯著優勢,最大橫向航程提高了52.23%,可達區的橫向航程覆蓋范圍提高了52.23%,縱向航程覆蓋范圍提高了39.1%,最大終端速度增大了705.6 m/s,最小飛行時間減小了6.88%。