周祥，張洪波,*，何睿智，湯國(guó)建，包為民

1. 國(guó)防科技大學(xué) 空天科學(xué)學(xué)院，長(zhǎng)沙 410073 2. 中國(guó)航天科技集團(tuán)有限公司科技委，北京 100048

可重復(fù)使用飛行器能夠?qū)崿F(xiàn)天地間的多次往返，從空間再入大氣層后，依靠大升阻比氣動(dòng)外形實(shí)現(xiàn)在臨近空間的遠(yuǎn)距離滑翔飛行，從而抵達(dá)預(yù)定目標(biāo)[1]。再入滑翔段軌跡規(guī)劃的準(zhǔn)確性和快速性對(duì)于改善該類飛行器的制導(dǎo)效果具有重要意義。

為了降低軌跡規(guī)劃的計(jì)算量，以航天飛機(jī)軌道器為代表的可重復(fù)使用飛行器采用固定攻角剖面，以傾側(cè)角為唯一控制量的方式進(jìn)行軌跡規(guī)劃，進(jìn)而實(shí)現(xiàn)制導(dǎo)任務(wù)。這種基于二維剖面的規(guī)劃方法具有計(jì)算負(fù)擔(dān)小、方法簡(jiǎn)單、容易實(shí)現(xiàn)等優(yōu)點(diǎn)，在航天飛機(jī)等飛行器上得到了成功應(yīng)用。然而，這種計(jì)算量小的規(guī)劃方法難以充分發(fā)揮可重復(fù)使用飛行器的大范圍機(jī)動(dòng)能力。飛行器固有的機(jī)動(dòng)能力主要由飛行器的總體布局和氣動(dòng)性能決定，上述因素確定后，如何選擇合適的控制策略從而充分發(fā)揮出其固有的機(jī)動(dòng)能力就顯得尤為重要。1998年，文獻(xiàn)[2]提出如果將傳統(tǒng)二維剖面規(guī)劃方法中固定攻角的條件放開，同時(shí)采用攻角和傾側(cè)角作為控制量，將有助于充分發(fā)揮飛行器的機(jī)動(dòng)能力。文獻(xiàn)[3]將同時(shí)考慮攻角和傾側(cè)角作為控制量的軌跡規(guī)劃問題描述為三維剖面規(guī)劃問題，使用阻力加速度和側(cè)向加速度作為中間控制變量(即剖面設(shè)計(jì)量)，引入降階動(dòng)力學(xué)系統(tǒng)降低問題求解的計(jì)算量，這里的三維剖面從控制量的角度講即是對(duì)攻角和傾側(cè)角變化規(guī)律的一種表征，而傳統(tǒng)的二維剖面規(guī)劃方法只能獲得對(duì)傾側(cè)角變化規(guī)律的表征。文獻(xiàn)[4]在加速度剖面規(guī)劃中同時(shí)考慮縱向和側(cè)向運(yùn)動(dòng)，首先基于二維剖面規(guī)劃一條可行軌跡，其次基于三維剖面獲得最優(yōu)軌跡。文獻(xiàn)[5]提出了一種分層策略的三維剖面規(guī)劃方法，該方法采用側(cè)向降階運(yùn)動(dòng)模型對(duì)縱向、側(cè)向剖面進(jìn)行反復(fù)迭代，直至滿足終端縱程和橫程要求。文獻(xiàn)[6]將基于分層策略的三維剖面規(guī)劃方法應(yīng)用到再入覆蓋區(qū)求解中，獲得了相比基于二維剖面規(guī)劃方法更大的覆蓋區(qū)，驗(yàn)證了機(jī)動(dòng)能力得到充分發(fā)揮后的效果。整體來看，目前對(duì)二維剖面規(guī)劃方法的研究較多[7-9]，而關(guān)于三維剖面規(guī)劃方法的研究還比較少，還有很多技術(shù)問題有待解決，其中，計(jì)算效率是限制三維剖面規(guī)劃技術(shù)的主要瓶頸之一。

近年來，內(nèi)點(diǎn)法的快速發(fā)展使得凸優(yōu)化在求解效率方面相比傳統(tǒng)方法如偽譜法的優(yōu)勢(shì)愈發(fā)明顯[10-11]。除此之外，凸優(yōu)化問題的局部極小點(diǎn)就是全局最小點(diǎn)，這一良好的理論性質(zhì)也為該方法的廣泛應(yīng)用奠定了基礎(chǔ)[12]。由于實(shí)際中大量的優(yōu)化問題都是非凸的，如何將軌跡規(guī)劃等這類具有多種約束的原始非凸問題轉(zhuǎn)化為凸問題是應(yīng)用凸優(yōu)化方法的關(guān)鍵，包括對(duì)非線性動(dòng)力學(xué)、過程約束、控制量約束、優(yōu)化目標(biāo)的凸化。截止目前，凸優(yōu)化已經(jīng)在再入軌跡規(guī)劃、行星著陸、無人機(jī)路徑規(guī)劃等方面得到了應(yīng)用，一系列凸化方法被廣泛使用在這些問題中，包括約束松弛和連續(xù)線性化等技術(shù)[13]。凸優(yōu)化早期在航天中的應(yīng)用主要是解決動(dòng)力下降段的軌跡規(guī)劃問題，在這類問題的凸化過程中提出了無損凸化的概念[14-15]。

