基于序列凸優化的高超聲速滑翔式再入軌跡快速優化

趙黨軍，梁步閣，楊德貴，時 偉

(中南大學航空航天學院，長沙 410083)

基于序列凸優化的高超聲速滑翔式再入軌跡快速優化

趙黨軍，梁步閣，楊德貴，時 偉

(中南大學航空航天學院，長沙 410083)

針對具有熱流、動壓、過載以及多個禁飛區約束的再入軌跡優化問題，提出采用序列凸優化方法快速求解。利用歸一化時間作為自變量解決終端時間自由問題，并引入輔助控制變量以減少序列優化結果中的高頻振蕩，在此基礎上，通過線性化、離散化和非凸約束的凸化處理，將非凸非線性優化問題轉化為二階錐規劃(Second Order Conic Programming, SOCP)問題，然后采用凸優化求解算法快速求解。數值優化結果與對比驗證表明該方法能快速高效求解多約束條件下的再入軌跡優化問題，且計算效率和性能均優于傳統的非線性規劃方法。

序列凸優化；再入軌跡優化；無損凸化；二階錐規劃

0 引言

升力式高超聲速飛行器具有速度高、突防能力強等優點，近年來受到了國內外各相關領域的密切關注，出現了大量的研究成果[1]。與現有的亞聲速/超聲速飛行器相比，高超聲速飛行器在近空間的稠密大氣層內進行長時間、遠距離飛行，其飛行軌跡必須滿足熱流、動壓、過載等方面的約束。另外，考慮到敵方反導系統的探測、追蹤以及攔截威脅，在飛行軌跡中還應考慮禁飛區約束，以避開敵方部署反導系統的區域。因此，飛行器軌跡規劃是一個典型的強約束條件的動態優化問題[2]。

到目前為止，求解該類問題的方法總體來說可以分為兩類：直接法和間接法[3]。所謂間接法主要源于龐特里亞金極大值原理，在狀態變量基礎上擴展協態變量獲得哈密頓邊值問題(Hamiltonian Boundary Value Problem, HBVP)，通過一系列復雜的解析推導計算獲得軌跡優化的解[4-5],這也使HBVP問題的求解異常復雜。為避免復雜的解析推導計算，間接法應運而生。間接法將無限維動態優化問題通過參數化技術轉化為有限維的非線性規劃(Nonlinear Programming, NLP)問題,然后通過序列二次規劃(Sequential Quadratic Programming, SQP)[6]方法求解,如廣泛使用的譜方法[7-10]就屬于典型的間接法。目前，已經有一些較成熟的基于間接法的軌跡優化軟件包可供使用，如GPOPS[9]、GPOCS[11]等。然而，SQP方法無法保證收斂速度，也無法保證必然獲得全局最優解。在實際使用過程中，基于SQP方法的軌跡優化，對于初值異常敏感，且優化速度較慢。

與一般NLP方法相比，凸優化方法具有多項式時間復雜度，收斂速度快，且具有獨特的理論優勢——定義域內必然收斂，因此受到青睞。特別是在20世紀80年代，隨著快速求解凸優化問題的內點法出現之后[12]，凸優化方法開始進入很多工程優化領域。最近凸優化方法也開始用于求解軌跡優化這類復雜的非線性動態優化問題。文獻[13-15]主要研究了帶動力軟著陸的軌跡優化問題，提出無損凸化處理方法，將非凸問題轉化為凸優化問題，從而利用對偶內點法快速求解。類似于序列二次規劃方法，Liu等給出了求解一般非凸非線性最優問題的序列凸優化算法[16]，該算法采用軌跡線性化、離散化措施將非線性動態約束轉化為凸約束，同時為保證線性化系統對原非線性系統的有效近似，增加了可信域約束。在此基礎上，Liu將該方法用于再入軌跡優化之中[17]：在速度-攻角剖面已知條件下，以能量為自變量對原問題進行轉化以避免終端時間自由問題；同時將熱流、動壓以及過載約束轉化為速度-高度邊界約束，在此基礎上利用序列凸優化算法進行迭代求解獲得再入最優軌跡。但以能量為自變量時，在飛行高度較高時，能量變化不顯著，導致在高度大于65km時的軌跡無法進行離散化并優化。

