高超聲速飛行器再入軌跡快速優(yōu)化研究

董春云，郭 志，趙培博，蔡遠(yuǎn)利，于振華

(1.西安交通大學(xué) 電子與信息工程學(xué)院,西安 710049；2.空軍工程大學(xué) 信息與導(dǎo)航學(xué)院,西安 710077)

?

高超聲速飛行器再入軌跡快速優(yōu)化研究

董春云1，郭 志1，趙培博1，蔡遠(yuǎn)利1，于振華2

(1.西安交通大學(xué) 電子與信息工程學(xué)院,西安 710049；2.空軍工程大學(xué) 信息與導(dǎo)航學(xué)院,西安 710077)

基于求解最優(yōu)控制問題的Chebyshev偽譜法(Chebyshev Pseudospectral Method, CPM)，研究了高超聲速飛行器再入軌跡快速優(yōu)化問題。針對遠(yuǎn)程多約束條件下再入軌跡優(yōu)化問題的難點(diǎn)，提出了一種線性初值與節(jié)點(diǎn)更新相結(jié)合的優(yōu)化策略，將攻角與傾側(cè)角同時(shí)作為控制變量，以再入飛行時(shí)間最短為優(yōu)化目標(biāo)，利用CPM將軌跡優(yōu)化問題轉(zhuǎn)化為非線性規(guī)劃問題，并使用SNOPT軟件包求解，使CPM成為一種再入軌跡快速優(yōu)化的通用算法。以某類高超聲速再入飛行器為對象進(jìn)行軌跡優(yōu)化計(jì)算，并對比相同仿真條件下粒子群(PSO)算法的優(yōu)化效果，仿真結(jié)果驗(yàn)證了該算法具有較高的求解效率和快速收斂性。

高超聲速飛行器；軌跡優(yōu)化；再入；Chebyshev偽譜方法

0 引言

近年來，以通用航空飛行器(Common Aero Vehicle，CAV)[1]為代表的一類高超聲速再入飛行器受到各軍事強(qiáng)國的高度重視。該類飛行器自身不攜帶發(fā)動(dòng)機(jī)，且具有較大的升阻比，從軌道或亞軌道高度再入，依靠氣動(dòng)力控制實(shí)現(xiàn)遠(yuǎn)距離滑翔再入飛行，突破了常規(guī)彈道式再入模式，具有機(jī)動(dòng)能力強(qiáng)、彈道靈活多變難以攔截等諸多優(yōu)點(diǎn)，在軍事上具有廣泛的應(yīng)用前景[1]。

軌跡優(yōu)化設(shè)計(jì)是高超聲速再入飛行器的關(guān)鍵技術(shù)之一，也是目前研究的熱點(diǎn)問題。高超聲速飛行器再入過程具有飛行速度快、空間跨度大、氣動(dòng)熱力環(huán)境惡劣、再入軌跡對控制變量高度敏感等特點(diǎn)，為了滿足安全飛行的需要，飛行器需滿足嚴(yán)格的過載、動(dòng)壓、駐點(diǎn)熱流密度等諸多非線性約束。同時(shí)，考慮到與末制導(dǎo)段交班的要求，再入終端狀態(tài)也需要滿足一定的約束。諸多約束條件使得再入軌跡的可行域被限制在較為狹窄的范圍內(nèi)，給軌跡優(yōu)化設(shè)計(jì)帶來了一定的困難。此外，由于升力式飛行器升阻比大，飛行時(shí)間過長極易造成再入過程的總吸熱量過高，給飛行器的熱防護(hù)系統(tǒng)設(shè)計(jì)帶來困難，也不利于提高有效載荷的質(zhì)量。因此，通過合適的優(yōu)化策略和方法得到滿足以上約束條件，且飛行時(shí)間最短的再入軌跡是非常有意義的。