文獻(xiàn)[16]利用序列凸優(yōu)化的方法求解再入軌跡規(guī)劃問題，相比偽譜法顯示出了很好的求解效率，文中所提的序列凸化方法是目前處理強(qiáng)非線性動(dòng)力學(xué)約束與非線性過程約束的主要方法。文獻(xiàn)[17]中討論了利用序列凸優(yōu)化求解再入軌跡規(guī)劃問題的收斂性問題，序列凸優(yōu)化的收斂性與解的等價(jià)性緊密相關(guān)。文獻(xiàn)[18]實(shí)現(xiàn)了基于凸優(yōu)化在線軌跡規(guī)劃的閉環(huán)跟蹤制導(dǎo)，將軌跡規(guī)劃子問題和軌跡跟蹤制導(dǎo)指令求解子問題均轉(zhuǎn)化為凸問題，從而可實(shí)現(xiàn)在線快速求解，取得了很好的制導(dǎo)效果。文獻(xiàn)[19-21]對(duì)于改進(jìn)序列凸優(yōu)化方法的性能提供了多種技巧。文獻(xiàn)[22]利用凸優(yōu)化求解再入軌跡規(guī)劃問題時(shí)將攻角與傾側(cè)角同時(shí)作為控制量，通過定義中間變量實(shí)現(xiàn)了對(duì)控制量約束的凸化處理。但并未給出規(guī)劃后的指令計(jì)算方法，限制了該方法在實(shí)際問題中的應(yīng)用。

總體來看，目前基于三維剖面的再入軌跡規(guī)劃方法研究還比較少，主要存在的問題包括：① 動(dòng) 力學(xué)模型簡(jiǎn)化較多，導(dǎo)致問題求解精度較低；② 為了減低計(jì)算量，將加速度剖面進(jìn)行簡(jiǎn)單參數(shù)化，人為限制了剖面求解的搜索空間，不利于充分發(fā)揮飛行器的能力?？偟膩砜?，已有三維剖面規(guī)劃方法并未很好解決軌跡規(guī)劃的求解精度與計(jì)算效率之間的矛盾。

基于上述分析，本文采用完整動(dòng)力學(xué)模型，利用凸優(yōu)化方法求解再入三維剖面規(guī)劃問題。主要貢獻(xiàn)是：① 將原始三維剖面規(guī)劃問題描述為最優(yōu)控制問題，采用合適的凸化技巧將考慮多種約束條件的原問題轉(zhuǎn)化為凸優(yōu)化問題，包括動(dòng)力學(xué)約束、路徑約束和控制量約束的凸化，其次對(duì)連續(xù)凸問題進(jìn)行離散化；② 引入指令反解步驟，通過優(yōu)化后的中間控制量計(jì)算攻角和傾側(cè)角指令，并建立了相應(yīng)的序列迭代算法，通過迭代求解凸優(yōu)化子問題，從而獲得原問題的可行解。通過數(shù)值仿真證明了本文方法的可行性和收斂性。

本文的框架結(jié)構(gòu)：引言部分介紹研究背景和研究現(xiàn)狀；第1部分介紹問題描述；第2部分介紹如何對(duì)原問題進(jìn)行凸化與離散處理，將其轉(zhuǎn)化為一個(gè)凸的數(shù)學(xué)規(guī)劃問題；第3部分介紹如何利用凸優(yōu)化的解(即中間變量和狀態(tài)量)反解攻角和傾側(cè)角指令，并給出了序列凸化算法的步驟；第4部分進(jìn)行了數(shù)值仿真；第5部分給出了結(jié)論。

1 問題描述

1.1 三維剖面規(guī)劃問題

首先建立軌跡規(guī)劃中使用的動(dòng)力學(xué)模型。在建模過程中主要引入以下假設(shè)： ① 考慮地球自轉(zhuǎn)；② 僅考慮球形引力；③ 認(rèn)為地球外形是一個(gè)標(biāo)準(zhǔn)圓球?？芍貜?fù)使用飛行器在再入過程中進(jìn)行無動(dòng)力飛行，由于大氣阻力的作用，飛行器能量是單調(diào)遞減的，為了降低優(yōu)化變量的數(shù)目，本文采用能量E作為自變量，定義為

E=V2/2-μ/r

(1)

式中：μ為地球引力常數(shù)；r為地心距；V為速度。

進(jìn)而得到以能量為自變量的再入質(zhì)心動(dòng)力學(xué)模型為

(2)

式中：λ為經(jīng)度；φ為緯度；θ為當(dāng)?shù)厮俣葍A角；σ為航跡偏航角；L、D分別為升力和阻力加速度；υ為傾側(cè)角；由地球自轉(zhuǎn)引起的哥氏加速度項(xiàng)Cθ和Cσ以及牽連加速度項(xiàng)C′θ和C′σ可分別表示為

(3)

上述所有變量均可采用如下歸一化方法進(jìn)行無量綱化處理。

(4)

(5)

式中:R0為地球平均半徑;g0為R0處的引力;ωe為地球旋轉(zhuǎn)角速度。升力和阻力加速度的定義為

(6)

式中：Mr為飛行器質(zhì)量；Sr為參考面積；CL、CD分別為升力系數(shù)和阻力系數(shù)，通常為攻角α和馬赫數(shù)Ma的函數(shù)；ρ為大氣密度，有

ρ=ρ0e-β h

(7)

其中：β=1/hs，hs=7 110 m；ρ0=1.225 kg/m3，h為高度。如無特殊說明，下文所有變量均為無量綱化處理后的值。

以E為自變量的動(dòng)力學(xué)方程包括5個(gè)狀態(tài)參數(shù)和2個(gè)控制指令。前者指地心距r、經(jīng)度λ、地心緯度φ、當(dāng)?shù)厮俣葍A角θ和航跡偏航角σ，而速度與自變量有關(guān)，可以通過式(8)確定

(8)

將獨(dú)立描述飛行器運(yùn)動(dòng)狀態(tài)的參數(shù)構(gòu)成狀態(tài)矢量x，可表示為

x=[r,λ,φ,θ,σ]T

(9)

2個(gè)控制指令是攻角和傾側(cè)角，利用它們實(shí)現(xiàn)對(duì)氣動(dòng)力大小和方向的控制，所以再入軌跡規(guī)劃中的控制矢量U′可表示為

U′=[α,υ]T

(10)

在軌跡規(guī)劃中需要對(duì)控制指令的幅值進(jìn)行約束，可以表示為