針對上述問題，提出以歸一化時間作為自變量，并選擇泛化升力系數和傾側角作為控制量，同時考慮多禁飛區約束，對高超聲速飛行器再入軌跡優化問題的序列凸優化方法做出進一步改進：同時優化泛化升力系數和傾側角兩個控制變量，避免能量為自變量的缺陷。數值優化結果表明，本文所提方法是行之有效的。

1 問題描述

1.1 CAV再入運動模型

基于上述假設和關于泛化升力系數的定義，建立再入飛行器無量綱運動方程為[2]：

(1)

(2)

1.2 約束條件

再入過程中，CAV軌跡需滿足一系列約束以保證飛行器安全。總體來說飛行約束分為等式約束和不等式約束。

(1)等式約束

根據事先確定的飛行器起點位置和目標位置，形成如下等式約束：

Φ[x(t0),x0]=0,Ψ[x(tf),xf]=0

(3)

值得注意的是，很多時候只指定了終端狀態中的部分狀態，如只規定了終端時的位置，也有可能只是規定了上下限，于是形成了不等式約束。

(2)不等式約束

不等式約束主要包括關于熱流、動壓和過載的約束

(4)

(5)

(6)

ρ=ρ0e-βReh

(7)

其中，h=(R-Re)/Re無量綱高度。為方便后續處理，對式(4)～式(6)進行歸一化處理，得

≤03×1

(8)

其中，

(9)

除上述不等式約束外，施加于控制量上的限制也應考慮：傾側角σ通常限制于某些特定區間(如[-π/3，π/3])以保證飛行穩定性；泛化升力系數λ因氣動特性有所限制，因此有如下關于控制變量的約束

(10)

1.3 優化問題描述

綜合考慮上述CAV運動模型以及相關約束條件，以時間最優為指標，CAV再入軌跡優化問題可以描述為

subject to: (2),(3),(8),(9),(10)

2 序列凸優化算法

2.1 問題重新描述

此外，若直接選擇傾側角σ和泛化升力系數λ作為控制變量，則會出現如文獻[17]指出的控制量高頻抖動現象，對飛行穩定不利，因此引入新的輔助控制變量，令

u1=λcosσ,u2=λsinσ,

(11)

控制約束可以改寫為

(12)

同時，應保證式(13)的等式約束成立。

(13)

至此，問題 P0′與下面的問題P1等價。

P1:mint(1)

(14)

(15)

(16)

(17)

(18)

(19)

其中，(·)′=d/dτ,τ∈[0,1]，約束(17)為不等式約束(8)和(9)的一般形式。問題P1中的指標函數J=t(1),表示為擴展狀態t在終端時刻τ=1時的大小，即tf。顯然，問題P0中的Lagrange項轉化為問題P1中的Mayer項，成為一個線性函數，從而使得問題P1滿足SOCP問題對于指標函數為凸函數的要求，且問題P1的解必然是問題P0的解。

2.2 線性化

(20)

為保證線性化系統式(20)合理逼近原系統式(15)，狀態及控制偏差應限制于可信域內，即

(21)

2.3 凸化處理

問題 P1中，路徑約束式(17)和控制約束式(18)和式(19)都是非凸的，為使問題能用凸優化方法求解，必須對這些非凸約束進行凸化處理。

2.3.1 路徑約束的凸化處理

路徑約束式(17)中包含m個不等式約束，即

(22)

(23)

2.3.2 控制約束的凸化處理

由于輔助控制變量的引入，增加了關于控制量的等式約束式(19)，顯然該約束非凸，對其進行松弛處理，即將等式約束放寬為不等式約束

(24)

2.4 離散化

(25)

其中，j=0,…,N。重新排列并合并相同項，得

(26)

經過上述線性化、凸化以及離散化處理后，問題P1的最優解可以通過求解下面序列優化問題P2來進行求解。

P2:mint(1)

(27)

(28)

(29)

(30)

(31)

0≤u3≤4

(32)

-tan(σmax)u1≤u2≤tan(σmax)u1

(33)

(34)

其中，τ∈[0,1],(i=1,…,m,k=1,2,…)。可以證明經松弛處理后問題P2的解是原問題P1的解(參見文獻[17])。