上述再入軌跡優(yōu)化問題是一類典型的具有路徑約束和終端約束的最優(yōu)控制問題，其數(shù)值求解方法主要分為間接法和直接法兩類。相對于間接法，直接法在解決實(shí)際問題的適應(yīng)性和收斂的魯棒性上更具優(yōu)勢。依據(jù)參數(shù)化方法的不同，直接法又分為僅離散控制變量的打靶法和同時(shí)離散狀態(tài)變量和控制變量的配點(diǎn)法。由于再入軌跡對控制變量高度敏感，僅離散控制變量的方法容易陷入局部解，甚至收斂不到可行解，而配點(diǎn)法以其求解精度較高、收斂性好、易實(shí)現(xiàn)等優(yōu)點(diǎn)，在軌跡優(yōu)化問題的求解中得到了更多的應(yīng)用。近年來，配點(diǎn)法中的偽譜法以其高精度和高效率等優(yōu)勢逐漸受到重視，其收斂性得到了理論證明[3-4]。根據(jù)離散節(jié)點(diǎn)和插值基函數(shù)選取的不同，偽譜法可分為Gauss偽譜法、Legendre偽譜法、Radau偽譜法和Chebyshev偽譜法(Chebyshev Pseudospectral Method, CPM)等。

Fahroo和Ross于2002年提出了一種用于求解最優(yōu)控制問題的Chebyeshev偽譜法[5]，但當(dāng)時(shí)并未引起廣泛關(guān)注。直到2008年，Trefethen證明了Clensshaw-Curtis數(shù)值求積公式與Gauss求積公式精度相當(dāng)[6]；Gong等于2009年～2010年進(jìn)一步完善了Chebyeshev偽譜法的理論體系，證明了其離散數(shù)值解會一致地趨近于原問題的最優(yōu)解[7]，并建立了協(xié)狀態(tài)估計(jì)理論[8]。此后，國內(nèi)外學(xué)者開始重新關(guān)注Chebyeshev偽譜法。

與Gauss偽譜法相比，Chebyeshev偽譜法的離散節(jié)點(diǎn)Chebyshev-Gauss-Lobatto(CGL)具有顯式表達(dá)式[9]，使用起來非常方便。而Gauss偽譜法只能通過大量的代數(shù)求根運(yùn)算來獲得節(jié)點(diǎn)的位置，且不能直接得到端點(diǎn)處的控制量，終端狀態(tài)必須通過積分獲得，這在一定程度上加劇了計(jì)算難度，同時(shí)也增加了計(jì)算耗時(shí)。近年來，在最短奔跑時(shí)間求解[10]、月球軟著陸軌跡優(yōu)化[11]、橋式吊車最優(yōu)控制[12]和制導(dǎo)炮彈彈道優(yōu)化[13]等領(lǐng)域，Chebyeshev偽譜法都以其高效的求解效率和良好的優(yōu)化性能得到了成功應(yīng)用。而對于高動(dòng)態(tài)、強(qiáng)約束和非線性的高超聲速飛行器再入過程軌跡優(yōu)化問題，Chebyeshev偽譜法是否適用，文中進(jìn)行了詳細(xì)的討論。

文中研究了基于Chebyshev偽譜法的高超聲速飛行器再入軌跡快速優(yōu)化問題。首先，給出了再入飛行器歸一化動(dòng)力學(xué)模型和相關(guān)約束條件，建立了以最短再入時(shí)間為優(yōu)化目標(biāo)的多約束軌跡優(yōu)化模型；在此基礎(chǔ)上，詳細(xì)闡述了利用Chebyshev偽譜法對該軌跡優(yōu)化問題進(jìn)行離散化求解的具體步驟，針對設(shè)計(jì)變量初值給定的難點(diǎn)，給出一種兼顧求解精度與實(shí)時(shí)性的串行求解策略；最后，通過仿真對本文算法的有效性與求解效率進(jìn)行了驗(yàn)證分析，并與基于另一種直接法-粒子群算法的軌跡優(yōu)化結(jié)果進(jìn)行了比較。

1 軌跡優(yōu)化模型

1.1 動(dòng)力學(xué)模型及歸一化處理

假設(shè)地球?yàn)樾D(zhuǎn)圓球，高超聲速飛行器無動(dòng)力再入的非線性動(dòng)力學(xué)模型可用6個(gè)狀態(tài)變量來描述。由于變量間存在不可公度性，取值范圍和量綱不盡相同，為了提高后續(xù)優(yōu)化計(jì)算效率與求解精度，需要對模型進(jìn)行歸一化處理。極坐標(biāo)系下無量綱運(yùn)動(dòng)方程[14]為