αmin≤α≤αmax

(11)

υmin≤|υ|≤υmax

(12)

式中：|·|表示絕對(duì)值，υmax>υmin≥0。

過程約束主要有：熱流密度、動(dòng)壓、法向過載和氣動(dòng)力，前3個(gè)約束用于避免飛行器結(jié)構(gòu)受到破壞，第4個(gè)約束用于確保飛行器有足夠氣動(dòng)力用于軌跡控制，分別表示為

(13)

(14)

(15)

(16)

軌跡規(guī)劃的邊界條件包括初始狀態(tài)約束和終端狀態(tài)約束，表示為

(17)

式中：x0和xf由具體任務(wù)給定;E0、Ef分別為初始狀態(tài)和終端狀態(tài)對(duì)應(yīng)的能量。

優(yōu)化目標(biāo)可以表示為一般形式

(18)

式中：Φ(·)為終端目標(biāo)函數(shù)；L(·)為過程目標(biāo)函數(shù)。

至此，將再入軌跡三維剖面規(guī)劃問題記為P0，表示如下

(19)

在問題P0中，動(dòng)力學(xué)方程和過程約束均是非線性的，因此這是一個(gè)連續(xù)的非線性最優(yōu)控制問題。同時(shí)，由于控制指令攻角并不顯含在方程式(2)中，這也增加了該問題的求解難度。

1.2 新控制量的選擇

再入軌跡三維剖面規(guī)劃問題的控制量是攻角和傾側(cè)角，但是，由于攻角在再入動(dòng)力學(xué)方程式(2) 中是隱含的，攻角通過影響氣動(dòng)系數(shù)從而改變氣動(dòng)力，傾側(cè)角則通過與三角函數(shù)進(jìn)行復(fù)合而顯含于式(2)中，因此，原動(dòng)力學(xué)方程式(2)關(guān)于攻角和傾側(cè)角是強(qiáng)非線性的。由于凸優(yōu)化問題中等式約束必須是線性的，若繼續(xù)使用攻角和傾側(cè)角作為軌跡規(guī)劃的控制量，則顯著增加了對(duì)原方程式(2) 進(jìn)行線性化的難度，并且線性化后的動(dòng)力學(xué)方程也容易引起再入軌跡的“抖振”現(xiàn)象[16]。

為了避免上述問題的出現(xiàn)，可以定義新控制量，從而獲得新的線性化方程。新控制量的選擇應(yīng)滿足以下2個(gè)條件：① 新控制量可以比較方便地將動(dòng)力學(xué)方程中與之相關(guān)的一項(xiàng)表示為關(guān)于它的線性表達(dá)式；② 新控制量與原始控制量之間應(yīng)該存在唯一映射關(guān)系。因此，本文定義如下新的控制變量:

(20)

(21)

這里引入u3的意義在于建立u1和u2之間的關(guān)系，方便后續(xù)進(jìn)行控制量約束的等價(jià)轉(zhuǎn)換。

定義了新的控制量后，式(2)所示的動(dòng)力學(xué)方程可以重新寫為

(22)

令U=[u1,u2,u3]Τ，將式(22)簡(jiǎn)寫為

(23)

式中：

(24)

(25)

(26)

在動(dòng)力學(xué)方程式(23)中，新控制量矢量U關(guān)于矩陣G(x)是線性的；g(x)為關(guān)于狀態(tài)量的非線性項(xiàng)，但其中包含阻力加速度項(xiàng)D,這與原始控制量攻角有關(guān)，關(guān)于該問題對(duì)求解精度的影響將在數(shù)值仿真中進(jìn)行分析；h(x)為與地球自轉(zhuǎn)相關(guān)的小量。

2 凸化與離散

2.1 凸 化

1)控制量約束的處理

定義了新控制量(u1,u2,u3)后，需將如式(11) 和式(12)所示的原始控制量約束轉(zhuǎn)化為對(duì)新控制量的約束。由于控制量u3僅與攻角和馬赫數(shù)有關(guān)，由參考軌跡對(duì)應(yīng)的高度和速度可得到參考馬赫數(shù)Mar?；诠ソ侨≈捣秶鶦α，采用一維優(yōu)化方法線搜索確定控制量u3的取值范圍，表示為

(27)

(28)

式中：Cα=[>αmin,αmax]，則攻角幅值約束可以轉(zhuǎn)化為控制量u3的幅值約束，表示為

(u3)min≤u3≤(u3)max

(29)

傾側(cè)角幅值約束可以轉(zhuǎn)化為對(duì)控制量u1的約束，表示為

u3cos(υmax)≤u1≤u3cos(υmin)

(30)

至此，實(shí)現(xiàn)了對(duì)控制量幅值約束的轉(zhuǎn)化處理，且轉(zhuǎn)化后的約束(29)和(30)均是凸的。

除了上述對(duì)控制量幅值約束進(jìn)行轉(zhuǎn)化外，對(duì)式(21)所示的非凸附加約束也要進(jìn)行凸化處理。這里采用松弛技術(shù)對(duì)其進(jìn)行凸化[16]，得到如下凸約束為

(31)

上述松弛過程的合理性證明將在后面進(jìn)行介紹。

2) 動(dòng)力學(xué)方程的線性化

動(dòng)力學(xué)方程作為等式約束，在凸優(yōu)化問題中必須是線性的。這里介紹如何對(duì)式(23)所示的動(dòng)力學(xué)方程進(jìn)行線性化處理。基于參考軌跡(x(k),α(k),υ(k))，其中上標(biāo)k表示參考軌跡序號(hào)，由一階泰勒展開可得

(32)

將式(32)化簡(jiǎn)為

(33)

其中：

(34)

B=G|x=x(k)

(35)

C=g(x(k))-Ax(k)+h(x(k))

(36)

根據(jù)推導(dǎo)，線性化矩陣A和B的表達(dá)式為

(37)

