李兆亭， 周 祥， 張洪波， 湯國建

(國防科技大學(xué) 空天科學(xué)學(xué)院，長(zhǎng)沙 410073)

再入可達(dá)域是指可重復(fù)使用飛行器在滿足約束條件下，到達(dá)著陸或者交班狀態(tài)下所能允許的區(qū)域范圍.它是飛行器縱向機(jī)動(dòng)能力與橫向機(jī)動(dòng)能力的綜合反映，體現(xiàn)了飛行器的整體機(jī)動(dòng)性能.可達(dá)域求解可為飛行器軌跡規(guī)劃與制導(dǎo)、應(yīng)急著陸點(diǎn)選擇、制導(dǎo)律性能評(píng)估提供依據(jù)，因此具有重要意義.

再入可達(dá)域的求解方法按照原理可以分為兩類：參數(shù)化方法和非參數(shù)化方法.參數(shù)化方法是指將再入軌跡規(guī)劃問題的狀態(tài)量或控制量等進(jìn)行參數(shù)化，采用某些優(yōu)化方法對(duì)參數(shù)化后的問題進(jìn)行優(yōu)化求解；非參數(shù)化方法則不包含此類參數(shù)化、問題優(yōu)化過程，常見的基于剖面的可達(dá)域求解方法就屬于此類.如曾夕娟等[1]在“r-V”空間內(nèi)引入約束的對(duì)數(shù)表達(dá)式簡(jiǎn)化再入走廊的表達(dá)，將再入軌跡分為3段設(shè)計(jì)，引入準(zhǔn)平衡滑翔條件獲取解析制導(dǎo)律以跟蹤所規(guī)劃軌跡，通過平移“邊界跟蹤段”軌跡使其布滿再入走廊，從而得到最終的可達(dá)域.常松濤等[2]在基于E剖面的再入走廊下，引入了傾側(cè)角反向策略，改進(jìn)了阻力加速度曲線插值方式和最大可行阻力加速度曲線的生成方法，使計(jì)算得到的可達(dá)區(qū)域更精確、邊界更均勻、彈道更平滑.He等[3]推導(dǎo)得到一個(gè)具有不確定性效應(yīng)的再入走廊，其中攻角剖面可調(diào)節(jié)，縱向阻力剖面設(shè)計(jì)為上、下擬合安全邊界的插值結(jié)果，橫向升力-阻力剖面則利用準(zhǔn)平衡滑翔條件得到，最后通過比例積分微分(PID)跟蹤器來跟蹤剖面，得到可達(dá)域.Zhang等[4]提出一種基于幾何預(yù)測(cè)軌跡的可達(dá)域計(jì)算方法，在D-E走廊下，對(duì)橫向運(yùn)動(dòng)狀態(tài)、轉(zhuǎn)彎半徑和角度進(jìn)行分析求解，得到D-E剖面的最遠(yuǎn)可達(dá)邊界，并提出幾何預(yù)測(cè)的軌跡；通過重復(fù)計(jì)算所有D-E剖面的可達(dá)邊界，生成可達(dá)域.上述基于剖面的可達(dá)域求解方法通常需要人工設(shè)計(jì)相關(guān)剖面，而剖面參數(shù)決定了其可達(dá)域范圍，因此得到的可達(dá)域通常只是實(shí)際可達(dá)域的一部分.此外，Liu等[5]提出一種使用多項(xiàng)式來描述可達(dá)域邊界的方法，使用反向傳播神經(jīng)網(wǎng)絡(luò)快速確定空對(duì)地飛行器的可達(dá)域.Zhang等[6]通過快速解析的方法求解3個(gè)不同的最優(yōu)控制問題，得到可達(dá)域的內(nèi)、外邊界.章吉力等[7]對(duì)考慮禁飛區(qū)條件下的空天飛機(jī)再入可達(dá)區(qū)域問題進(jìn)行研究，并基于極限繞飛軌跡提出一種不限禁飛區(qū)位置的可達(dá)區(qū)域求解方法.吳楠等[8]基于平衡滑翔假設(shè)和最優(yōu)化飛行準(zhǔn)則，通過數(shù)值積分獲得最大縱程和橫程彈道，基于橢圓近似法利用3個(gè)末端點(diǎn)構(gòu)建可達(dá)域橢圓邊界.