(1)

(2)

(3)

(sinγcosφ-cosγsinφcosΨ)

(4)

Ω2rcosφ(cosγcosφ+sinγsinφcosΨ)]

(5)

2ΩV(tanγcosφcosΨ-sinφ)+

(6)

式中m、S分別為再入飛行器質(zhì)量和氣動(dòng)參考面積；ρ為大氣密度；CL、CD分別為升力和阻力系數(shù)。

氣動(dòng)力的計(jì)算采用1976年美國標(biāo)準(zhǔn)大氣模型USSA76[15]，在0～100 km高度范圍內(nèi)大氣劃分為8層，根據(jù)飛行高度計(jì)算空氣密度ρ。CL、CD通過氣動(dòng)力數(shù)據(jù)擬合得到[16]。

1.2 控制變量

在高超聲速飛行條件下，氣動(dòng)力系數(shù)可近似表示為攻角的函數(shù)，同時(shí)將攻角α和傾側(cè)角σ作為優(yōu)化控制變量[17]，能夠充分體現(xiàn)出攻角在再入飛行中的控制調(diào)節(jié)作用。因此，控制變量為

u=(ασ)T

1.3 約束條件

(1)過程約束

飛行器再入飛行是一個(gè)高動(dòng)態(tài)過程，要求嚴(yán)格滿足過程約束，以保證飛行器在結(jié)構(gòu)和熱防護(hù)上的可靠性，主要包括熱流密度、動(dòng)壓、過載和控制變量約束，即

(7)

(8)

(9)

(10)

式中R為飛行器頭部曲率半徑；ρs為海平面大氣密度；C1為常數(shù)。

(2)終端約束

飛行器再入終端狀態(tài)對末制導(dǎo)段的飛行有重要的影響，終端約束與再入飛行的任務(wù)相關(guān)，一般包括速度、飛行高度和地理位置約束等：

|V(τf)-Vf|≤εv，|r(τf)-rf|≤εr

|θ(τf)-θf|≤εθ，|φ(τf)-φf|≤εφ

(11)

其中，εV，εr，εθ，εφ分別是終端狀態(tài)的誤差界，考慮精確制導(dǎo)時(shí)，可設(shè)置εV=0，εr=0，εθ=0，εφ=0。

1.4 性能指標(biāo)

優(yōu)化目標(biāo)按照設(shè)計(jì)要求可以有不同的形式。由于再入環(huán)境惡劣，最短的飛行時(shí)間能夠減小氣動(dòng)加熱，降低能量損耗，縮短敵方響應(yīng)時(shí)間。因此，以最短飛行時(shí)間為性能指標(biāo)：

J=min(τf)

(12)

1.5 優(yōu)化模型

高超聲速飛行器無動(dòng)力再入軌跡優(yōu)化問題可具體描述為：在[τ0,τf]時(shí)間內(nèi)(τf未知)，確定控制變量u=(ασ)T和終端時(shí)刻τf，使得目標(biāo)函數(shù)式(12)最小，并滿足系統(tǒng)狀態(tài)方程約束式(1)～式(6)、過程約束式(7)～式(10)以及終端約束式(11)。

將Chebyshev偽譜法用于該軌跡優(yōu)化問題的求解，力求使其成為一種再入軌跡快速優(yōu)化的通用算法。

2 Chebyshev偽譜法與軌跡優(yōu)化

2.1 方法描述

對于上述連續(xù)時(shí)間最優(yōu)控制問題，Chebyshev偽譜法的基本求解思路為：將連續(xù)時(shí)間狀態(tài)變量和控制變量在一系列CGL點(diǎn)上離散，并以這些離散點(diǎn)為節(jié)點(diǎn)構(gòu)造Largrange插值多項(xiàng)式來逼近真實(shí)狀態(tài)和控制；通過對全局插值多項(xiàng)式求導(dǎo)來近似狀態(tài)變量對時(shí)間的導(dǎo)數(shù)，將微分方程約束轉(zhuǎn)換為代數(shù)約束；性能指標(biāo)中的積分項(xiàng)由Clenshaw-Curtis數(shù)值積分計(jì)算。通過上述方法，可將最優(yōu)控制問題轉(zhuǎn)化為具有一系列代數(shù)約束的NLP問題[5,8]。

