基于序列凸優(yōu)化的多約束軌跡快速優(yōu)化

楊 奔 李天任 馬曉媛

中國運(yùn)載火箭技術(shù)研究院，北京100071

0 引言

再入飛行器利用氣動力在大氣層內(nèi)進(jìn)行長時(shí)間、遠(yuǎn)距離、高馬赫數(shù)的飛行，期間會受到極其復(fù)雜的非線性熱流、動壓、過載約束以及終端的高精度約束，使其軌跡規(guī)劃問題成為一個極富挑戰(zhàn)性的問題，引起了國內(nèi)外許多學(xué)者的研究興趣[1]。傳統(tǒng)的再入彈道優(yōu)化設(shè)計(jì)方法主要有基于極大值原理的間接法和基于非線性規(guī)劃理論的直接法[2]等優(yōu)化算法，兩者或因巨大的計(jì)算量，或因不易收斂，都只適合離線的軌跡優(yōu)化[3-6]。雖然航天飛機(jī)所采用的基于阻力加速度剖面的軌跡規(guī)劃[7]，曾取得過輝煌的成績，但該方法的簡化過程限制了飛行器的機(jī)動能力?？偟膩碚f，基于極大值原理的直接法推導(dǎo)最優(yōu)解的過程較為繁瑣，對于復(fù)雜問題幾乎無法得到解析解；以偽譜法為代表的直接法對多階段等復(fù)雜問題的建模過程比較麻煩；而以粒子群法等現(xiàn)代啟發(fā)式算法對于計(jì)算量大、優(yōu)化參數(shù)多的復(fù)雜模型優(yōu)化問題，需反復(fù)交叉、變異、迭代，計(jì)算效率較低。

強(qiáng)耦合性、強(qiáng)非線性使得再入軌跡規(guī)劃問題難以直接求解，通常需要采用一系列易于求解的子問題逼近原問題。凸優(yōu)化具有多項(xiàng)式復(fù)雜度，局部最優(yōu)解即全局最優(yōu)解等特點(diǎn)，加上獨(dú)特的內(nèi)點(diǎn)法，使得數(shù)學(xué)上所有易解的規(guī)劃問題幾乎都指向了凸優(yōu)化[8-9]。美國加州理工大學(xué)的文獻(xiàn)[10-15]最早開展了基于凸優(yōu)化的軌跡優(yōu)化研究。文獻(xiàn)[10-11]針對火星軟著陸問題最先提出了無損凸化(Lossless Convexification)的概念，使用松弛的凸形式約束代替原非凸約束，將原問題轉(zhuǎn)化為凸問題，并利用極大值原理對二者的等價(jià)性進(jìn)行了證明，保證轉(zhuǎn)化過程的“無損性”。針對具有非凸約束的一般形式線性系統(tǒng)最優(yōu)控制問題，文獻(xiàn)[12-13]討論了無損凸化的方法，對原問題與凸優(yōu)化問題的等價(jià)性進(jìn)行了證明。相關(guān)理論研究成果已經(jīng)在NASA的G-FOLD飛行試驗(yàn)中得到了成功應(yīng)用[14]。

近年來，許多學(xué)者都在嘗試將凸優(yōu)化引入到再入軌跡規(guī)劃問題中來。文獻(xiàn)[15]在給定攻角剖面的條件下，通過變量替換得到以能量為自變量的動力學(xué)方程，采用序列凸優(yōu)化解決了再入軌跡優(yōu)化問題。文獻(xiàn)[16]則以時(shí)間為自變量，以傾側(cè)角的導(dǎo)數(shù)為控制量，求解了終端時(shí)間固定的軌跡規(guī)劃?，F(xiàn)有的基于凸優(yōu)化的滑翔再入軌跡規(guī)劃方法基本上都直接使用了小擾動線性化來處理微分方程的非線性，同時(shí)添加了信賴域約束來保證子問題向原問題逼近。雖然目前的序列凸優(yōu)化技術(shù)簡單可行，尋優(yōu)時(shí)間較傳統(tǒng)直接法有一定優(yōu)勢，但是在實(shí)際仿真計(jì)算中還會遇到以下問題：1)反復(fù)迭代的序列凸優(yōu)化技術(shù)本質(zhì)上還是牛頓法，其收斂范圍小，并且在迭代后期可能會出現(xiàn)類似與牛頓法的Maratos效應(yīng)，導(dǎo)致迭代尋優(yōu)失敗。2)當(dāng)初值給的比較粗糙時(shí)，子問題的可行性將是一個嚴(yán)峻的問題，可能出現(xiàn)“偽不可行”[17](artificial infeasibility)的情況使得迭代失敗。

針對上述序列凸優(yōu)化的不足之處，本文提出了幾種改善其收斂性的方法。首先通過增加虛擬控制量，可避免迭代初期由于控制能力不足帶來的“偽不可行”問題，降低該算法對初始參考軌跡的依賴。其次通過B樣條曲線離散控制量，抑制由于離散帶來的鋸齒化現(xiàn)象。最后，引入一種“回溯搜索”的方法，改善迭代過程中狀態(tài)量的振蕩，從而加快算法的收斂。

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