B(xk)=

(38)

其中：

(39)

(40)

(41)

(42)

為了降低線性化誤差，需添加信賴域約束，表示為

|x-x(k)|≤δ

(43)

式中：δ表示信賴域半徑，與x維數(shù)相同。

3) 過程約束的線性化

在式(13)～式(16)所示的4項(xiàng)過程約束中，對(duì)應(yīng)函數(shù)均具有強(qiáng)非線性，且既不是凸函數(shù)也不是凹函數(shù)，因此仍需要對(duì)這些函數(shù)進(jìn)行線性化。本文采用文獻(xiàn)[16]提出的處理方法，利用一階泰勒展開，將4項(xiàng)過程約束轉(zhuǎn)化為關(guān)于地心距的線性約束，統(tǒng)一記為

(44)

4) 目標(biāo)函數(shù)的選擇

在凸優(yōu)化問題中，目標(biāo)函數(shù)必須是凸的或線性的。為了簡(jiǎn)化處理，這里以終端位置誤差最小為主要指標(biāo)。由于終端位置由經(jīng)度和緯度共同描述，這也是一個(gè)終端性能指標(biāo)，可以表示為

(45)

(46)

式中：C1和C2均為權(quán)重系數(shù)，用于平衡不同項(xiàng)之間的量級(jí)差異。

由于引入的附加項(xiàng)是線性的，因此式(46)所示的目標(biāo)函數(shù)仍然是凸的。

基于上述凸化處理，可以得到具有凸目標(biāo)函數(shù)、凸或線性約束的連續(xù)優(yōu)化問題，記為問題P1，表示為

(47)

由于問題P1中優(yōu)化目標(biāo)是終端位置誤差最小，因此相應(yīng)的終端約束應(yīng)去除對(duì)經(jīng)緯度的限制，將終端約束調(diào)整為

xf=[rf, ～, ～,θf,σf]Τ

(48)

式中：～表示該項(xiàng)無約束。從工程實(shí)現(xiàn)的角度考慮，在終端約束中對(duì)速度傾角和航跡偏航角施加了嚴(yán)格等式約束。

為了證明問題P1中松弛過程的合理性，這里引入如下2個(gè)假設(shè)：

假設(shè)1式(43)所示的信賴域約束|x-x(k)|≤δ在定義域的大部分區(qū)間內(nèi)不會(huì)取到邊界，即除了少數(shù)有限個(gè)連接點(diǎn)外，在[>E0,Ef]上|x-x(k)|<δ是恒成立的。

假設(shè)2式(44)所示的過程約束rmin≤r≤rmax在定義域的大部分區(qū)間內(nèi)不會(huì)取到邊界，即除了少數(shù)有限個(gè)連接點(diǎn)外，在[>E0,Ef]上rmin

對(duì)于假設(shè)1，只要在求解過程中設(shè)置一個(gè)充分大的δ，即可保證在[>E0,Ef]上|x-x(k)|<δ是恒成立的。對(duì)于假設(shè)2，由于求解過程采用了梯度下降算法，當(dāng)?shù)饨咏^程約束邊界時(shí)，控制量將迫使當(dāng)前解遠(yuǎn)離邊界，所以能夠保證在大部分連續(xù)定義域內(nèi)rmin

2.2 離散

問題P1是一個(gè)連續(xù)凸優(yōu)化問題，其目標(biāo)函數(shù)是凸的，約束是凸或線性的，但狀態(tài)量和控制量是連續(xù)的，這為直接求解該問題帶來了困難。因此，將其轉(zhuǎn)化為對(duì)有限參數(shù)進(jìn)行尋優(yōu)的數(shù)學(xué)規(guī)劃問題，應(yīng)用數(shù)值方法進(jìn)行求解是一種可行的方法。這里介紹對(duì)問題P1進(jìn)行離散化，從而便于求解。

將自變量變化域[>E0,Ef]等間距離散為N個(gè)間隔，則自變量離散為E0,E1,E2,…,EN-1,EN，且EN=Ef；狀態(tài)量離散為x0,x1,x2,…,xN-1,xN，且xN=xf；控制量離散為U0,U1,U2,…,UN-1,UN。

采用梯形法則對(duì)問題P1中的線性動(dòng)力學(xué)方程進(jìn)行近似離散，則離散后的動(dòng)力學(xué)方程為

Bj+1Uj+1+Cj+1)

(49)

式中：ΔE=(Ef-E0)/N，下標(biāo)j表示離散點(diǎn)編號(hào)，j=0,1,…,N-1。

整理式(49)，可得到更簡(jiǎn)潔的動(dòng)力學(xué)方程形式，表示為

LjXj-Lj+1Xj+1+MjUj+Mj+1Uj+1=Gj

(50)

完成動(dòng)力學(xué)方程離散化后，再對(duì)問題P1中的其他約束和目標(biāo)函數(shù)積分項(xiàng)進(jìn)行離散化，進(jìn)而得到一個(gè)凸數(shù)學(xué)規(guī)劃問題，記為問題P2，表示為

(51)

問題P2中的優(yōu)化變量包括每個(gè)離散點(diǎn)處的狀態(tài)量xj和控制量Uj。由于N個(gè)離散間隔對(duì)應(yīng)N+1個(gè)離散點(diǎn)，因此優(yōu)化變量數(shù)目為8(N+1)。

3 指令反解與迭代算法

3.1 指令反解

由于問題P2中采用狀態(tài)量x和新控制量U作為優(yōu)化變量，經(jīng)過凸優(yōu)化求解后不能直接獲得用于制導(dǎo)的攻角和傾側(cè)角指令，因此需要利用凸優(yōu)化解進(jìn)行指令反解，獲得每個(gè)離散點(diǎn)處相應(yīng)的攻角和傾側(cè)角指令。