CPM的具體計(jì)算步驟簡述如下：

(1)時(shí)間區(qū)間轉(zhuǎn)換

CPM的CGL離散點(diǎn)定義在[-1,1]區(qū)間上，需要將原最優(yōu)控制問題的時(shí)間域τ從[τ0,τf]線性轉(zhuǎn)換到t∈[-1,1]區(qū)間，轉(zhuǎn)換式為

τ=[(τf-τ0)t+(τf+τ0)]/2

(2)離散節(jié)點(diǎn)計(jì)算

CPM離散節(jié)點(diǎn)選取為N階Chebyshev多項(xiàng)式TN(t)=cos(Ncos-1t)的極值點(diǎn)，即CGL點(diǎn)，它們不均勻的分布在[-1,1]區(qū)間上，顯式計(jì)算式為

顯然，離散節(jié)點(diǎn)tk滿足t0=-1，tN=1。

(3)狀態(tài)變量和控制變量插值近似

取上述N+1個(gè)離散點(diǎn)處的狀態(tài)變量和控制變量，分別構(gòu)造Lagrange插值多項(xiàng)式作為連續(xù)狀態(tài)和控制的近似。真實(shí)狀態(tài)變量x(t)與控制變量u(t)的近似表達(dá)式為

(13)

(14)

其中，Lagrange插值基函數(shù)：

由Lagrange插值的性質(zhì)可知，離散節(jié)點(diǎn)處的狀態(tài)近似值與實(shí)際狀態(tài)相等，控制近似值與實(shí)際控制相等。

(4)動(dòng)態(tài)約束處理

對式(13)求導(dǎo)，得到狀態(tài)向量在tk點(diǎn)處的導(dǎo)數(shù)近似表達(dá)式為

(15)

式中Dkj為(N+1)×(N+1)微分矩陣D的第k行第j列元素；D的計(jì)算公式見文獻(xiàn)[5]。

(16)

式中f為式(1)～式(6)所示的狀態(tài)方程。

對于式(7)～式(10)所示的過程約束，需在離散節(jié)點(diǎn)處嚴(yán)格滿足。

(5)性能指標(biāo)積分近似

對于優(yōu)化性能指標(biāo)中存在積分項(xiàng)的情況，利用Clenshaw-Curtis數(shù)值積分對其進(jìn)行近似。對于 [-1,1]區(qū)間上的任一連續(xù)函數(shù)p(t)，其積分可用N+1個(gè)CGL離散點(diǎn)處的函數(shù)值累加和近似，即

(17)

式中，ωk(k=0,1,…,N)為Clenshaw-Curtis加權(quán)，其計(jì)算式參見文獻(xiàn)[5]。

通過以上對最優(yōu)控制問題在CGL節(jié)點(diǎn)的近似，將原連續(xù)時(shí)間最優(yōu)控制問題轉(zhuǎn)化為以下離散化的非線性規(guī)劃問題，即確定CGL節(jié)點(diǎn)處的離散狀態(tài)變量X、離散控制變量U和終端時(shí)間τf，使得性能指標(biāo)J最小，并滿足離散化后的狀態(tài)方程約束、過程約束與終端約束等。

2.2 軌跡優(yōu)化求解策略

利用SNOPT軟件包[18]求解經(jīng)Chebyshev偽譜法轉(zhuǎn)化獲得的非線性規(guī)劃問題，該軟件包基于序列二次規(guī)劃算法，在處理大規(guī)模NLP的求解上具有良好的性能，目前已得到廣泛應(yīng)用。