1.1 再入動力學(xué)模型

考慮地球自轉(zhuǎn)的影響，建立如下再入飛行器無量綱運(yùn)動方程：

(1)

其中，r,θ,φ,V,γ,ψ分別為無量綱化的地心距、經(jīng)度、緯度、對地速度、航跡角、航向角，地心距按地球平均半徑歸一化，對地速度按第一宇宙速度歸一化。速度σ為傾側(cè)角，無量綱化升力L、阻力D的表達(dá)式為

(2)

其中，ρ是大氣密度，為簡化問題，采用指數(shù)形式的大氣模型ρ=ρ0exp(-βh)，S為飛行器的氣動參考面積，升力系數(shù)CL、阻力系數(shù)CD均認(rèn)為是攻角和馬赫數(shù)的函數(shù)。

1.2 基本約束

飛行器再入過程中，為保證飛行任務(wù)的成功，需滿足熱流、動壓、過載等路徑約束，可表示為

(3)

q=0.5ρV2≤qmax

(4)

(5)

至此，得到該最優(yōu)控制問題(P1)如下

(6)

(7)

(8)

(9)

式(7)中動力學(xué)方程可表示為

(10)

其中，第1項(xiàng)是與狀態(tài)變量相關(guān)的項(xiàng)，第2項(xiàng)是與控制量相關(guān)的項(xiàng)，第3項(xiàng)是與地球自旋相關(guān)項(xiàng)。可以看出，選擇傾側(cè)角變化率為控制量，使得控制量與狀態(tài)量自然解耦，更加有利于優(yōu)化收斂。

1.3 凸化處理

由1.2節(jié)可知，P1是一個高度的非線性最優(yōu)控制問題，無法直接采用凸優(yōu)化方法求解。本文采用基于信賴域約束的序列凸優(yōu)化方法。通過求解一系列的凸優(yōu)化子問題，逐漸逼近原有問題的解。

將式(10)一階泰勒展開，由于地球自旋速度較慢，可認(rèn)為fΩ(x)≈fΩ(xk) ，則線性化的動力學(xué)方程為

(11)

其中，

(12)

(13)

為避免迭代初期“偽不可行”的問題，降低對初值的依賴程度，同時(shí)彌補(bǔ)線性化帶來的誤差，在式(11)中增加一個虛擬控制量w，可得

(14)

可在迭代時(shí)，在目標(biāo)函數(shù)中控制虛擬控制量的范圍，使其在迭代后期為小量，保證問題的準(zhǔn)確性。

J′=J+λCw

(15)

2 算法表述

2.1 問題離散

為求解無限維最優(yōu)控制問題的數(shù)值解，必須對其進(jìn)行離散。為了簡單起見，狀態(tài)變量像通常一樣離散成等距的數(shù)值點(diǎn)，各個數(shù)值點(diǎn)處的控制量則通過B樣條曲線控制點(diǎn)的表達(dá)式表示，控制點(diǎn)個數(shù)的選取可在實(shí)際仿真中調(diào)整，優(yōu)化過程中的積分項(xiàng)均采用梯形積分表示。仿真結(jié)果顯示這種方法可以十分有效地抑制數(shù)值求解中控制量的鋸齒化現(xiàn)象，更加有利于迭代收斂，并且減少了需要求解的變量個數(shù)。將時(shí)間N等分，可得離散后的動力學(xué)方程。

(16)

(17)

2.2 回溯搜索

大氣層內(nèi)各個狀態(tài)變量與控制量相互耦合，具有極強(qiáng)的非線性，且各個約束集也均為非凸，導(dǎo)致無法直接運(yùn)用內(nèi)點(diǎn)法進(jìn)行求解。受序列二次規(guī)劃算法的啟發(fā)，國內(nèi)外學(xué)者提出了基于信賴域約束的序列凸優(yōu)化技術(shù)。其信賴域約束是為了保障線性化的可靠程度。當(dāng)然，每次外層迭代所使用的參考軌跡的不同使得每次迭代過程并不是同一個問題，也就是說并不是嚴(yán)格的最優(yōu)解。若取固定信賴域，且每次迭代間均采用全步長更新尋優(yōu)結(jié)果，在迭代后期可能會出現(xiàn)如圖1的情況，各個狀態(tài)變量在迭代后期產(chǎn)生劇烈振蕩，難以符合停機(jī)準(zhǔn)則，導(dǎo)致尋優(yōu)失敗。

圖1 振蕩現(xiàn)象

針對這種具有兩個局部極小值，得益于牛頓法的啟發(fā)，本文提出一種“回溯搜索”的收斂辦法。通過設(shè)置阻尼步長和回歸搜索來快速處理這種局部極小值問題。具體過程如下：

(1)初值選取，并線性化，取兩個回溯參數(shù)α∈(0,1],β∈(0,1]以及信賴域Δ0