假設(shè)優(yōu)化后第j個(gè)離散點(diǎn)處的狀態(tài)量為xj，新控制量為Uj=[u1,j,u2,j,u3,j]Τ，則對(duì)應(yīng)攻角指令αj可利用牛頓法或反插值方法求解下列非線性方程確定。

F(αj,Maj)=u3,j

(52)

式中：Maj為由xj所確定的當(dāng)前馬赫數(shù)。

傾側(cè)角指令可利用式(53)確定

(53)

在問題P2的求解過程中，控制量附加約束很難保證在每個(gè)離散點(diǎn)處都滿足等式成立。這里定義附加約束違反程度函數(shù)，用于定量描述該約束的違反情況，表示為

(54)

當(dāng)違反程度VU超出給定容許值εU時(shí)，則由式(52)和式(53)所得的反解指令(αj,υj)是不可信的，反之，則是可信的。本文將在仿真中分析這種情況對(duì)解精度的影響。

違反程度容許值的取值可以根據(jù)如下經(jīng)驗(yàn)公式確定:

(55)

式中：Z為經(jīng)驗(yàn)常數(shù)，一般可以取100～500之間。如果Z取值過大，則容許值εU較小，使得較多反解指令是不可信的；反之，若容許值較大，則大部分反解指令是可信的，則失去了判斷反解指令可信度的意義。常數(shù)Ζ的選擇應(yīng)根據(jù)具體問題進(jìn)行權(quán)衡確定。

3.2 序列凸化算法

由于原最優(yōu)控制問題P0具有強(qiáng)非線性，如果僅采用一次線性化來處理非線性動(dòng)力學(xué)與過程約束，即僅求解一次問題P2，則所得凸優(yōu)化解對(duì)于原問題P0基本是不可行的。因此需采用連續(xù)線性化方法，對(duì)動(dòng)力學(xué)與過程約束進(jìn)行連續(xù)線性化，從而產(chǎn)生一個(gè)迭代求解序列。

在這個(gè)序列中，每次迭代都會(huì)求解一個(gè)凸優(yōu)化子問題(即問題P2)，利用該子問題的解再進(jìn)行下一次的線性化，直至這個(gè)解序列最終滿足收斂條件。為了保證該求解序列的收斂性[15]，在信賴域約束中引入松弛系數(shù)κ，并將其作為一個(gè)附加項(xiàng)添加至目標(biāo)函數(shù)中，從而得到新問題P3，表示如下

(56)

式中：C3為松弛系數(shù)項(xiàng)的權(quán)重系數(shù)。

動(dòng)力學(xué)約束是軌跡規(guī)劃問題中必須滿足的一個(gè)約束，問題P3所得的解僅滿足線性動(dòng)力學(xué)約束，相比問題P0中的非線性約束必然存在一定誤差，定量描述這個(gè)誤差可用于表示凸優(yōu)化解相對(duì)原問題P0的可行性。因此，下面定義非線性約束違反程度，用于定量描述凸優(yōu)化解對(duì)原非線性動(dòng)力學(xué)約束的違反情況。

求解問題P3得到的凸優(yōu)化解可表示為(x(k),U(k))，通過指令反解得到的控制指令為(α(k),υ(k))。以該控制指令為輸入，對(duì)式(2)所示的原始動(dòng)力學(xué)方程進(jìn)行積分，得到相對(duì)于凸優(yōu)化解的積分驗(yàn)證解，表示為

(57)

由此可定義非線性約束違反程度，也可稱為凸優(yōu)化解的近似誤差，用于表征優(yōu)化后的狀態(tài)量x(k)與控制指令(α(k),υ(k))對(duì)原始非線性動(dòng)力學(xué)約束的違反程度，計(jì)算公式為

(58)

該迭代算法以相鄰兩次凸優(yōu)化解對(duì)應(yīng)的狀態(tài)量最大偏差作為收斂準(zhǔn)則，表示為

(59)

式中：收斂誤差限定義為ε=[εr,ελ,εφ,εθ,εσ]Τ，上述不等式左邊的表達(dá)式可定義為每個(gè)狀態(tài)量的相鄰2次凸優(yōu)化解最大偏差。

基于上述分析，這里給出基于凸優(yōu)化的序列凸化算法求解原最優(yōu)控制問題P0的步驟：

步驟1給出一條初始參考軌跡x(0)，初始參考指令(α(0),υ(0))，令k=0。

步驟2基于參考軌跡x(k)和參考指令(α(k),υ(k))計(jì)算控制量約束邊界、過程約束邊界以及線性化矩陣。

步驟3求解凸優(yōu)化子問題P3，獲得凸優(yōu)化解(x(k+1),U(k+1))。

步驟4利用凸優(yōu)化解進(jìn)行指令反解，得到新的指令序列(α(k+1),υ(k+1))，以此作為積分輸入，計(jì)算非線性動(dòng)力學(xué)違反程度函數(shù)Rnonlinear。

步驟5判斷是否滿足如式(59)所示的收斂條件，若滿足，則算法結(jié)束；否則，令k=k+1，轉(zhuǎn)入步驟2。

4 數(shù)值仿真

本文采用CAV-H模型進(jìn)行仿真分析[23]，初始化軌跡生成方式可參考文獻(xiàn)[16]，初始攻角指令取為20°常值，初始傾側(cè)角指令取為0°常值。邊界條件設(shè)置如表1所示，主要約束設(shè)置如表2所示。

表1 邊界條件取值Table 1 Values of boundary constraints

表2 其他約束條件取值Table 2 Values of other constraints

仿真中采用MATLAB平臺(tái)下CVX作為建模語言[24]，計(jì)算精度設(shè)為1×10-6，即當(dāng)相鄰2次迭代的目標(biāo)函數(shù)相對(duì)變化量小于預(yù)設(shè)計(jì)算精度時(shí)，求解器迭代終止。默認(rèn)SDPT3求解器求解凸優(yōu)化子問題P3，運(yùn)行計(jì)算機(jī)搭載Intel Core i5-3470 3.20 GHz處理器。