可達(dá)域的參數(shù)化求解方法包含的種類較多，常見方法主要包括基于虛擬目標(biāo)點(diǎn)[9-11]、基于偽譜法[12-14]和基于凸優(yōu)化方法[15]等.如趙吉松等[16]提出一種基于稀疏差分法和網(wǎng)格細(xì)化技術(shù)的快速、高精度軌跡求解方法，并利用可達(dá)域的求解驗(yàn)證了方法的有效性.藺君等[17]對(duì)攻角、傾側(cè)角進(jìn)行參數(shù)化，利用帶約束的差分進(jìn)化算法求解滿足再入過程約束和終端約束的再入軌跡.梁巨平等[18]選取傾側(cè)角為時(shí)間的分段常值函數(shù)，采用遺傳算法來求解一系列的軌跡優(yōu)化問題.趙江等[19]設(shè)計(jì)了一種參數(shù)化的傾側(cè)角剖面，利用約束粒子群優(yōu)化算法求解滿足再入過程約束和末端約束的最優(yōu)滑翔軌跡.相比基于剖面的求解方法，可達(dá)域的參數(shù)化求解方法不需要設(shè)計(jì)相關(guān)剖面，其可達(dá)域范圍更大且更接近實(shí)際值，但受限于具體算法的性能，其全局最優(yōu)性仍難以保證.孫勇等[20]基于擬能量概念提出了多段軌跡優(yōu)化方法，將軌跡按照擬能量進(jìn)行分段，對(duì)各段分別進(jìn)行優(yōu)化，減少了每次優(yōu)化的變量個(gè)數(shù)，計(jì)算效率得到較大提高.

偽譜法具有良好的適用性和強(qiáng)大的求解能力，基于偽譜法求解可達(dá)域的方案較多.如王濤等[21]提出一種基于Gauss偽譜法的再入可達(dá)域計(jì)算方法，采用固定的攻角剖面，僅對(duì)傾側(cè)角進(jìn)行單變量尋優(yōu).李柯等[22]基于Radau偽譜法求解得到飛行器的可達(dá)域，著重分析了升阻比、終端速度以及終端傾角對(duì)可達(dá)域的影響.Wang等[23]在準(zhǔn)平衡滑行條件下，使用高斯偽譜法計(jì)算得到了返回航天器的可達(dá)域.解永鋒等[24]通過定義加權(quán)的橫程、縱程組合性能指標(biāo)函數(shù)，將可達(dá)域求解問題轉(zhuǎn)化為組合性能指標(biāo)最優(yōu)的控制問題，采用勒讓德偽譜法快速計(jì)算可達(dá)域.

本文繼承發(fā)展了文獻(xiàn)[21]的求解思路，將攻角、傾側(cè)角同時(shí)作為控制量代入到優(yōu)化問題中，采用偽譜法進(jìn)行軌跡規(guī)劃，得到的可達(dá)域范圍更廣且更接近實(shí)際值.此外，描述了基于偽譜法的可達(dá)域快速求解方法，并對(duì)偽譜法進(jìn)行了簡(jiǎn)要闡述；然后基于上述方法對(duì)影響可達(dá)域的相關(guān)因素進(jìn)行仿真研究.仿真結(jié)果表明，飛行器質(zhì)量、氣動(dòng)參考面積、大氣密度等會(huì)對(duì)短縱程軌跡產(chǎn)生明顯影響；同時(shí)，升阻比大小與可達(dá)域范圍成正相關(guān).

1 再入飛行器及其可達(dá)域描述

介紹再入飛行器的動(dòng)力學(xué)模型，并在換極坐標(biāo)系下給出利用偽譜法求解可達(dá)域的思路.

1.1 動(dòng)力學(xué)方程及約束

在求解可達(dá)域時(shí)，通常將飛行器的運(yùn)動(dòng)轉(zhuǎn)換到換極坐標(biāo)系[25]中.在換極坐標(biāo)系中，飛行器起點(diǎn)的經(jīng)緯度均為0°，終點(diǎn)的經(jīng)度和緯度分別描述了再入縱程和再入橫程.利用這些特性，可方便地標(biāo)定偽譜法中的參數(shù)范圍，簡(jiǎn)化彈道規(guī)劃算法.

在換極坐標(biāo)系下，可以推導(dǎo)得到飛行器再入運(yùn)動(dòng)方程：

