基于自適應偽譜法的高超聲速飛行器再入軌跡優化

任鵬飛， 王洪波， 周國峰

(中國運載火箭技術研究院, 北京 100076)

高超聲速飛行器飛行速度高，機動范圍大，具備入軌和再入大氣的能力，可實現遠程快速物資運送或精確打擊，具有較高的經濟軍事價值，近年來已成為各國的研究熱點[1]，其中具有代表性的為美國通用空天飛行器(Common Aero Vehicle, CAV)。再入飛行段作為此類飛行器重要的飛行階段，由于飛行器的力熱載荷環境嚴酷，且對于控制量高度敏感，因此再入軌跡優化設計也面臨著一定挑戰。

此類問題可以描述為具有多種復雜約束的非線性連續時間最優控制問題。在近20年里，直接配點法已經成為數值求解非線性最優控制問題的有力工具之一。直接配點法的基本思想是在所研究時域內合理的選取一系列不同的點，對狀態變量和控制變量進行離散并建立一定的約束關系，將連續時間最優控制問題轉化為有限維的非線性規劃 (Nonlinear Programming,NLP)問題，從而通過著名的NLP求解器，例如SNOPT[2]、IPOPT[3]、 KNITRO[4]等進行有效求解。相比于間接法，直接配點法無需推導Pontryagin極小值原理的最優控制一階必要條件[5]，降低了初值敏感性，減小了求解的難度，在飛行器軌跡優化等領域得到了廣泛應用。

起初，直接配點法中發展了h方法，例如Runge-Kutta方法等，其通過將時間區間劃分為多個子區間，在每個子區間上利用固定階數的多項式來近似狀態變量和控制，通過增加子區間數量來實現h方法的收斂[6]。

近年來，大量的研究都聚焦于采用一系列Gauss正交配點的方法[7-10]，稱為偽譜法或正交配點法。由于高精度Gauss積分規則，LG(Legendre- Gauss)、LGR(Legendre-Gauss-Radau)以及LGL(Legendre-Gauss-Lobatto)這3種配點類型被廣泛采用，它們分別是Legendre多項式的根、Legendre多項式線性組合的根以及Legendre多項式導數的根。偽譜法屬于p方法，通常基于Gauss正交規則，在整個時間區間選取一系列配點，通常采用Lagrange全局插值多項式近似狀態變量和控制變量，通過配置微分代數方程來近似微分約束，將連續時間最優控制問題轉換為NLP問題。通過增加多項式的階數來實現p方法的收斂[11]。

自適應偽譜法[12-14]結合了h方法和p方法的特點，通過估計當前解的相對誤差，按照一定規則調整區間數量和多項式階數，相比單獨使用h方法或p方法，提高了求解的精度及效率，并能夠有效應對狀態控制等不連續的問題。

本文針對高超聲速飛行器再入軌跡優化問題，建立了三自由度再入運動方程，以CAV-H飛行器為對象，考慮主要約束條件，基于Legendre多項式近似理論，建立了LGR偽譜法誤差估計關系式，采用自適應網格重構技術，實現了3種典型性能指標再入軌跡優化問題的高效求解。本文算法相對傳統自適應偽譜法具有配點和區間分配合理，收斂速度快及對人工參數不敏感等優點。

1 再入軌跡優化問題

1.1 三自由度再入運動方程

本文考慮地球為旋轉圓球，飛行器無側滑情況下，建立無動力高超聲速飛行器再入運動方程：

(1)

式中：h、V、θ、φ、γ、ψ、m、ωe、σ和r分別為高度、速度、經度、緯度、航跡角、航向角、飛行器質量、地球自轉角速度、傾側角和地心距；引力常數μ=3.986×1014m3/s2；地球半徑Re=6 378 137 m；L和D分別為升力和阻力，其表達式為

(2)

其中：ρ為大氣密度；海平面大氣密度ρ0=1.2 2 56kg/m3；歸一化高度H0=7 110 m；升力系數CL和阻力系數CD由飛行器外形和飛行狀態決定；S為參考面積。本文以CAV-H為例，相關總體和氣動參數來自文獻[15]，將氣動系數擬合為馬赫數Ma和迎角α的函數。

(3)

1.2 約束條件

1) 端點約束

為滿足再入段末端能量管理要求，對末端高度hf、速度Vf及航跡角γf提出限制：

(4)

2) 過程約束

(5)

式中：KQ為與飛行器外形和材料相關的常數，本文取KQ=2.582 0×10-7;g0為海平面重力加速度。

3) 控制約束

為減輕姿態控制系統壓力，要求迎角和傾側角變化較為平緩，同時對其大小提出限制：

(6)

1.3 目標函數

考慮3種典型再入段軌跡優化目標，即最大縱程，最大橫程及最小航跡角變化。通過引入權重系數表征設計偏好，則一般性目標函數可以表示為

(7)

式中：w1、w2、w3為權重系數；t0為初始時刻；tf為末端時刻。

2 Bolza型連續時間最優控制問題

考慮如下Bolza型連續時間最優控制問題：