由于整個(gè)滑翔飛行過程中能量E一直為負(fù)值，為了顯示方便，將其轉(zhuǎn)換至[0,1]范圍內(nèi)，公式為

(60)

此時(shí)由再入點(diǎn)開始，能量E′從0到1單調(diào)遞增。

4.1 可行性分析

本節(jié)主要分析所提三維剖面規(guī)劃方法優(yōu)化解的可行性，并驗(yàn)證該方法相比二維剖面規(guī)劃方法在發(fā)揮飛行器固有機(jī)動(dòng)能力方面的優(yōu)勢(shì)。這里的可行主要包括2個(gè)方面，一是優(yōu)化解的狀態(tài)軌跡與控制指令(包括攻角和傾側(cè)角)之間是否滿足非線性動(dòng)力學(xué)約束，二是由控制指令積分獲得的過程量是否滿足非線性約束要求。

圖1～圖5給出了本文仿真條件下利用所提方法得到的再入軌跡。其中，藍(lán)色粗實(shí)線代表的“優(yōu)化解”為迭代過程收斂時(shí)最后一次求解問題P3獲得的解，以紅色點(diǎn)劃線代表的“積分解”由式(57) 定義可得。這里以積分解的結(jié)果作為精確值，驗(yàn)證“優(yōu)化解”的可行性。

圖1中對(duì)軌跡規(guī)劃的起始點(diǎn)進(jìn)行了標(biāo)注，以該點(diǎn)對(duì)應(yīng)的狀態(tài)作為軌跡規(guī)劃的起始狀態(tài)，在此之前的軌跡狀態(tài)按照常值控制量進(jìn)行積分獲得。由圖1～圖3可以看出，優(yōu)化解對(duì)應(yīng)的狀態(tài)量與積分解的結(jié)果基本一致，說明優(yōu)化解中的狀態(tài)量與控制指令之間基本滿足原始非線性動(dòng)力學(xué)約束，沒有出現(xiàn)較大偏差。由圖3可知，優(yōu)化解對(duì)應(yīng)的軌跡滿足期望目標(biāo)點(diǎn)位置約束。

圖4給出了熱流密度、動(dòng)壓和過載3項(xiàng)過程約束的變化規(guī)律。可以看出，優(yōu)化解與積分解對(duì)應(yīng)的軌跡均滿足給定的過程約束，說明相應(yīng)的狀態(tài)量基本一致。同時(shí)，由凸優(yōu)化模型得到的狀態(tài)量仍然滿足原始過程約束，即相對(duì)于原問題是可行的，說明應(yīng)用線性化方法處理非線性過程約束在該問題中是可行的。

圖1 高度和速度變化規(guī)律Fig.1 History of altitude and velocity

圖2 速度傾角和航跡偏航角變化規(guī)律Fig.2 History of flight path angle and heading angle

圖3 地面軌跡對(duì)比示意圖Fig.3 Comparison of footprints

圖5給出了利用優(yōu)化解進(jìn)行指令反解得到的攻角和傾側(cè)角曲線，從控制量的角度講它們即是對(duì)三維剖面的一種表征。在規(guī)劃起點(diǎn)之前，攻角設(shè)置為20°常值，傾側(cè)角設(shè)置為0°常值。圖1～圖4 中積分解均是以圖5所示的控制指令為輸入進(jìn)行積分獲得的。

下面定量分析優(yōu)化解的可行性，主要指優(yōu)化解相對(duì)于原最優(yōu)控制問題的非線性約束違反程度。圖6給出了優(yōu)化解對(duì)應(yīng)的控制量附加約束違反程度V在所有離散點(diǎn)上的分布情況。可以看出，在第153～157個(gè)離散點(diǎn)之間，違反程度異常增大，說明這幾個(gè)離散點(diǎn)處的反解指令是不可信的，而其他離散點(diǎn)處的反解指令是可信的。

圖7給出了優(yōu)化解在每個(gè)離散點(diǎn)處的近似誤差分布情況?？梢悦黠@發(fā)現(xiàn)，在反解指令異常情況出現(xiàn)前后，單個(gè)離散點(diǎn)的近似誤差出現(xiàn)明顯變化，由此也導(dǎo)致后續(xù)離散點(diǎn)處的近似誤差都顯著增大。這個(gè)現(xiàn)象與圖2中，在歸一化能量取值0.8以后，優(yōu)化解與積分解開始出現(xiàn)微小差異是一致的。這些誤差的產(chǎn)生都是由于優(yōu)化解在極個(gè)別離散點(diǎn)處不嚴(yán)格滿足如式(21)所示的控制量附加約束引起的，而這種個(gè)別離散點(diǎn)處的現(xiàn)象一般是由于數(shù)值求解不穩(wěn)定造成的。

圖4 三項(xiàng)過程約束變化規(guī)律Fig.4 History of three path constraints

圖5 控制指令曲線Fig.5 History of control command

表3給出了優(yōu)化解與積分解的終端狀態(tài)對(duì)比結(jié)果?？梢钥闯鼋K端狀態(tài)偏差是比較小的，即解精度在一定范圍內(nèi)是可以保證的。這也說明根據(jù)本文定義的新控制量得到的線性動(dòng)力學(xué)相對(duì)于原始非線性動(dòng)力學(xué)的約束誤差是有限的。所以，在一定精度條件下，由本文方法所得的凸優(yōu)化解可以作為原問題的一個(gè)可行解。

圖6 控制量附加約束違反程度分布Fig.6 Violation distribution of additional constraint

圖7 逐點(diǎn)近似誤差分布Fig.7 Distribution of approximation error on each discretization point

表3 終端狀態(tài)對(duì)比Table 3 Comparison of terminal states