(1)

(2)

式中：ρ為大氣密度；M為飛行器質(zhì)量；S為參考面積；CL和CD分別為升力系數(shù)和阻力系數(shù)，均與攻角α有關(guān).同時(shí)，再入飛行過程中還需考慮多種過程約束，如熱流密度約束、過載約束、動(dòng)壓約束等，具體為

(3)

在本文中，目標(biāo)飛行器為CAV-H，其相關(guān)參數(shù)和氣動(dòng)數(shù)據(jù)可參考文獻(xiàn)[26].

1.2 換極坐標(biāo)系下的可達(dá)域求解

為了加快優(yōu)化速度和提高優(yōu)化質(zhì)量，文獻(xiàn)[21]提出了一種基于Gauss偽譜法的再入可達(dá)域計(jì)算方法，采用固定的攻角剖面，將攻角作為狀態(tài)變量，僅對(duì)傾側(cè)角進(jìn)行單變量尋優(yōu).本文延續(xù)該求解思路，但與之不同的是，攻角剖面并不是固定的，即攻角、傾側(cè)角將同時(shí)作為待優(yōu)化變量.

再入可達(dá)域問題可表述為求解不同縱程條件下的最大橫程.在換極坐標(biāo)系內(nèi)，落點(diǎn)的橫程可用緯度表征:

J=min(cosφpf)

(4)

式中：φpf為換極坐標(biāo)系中終點(diǎn)的緯度，如圖1所示.

圖1中，Lset為設(shè)定的縱程，λpf為換極坐標(biāo)系中終點(diǎn)的經(jīng)度，二者數(shù)值相等.不妨將Lset視為終端約束，通過改變Lset值得到不同縱程條件下的最大橫程軌跡.同時(shí)將換極坐標(biāo)系下，φpf>0°的所有軌跡稱之為北半球軌跡，反之稱之為南半球軌跡.在進(jìn)行偽譜法的參數(shù)設(shè)置時(shí)，考慮到地球自轉(zhuǎn)的影響，即微分方程(1)中的科氏加速度項(xiàng)和牽連加速度項(xiàng)均不為0，不妨將把北半球軌跡的緯度變化范圍設(shè)置為 [-2°, 90°]，將南半球軌跡的緯度變化范圍設(shè)置為 [-90°, 2°].

圖1 不同縱程下的最大橫程Fig.1 Maximum cross-ranges of different longitudinal ranges

求解上述不同縱程下的最大橫程問題，即可得到有效的可達(dá)域.但實(shí)際仿真中，同樣需要確定最大、最小縱程(Lmax和Lmin)，以便確定Lset的取值范圍.因此，最終的規(guī)劃流程如圖2所示.首先求解最大、最小縱程，標(biāo)定Lmax和Lmin的具體數(shù)值；然后從縱程Lmax-ΔL開始，分別求解兩個(gè)最大橫程問題，依次減小縱程ΔL，直到達(dá)到Lmin；最后即可得到飛行器的再入可達(dá)域.

圖2 數(shù)值求解可達(dá)域的一般流程Fig.2 General procedure for numerically solving footprint

2 偽譜軌跡規(guī)劃方法

介紹偽譜法的原理，描述將原問題轉(zhuǎn)化為非線性規(guī)劃(NLP)問題的過程，給出偽譜法求解問題的一般流程.

2.1 偽譜離散化方法

偽譜法作為直接法中的典型代表，其原理與直接法一致：基于離散化一組節(jié)點(diǎn)上的控制變量或狀態(tài)變量，將最優(yōu)控制問題轉(zhuǎn)換為NLP問題.不同的是，偽譜法使用全局多項(xiàng)式同時(shí)對(duì)狀態(tài)量和控制量進(jìn)行參數(shù)化，并通過正交配點(diǎn)來近似微分代數(shù)方程.其一般流程如下所示.

首先引入新的時(shí)間變量τ，對(duì)原問題做如下轉(zhuǎn)換：

(5)

式中：t0為初始時(shí)間；tf為終端時(shí)間.轉(zhuǎn)換后得到的狀態(tài)量、控制量均采用拉格朗日多項(xiàng)式進(jìn)行逼近.

(6)

