基于自適應(yīng)偽譜法的高超聲速飛行器再入軌跡優(yōu)化

任鵬飛， 王洪波， 周?chē)?guó)峰

(中國(guó)運(yùn)載火箭技術(shù)研究院, 北京 100076)

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

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

起初，直接配點(diǎn)法中發(fā)展了h方法，例如Runge-Kutta方法等，其通過(guò)將時(shí)間區(qū)間劃分為多個(gè)子區(qū)間，在每個(gè)子區(qū)間上利用固定階數(shù)的多項(xiàng)式來(lái)近似狀態(tài)變量和控制，通過(guò)增加子區(qū)間數(shù)量來(lái)實(shí)現(xiàn)h方法的收斂[6]。

近年來(lái)，大量的研究都聚焦于采用一系列Gauss正交配點(diǎn)的方法[7-10]，稱(chēng)為偽譜法或正交配點(diǎn)法。由于高精度Gauss積分規(guī)則，LG(Legendre- Gauss)、LGR(Legendre-Gauss-Radau)以及LGL(Legendre-Gauss-Lobatto)這3種配點(diǎn)類(lèi)型被廣泛采用，它們分別是Legendre多項(xiàng)式的根、Legendre多項(xiàng)式線性組合的根以及Legendre多項(xiàng)式導(dǎo)數(shù)的根。偽譜法屬于p方法，通?；贕auss正交規(guī)則，在整個(gè)時(shí)間區(qū)間選取一系列配點(diǎn)，通常采用Lagrange全局插值多項(xiàng)式近似狀態(tài)變量和控制變量，通過(guò)配置微分代數(shù)方程來(lái)近似微分約束，將連續(xù)時(shí)間最優(yōu)控制問(wèn)題轉(zhuǎn)換為NLP問(wèn)題。通過(guò)增加多項(xiàng)式的階數(shù)來(lái)實(shí)現(xiàn)p方法的收斂[11]。

自適應(yīng)偽譜法[12-14]結(jié)合了h方法和p方法的特點(diǎn)，通過(guò)估計(jì)當(dāng)前解的相對(duì)誤差，按照一定規(guī)則調(diào)整區(qū)間數(shù)量和多項(xiàng)式階數(shù)，相比單獨(dú)使用h方法或p方法，提高了求解的精度及效率，并能夠有效應(yīng)對(duì)狀態(tài)控制等不連續(xù)的問(wèn)題。

本文針對(duì)高超聲速飛行器再入軌跡優(yōu)化問(wèn)題，建立了三自由度再入運(yùn)動(dòng)方程，以CAV-H飛行器為對(duì)象，考慮主要約束條件，基于Legendre多項(xiàng)式近似理論，建立了LGR偽譜法誤差估計(jì)關(guān)系式，采用自適應(yīng)網(wǎng)格重構(gòu)技術(shù)，實(shí)現(xiàn)了3種典型性能指標(biāo)再入軌跡優(yōu)化問(wèn)題的高效求解。本文算法相對(duì)傳統(tǒng)自適應(yīng)偽譜法具有配點(diǎn)和區(qū)間分配合理，收斂速度快及對(duì)人工參數(shù)不敏感等優(yōu)點(diǎn)。

1 再入軌跡優(yōu)化問(wèn)題

1.1 三自由度再入運(yùn)動(dòng)方程

本文考慮地球?yàn)樾D(zhuǎn)圓球，飛行器無(wú)側(cè)滑情況下，建立無(wú)動(dòng)力高超聲速飛行器再入運(yùn)動(dòng)方程：

(1)

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

(2)

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

(3)

1.2 約束條件

1) 端點(diǎn)約束

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

(4)

2) 過(guò)程約束

(5)

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

3) 控制約束

為減輕姿態(tài)控制系統(tǒng)壓力，要求迎角和傾側(cè)角變化較為平緩，同時(shí)對(duì)其大小提出限制：

(6)

1.3 目標(biāo)函數(shù)

考慮3種典型再入段軌跡優(yōu)化目標(biāo)，即最大縱程，最大橫程及最小航跡角變化。通過(guò)引入權(quán)重系數(shù)表征設(shè)計(jì)偏好，則一般性目標(biāo)函數(shù)可以表示為