為了比較三維剖面規(guī)劃方法相較于二維剖面規(guī)劃方法的優(yōu)勢(shì)，以文獻(xiàn)[16]中的二維剖面規(guī)劃方法為參考，仿真條件與本文相同，2種方法所得的地面軌跡對(duì)比結(jié)果如圖8所示，終端位置偏差對(duì)比結(jié)果如表4所示?？梢钥闯?，在相同條件下，三維剖面規(guī)劃方法的優(yōu)化解對(duì)應(yīng)的終端位置可以到達(dá)給定目標(biāo)點(diǎn)，距離偏差為0，但二維剖面規(guī)劃方法的優(yōu)化解對(duì)應(yīng)的終端位置與給定目標(biāo)點(diǎn)相距1 678.57 km。由于仿真設(shè)置的優(yōu)化目標(biāo)是終端位置偏差最小，此時(shí)二維剖面規(guī)劃方法的優(yōu)化結(jié)果代表了其所能發(fā)揮的最大機(jī)動(dòng)能力，對(duì)比而言，三維剖面規(guī)劃方法所能發(fā)揮的機(jī)動(dòng)能力更大。由此可以得出，三維剖面規(guī)劃方法在再入滑翔軌跡規(guī)劃時(shí)可以充分發(fā)揮飛行器固有的機(jī)動(dòng)能力。

4.2 收斂性分析

序列凸化算法的收斂性對(duì)于該方法的可行性至關(guān)重要，本節(jié)對(duì)4.1節(jié)算例求解過程的收斂性進(jìn)行分析。

根據(jù)式(59)所示的收斂準(zhǔn)則，當(dāng)序列凸化算法迭代至第15次時(shí)，各個(gè)狀態(tài)量的相鄰兩次凸優(yōu)化解最大偏差分別為

|Δrj|max=78.286 6 m，|Δλj|max=0.044 5°，|Δφj|max=0.020 0°，|Δθj|max=0.045 4°，|Δσj|max=0.198 4°，滿足仿真中所給的收斂誤差限，迭代終止。

為了體現(xiàn)凸優(yōu)化的求解效率，在相同計(jì)算條件和計(jì)算精度下，利用GPOPS 5.0偽譜法工具包求解該問題，則本文所提方法與偽譜法的比較結(jié)果如表5所示，其中，凸優(yōu)化方法的求解時(shí)間由每次子問題求解輸出的CPU時(shí)間累加獲得，偽譜法的求解時(shí)間由每次迭代輸出的實(shí)際求解時(shí)間累加獲得。

由表5可以看出，凸優(yōu)化和偽譜法均可獲得滿足各項(xiàng)約束條件的再入軌跡，實(shí)際計(jì)算精度均小于預(yù)設(shè)計(jì)算精度，2種方法獲得的性能指標(biāo)差異小于10%。相比偽譜法，凸優(yōu)化方法迭代次數(shù)更少，在計(jì)算效率方面具有明顯優(yōu)勢(shì)。

圖9給出了序列凸化迭代過程中得到的歷次軌跡與初始化軌跡的對(duì)比結(jié)果，其中，藍(lán)色實(shí)線表示由線性假設(shè)得到的初始化軌跡，紅色點(diǎn)劃線表示迭代收斂后輸出的優(yōu)化軌跡，其他細(xì)實(shí)線表示中間迭代解對(duì)應(yīng)的軌跡?？梢钥闯?，隨著迭代次數(shù)的增加，相鄰2次軌跡的偏差逐漸減小，優(yōu)化軌跡逐漸收斂。同時(shí)，收斂后的優(yōu)化軌跡與初始化軌跡相差較大，說明本文所提方法的收斂性對(duì)于初始化軌跡的選擇具有一定的適應(yīng)能力。

圖10給出了5個(gè)狀態(tài)量的收斂情況，所有子圖的橫軸取值從1開始。圖10(a)中高度收斂圖的縱坐標(biāo)是對(duì)數(shù)坐標(biāo)。本文中狀態(tài)量的收斂過程主要由相鄰2次凸優(yōu)化解的最大偏差來衡量，由式(59)計(jì)算可得?？偟膩砜?，5個(gè)狀態(tài)量的收斂趨勢(shì)是比較一致的，最大高度偏差在最后2次迭代時(shí)趨于收斂，而最大速度傾角偏差在第10次迭代后趨于收斂，剩余3個(gè)最大狀態(tài)偏差則在第6次迭代后趨于收斂。這說明不同狀態(tài)量的收斂速率是不一致的。

圖11給出了信賴域松弛系數(shù)的收斂情況。可以看出，除了前2次迭代時(shí)得到的松弛系數(shù)比較大以外，后續(xù)迭代過程中松弛系數(shù)始終比較小，且最后達(dá)到收斂條件時(shí)松弛系數(shù)取值為κ*=0.002 7，說明隨著信賴域松弛系數(shù)的減小，序列迭代算法逐漸趨近于收斂。

圖10 5個(gè)狀態(tài)量收斂情況Fig.10 Convergence process of five state variables

圖12給出了最優(yōu)性能指標(biāo)的收斂情況，該圖中縱坐標(biāo)為對(duì)數(shù)坐標(biāo)?？梢钥闯觯瑑?yōu)化目標(biāo)的收斂是比較快的，基本在第6次迭代以后就趨于穩(wěn)定。當(dāng)序列凸化算法收斂時(shí)，目標(biāo)函數(shù)中的終端位置偏差和信賴域松弛變量都將是一個(gè)小量，此時(shí)影響目標(biāo)函數(shù)取值的主要是航跡偏航角的積分項(xiàng)，而該項(xiàng)的變化僅與側(cè)向運(yùn)動(dòng)相關(guān)。所以，經(jīng)度、緯度和航跡偏航角這3個(gè)與側(cè)向運(yùn)動(dòng)相關(guān)的狀態(tài)量，其收斂過程與優(yōu)化目標(biāo)基本一致。由于目標(biāo)函數(shù)中缺少縱向運(yùn)動(dòng)參數(shù)，所以與縱向運(yùn)動(dòng)相關(guān)的高度和速度傾角這2個(gè)狀態(tài)量在迭代過程中出現(xiàn)了較多震蕩，一定程度上降低了算法的整體收斂速率。