式中：τ∈[-1, +1]；x為實(shí)際的狀態(tài)量；X為多項(xiàng)式擬合的狀態(tài)量；lj(τ),j=1, 2, …,Nk+1為拉格朗日多項(xiàng)式；τ1,τ2, …,τNk為配點(diǎn)，通?？梢赃x擇Legendre-Gauss(LG)、Legendre-Gauss-Radau(LGR)、Legendre-Gauss-Lobatto(LGL)等正交配點(diǎn)[27].對(duì)于拉格朗日多項(xiàng)式，顯然有：

(7)

對(duì)拉格朗日多項(xiàng)式進(jìn)行微分，可以得到微分矩陣算子：

(8)

DX-f(X)=0

(9)

偽譜法中對(duì)于動(dòng)力學(xué)微分方程的處理如式(9)所示，其作為等式約束與其他約束一起組成最終的NLP問題.可以看到，由于引入了微分方程的右側(cè)項(xiàng)f(X)，此約束一般而言是非線性約束，所以最終形成的離散問題，也就是NLP問題.關(guān)于再入過程中的其他約束，如式(3)所示，可以通過狀態(tài)量、控制量的擬合式(6)進(jìn)行顯示表達(dá).

最終，原再入軌跡規(guī)劃問題就被轉(zhuǎn)化為NLP問題.實(shí)際上，考慮到軌跡規(guī)劃問題的復(fù)雜性，求解時(shí)不僅需要增加多項(xiàng)式的階數(shù)，通常也需要增加多項(xiàng)式的個(gè)數(shù)，即為常用的ph網(wǎng)格細(xì)化方法，將歸一化后的總時(shí)間[-1, 1]劃分為多個(gè)區(qū)間，每一個(gè)區(qū)間對(duì)應(yīng)一個(gè)拉格朗日多項(xiàng)式.

2.2 NLP問題及其求解流程

利用上述偽譜離散化方法，可以將原問題轉(zhuǎn)化為NLP問題.通常可以將NLP問題描述為

minΦ(Z)

(10)

s.t.Fmin≤F(Z)≤Fmax

Zmin≤Z≤Zmax

式中：Φ為性能指標(biāo)函數(shù)；Z為自變量，通常包括狀態(tài)量、控制量、時(shí)間變量、參數(shù)變量等；Zmax和Zmin分別為上述各變量的上、下限；F為等式和不等式約束，通常包含關(guān)于動(dòng)力學(xué)微分方程的等式約束、動(dòng)壓約束、過載約束、熱流約束等；Fmax和Fmin分別對(duì)應(yīng)約束的上、下限，特別地，等式約束對(duì)應(yīng)的上、下限均為0.

得到如式(10)所示的NLP問題后，可以采用通用的非線性規(guī)劃求解器對(duì)問題進(jìn)行求解，如稀疏非線性優(yōu)化器(SNOPT)[28]或內(nèi)點(diǎn)優(yōu)化器(IPOPT)[29]等.為使NLP問題能夠被求解，需要即時(shí)計(jì)算一些關(guān)鍵信息，如梯度、雅可比矩陣、海森矩陣等.此類信息的解算速度嚴(yán)重影響了最終的計(jì)算時(shí)間，為此，文獻(xiàn)[30]推導(dǎo)出了NLP問題的梯度、雅可比矩陣和拉格朗日海森矩陣的明確表達(dá)式，確定了NLP問題的稀疏結(jié)構(gòu)，可以顯著提高NLP問題的計(jì)算效率和可靠性.

在NLP問題被求解后，通常來說此時(shí)的解很難滿足期望的精度要求.因此，需要對(duì)網(wǎng)格進(jìn)行細(xì)化，在更細(xì)密的網(wǎng)格下離散原問題，并求解得到更大規(guī)模的NLP問題.如此，上述網(wǎng)格被不斷細(xì)化，直到滿足迭代截止條件.

綜上所述，偽譜法求解優(yōu)化問題的流程如圖3所示.圖中，原問題首先在初始網(wǎng)格條件下進(jìn)行離散化得到NLP問題，然后調(diào)用相關(guān)求解器求解此類問題，最后不斷加密網(wǎng)格并求解離散化得到的NLP問題，直到滿足截止條件.需要注意的是，由于稀疏網(wǎng)格條件下NLP問題與原問題不等價(jià)，NLP問題的收斂容限可以適當(dāng)增加以減少計(jì)算時(shí)間、提高求解效率.