根據(jù)數(shù)值驗(yàn)證結(jié)果，CGL節(jié)點(diǎn)的數(shù)量越大，求解精度越高，但計(jì)算量與計(jì)算耗時(shí)也會增加。另一方面，高超聲速飛行器再入軌跡優(yōu)化模型復(fù)雜，考慮的約束條件較多，根據(jù)本文模型，當(dāng)CGL離散節(jié)點(diǎn)數(shù)為N+1時(shí)，設(shè)計(jì)變量數(shù)目為6(N+1)+2(N+1)+1，約束總數(shù)為6(N+1)+4(N+1)+4，N的數(shù)值較大時(shí)，設(shè)計(jì)變量初值的設(shè)置會非常繁雜，不恰當(dāng)?shù)某踔禃黾佑?jì)算復(fù)雜度與程序耗時(shí)，甚至使問題無法收斂到可行解。因此，如何確定CGL節(jié)點(diǎn)的數(shù)目與設(shè)計(jì)變量的初值，是優(yōu)化過程中首先需要解決的問題。為此，提出了一種線性初值與節(jié)點(diǎn)更新相結(jié)合的優(yōu)化策略。

如圖1所示，首先采用線性化方法獲得設(shè)計(jì)變量初值，求解較少CGL節(jié)點(diǎn)下的初步優(yōu)化數(shù)值解；然后，增加CGL節(jié)點(diǎn)，利用初步優(yōu)化數(shù)值解插值獲得需要的設(shè)計(jì)變量初值，從而提高求解效率。

所謂線性化方法即根據(jù)設(shè)計(jì)變量初始位置x0與終端約束xf構(gòu)造簡單線性函數(shù)，具體可表述為

節(jié)點(diǎn)更新過程中，通過循環(huán)依次增加節(jié)點(diǎn)的數(shù)目，計(jì)算優(yōu)化結(jié)果的數(shù)值積分終端約束偏差，當(dāng)最大相對誤差小δ時(shí)，認(rèn)為精度符合要求，停止計(jì)算。圖1中，變量nx、ny可根據(jù)仿真效果靈活調(diào)整。

2.3 CPM數(shù)值驗(yàn)證

在應(yīng)用CPM求解再入軌跡優(yōu)化問題之前，首先選擇標(biāo)準(zhǔn)最速降線問題對本文算法的有效性進(jìn)行驗(yàn)證，并分析其求解精度。

圖1 軌跡優(yōu)化求解策略Fig.1 Combination strategy for trajectory optimization

最速降線問題即尋找一條最優(yōu)曲線，使得質(zhì)點(diǎn)在重力作用下，沿該曲線從給定點(diǎn)到不在它垂直下方的另一點(diǎn)的時(shí)間最短。其數(shù)學(xué)模型可以概括為

式中，控制量θ是曲線斜率隨時(shí)間變化的函數(shù)。

該問題的解析形式為

x(τ)=(gτf/π){τ-(τf/π)sin[π(1-τ/τf)]}

將CPM得到的控制變量插值后代入原運(yùn)動(dòng)方程，采用四階龍格庫塔法積分求解運(yùn)動(dòng)軌跡，驗(yàn)證方法可行性。CGL節(jié)點(diǎn)數(shù)目K=N+1，設(shè)N=10時(shí)對比CPM數(shù)值解、積分求解結(jié)果及原問題解析解曲線圖，如圖2與圖3所示。結(jié)果表明，三者的位置坐標(biāo)x、y軌跡均能較好吻合，CPM和SNOPT軟件包的計(jì)算時(shí)間短，求解精度高。

圖2 最速降線控制變量對比圖Fig.2 Time history of control for the brachistochrone problem

為了分析CGL節(jié)點(diǎn)數(shù)目對CPM求解精度與收斂速度的影響，分別取N=5,10,20，對比CPM數(shù)值解與原問題解析解在離散節(jié)點(diǎn)處的求解誤差以及程序仿真運(yùn)行時(shí)間，結(jié)果如表1所示。

從表1看出，3種情況下CPM獲得的數(shù)值解精度均較高，且隨著離散點(diǎn)數(shù)的增加，求解精度越來越高，一致趨向于原問題最優(yōu)解，與文獻(xiàn)[7]結(jié)論一致。需注意的是，仿真程序耗時(shí)隨著離散點(diǎn)數(shù)增加而增加。

圖3 最速降線狀態(tài)變量對比圖Fig.3 Time history of states for the brachistochrone problem

表1 CPM求解精度與仿真時(shí)間比較Table 1 Comparison of solution accuracy and computation cost of CPM

3 仿真校驗(yàn)與結(jié)果分析