minJ=Φ(x(t0),t0,x(tf),tf)+

(8)

式中：x(t)∈RNx為狀態變量，Nx為狀態變量維數；u(t)∈RNu為控制變量，Nu為控制變量維數。Mayer型性能指標Φ，Lagrange型性能指標L，動力學微分約束f，路徑約束C及端點約束φ定義為

(9)

其中：NC和Nφ分別為路徑約束和端點約束維數。

3 多區間LGR偽譜法

多區間偽譜法的思想是將整個時間區間劃分為多個小的子區間，在每個子區間內采用較少數量的配點來離散最優控制問題。全局偽譜法則可認為是只有一個區間的特例。不失一般性，將t∈[t0,tf]分為s(s≥1)個子區間Sk=[tk-1,tk]，其中:k=1,2,…,s;t0

(10)

LGR偽譜法的配點采用τ∈[-1,1]上的N個LGR點或反轉的LGR點，本文選用前者。LGR點為(τ+1)(PN(τ)+PN-1(τ))的根，其中:PN(τ)為N階Legendre多項式。因此對t進行區間轉化，可使每個子區間滿足Legendre正交多項式的定義區間。

LGR配點和Gauss-Radau積分權重的計算方法主要為特征值方法和牛頓迭代法[16]。本文采用特征值方法進行求解，N-1階Jacobi矩陣為

(11)

式中：

(12)

AN-1的第i個特征值為λi，對應單位特征向量的首個元素為v1i。根據Radau點與Gauss點的對應關系[16]，可得到LGR配點τi及積分權值ωi為

(13)

式中：i=1,2,…,Nk-1。

在每個時間子區間內采用Lagrange基函數近似狀態變量和控制變量。狀態變量的插值節點為Nk個LGR點和τNk=1點，即插值節點數比配點數多1個；控制變量的插值節點為Nk個LGR點，末端控制點通過外插得到。因此第k個子區間的狀態變量和控制變量可近似表示為

(14)

式中：

(15)

狀態變量在配點處對τ的微分可表示為

j= 0,1,…,Nk-1

(16)

式中：DNk×(Nk+1)為Radau微分矩陣，表示各Lagrange基函數在各LGR配點處的微分值，其表達式為[16]

(17)

Z(k)(τ)=(1-τ)(PNk(τ)+PNk-1(τ))

(18)

至此最優控制問題(8)可以離散為微分形式的NLP問題：

minJ=Φ(X(1)(-1),t0,X(s)(1),tf)+

(19)

4 自適應網格重構技術

自適應偽譜法的基本思想：基于當前解的最大相對誤差來評估解的質量，當不滿足求解精度要求時，按照一定規則調整子區間內多項式階數或調整子區間數量，以達到改善解精度的目的。

4.1 誤差評估

文獻[13]給出了一種近似相對誤差的算法，其基本思想是通過對比2種狀態近似解，從而構建相對誤差的近似表達式。

j=1,2,…,Nk+1

(20)

(21)

此時可以構建第i維狀態變量絕對誤差和相對誤差的近似表達式為

(22)

則Sk內最大相對誤差為

(23)

4.2 Legendre多項式近似

根據收斂性理論[17]，對于光滑問題，全局偽譜法以指數形式收斂。文獻[17]指出，對于分段光滑函數，采用Legendre多項式的近似誤差與多項式系數具有相同的衰減率，且衰減率與多項式的階數相關。因此可以基于Legendre多項式系數構建多項式階數與近似誤差的關系。

對于分段光滑有界函數y(τ)可以展開為Legendre級數：

(24)

式中：pi為Legendre系數；Pi(τ)為Legendre多項式。取N階近似為

(25)

則截斷誤差為

(26)

根據Legendre多項式的正交性：

(27)

式中：δmn為Kronecker函數。則式(26)化簡為

(28)

根據文獻[17]分段光滑有界函數y(τ)的Legendre系數近似為

(29)

(30)

因此可根據離散Legendre近似形式：

(31)

采用最小二乘法求解c和υ。

同時考慮到級數關系式：

(32)

因此可構建近似誤差關系式：

(33)

4.3 光滑性判定

4.4 多項式階數更新

當υ>υε時，根據式(34)可計算當前網格的求解誤差：

(34)

(35)

(36)

為了嚴格增加配點數量，采用向上取整算子:

(37)

4.5 區間更新

當υ≤υε時，考慮到υ接近于0時，由式(33)確定區間分段數會產生很大的區間數量，從而導致NLP問題的規模過大，使得求解效率降低。因此采用υε代替υ，根據式(36)估計所需的多項式階數：

(38)

令每個新子區間與當前子區間的配點數相同，則區間分段數Hk為

(39)

4.6 網格縮減

減少區間配點數的基本思想是將Lagrange插值多項式展開為τ的冪級數形式，并計算各階次系數。為了減少該多項式階數同時保證相對誤差仍滿足要求，考慮到τ∈[-1,1]，令τ=τmax=1，因此僅需從最高階次系數向最低階次系數依次累加，直至累加值超過容許誤差ε，則可減少的配點數不超過累加次數。