(2)解決子SOCP問題，得到Δx。

(3)判斷停機(jī)準(zhǔn)則。

(4)判斷回溯條件：IFJ(k+1)>J(k)，x(k+1)=x(k)+αΔx，Δ0=βΔ0，進(jìn)入步驟(1)。

該回溯搜索的關(guān)鍵點(diǎn)在于信賴域的收縮權(quán)值以及阻尼步長權(quán)值的選取，調(diào)試初期可先設(shè)為黃金分割數(shù)，后期在具體的仿真過程中適當(dāng)調(diào)整。綜上所述，基于序列凸優(yōu)化的軌跡優(yōu)化流程如圖2。

圖2 序列凸優(yōu)化的求解流程

3 仿真分析

3.1 問題描述

本節(jié)以美國某一可重復(fù)使用飛行器的再入軌跡優(yōu)化為例，對第2節(jié)提出的軌跡優(yōu)化方法進(jìn)行仿真分析驗(yàn)證。其質(zhì)量m=104305kg，氣動參考面積S=391.22m2，氣動模型如下

CL=-0.041065+0.016292α+0.0002602α2

(18)

(19)

借鑒航天飛機(jī)程序攻角的設(shè)計(jì)方法，采用典型的二次分段標(biāo)稱攻角剖面，其具體表達(dá)形式為

(20)

仿真中暫不考慮飛行時(shí)間對軌跡的影響，規(guī)定總再入時(shí)間tf=1600s，離散點(diǎn)總數(shù)為200，為保證計(jì)算的效率，設(shè)定最高迭代次數(shù)為100次。仿真所需參數(shù)見表2。仿真中，序列凸優(yōu)化算法中信賴域約束和收斂條件如下

(21)

(22)

3.2 結(jié)果分析

本文的仿真均是基于Matlab環(huán)境，調(diào)用Mosek8.0[18]的求解子凸問題的求解器，以初始狀態(tài)和末端狀態(tài)的直線為初始參考軌跡，由于凸優(yōu)化具有多項(xiàng)式時(shí)間復(fù)雜度，且局部最優(yōu)解即全局最優(yōu)解等特點(diǎn)，因此該尋優(yōu)時(shí)間比較短，大約迭代5～7次就基本收斂，整個尋優(yōu)過程耗時(shí)大約為7s，具有在線軌跡計(jì)算的潛力，具體結(jié)果如圖3。

表1 再入飛行參數(shù)

圖3 序列凸優(yōu)化優(yōu)化結(jié)果

圖3(a)是優(yōu)化結(jié)果中的速度-高度曲線和熱流、動壓、過載等路徑約束的示意圖，在提供再入點(diǎn)與終點(diǎn)連線的軌跡作為初始參考值的情況下，序列凸優(yōu)化技術(shù)也很好地求解出滿足各非線性約束的最優(yōu)軌跡，表明該方法對初值的依賴度較低，魯棒性較強(qiáng)。圖3(c)是對控制量分別進(jìn)行均勻離散和B樣條曲線離散的尋優(yōu)結(jié)果，可以發(fā)現(xiàn)，本文采取的方法有效避免了數(shù)值計(jì)算中鋸齒化現(xiàn)象，使得尋優(yōu)結(jié)果更為光滑可行。

為考察序列凸優(yōu)化方法對禁飛區(qū)的規(guī)避機(jī)動能力，在問題P1的基礎(chǔ)上加入禁飛區(qū)約束，該約束模型具體為

(23)

其中，θc和φc分別為禁飛區(qū)中心經(jīng)度與緯度，RC為禁飛區(qū)半徑。很顯然，該約束也是非凸的，也需要進(jìn)行凸化處理，這里依然采用線性化的方式進(jìn)行松弛凸化。

(24)

現(xiàn)分別在仿真過程中設(shè)置2個禁飛區(qū)，其禁飛區(qū)半徑分別為200km和400km，最終經(jīng)過數(shù)次迭代，得到如下結(jié)果：

圖4 考慮禁飛區(qū)的優(yōu)化結(jié)果

由圖4可知，本文所提的序列凸優(yōu)化算法對于禁飛區(qū)有比較強(qiáng)的機(jī)動繞飛能力，即使在初始參考軌跡違背約束的情況下，也能通過虛擬控制量以及數(shù)次迭代的方法將其拉回安全區(qū)域，保證了飛行器能夠安全、可靠地完成再入任務(wù)。

4 結(jié)論

本文針對傳統(tǒng)軌跡優(yōu)化算法效率低、對初值敏感程度高等問題，提出了一種改進(jìn)的序列凸優(yōu)化算法，通過b樣條離散、回溯搜索等策略，提高了算法的穩(wěn)定性、快速性和尋優(yōu)結(jié)果的光滑性。仿真結(jié)果表明，即使是在不符合動力學(xué)約束的初末狀態(tài)連線作為初值的情況下，該算法依然快速、有效地得到了滿足各項(xiàng)約束的最優(yōu)軌跡，具有一定的工程應(yīng)用價(jià)值。