2.5 序列凸優化算法

顯然，問題P2是SOCP問題，其目標函數為線性的，且各約束或者為線性函數，或者為二階錐約束形式。如果選擇足夠小的離散化步長，SOCP問題P2的解與優化問題P1充分接近。求解問題P2的序列凸優化算法流程如下：

Step3:檢查下面的收斂條件是否滿足

(35)

3 數值優化結果

本文以CAV為例在MATLAB環境中進行數值優化求解以驗證上述算法的有效性。為方便算法實施，采用Yalmip建模工具箱[19]進行優化問題建模，利用輕量級凸優化求解算法包ECOS[20]進行SOCP問題求解。所有優化算法運行的計算機配置為： CPUi7-4790@3.6GHz, 8G內存。

表1給出了CAV任務描述，給出了任務起點和目標點的水平位置、高度，并給出了禁飛區的中心點坐標以及半徑。優化過程中禁飛區次序未知，但根據算法流程步驟2中的逐點檢測可以方便地施加禁飛區約束。第一個禁飛區的半徑遠小于CAV的轉彎能力，第二個禁飛區半徑則足夠大。

表1 CAV任務描述

式(35)中的迭代收斂停止條件設置為δ=0.03，離散化區間Δτ=1/300，即N=300。

為對比說明本文所提方法的有效性，將序列凸優化求解結果(記為SOCP)與偽譜法求解結果(記做GPOPS)進行了對比，其中偽譜法求解主要基于開源Matlab工具箱GPOPS[21]進行。

數值求解過程中，序列凸優化SOCP經過7次迭代獲得最優解，整個優化計算的CPU時間為8.918s，而GPOPS優化則需要21.109s；SOCP和GPOPS求得的最優飛行時間分別是2869.2s和3061.7s，顯然，序列凸優化方法無論是效率還是性能上都要明顯優于GPOPS。兩種方法的優化結果如圖2～圖6所示。

圖2為地面航跡，可以看出SOCP解和GPOPS解都滿足兩個禁飛區約束，但SOCP的軌跡更逼近禁飛區的邊界。圖3是高度曲線和速度曲線，兩種方法的解均嚴格滿足終端高度和速度約束，其中SOCP解的高度曲線變化更加平滑。兩種優化方法得到的彈道傾角和航向角曲線變化趨勢與幅度(如圖4所示)基本保持一致，彈道傾角變化較小，航向角平滑。圖5給出泛化升力系數和傾側角曲線，整個過程均滿足控制變量約束條件，但GPOPS得到的控制量抖動幅度較大，而SOCP由于引入新的控制變量進行松弛處理，抖動幅度較小。圖6給出了兩種方法得到的熱流、動壓以及過載曲線，從圖中可以看出，盡管SOCP方法中過載曲線末端較大，但所有約束在整個過程中均滿足給定條件。

4 結論

具有熱流、動壓、過載以及禁飛區約束的再入軌跡優化問題是一個強約束的非線性動態規劃問題，本文通過引入歸一化時間和輔助控制變量，結合軌跡線性化和無損凸化處理技術，將原問題轉化為標準的SOCP問題，利用開源凸優化求解器ECOS求解，快速獲得最優軌跡。歸一化時間盡管增加了控制變量，但有效避免了終端時間自由所帶來的離散化問題，同時也避免了以能量為自變量時無法從起始點進行離散化的缺陷。輔助控制變量的引入則在一定程度上避免了序列優化過程中容易引起的控制量高頻振蕩的不利因素，提高了優化解的工程適用性。數值優化結果以及對比實驗結果表明，本文所提方法是一種行之有效的再入軌跡快速優化方法，計算效率和性能均優于傳統的非線性規劃方法。

[1] Lu P. Entry guidance: a unified method[J]. Journal of Guidance Control and Dynamics, 2014, 37(3):713-729.

[2] Jorris T R, Cobb R G.Three-dimensional trajectory optimization satisfying waypoint and no-fly zone constraints[J]. Journal of Guidance Control and Dynamics, 2009, 32(2): 551-572.

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

[4] Poustini M J, Esmaelzadeh R, Adami A. A new approach to trajectory optimization based on direct transcription and differential flatness[J]. Acta Astronautica, 2015,107:1-13.