以遠(yuǎn)程高超聲速滑翔式再入飛行為例，研究Chebyshev偽譜法在再入軌跡快速優(yōu)化問題中的應(yīng)用。相同仿真條件下，將本文算法的優(yōu)化結(jié)果與粒子群(Particle Swarm Optimization, PSO)算法的優(yōu)化結(jié)果進(jìn)行對比分析。

PSO算法的詳細(xì)介紹與具體步驟可參考文獻(xiàn)[19-20]，在此不展開討論。

3.1 參數(shù)設(shè)定

飛行器結(jié)構(gòu)參數(shù)和氣動(dòng)參數(shù)參考1998年美國洛克希勒-馬丁公司設(shè)計(jì)的高超聲速飛行器CAV-H[16]，該氣動(dòng)參考面積為0.483 9 m2，質(zhì)量為907.2 kg，長度為2.717 8～3.627 6 m。

3.2 結(jié)果及分析

仿真平臺是2.33G主頻CPU、2G內(nèi)存的PC機(jī)，軟件環(huán)境是Matlab 2009。初始節(jié)點(diǎn)個(gè)數(shù)取為K=5，采用線性化方法設(shè)置初值。節(jié)點(diǎn)更新過程中，將獲得的控制變量帶入動(dòng)力學(xué)方程，利用四階龍格庫塔法進(jìn)行數(shù)值積分，積分步長取為0.05，積分軌跡終端約束最大相對誤差閾值取為0.05。

初始優(yōu)化結(jié)果獲得的積分軌跡中，終端最大相對誤差0.496。更新CGL節(jié)點(diǎn)數(shù)目，重新優(yōu)化并積分，當(dāng)節(jié)點(diǎn)個(gè)數(shù)K=15時(shí)，終端最大相對誤差為0.117。繼續(xù)更新節(jié)點(diǎn)數(shù)目，當(dāng)K=28時(shí)，獲得終端高度20.54 km，速度1 018.21 m/s，經(jīng)度236.36°，緯度36.95°，最大相對誤差為0.027，小于最大誤差閾值，停止更新。

至此，CPM優(yōu)化獲得的再入飛行時(shí)間為1 486.46 s，初始優(yōu)化耗時(shí)約1～2 s，末次優(yōu)化耗時(shí)約為75～80 s。相同仿真條件下，PSO優(yōu)化程序耗時(shí)約16 min，優(yōu)化得到的再入飛行時(shí)間為1 485.11 s。Chebyshev偽譜法較一般直接法(如本文PSO算法)計(jì)算效率更高，另一方面，作者在實(shí)際應(yīng)用過程中發(fā)現(xiàn)，良好的初值條件能夠有效改善優(yōu)化效率，減少程序耗時(shí)。

飛行器的主要軌跡參數(shù)，即高度、速度及航跡的優(yōu)化結(jié)果如圖4～圖6所示。圖中，實(shí)線為Chebyeshev偽譜法優(yōu)化結(jié)果，雙劃線為其數(shù)值積分結(jié)果，虛線表示粒子群算法。可以發(fā)現(xiàn)，CPM及PSO 2種方法獲得的飛行器高度、速度及航跡曲線變化趨勢較為一致，均能較好地收斂到指定終端狀態(tài)，且優(yōu)化獲得的再入飛行時(shí)間非常接近。從程序耗時(shí)角度考慮，本文算法運(yùn)行速度更快，能夠滿足再入軌跡快速優(yōu)化的要求。

圖4 飛行高度對比Fig.4 Comparison of altitude profile

圖5 飛行速度對比Fig.5 Comparison of velocity profile

圖6 地面航跡對比Fig.6 Comparison of ground track

從圖4～圖6的積分曲線可看出，CPM優(yōu)化結(jié)果與數(shù)值積分結(jié)果基本一致，具有較高的精度，說明本文給出的節(jié)點(diǎn)更新策略有效保證了算法優(yōu)化結(jié)果的可行性與有效性。其中，飛行高度數(shù)值積分結(jié)果在約500～1 000 s的中間階段略有偏差，而在初始與末尾階段均保持了較高的積分精度。這是因?yàn)镃GL節(jié)點(diǎn)分布不等距，兩端較密集中間較稀疏，且實(shí)際飛行中彈道中段高度變化劇烈。此外，由高度軌跡可看出，為了增大射程，從而減小飛行時(shí)間，飛行器進(jìn)行了若干次跳躍滑翔，開始跳躍幅度較大，但在滑翔末段，軌跡趨于平緩，為與末制導(dǎo)的交班提供了良好的條件。