圖3 偽譜法求解優(yōu)化問題的一般流程Fig.3 General flow of pseudospectral method for solving optimization problems

3 仿真分析

將上述偽譜法應(yīng)用于圖3所示的可達(dá)域求解方案中，即可快速求解再入飛行器的可達(dá)域.

3.1 標(biāo)準(zhǔn)情況下的可達(dá)域求解

飛行器在初始時(shí)刻的經(jīng)度、緯度被設(shè)定為(0°，0°)，仿真采用的大氣環(huán)境為指數(shù)型大氣模型，其他參數(shù)如表1所示.仿真計(jì)算機(jī)硬件條件為AMD Ryzen 7 5800H@3.2 GHz，偽譜法求解器SWIFT由本小組獨(dú)立開發(fā).

表1中給出了偽譜法中一些參數(shù)的設(shè)置范圍.其中，h為飛行器的滑翔高度，φS為求解換極坐標(biāo)系下南半球的飛行軌跡時(shí)的緯度，φN為求解換極坐標(biāo)系下北半球的飛行軌跡時(shí)的緯度.為保證問題的全局最優(yōu)性，部分參數(shù)范圍取得較寬泛，如速度方位角、速度傾角等.

表1 偽譜法求解可達(dá)域時(shí)的部分參數(shù)設(shè)置Tab.1 Parameters for solving footprint by pseudospectral method

此外，一些其他參數(shù)設(shè)置如下：初始網(wǎng)格設(shè)置為2×2，即將總時(shí)間分為2段，每一段的軌跡采用2階多項(xiàng)式進(jìn)行逼近；最大迭代次數(shù)設(shè)置為5；原問題期望的收斂容限設(shè)置為10-4；最大飛行時(shí)間設(shè)置為 4 500 s；過載不超過3 g，動(dòng)壓不超過100 kPa，熱流不超過150 kW/m2.最終得到的可達(dá)域仿真結(jié)果如圖4所示.

圖4(a)～(c)中， 紅色、綠色、藍(lán)色曲線分別代表最大縱程軌跡、北半球最小縱程軌跡、南半球最小縱程軌跡，其中最大縱程約為136.4°，南北半球最小縱程約為3.62°；32條黑色曲線代表不同縱程條件下的最值橫程軌跡，其中所有軌跡的最值橫程約為±54.52°.將所有軌跡的終點(diǎn)連接在一起，即可形成最終的飛行器可達(dá)域.可見，標(biāo)準(zhǔn)情況下此類飛行器的可達(dá)域接近于橢圓形.同時(shí)給出文獻(xiàn)[21]中帶有攻角剖面的可達(dá)域計(jì)算結(jié)果，如圖4(d)中黃色實(shí)線所示，可見受攻角剖面的限制，這種情況下得到的可達(dá)域相對(duì)較小.此時(shí)，最大縱程約為112.6°，相比136.4°，損失約17.45%.最值橫程約為±41.53°，相比±54.52°，損失23.83%.

整體而言，對(duì)于構(gòu)成可達(dá)域的圖4(a)～(c)中的35條軌跡，其速度、高度、速度方位角等隨時(shí)間變化較為規(guī)律.對(duì)于速度而言，縱程越大，速度減小得越慢.同時(shí)從4(b)也可以看到，速度減小變慢的原因在于飛行器的滑翔高度較高；從4(c)可知，整個(gè)可達(dá)域條件下的速度方位角在-90°～270°之間變化.

圖4 各條軌跡的狀態(tài)量變化Fig.4 State change for each trajectory

給出上述35條規(guī)劃軌跡的控制量變化，如圖5所示.其中，攻角的變化總是以一個(gè)較大的值開始，逐步過渡到最大升阻比攻角附近，最后階段受表1中終端速度傾角范圍設(shè)置的影響，攻角快速拉升到30°附近，以抬高速度傾角至設(shè)定范圍.傾側(cè)角的變化則明顯分為南、北半球兩部分.當(dāng)規(guī)劃北半球軌跡時(shí)，傾側(cè)角多數(shù)情況下小于0°；南半球軌跡則相反.同時(shí)可以看到，縱程越大的軌跡，起點(diǎn)處的傾側(cè)角絕對(duì)值越小，相對(duì)時(shí)間的平均斜率的絕對(duì)值也越?。欢v程越小的軌跡，起點(diǎn)處的傾側(cè)角絕對(duì)值越大，相對(duì)時(shí)間的平均斜率的絕對(duì)值也越大.