[5] Petropoulos A E, Sims J A. A review of some exact solutions to the planar equations of motion of a thrusting spacecraft[C]. NASA Jet Propolsion Laboratory Technical Reports,USA,2002.

[6] Betts J T. Very low-thrust trajectory optimization using a direct SQP method[J].Journal of Computational and Applied Mathematics, 2000, 120:27-40.

[7] Williams P.Hermite-Legendre-Gauss-Lobatto direct transcription in trajectory optimization[J]. Journal of Guidance Control and Dynamics,2009, 32(4): 1392-1395.

[8] Zhao J, Zhou R， Jin X. Reentry trajectory optimization based on a multistage pseudospectral method[J]. Scientific World Journal, 2014, 2014:1-13.

[9] Garg D, Patterson M， Darby C, et al. Direct trajectory optimization and costate estimation of general optimal control problems using a Radau pseudospectral method[J]. Computational Optimization and Applications,2011, 49(2):335-358.

[10] Chai D, Fang Y W，Wu Y L，et al. Boost-skipping trajectory optimization for air-breathing hypersonic missile[D]. Astronautical Science and Technology, 2007.

[11] Huntington G T. Advancement and analysis of a Gauss pseudospectral transcription for optimal control problems[J]. Massachusetts Institute of Technology,2007.

[12] Andersen E D， Roos C， Terlaky T. On implementing a primal-dual interior-point method for conic quadratic optimization[J]. Mathematical Programming, 2003, 95(2): 249-277.

[13] Blackmore L,Acikmese B，Scharf D P .Minimum-landing-error powered-descent guidance for Mars landing using convex optimization[J]. Journal of Guidance Control, and Dynamics, 2010, 33(4): 1161-1171.

[16] Liu X F, Lu P. Solving nonconvex optimal control problems by convex optimization[J]. Journal of Guidance Control and Dynamics, 2014, 37(3): 750-765.

[17] Liu X F, Shen Z J，Lu P. Entry trajectory optimization by second-order cone programming[J]. Journal of Guidance Control and Dynamics (Articles in Advance),2015,39(2): 1-15.

[18] Jorris T R. Common aero vehicle autonomous reentry trajectory optimization satisfying waypoint and no-fly zone constraints[D]. Doctor, Engineering and Management, Air University, Ohio, US, 2007.

[19] L?fberg J.YALMIP : a toolbox for modeling and optimization in MATLAB2004[C]. Proceedings of the CACSD Conference, Taipei, Taiwan, 2004.

[20] Domahidi A， Chu E， Boyd S.ECOS: an SOCP solver for embedded systems[C].European Control Conference, Zurich, Switzerland, 2013.

[21] Garg D, Hager W W，Rao A V .Pseudospectral methods for solving infinite-horizon optimal control problems[J]. Automatica,2011,47: 829-837.

Rapid Planning of Reentry Trajectory viaSequential Convex Optimization

ZHAO Dang-Jun，LIANG Bu-Ge，YANG De-Gui，SHI Wei

(School of Aeronautics and Astronautics, Central South University, Changsha 410083， China)

A sequential convex optimization scheme is proposed for rapid solving the reentry trajectory optimization problem constrained by heat flux, dynamic pressure, normal load, and multiple no-fly zones. The normalized time is used as the independent variable to accommodate the problem of free terminal time, and auxiliary control variables are introduced to alleviate the high frequent chatter in the sequential convex optimization solution. Further, the linearization, discretization, and convexification techniques are used to convert the original concave and nonlinear optimization problem into a standard second order conic programming problem, which can be rapid solved by convex algorithms. Numerical results and comparison study reveal that the proposed method is efficient and effective to solve the problem of reentry trajectory optimization with multiple constraints, and the computational efficiency and performance of the proposed method is superior to that of the classical nonlinear programming.

Sequential convex optimization; Reentry trajectory optimization; Lossless convexification; Second order conic programming

2017-03-20；

2017-04-13基金項目：湖南省自然科學基金(14JJ3024)

趙黨軍(1978-)，男，博士,副教授，主要從事飛行器制導與控制研究。E-mail:zhao_dj@esu.edu.cn

V411

A

2096-4080(2017)01-0034-07