(7)

式中：w1、w2、w3為權(quán)重系數(shù)；t0為初始時(shí)刻；tf為末端時(shí)刻。

2 Bolza型連續(xù)時(shí)間最優(yōu)控制問(wèn)題

考慮如下Bolza型連續(xù)時(shí)間最優(yōu)控制問(wèn)題：

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

(8)

式中：x(t)∈RNx為狀態(tài)變量，Nx為狀態(tài)變量維數(shù)；u(t)∈RNu為控制變量，Nu為控制變量維數(shù)。Mayer型性能指標(biāo)Φ，Lagrange型性能指標(biāo)L，動(dòng)力學(xué)微分約束f，路徑約束C及端點(diǎn)約束φ定義為

(9)

其中：NC和Nφ分別為路徑約束和端點(diǎn)約束維數(shù)。

3 多區(qū)間LGR偽譜法

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

(10)

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

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

(11)

式中：

(12)

AN-1的第i個(gè)特征值為λi，對(duì)應(yīng)單位特征向量的首個(gè)元素為v1i。根據(jù)Radau點(diǎn)與Gauss點(diǎn)的對(duì)應(yīng)關(guān)系[16]，可得到LGR配點(diǎn)τi及積分權(quán)值ωi為

(13)

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

在每個(gè)時(shí)間子區(qū)間內(nèi)采用Lagrange基函數(shù)近似狀態(tài)變量和控制變量。狀態(tài)變量的插值節(jié)點(diǎn)為Nk個(gè)LGR點(diǎn)和τNk=1點(diǎn)，即插值節(jié)點(diǎn)數(shù)比配點(diǎn)數(shù)多1個(gè)；控制變量的插值節(jié)點(diǎn)為Nk個(gè)LGR點(diǎn)，末端控制點(diǎn)通過(guò)外插得到。因此第k個(gè)子區(qū)間的狀態(tài)變量和控制變量可近似表示為

(14)

式中：

(15)

狀態(tài)變量在配點(diǎn)處對(duì)τ的微分可表示為

j= 0,1,…,Nk-1

(16)

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

(17)

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

(18)

至此最優(yōu)控制問(wèn)題(8)可以離散為微分形式的NLP問(wèn)題：

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

(19)

4 自適應(yīng)網(wǎng)格重構(gòu)技術(shù)

自適應(yīng)偽譜法的基本思想：基于當(dāng)前解的最大相對(duì)誤差來(lái)評(píng)估解的質(zhì)量，當(dāng)不滿足求解精度要求時(shí)，按照一定規(guī)則調(diào)整子區(qū)間內(nèi)多項(xiàng)式階數(shù)或調(diào)整子區(qū)間數(shù)量，以達(dá)到改善解精度的目的。

4.1 誤差評(píng)估

文獻(xiàn)[13]給出了一種近似相對(duì)誤差的算法，其基本思想是通過(guò)對(duì)比2種狀態(tài)近似解，從而構(gòu)建相對(duì)誤差的近似表達(dá)式。

j=1,2,…,Nk+1

(20)

(21)

此時(shí)可以構(gòu)建第i維狀態(tài)變量絕對(duì)誤差和相對(duì)誤差的近似表達(dá)式為

(22)

則Sk內(nèi)最大相對(duì)誤差為

(23)

4.2 Legendre多項(xiàng)式近似

根據(jù)收斂性理論[17]，對(duì)于光滑問(wèn)題，全局偽譜法以指數(shù)形式收斂。文獻(xiàn)[17]指出，對(duì)于分段光滑函數(shù)，采用Legendre多項(xiàng)式的近似誤差與多項(xiàng)式系數(shù)具有相同的衰減率，且衰減率與多項(xiàng)式的階數(shù)相關(guān)。因此可以基于Legendre多項(xiàng)式系數(shù)構(gòu)建多項(xiàng)式階數(shù)與近似誤差的關(guān)系。

對(duì)于分段光滑有界函數(shù)y(τ)可以展開(kāi)為L(zhǎng)egendre級(jí)數(shù)：

(24)

式中：pi為L(zhǎng)egendre系數(shù)；Pi(τ)為L(zhǎng)egendre多項(xiàng)式。取N階近似為