圖5 各條軌跡的控制量變化Fig.5 Control changes for each trajectory

3.2 參數(shù)拉偏條件下的可達(dá)域

研究某些參數(shù)對(duì)可達(dá)域的影響，將某些參數(shù)進(jìn)行拉偏，在其他參數(shù)不變的條件下求解可達(dá)域.一些拉偏情形如表2所示；仿真得到的可達(dá)域隨拉偏參數(shù)的變化如圖6所示.

表2 5種參數(shù)拉偏情形Tab.2 Five cases of parameter deviation

圖6(a)中，在質(zhì)量拉偏±5%之間，飛行器的可達(dá)域基本沒有變化；而當(dāng)質(zhì)量拉偏+10%時(shí)，短縱程軌跡的某些約束達(dá)到上限，其最小縱程有了明顯變化，但最大縱程基本沒有變化.同樣的現(xiàn)象出現(xiàn)在圖6(b)和6(c)中，氣動(dòng)參考面積和大氣密度在±5%之間波動(dòng)時(shí)，可達(dá)域不發(fā)生變化；當(dāng)拉偏在 -8% 時(shí)，其最小縱程有了明顯的損失，而最大縱程不變.

圖6 可達(dá)域隨參數(shù)拉偏情況的變化Fig.6 Variation of footprint with parameter deviation

結(jié)合式(1)～(3)可知，飛行器的質(zhì)量、氣動(dòng)參考面積、大氣密度僅出現(xiàn)在式(2)和式(3)中，并不顯含在動(dòng)力學(xué)微分方程(1)中.并且質(zhì)量、面積、大氣密度的變化僅會(huì)導(dǎo)致升力、阻力的同步變化，因此對(duì)于飛行器某一時(shí)刻的飛行狀態(tài)而言，通過適當(dāng)改變攻角、傾側(cè)角，就可以使得式(1)中的左側(cè)微分項(xiàng)完全一致，則在一定范圍內(nèi)，可達(dá)域是不變的.而超出此范圍后，受到熱流、過載、動(dòng)壓等約束的影響，部分軌跡，特別是短縱程軌跡，可能存在達(dá)到約束上界的情況，因此必然要通過增加縱程以降低過程約束.

此外，圖6(d)為升力和阻力系數(shù)拉偏條件下的可達(dá)域計(jì)算結(jié)果.其中，實(shí)線代表升力系數(shù)拉偏的可達(dá)域，虛線代表阻力系數(shù)拉偏的可達(dá)域.可見，升力系數(shù)的增加和阻力系數(shù)的減小均會(huì)導(dǎo)致可達(dá)域范圍增加，反之亦然.同時(shí)，升力系數(shù)拉偏+10%和阻力系數(shù)拉偏-10% 的結(jié)果一致，升力系數(shù)拉偏-10% 和阻力系數(shù)拉偏 +10% 的結(jié)果也一致，說明升力、阻力系數(shù)的變化對(duì)可達(dá)域的影響是等價(jià)的.更準(zhǔn)確地，升阻比的大小影響了可達(dá)域的大?。荷璞仍酱螅蛇_(dá)域越大；升阻比越小，可達(dá)域越小.

4 結(jié)語

本文描述了基于偽譜法的可達(dá)域快速求解方法，并基于此方法對(duì)影響可達(dá)域的相關(guān)因素進(jìn)行了研究與仿真.仿真結(jié)果表明，在給定參數(shù)條件下，飛行器質(zhì)量、氣動(dòng)參考面積、大氣密度等在±5%范圍內(nèi)不會(huì)導(dǎo)致可達(dá)域的改變；超出一定范圍后，如飛行器質(zhì)量增加10%、氣動(dòng)參考面積減少8%、大氣密度減少8%等，均會(huì)對(duì)短縱程軌跡產(chǎn)生明顯影響，且體現(xiàn)在可達(dá)域的左半部分，其右半部分不受影響；升阻比對(duì)可達(dá)域的影響較大，其大小與可達(dá)域范圍成正相關(guān).