圖7給出了控制變量優(yōu)化結(jié)果的對比曲線，包括攻角與傾側(cè)角曲線。從圖7可看出，Chebyshev偽譜法與粒子群算法獲得的控制變量變化規(guī)律有一定的差別。這是由于高超聲速飛行器再入環(huán)境惡劣、再入過程復(fù)雜，不同優(yōu)化方法得到的飛行軌跡并不一定是嚴(yán)格意義上的最優(yōu)，而有可能是“次優(yōu)控制”[21]，這也導(dǎo)致了2種方法在飛行軌跡上存在差異。

Chebyshev偽譜法獲得的飛行器再入熱流密度、動(dòng)壓和過載嚴(yán)格滿足過程約束且數(shù)值較小，航跡角與航向角曲線變化均較為平緩，滿足設(shè)計(jì)指標(biāo)要求，結(jié)果如圖8～圖12所示。從而可看出，再入初期，飛行器速度較高，熱流密度約束起主要作用；在隨后的滑翔階段，隨著再入高度的降低以及大氣密度的增加，動(dòng)壓和過載逐漸起主要作用。

圖7 控制變量對比Fig.7 Comparison of control profile

圖8 CPM熱流密度變化圖Fig.8 Heating rate profile

圖9 CPM動(dòng)壓變化圖Fig.9 Dynamic pressure profile

圖10 CPM過載變化圖Fig.10 Overload profile

圖11 CPM航跡角變化圖Fig.11 Flight path angle profile

圖12 CPM航向角變化圖Fig.12 Heading angle profile

4 結(jié)論

將Chebyshev偽譜法用于高超聲速無動(dòng)力再入飛行器的軌跡優(yōu)化領(lǐng)域，提出了一種線性初值與節(jié)點(diǎn)更新相結(jié)合的優(yōu)化策略。仿真結(jié)果表明，與一般直接法比較，該方法能夠利用較少的離散節(jié)點(diǎn)，在短時(shí)間內(nèi)生成滿足約束條件的再入飛行軌跡，求解效率高，收斂速度快，且步驟簡單，使用方便，是一種有效的高超聲速飛行器再入軌跡快速優(yōu)化算法。

[1] George Richie.The common aero vehicle:space delivery system of the future[C] //AIAA Space Technology Conference and Exposition.Albuquerque,New Mexico,1999.

[2] 謝愈,劉魯華,湯國建,等.多約束條件下高超聲速滑翔飛行器軌跡優(yōu)化[J].宇航學(xué)報(bào),2011,32(12):2499- 2504.

[3] Rao A V.A survey of numerical methods for optimal control[J].Advances in the Astronautical Sciences,2009,135(1):497-528.

[4] Begum Senses,Rao A V.A preliminary analysis of small spacecraft finite-thrust aeroassisted orbital transfer[C] //AIAA/AAS Astrodynamics Specialist Conference.Minneapolis,Minnesota,2012.

[5] Fahroo Fariba,Ross I Michael.Direct trajectory optimization by a Chebyshev pseudospectral method[J].Journal of Guidance,Control,and Dynamics,2002,25(1):160-166.

[6] Trefethen Lloyd N.Is Gauss quadrature better than Clenshaw-Curtis[J].SIAM review,2008,50(1):67-87.

[7] Qi Gong,Ross I Michael,Fahroo Fariba.A Chebyshev pseudospectral method for nonlinear constrained optimal control problems[C]//The 48th IEEE Conference on Decision and Control.Shanghai,China:IEEE,2009:16-18.

[8] Gong Qi,Michael Ross I,Fahroo Fariba.Costate computation by a Chebyshev pseudospectral method[J].Journal of Guidance,Control,and Dynamics,2010,33(2):623-628.

[9] Vlassenbroeck J,Van Dooren R.A Chebyshev technique for solving nonlinear optimal control problems[J].IEEE Transactions on Automatic Control,1988,33:333-340.