以上分析從數(shù)值上驗(yàn)證了本文所提方法的收斂性，且收斂后的優(yōu)化解僅可能是原問題的一個(gè)局部最優(yōu)解，關(guān)于這種方法最優(yōu)性的理論證明目前還并不充分[17]。

圖12 最優(yōu)性能指標(biāo)收斂情況Fig.12 Convergence process of optimal objective function

由于控制量附加約束的違反程度對(duì)于優(yōu)化解精度有一定影響，因此在目標(biāo)函數(shù)中引入關(guān)于航跡偏航角積分的附加項(xiàng)，如式(46)所示，下面分析該附加項(xiàng)的引入對(duì)解精度的影響。圖13給出了迭代過程中，每次迭代的解中，嚴(yán)格滿足附加約束的離散點(diǎn)占總離散點(diǎn)數(shù)目的比例變化情況?？梢钥闯?，引入附加項(xiàng)后，除了第8次和第11次的比例很低外，其他每次迭代解對(duì)應(yīng)的比例都比較高，迭代后期這個(gè)比例值基本趨于穩(wěn)定，而在不引入附加項(xiàng)的迭代過程中，滿足附加約束的離散點(diǎn)占比在后期一直較低，由此說明該附加項(xiàng)的引入可以有效確保附加控制量約束松弛過程的有效性。

圖14給出了迭代過程中近似誤差的收斂情況?？梢钥闯觯敫郊禹?xiàng)后，由式(58)所定義的近似誤差在迭代過程中的收斂性是比較好的，在迭代前期很快能夠減小到一個(gè)很低的水平，說明隨著迭代算法逐步趨于收斂，凸優(yōu)化解對(duì)原始非線性約束的違反程度逐漸降低，即解的精度逐步提高。但在不引入附加項(xiàng)的迭代過程中，由于大部分離散點(diǎn)處的優(yōu)化值均不滿足附加約束，所以導(dǎo)致解精度較低。由此說明附加項(xiàng)的引入對(duì)于提高優(yōu)化解精度是有顯著意義的。

圖13 附加項(xiàng)對(duì)離散點(diǎn)占比的影響Fig.13 Effect of additional term on proportions of discretization points

圖14 附加項(xiàng)對(duì)凸優(yōu)化解的近似誤差影響Fig.14 Effect of additional term on approximation errors of solution by convex optimization

序列凸化算法的啟動(dòng)必須依賴于一條初始軌跡。本文采用線性假設(shè)進(jìn)行初始軌跡的生成，即認(rèn)為運(yùn)動(dòng)狀態(tài)的滑翔段均是線性變化；作為對(duì)比的是以表5中提到的偽譜法所得優(yōu)化軌跡作為初始軌跡。下面分析不同初始軌跡對(duì)本文所提算法收斂性的影響。圖15～圖17分別給出了信賴域松弛系數(shù)、最優(yōu)性能指標(biāo)和凸優(yōu)化解近似誤差的收斂對(duì)比情況。

由圖15～圖17可以看出，采用偽譜法得到>的優(yōu)化軌跡作為初始軌跡后，序列凸化算法在迭代12次即可達(dá)到收斂條件，而同樣條件下線性假設(shè)得到的初始軌跡需要15次迭代。由于收斂條件的參數(shù)不變，所以2組優(yōu)化解的計(jì)算精度和近似誤差也基本一致。在相同計(jì)算精度條件下，由于迭代次數(shù)的減少，以偽譜法優(yōu)化軌跡作為初始軌跡的迭代過程耗時(shí)也更短，由于每次迭代的CPU時(shí)間基本一致，所以節(jié)約時(shí)間比例與迭代次數(shù)的減少比例也是基本一致的。

圖15 初始軌跡對(duì)信賴域松弛系數(shù)的影響Fig.15 Effect of initial trajectory on trust region relaxation coefficient

圖16 初始軌跡對(duì)最優(yōu)性能指標(biāo)的影響Fig.16 Effect of initial trajectory on performance index

圖17 初始軌跡對(duì)凸優(yōu)化解的近似誤差影響Fig.17 Effect of initial trajectory on approximation errors of solution by convex optimization

雖然采用偽譜法優(yōu)化軌跡作為初始軌跡可以減少序列凸化算法的求解時(shí)間，但如果后者依賴于偽譜法提供初始軌跡，則凸優(yōu)化的計(jì)算時(shí)間優(yōu)勢(shì)將無法體現(xiàn)。因此，有必要尋找更快速合適的初始軌跡生成方法，從而進(jìn)一步提高序列凸化算法的求解效率。

5 結(jié) 論

1) 建立了再入軌跡三維剖面凸優(yōu)化模型，相比二維剖面可充分發(fā)揮飛行器固有的機(jī)動(dòng)能力；相比偽譜法，凸優(yōu)化在軌跡規(guī)劃問題上具有更高的求解效率，具備實(shí)現(xiàn)軌跡在線規(guī)劃的潛力。

2) 為了實(shí)現(xiàn)問題凸化，引入了新的控制量和附加約束。迭代求解中，在大部分離散點(diǎn)處該附加約束都是嚴(yán)格滿足的，個(gè)別點(diǎn)處的不滿足對(duì)收斂后優(yōu)化解的精度影響是有限的。

3) 序列凸化方法在數(shù)值求解中具有穩(wěn)定的收斂性，隨著迭代過程趨于收斂，凸優(yōu)化解的精度逐漸提高。由于目標(biāo)函數(shù)中主要包含側(cè)向運(yùn)動(dòng)項(xiàng)，所以側(cè)向運(yùn)動(dòng)狀態(tài)量的收斂速率快于縱向運(yùn)動(dòng)狀態(tài)量。