(25)

則截?cái)嗾`差為

(26)

根據(jù)Legendre多項(xiàng)式的正交性：

(27)

式中：δmn為Kronecker函數(shù)。則式(26)化簡(jiǎn)為

(28)

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

(29)

(30)

因此可根據(jù)離散Legendre近似形式：

(31)

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

同時(shí)考慮到級(jí)數(shù)關(guān)系式：

(32)

因此可構(gòu)建近似誤差關(guān)系式：

(33)

4.3 光滑性判定

4.4 多項(xiàng)式階數(shù)更新

當(dāng)υ>υε時(shí)，根據(jù)式(34)可計(jì)算當(dāng)前網(wǎng)格的求解誤差：

(34)

(35)

(36)

為了嚴(yán)格增加配點(diǎn)數(shù)量，采用向上取整算子:

(37)

4.5 區(qū)間更新

當(dāng)υ≤υε時(shí)，考慮到υ接近于0時(shí)，由式(33)確定區(qū)間分段數(shù)會(huì)產(chǎn)生很大的區(qū)間數(shù)量，從而導(dǎo)致NLP問(wèn)題的規(guī)模過(guò)大，使得求解效率降低。因此采用υε代替υ，根據(jù)式(36)估計(jì)所需的多項(xiàng)式階數(shù)：

(38)

令每個(gè)新子區(qū)間與當(dāng)前子區(qū)間的配點(diǎn)數(shù)相同，則區(qū)間分段數(shù)Hk為

(39)

4.6 網(wǎng)格縮減

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

合并相鄰子區(qū)間的基本思想與上述算法類(lèi)似，只是將相鄰子區(qū)間的值在當(dāng)前區(qū)間對(duì)應(yīng)多項(xiàng)式上進(jìn)行外插，然后類(lèi)似地判斷容許誤差是否滿足。

網(wǎng)格縮減算法的具體推導(dǎo)過(guò)程見(jiàn)文獻(xiàn)[14]，本文不再贅述。

4.7 算法流程

本文自適應(yīng)偽譜法的計(jì)算流程如下：

步驟1給定求解容許誤差ε和初始網(wǎng)格。

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

步驟5重復(fù)步驟2。

綜上，本文算法的計(jì)算流程圖見(jiàn)圖1。

圖1 自適應(yīng)偽譜法流程圖Fig.1 Flowchart of adaptive pseudospectral method

5 算例分析

再入軌跡優(yōu)化算例設(shè)置見(jiàn)表1，邊界條件設(shè)置見(jiàn)表2。計(jì)算環(huán)境為Windows XP 2.83 GHz，3.25 GB內(nèi)存，軟件環(huán)境為MATLAB R2014a，NLP求解器采用SNOPT v7.0，梯度計(jì)算采用有限差分算法。

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

表1 算例設(shè)置Table 1 Setup of calculation examples

表2 邊界條件設(shè)置Table 2 Setup of boundary conditions

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

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

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

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

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

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

表3 不同算法軌跡優(yōu)化結(jié)果Table 3 Trajectory optimization results by different methods

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

6 結(jié) 論

本文針對(duì)高超聲速飛行器再入段軌跡優(yōu)化問(wèn)題，以CAV-H飛行器對(duì)象，基于多區(qū)間Radau偽譜法，對(duì)3種典型性能指標(biāo)的連續(xù)時(shí)間最優(yōu)控制問(wèn)題進(jìn)行離散，基于Legendre多項(xiàng)式近似理論構(gòu)建相對(duì)誤差估計(jì)關(guān)系式，采用自適應(yīng)網(wǎng)格重構(gòu)技術(shù)實(shí)現(xiàn)了再入軌跡優(yōu)化的高效求解。仿真結(jié)果表明：

1) 本文算法結(jié)果與變步長(zhǎng)Runge-Kutta-Fehlberg法積分結(jié)果一致。

2) 本文算法較2種傳統(tǒng)自適應(yīng)偽譜法，配點(diǎn)與區(qū)間分配更合理，計(jì)算效率更高，且對(duì)于人工參數(shù)不敏感，具有較好的工程應(yīng)用價(jià)值。