高超聲速飛行器再入軌跡快速優化研究

董春云，郭 志，趙培博，蔡遠利，于振華

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

?

高超聲速飛行器再入軌跡快速優化研究

董春云1，郭 志1，趙培博1，蔡遠利1，于振華2

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

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

高超聲速飛行器；軌跡優化；再入；Chebyshev偽譜方法

0 引言

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

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

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

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

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

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

1 軌跡優化模型

1.1 動力學模型及歸一化處理

假設地球為旋轉圓球，高超聲速飛行器無動力再入的非線性動力學模型可用6個狀態變量來描述。由于變量間存在不可公度性，取值范圍和量綱不盡相同，為了提高后續優化計算效率與求解精度，需要對模型進行歸一化處理。極坐標系下無量綱運動方程[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分別為再入飛行器質量和氣動參考面積；ρ為大氣密度；CL、CD分別為升力和阻力系數。

氣動力的計算采用1976年美國標準大氣模型USSA76[15]，在0～100 km高度范圍內大氣劃分為8層，根據飛行高度計算空氣密度ρ。CL、CD通過氣動力數據擬合得到[16]。

1.2 控制變量

在高超聲速飛行條件下，氣動力系數可近似表示為攻角的函數，同時將攻角α和傾側角σ作為優化控制變量[17]，能夠充分體現出攻角在再入飛行中的控制調節作用。因此，控制變量為

u=(ασ)T

1.3 約束條件

(1)過程約束

飛行器再入飛行是一個高動態過程，要求嚴格滿足過程約束，以保證飛行器在結構和熱防護上的可靠性，主要包括熱流密度、動壓、過載和控制變量約束，即

(7)

(8)

(9)

(10)

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

(2)終端約束

飛行器再入終端狀態對末制導段的飛行有重要的影響，終端約束與再入飛行的任務相關，一般包括速度、飛行高度和地理位置約束等：

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

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

(11)

其中，εV，εr，εθ，εφ分別是終端狀態的誤差界，考慮精確制導時，可設置εV=0，εr=0，εθ=0，εφ=0。

1.4 性能指標

優化目標按照設計要求可以有不同的形式。由于再入環境惡劣，最短的飛行時間能夠減小氣動加熱，降低能量損耗，縮短敵方響應時間。因此，以最短飛行時間為性能指標：

J=min(τf)

(12)

1.5 優化模型

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

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

2 Chebyshev偽譜法與軌跡優化

2.1 方法描述

對于上述連續時間最優控制問題，Chebyshev偽譜法的基本求解思路為：將連續時間狀態變量和控制變量在一系列CGL點上離散，并以這些離散點為節點構造Largrange插值多項式來逼近真實狀態和控制；通過對全局插值多項式求導來近似狀態變量對時間的導數，將微分方程約束轉換為代數約束；性能指標中的積分項由Clenshaw-Curtis數值積分計算。通過上述方法，可將最優控制問題轉化為具有一系列代數約束的NLP問題[5,8]。

CPM的具體計算步驟簡述如下：

(1)時間區間轉換

CPM的CGL離散點定義在[-1,1]區間上，需要將原最優控制問題的時間域τ從[τ0,τf]線性轉換到t∈[-1,1]區間，轉換式為

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

(2)離散節點計算

CPM離散節點選取為N階Chebyshev多項式TN(t)=cos(Ncos-1t)的極值點，即CGL點，它們不均勻的分布在[-1,1]區間上，顯式計算式為

顯然，離散節點tk滿足t0=-1，tN=1。

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

取上述N+1個離散點處的狀態變量和控制變量，分別構造Lagrange插值多項式作為連續狀態和控制的近似。真實狀態變量x(t)與控制變量u(t)的近似表達式為

(13)

(14)

其中，Lagrange插值基函數：

由Lagrange插值的性質可知，離散節點處的狀態近似值與實際狀態相等，控制近似值與實際控制相等。

(4)動態約束處理

對式(13)求導，得到狀態向量在tk點處的導數近似表達式為

(15)

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

(16)

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

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

(5)性能指標積分近似

對于優化性能指標中存在積分項的情況，利用Clenshaw-Curtis數值積分對其進行近似。對于 [-1,1]區間上的任一連續函數p(t)，其積分可用N+1個CGL離散點處的函數值累加和近似，即

(17)

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

通過以上對最優控制問題在CGL節點的近似，將原連續時間最優控制問題轉化為以下離散化的非線性規劃問題，即確定CGL節點處的離散狀態變量X、離散控制變量U和終端時間τf，使得性能指標J最小，并滿足離散化后的狀態方程約束、過程約束與終端約束等。

2.2 軌跡優化求解策略

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

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

如圖1所示，首先采用線性化方法獲得設計變量初值，求解較少CGL節點下的初步優化數值解；然后，增加CGL節點，利用初步優化數值解插值獲得需要的設計變量初值，從而提高求解效率。

所謂線性化方法即根據設計變量初始位置x0與終端約束xf構造簡單線性函數，具體可表述為

節點更新過程中，通過循環依次增加節點的數目，計算優化結果的數值積分終端約束偏差，當最大相對誤差小δ時，認為精度符合要求，停止計算。圖1中，變量nx、ny可根據仿真效果靈活調整。

2.3 CPM數值驗證

在應用CPM求解再入軌跡優化問題之前，首先選擇標準最速降線問題對本文算法的有效性進行驗證，并分析其求解精度。

圖1 軌跡優化求解策略Fig.1 Combination strategy for trajectory optimization

最速降線問題即尋找一條最優曲線，使得質點在重力作用下，沿該曲線從給定點到不在它垂直下方的另一點的時間最短。其數學模型可以概括為

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

該問題的解析形式為

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

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

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

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

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

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

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

3 仿真校驗與結果分析

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

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

3.1 參數設定

飛行器結構參數和氣動參數參考1998年美國洛克希勒-馬丁公司設計的高超聲速飛行器CAV-H[16]，該氣動參考面積為0.483 9 m2，質量為907.2 kg，長度為2.717 8～3.627 6 m。

3.2 結果及分析

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

初始優化結果獲得的積分軌跡中，終端最大相對誤差0.496。更新CGL節點數目，重新優化并積分，當節點個數K=15時，終端最大相對誤差為0.117。繼續更新節點數目，當K=28時，獲得終端高度20.54 km，速度1 018.21 m/s，經度236.36°，緯度36.95°，最大相對誤差為0.027，小于最大誤差閾值，停止更新。

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

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

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

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

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

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

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

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

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

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

圖9 CPM動壓變化圖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 結論

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

[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] 謝愈,劉魯華,湯國建,等.多約束條件下高超聲速滑翔飛行器軌跡優化[J].宇航學報,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] 劉熔潔,李世華.橋式吊車系統的偽譜最優控制設計[J].控制理論與應用,2013,30(8):981-989.

[13] 陳琦,王中原,常思江,等.不確定飛行環境下的滑翔制導炮彈方案彈道優化[J].航空學報,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] 張毅,肖龍旭,等.彈道導彈彈道學[M].長沙:國防科技大學出版社,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] 趙欣,閆循良,張金生,等.助推-滑翔導彈再入彈道快速優化[J].固體火箭技術,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 偽譜法的航天器氣動力輔助平面變軌問題研究[J].國防科技大學學報,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。

國家自然科學基金(61202128)；宇航動力學國家重點實驗室開放基金(2011ADL-JD0202)。

董春云(1989—)，女，博士生，研究方向為飛行器軌跡優化與方法評估。E-mail:dongdong_2007y@163.com

V412.4

A

1006-2793(2015)06-0757-07

10.7673/j.issn.1006-2793.2015.06.002