合并相鄰子區間的基本思想與上述算法類似，只是將相鄰子區間的值在當前區間對應多項式上進行外插，然后類似地判斷容許誤差是否滿足。

網格縮減算法的具體推導過程見文獻[14]，本文不再贅述。

4.7 算法流程

本文自適應偽譜法的計算流程如下：

步驟1給定求解容許誤差ε和初始網格。

步驟2求解NLP問題(19)。

步驟5重復步驟2。

綜上，本文算法的計算流程圖見圖1。

圖1 自適應偽譜法流程圖Fig.1 Flowchart of adaptive pseudospectral method

5 算例分析

再入軌跡優化算例設置見表1，邊界條件設置見表2。計算環境為Windows XP 2.83 GHz，3.25 GB內存，軟件環境為MATLAB R2014a，NLP求解器采用SNOPT v7.0，梯度計算采用有限差分算法。

為了驗證本文自適應偽譜法求解的有效性與快速性，采用本文算法(記為hpL)求解最優控制變量，并采用變步長Runge-Kutta-Fehlberg法(記為RKF45)進行積分驗證，積分相對誤差為10-13。對比算法分別采用文獻[12]算法(記為hp)以及文獻[13]算法(記為ph)。本文算法軌跡優化結果見圖2～圖6，不同算法對比結果見表3。

表1 算例設置Table 1 Setup of calculation examples

表2 邊界條件設置Table 2 Setup of boundary conditions

由圖2～圖4可以發現，對于3種典型算例，采用hpL算法的求解結果與采用變步長RKF45算法積分的結果均保持一致。優化獲得的飛行器再入軌跡末端均嚴格滿足末端高度約束，速度約束以及航跡角約束。由圖4可以發現，h-V平面內,3種算例的最優軌跡均位于再入走廊下邊界的上方，即全程滿足飛行過程中的熱流率約束、法向過載約束以及動壓約束要求。可以注意到，對于本文CAV-H飛行器，法向過載在再入段過程中均不是主導約束，再入初始階段，飛行器飛行高度高，大氣非常稀薄，即使飛行速度較高，相比于動壓，熱流率為主導約束，由圖4可以發現，此時飛行器采用較大迎角以獲得較大升力，從而有效減緩飛行高度的快速下降，以滿足熱流率要求。當飛行器能量下降至一定時，動壓為主導約束。對于最小航跡角變化率問題，整個飛行軌跡平滑無跳躍，初始以較大迎角飛行，飛行速度快速下降以匹配平穩飛行的能量要求。由圖5和圖6可以發現，整個飛行過程中，3種算例的迎角和傾側角均變化緩慢，其變化率均不超過1(°)/s。從而驗證了本文算法的有效性。

圖2 3D軌跡Fig.2 Three-dimensional trajectory

圖3 航跡角變化Fig.3 Path angle versus time

圖4 再入走廊邊界Fig.4 Boundary of reentry corridor

圖5 迎角與傾側角變化Fig.5 Angle of attack and heeling angle versus time

圖6 迎角與傾側角變化率變化Fig.6 Change rate of angle of attack and heeling angle versus time

表3 不同算法軌跡優化結果Table 3 Trajectory optimization results by different methods

由表3可以發現，本文hpL算法具有較好的穩健性，人工參數衰減容限υε一般取值為(0,1]之間，本文設置為0.25和0.5時，3種算例的求解時間和網格迭代次數差別不大，原因是采用Legendre衰減系數算法進行光滑性判定時，衰減系數υ的區分度較大，使得人工參數在一定區間內取值時，對增加配點數和分割區間數的選擇影響較小，因此軌跡優化效率對于人工參數的選擇不敏感。對于hpL算法，增加的配點數取決于衰減系數υ，即與狀態平滑性相關，而ph算法可認為是υ≡lgN的算法，當N<10時，ph算法估計的多項式階數偏大，其求解效率受到區間最大配點數的影響較大，對于最大射程問題，最大配點數為6的求解時間僅為最大配點數為10的0.1倍，此時ph算法更傾向于分段，從而避免稠密的約束Jacobian矩陣。類似地，對于hp算法可認為是υ≡1的算法。本文phL算法對于3種算例均展現了較好的求解快速性，并且相比于hp算法和ph算法，hpL算法的網格迭代次數更少，收斂速度更快，同時采用網格縮減技術，使得配點和區間總數均相對更少。驗證了基于Legendre衰減系數的誤差估計的有效性以及本文算法的快速性。

6 結 論

本文針對高超聲速飛行器再入段軌跡優化問題，以CAV-H飛行器對象，基于多區間Radau偽譜法，對3種典型性能指標的連續時間最優控制問題進行離散，基于Legendre多項式近似理論構建相對誤差估計關系式，采用自適應網格重構技術實現了再入軌跡優化的高效求解。仿真結果表明：

1) 本文算法結果與變步長Runge-Kutta-Fehlberg法積分結果一致。

2) 本文算法較2種傳統自適應偽譜法，配點與區間分配更合理，計算效率更高，且對于人工參數不敏感，具有較好的工程應用價值。