[10] Maroński R,Rogowski K.Minimum-time running:a numerical approach[J].Acta of Bioengineering and Biomechanics/Wroclaw University of Technology,2011,13(2):83-86.

[11] Jianhui Z,Peng Q.Lunar soft landing trajectory optimization by a Chebyshev pseudospectral method[C]//Computer Science and Automation Engineering (CSAE),2011 IEEE International Conference on IEEE,2011,2:425-430.

[12] 劉熔潔,李世華.橋式吊車系統(tǒng)的偽譜最優(yōu)控制設(shè)計(jì)[J].控制理論與應(yīng)用,2013,30(8):981-989.

[13] 陳琦,王中原,常思江,等.不確定飛行環(huán)境下的滑翔制導(dǎo)炮彈方案彈道優(yōu)化[J].航空學(xué)報(bào),2014,35(9):2593-2604.

[14] Vinh Nguyen X,Adolf Busemann,Robert D.Hypersonic and planetary entry flight mechanics[J].NASA STI/Recon Technical Report A,1980,81:16245.

[15] 張毅,肖龍旭,等.彈道導(dǎo)彈彈道學(xué)[M].長沙:國防科技大學(xué)出版社,2005:82-85.

[16] Jorris Timothy R.Common aero vehicle autonomous reentry trajectory optimization satisfying waypoint and no-fly zone constraints[D].Wright-Patterson:Air Force Institute of Technology,2007:88-94.

[17] 趙欣,閆循良,張金生,等.助推-滑翔導(dǎo)彈再入彈道快速優(yōu)化[J].固體火箭技術(shù),2012,35(4):427-433.

[18] Gill Philip E,Walter Murray,Michael A.SNOPT:An SQP algorithm for large-scale constrained optimization[J].SIAM Journal on Optimization,2002,12(4):979-1006.

[19] Rahimi A,Dev Kumar K,Alighanbari H.Particle swarm optimization applied to spacecraft reentry trajectory[J].Journal of Guidance,Control,and Dynamics,2012,36(1):307-310.

[20] Gao,Xiao-zhi,Wu Ying,Wang Xiao-lei,et al.Trajectory optimization in reentry phase for hypersonic gliding vehicles using swarm intelligence algorithms[M].Practical Applications of Intelligent Systems.Springer Berlin Heidelberg,2012:361-371.

[21] 符俊,蔡洪,李安梁.基于 Gauss 偽譜法的航天器氣動(dòng)力輔助平面變軌問題研究[J].國防科技大學(xué)學(xué)報(bào),2011,33(6):95-99.

(編輯：薛永利)

Rapid trajectory optimization for hypersonic reentry vehicle

DONG Chun-yun1, GUO Zhi1, ZHAO Pei-bo1, CAI Yuan-li1, YU Zhen-hua2

(1.School of Electronic and Information Engineering,Xi'an Jiaotong University,Xi'an 710049,China;2.School of Information and Navigation, Air Force Engineering University,Xi'an 710077,China)

Trajectory optimization of a hypersonic reentry vehicle was investigated via the optimal control method-Chebyshev Pseudospectral Method (CPM). Upon the difficulties of long-range reentry process with many constraints, a combination strategy of the linear initial guess and the nodes update for optimization was adopted. Both the angle of attack and the bank angle were chosen as the control variables to minimize the flight time. The trajectory optimization problem was then translated to a nonlinear programming problem (NLP) via CPM and optimized by the SNOPT software package to obtain a general algorithm for rapid trajectory optimization. In the trajectory optimization of a hypersonic reentry vehicle，simulation results show high efficiency and rapid convergence of the presented strategy compared to the Particle Swarm Optimization (PSO) algorithm.

hypersonic vehicle; trajectory optimization; reentry; Chebyshev pseudospectral method

2014-08-26；

：2014-11-04。

國家自然科學(xué)基金(61202128)；宇航動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室開放基金(2011ADL-JD0202)。

董春云(1989—)，女，博士生，研究方向?yàn)轱w行器軌跡優(yōu)化與方法評估。E-mail:dongdong_2007y@163.com

V412.4

A

1006-2793(2015)06-0757-07

10.7673/j.issn.1006-2793.2015.06.002