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

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

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

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

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

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

序列凸優(yōu)化；再入軌跡優(yōu)化；無損凸化；二階錐規(guī)劃

0 引言

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

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

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

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

1 問題描述

1.1 CAV再入運動模型

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

(1)

(2)

1.2 約束條件

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

(1)等式約束

根據(jù)事先確定的飛行器起點位置和目標(biāo)位置，形成如下等式約束：

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

(3)

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

(2)不等式約束

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

(4)

(5)

(6)

ρ=ρ0e-βReh

(7)

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

≤03×1

(8)

其中，

(9)

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

(10)

1.3 優(yōu)化問題描述

綜合考慮上述CAV運動模型以及相關(guān)約束條件，以時間最優(yōu)為指標(biāo)，CAV再入軌跡優(yōu)化問題可以描述為

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

2 序列凸優(yōu)化算法

2.1 問題重新描述

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

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

(11)

控制約束可以改寫為

(12)

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

(13)

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

P1:mint(1)

(14)

(15)

(16)

(17)

(18)

(19)

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

2.2 線性化

(20)

為保證線性化系統(tǒng)式(20)合理逼近原系統(tǒng)式(15)，狀態(tài)及控制偏差應(yīng)限制于可信域內(nèi)，即

(21)

2.3 凸化處理

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

2.3.1 路徑約束的凸化處理

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

(22)

(23)

2.3.2 控制約束的凸化處理

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

(24)

2.4 離散化

(25)

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

(26)

經(jīng)過上述線性化、凸化以及離散化處理后，問題P1的最優(yōu)解可以通過求解下面序列優(yōu)化問題P2來進(jìn)行求解。

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,…)?？梢宰C明經(jīng)松弛處理后問題P2的解是原問題P1的解(參見文獻(xiàn)[17])。

2.5 序列凸優(yōu)化算法

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

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

(35)

3 數(shù)值優(yōu)化結(jié)果

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

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

表1 CAV任務(wù)描述

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

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

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

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

4 結(jié)論

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

[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基金項目：湖南省自然科學(xué)基金(14JJ3024)

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

V411

A

2096-4080(2017)01-0034-07