基于混沌多項式的再入滑翔魯棒軌跡凸優(yōu)化

閆循良, 王培臣, 夏文杰, 王寬, 楊春偉

(1.西北工業(yè)大學 航天學院 陜西省空天飛行器設(shè)計重點實驗室, 陜西 西安 710072;2.中國人民解放軍96901部隊,北京 100094)

軌跡優(yōu)化設(shè)計技術(shù)是高超聲速滑翔飛行器的關(guān)鍵技術(shù)之一,也是近年來研究的熱點問題[1-3]。多種復雜約束和典型飛行任務(wù)條件給高超聲速滑翔彈道優(yōu)化設(shè)計帶來了嚴峻挑戰(zhàn),這使得以直接法[3-4]為代表的數(shù)值優(yōu)化方法成為了解決該問題的主要途徑。其中,偽譜法以較少的計算量和較高的精度優(yōu)勢被普遍使用[3],但仍難以實現(xiàn)多約束復雜條件下的彈道快速生成,而凸優(yōu)化方法[4]因具有多項式復雜度、收斂速度快等特點,已被廣泛應(yīng)用于高超聲速等航空航天領(lǐng)域的軌跡優(yōu)化任務(wù)[4-7]。然而,現(xiàn)有的高超聲速軌跡優(yōu)化研究大都未考慮不確定性的影響,即軌跡優(yōu)化研究通常僅注重標稱情況下軌跡性能提升,而未考慮真實情況下軌跡的魯棒性和可靠性。因此,部分學者逐步將研究重點聚焦于不確定性量化傳播技術(shù)[8-10]以及考慮參數(shù)不確定性的魯棒軌跡優(yōu)化方法[10-15]。

典型的不確定性量化及傳播方法包括蒙特卡洛法(Monte Carlo,MC)、線性協(xié)方差分析法、無跡變換法以及混沌多項式展開法等[8-12],這些方法在不確定性量化傳播及魯棒優(yōu)化中得到了較為廣泛的應(yīng)用。其中,混沌多項式(polynomial chaos,PC)方法將隨機變量表示為隨機正交多項式的加權(quán)求和,能夠以較低的計算代價獲得與MC法相當?shù)木萚12]。鑒于PC與直接法的優(yōu)勢,諸多研究將二者結(jié)合,設(shè)計了基于嵌入式[10-11]與非嵌入式[12-15]的問題轉(zhuǎn)化策略。Fisher等[10]首次基于嵌入式策略將隨機最優(yōu)控制問題轉(zhuǎn)化為高維確定性問題,并采用配點法對該問題進行求解。Jiang等[13]則將非嵌入式策略與hp偽譜法相結(jié)合,解決了火星再入魯棒軌跡優(yōu)化設(shè)計問題。Wang等[11]提出了一種基于嵌入式策略與凸優(yōu)化的魯棒軌跡優(yōu)化方法,并與基于高斯偽譜法的方法進行對比,證明了所提方法在精度相當?shù)那闆r下,計算效率顯著提高。楊奔等[15]則將非嵌入式策略與凸優(yōu)化結(jié)合,提升了氣動參數(shù)擾動下再入軌跡設(shè)計的可靠性。

基于PC和直接法的隨機問題求解技術(shù)雖然已被逐步用于魯棒軌跡優(yōu)化,但仍存在計算效率低、通用性較差以及難于處理多維不確定性等問題,因而限制了其在復雜約束及強非線性魯棒優(yōu)化問題中的應(yīng)用。此外,現(xiàn)有的公開成果亦存在未綜合考慮軌跡魯棒性和可靠性等不足。因此,本文提出了一種兼顧過程安全性與終端魯棒性的再入滑翔軌跡優(yōu)化設(shè)計方法,該方法主要包含兩部分關(guān)鍵技術(shù):①構(gòu)建基于非嵌入式混沌多項式和高斯求積策略的不確定性量化傳播模型,實現(xiàn)隨機問題的轉(zhuǎn)化;②設(shè)計基于序列凸優(yōu)化的魯棒軌跡快速優(yōu)化求解策略及算法。最終,以X-33再入滑翔飛行為例進行仿真,驗證了方法的有效性與快速收斂性。

1 軌跡優(yōu)化問題描述

1.1 確定性再入軌跡優(yōu)化建模

假設(shè)地球為旋轉(zhuǎn)圓球,建立以無量綱能量e為自變量的三自由度再入滑翔無量綱運動方程[4]

(1)

式中,自變量為無量綱能量e=1/r-V2/2,無量綱狀態(tài)量X=[r,θ,φ,γ,ψ]T;r為無量綱地心距,V為無量綱速度,θ為地心經(jīng)度,φ為地心緯度,γ為當?shù)厮俣葍A角,ψ為航跡偏航角,υ為傾側(cè)角,控制量取為攻角α和傾側(cè)角υ,即u=[α,υ]T。

無量綱升力和阻力加速度L,D計算公式為

(2)

式中:m0為常值質(zhì)量;Sref參考面積;Re為地球半徑;升力及阻力系數(shù)CL,CD均為攻角α和馬赫數(shù)Ma的函數(shù);大氣密度ρ采用指數(shù)模型ρ=ρ0e-βh進行計算;h為飛行高度。

攻角控制指令由預設(shè)剖面α=α(Ma)直接給出,故待設(shè)計控制量僅為傾側(cè)角υ,且需要滿足以下約束

υmin≤|υ|≤υmax

(3)

軌跡優(yōu)化的邊界條件包括初始狀態(tài)約束和終端狀態(tài)約束

(4)

式中:X0與Xf由具體任務(wù)給定;e0與ef分別為初始狀態(tài)和終端狀態(tài)對應(yīng)的能量。

(5)

考慮到約束施加的便捷性,(5)式的過程約束可轉(zhuǎn)化為地心距約束[4]

(6)

式中,lQ(e),lq(e),ln(e)可以通過反解(5)式得到。

定義L(e)=max{lQ(e),lq(e),ln(e)},(6)式可以表示為

r(e)≥L(e)

(7)

考慮飛行時間最短問題,即以飛行時間作為目標函數(shù),則有

(8)

式中,τ為無量綱時間。

綜合(1)、(3)、(4)、(7)、(8)式,即可構(gòu)建以飛行時間最短作為優(yōu)化目標的確定性再入軌跡優(yōu)化問題P0。

1.2 不確定性因素建模

再入滑翔飛行所面臨的不確定性主要包括再入初始狀態(tài),以及大氣密度、氣動力系數(shù)等動力學參數(shù)的不確定性,具體可建模描述為

(9)

且有

(10)

(11)

1.3 魯棒再入軌跡優(yōu)化問題描述

在問題P0中引入?yún)?shù)不確定性,則其狀態(tài)量、目標函數(shù)、約束條件均為關(guān)于隨機變量S的隨機函數(shù),為盡可能降低隨機干擾對尋優(yōu)結(jié)果的影響,即當存在不確定性時,預期軌跡離散度在統(tǒng)計意義上盡可能小,需要在目標函數(shù)、約束函數(shù)中引入隨機量的統(tǒng)計矩,故可將問題P0轉(zhuǎn)化為P1進行描述

(12)

式中:μ(·)和σ(·)分別表示隨機函數(shù)(·)的均值和標準差;X(e,S)表示考慮隨機干擾的狀態(tài)量;kJ,kC為非負權(quán)重系數(shù)。

問題P1中,狀態(tài)、約束、目標函數(shù)均為非線性隨機函數(shù),傳統(tǒng)的確定性優(yōu)化算法無法直接對其進行求解,一般需要引入不確定性量化傳播技術(shù)計算相關(guān)統(tǒng)計量并處理隨機動力學方程,進而采用高效優(yōu)化算法求解轉(zhuǎn)化所得到的高維確定性問題。

2 再入滑翔隨機最優(yōu)控制問題轉(zhuǎn)化

本節(jié)將非嵌入式混沌多項式與高斯求解策略相結(jié)合,建立了不確定性量化傳播模型,從而將再入滑翔魯棒軌跡優(yōu)化問題轉(zhuǎn)化為高維狀態(tài)空間中的等價確定性優(yōu)化問題。以下給出算法原理及實現(xiàn)步驟。

2.1 非嵌入式混沌多項式算法原理

非嵌入式混沌多項式方法的核心思路在于對彈道模型的轉(zhuǎn)換和樣本點的選取,即通過隨機變量取值空間上的積分點代替隨機樣本點,極大減少了模型的求解次數(shù)[16]。

首先,考慮形如(12)式的隨機微分方程,其解X=X(e,S)采用以下混沌多項式進行逼近[12]

(13)

多項式展開的項數(shù)P+1可由(14)式確定

(14)

式中:p為多項式的階數(shù),其取值決定了多項式逼近的程度;d為隨機向量S的維數(shù)。

(15)

式中:Q=nd為隨機空間的配點總數(shù);n為S中單維隨機變量的配點數(shù);τm為配點sm對應(yīng)的權(quán)重系數(shù);X(e,sm)可在配點sm處積分確定性運動方程得到。

類似地,考慮(12)式中的目標函數(shù)J(X,e,S)、約束函數(shù)L(X,e,S)一般亦為狀態(tài)X(e,S)的函數(shù),故其混沌多項式展開形式為

式中:Jm=J(X(e,sm),e,sm);Lm=L(X(e,sm),e,sm)可由確定性狀態(tài)X(e,sm)求解得到。

2.2 再入滑翔隨機問題轉(zhuǎn)化

(20)

類似地,目標函數(shù)J(X,e,S)與約束函數(shù)L(X,e,S)的統(tǒng)計量亦可表征為相似形式。

結(jié)合2.1節(jié)的非嵌入式混沌多項式算法和上述統(tǒng)計量計算公式,可將問題P1轉(zhuǎn)化為狀態(tài)擴展的確定性軌跡優(yōu)化問題P2,即

(21)

式中,Xm=[rm,θm,φm,γm,ψm]T。該問題為一高維確定性優(yōu)化問題,其狀態(tài)維數(shù)擴展為P1問題的Q倍。因此,下文采用凸優(yōu)化方法對該高維確定性問題進行求解。

3 基于序列凸優(yōu)化的軌跡優(yōu)化求解

3.1 凸化處理

為避免抖振問題[4]出現(xiàn),此處定義新的控制變量u=[u1,u2]T,有

u1=cosυ,u2=sinυ

(22)

且滿足約束

(23)

(24)

定義Dm=D(rm,e,sm),Lm=L(rm,e,sm),則有

(25)

(26)

式中:fΩ(Xm,sm)為與地球自轉(zhuǎn)相關(guān)的小量,且f0(Xm,sm)關(guān)于Xm是非線性的。以下采用逐次線性化方法,對相關(guān)非線性方程進行凸化處理。

1) 動力學方程凸化

(27)

且有

(28)

(29)

式中,δ表示信賴域半徑,m=1,2,…,Q。

2) 目標函數(shù)凸化

為使班會的內(nèi)容更為飽滿，形式更為生動，有效實現(xiàn)主題，我嘗試著將寓理式微型班會與體驗式微型班會的結(jié)構(gòu)模式相結(jié)合。下面我將從主題的導入、展開及深化三個方面來談?wù)勥@堂課的環(huán)節(jié)設(shè)計：

(30)

定義a0=2m/(ReSrefCD(α,Ma,sm)V3),則

(31)

可得凸函數(shù)

(32)

此外,目標函數(shù)第二項為凸二次函數(shù),故其整體即為凸函數(shù)。

3) 控制約束凸化

將(21)式中傾側(cè)角幅值約束轉(zhuǎn)化為對控制量u1的約束,即

cos(υmax)≤u1≤cos(υmin)

(33)

此外,利用松弛技術(shù)對(23)式的非凸附加約束進行凸化處理

(34)

4) 過程約束凸化

(35)

且有

(36)

至此,問題P2即轉(zhuǎn)化為等價的高維確定性軌跡凸優(yōu)化問題P3,限于篇幅,此處不再重復列出。

3.2 離散化處理

為了應(yīng)用數(shù)值方法對凸化后的問題P3進行求解,需要對狀態(tài)量和控制量進行離散化處理。將自變量變化域[e0,ef]等間距離散為N個間隔,自變量離散為[e0,e1,…,eN],狀態(tài)量Xm離散為[Xm(0),Xm(1),…,Xm(N)],控制量離散為[u0,u1,…,uN]。

采用梯形法則對(27)式進行離散,可得

(37)

類似地,對問題P3的其他約束及目標函數(shù)進行離散化處理,得到

(38)

最終可以得到離散化的問題P4,可通過經(jīng)典內(nèi)點法求解,其優(yōu)化變量包括每個離散點處的Q組狀態(tài)量xm(h)和控制量uh,待優(yōu)化變量數(shù)目為(5Q+2)N。

3.3 序列凸優(yōu)化求解策略

問題P4將動力學約束線性化處理,故必然存在模型偏差,序列凸優(yōu)化算法通過迭代求解P4問題對初始問題進行逼近,即利用外層迭代更新求解變量序列,內(nèi)層求解問題P4,逐步逼近原問題P2的解,以保證凸優(yōu)化解的收斂性和精度,即令

(39)

以相鄰2次迭代的凸優(yōu)化解對應(yīng)的狀態(tài)量最大偏差作為收斂準則,即

(40)

式中,收斂誤差限定義為ε=[εr,εθ,εφ,εγ,εψ]T。

以下給出基于序列凸優(yōu)化求解原最優(yōu)控制問題P2的步驟和策略:

1) 根據(jù)隨機向量S服從的分布確定配點sm和對應(yīng)積分權(quán)重τm,其中m=1,2,…,Q;

5) 判斷(40)式是否成立,若成立,則算法結(jié)束;否則,令k=k+1,轉(zhuǎn)入步驟3)。

4 仿真分析

4.1 仿真條件設(shè)置

表1 邊界條件取值

(41)

所有仿真均在搭載Intel Core i7-8700 3.20 GHz Intel處理器的臺式機完成,仿真環(huán)境為MATLAB 2016b平臺?；谧孕虚_發(fā)的代碼進行不確定性量化傳播,同時基于CVX工具包進行軌跡優(yōu)化算法開發(fā),并調(diào)用SDPT3求解器求解凸優(yōu)化子問題P4。

仿真中采用固定數(shù)目等距離散點,取值N=300,序列凸優(yōu)化算法的信賴域半徑和收斂誤差限設(shè)置為:

(42)

4.2 確定性軌跡優(yōu)化及不確定性量化傳播仿真分析

首先,利用序列凸優(yōu)化算法對問題P0進行仿真求解,得到標稱條件下的再入滑翔軌跡優(yōu)化結(jié)果,如圖1所示。圖1a)給出了速度-高度初始參考及優(yōu)化結(jié)果曲線,可以看出,優(yōu)化軌跡連續(xù)光滑,滿足過程約束;由于性能指標為時間最短,故高度曲線始終處于下邊界附近。由圖1b)可知,優(yōu)化結(jié)果滿足過程約束要求,熱流密度和過載曲線存在一定時間區(qū)間內(nèi)達到約束峰值的情況。

圖1 標稱情況下滑翔軌跡優(yōu)化結(jié)果

以下基于標稱條件下的優(yōu)化結(jié)果分別開展蒙特卡羅打靶法(MC)與混沌多項式方法(PC)的不確定性傳播量化仿真,以驗證所提不確定量化方法的有效性。設(shè)定MC打靶次數(shù)為5 000次,所得結(jié)果如圖2所示。由圖2可知,再入飛行存在超出熱流約束邊界的情況,統(tǒng)計結(jié)果顯示超出過程約束邊界的概率達到了78.5%,可見,標稱情況下的優(yōu)化結(jié)果易受到參數(shù)不確定性的影響,從而出現(xiàn)終端精度下降和過程約束超限等問題,在實際飛行過程中,上述問題會顯著增大任務(wù)失敗的風險,因此,需要進一步開展考慮不確定性的魯棒軌跡優(yōu)化研究,以提高優(yōu)化軌跡的抗干擾能力和可靠性。

圖2 標稱情況下優(yōu)化控制量打靶結(jié)果

以下對MC與PC方法的不確定性量化傳播性能進行比較。PC方法仿真時,設(shè)置一維配點數(shù)n=3,多項式的階數(shù)p=2,隨機變量維數(shù)d=3。

圖3～4給出了MC與PC方法的部分狀態(tài)量均值與標準差的對比曲線。由結(jié)果可知,2種方法的均值結(jié)果幾乎完全一致,而PC方法計算的標準差有一定偏差但在可接受范圍內(nèi),如經(jīng)緯度標準差的相對誤差僅約3%,即PC方法可替代MC方法進行不確定性量化傳播。此外,在計算結(jié)果相當?shù)那闆r下,PC方法只需對動力學方程進行27次積分,而本文的MC方法打靶次數(shù)為5 000次,充分體現(xiàn)了PC在計算效率方面的顯著優(yōu)勢。

圖3 不確定性量化傳播的經(jīng)度均值對比

圖4 不確定性量化傳播的經(jīng)度標準差對比

4.3 考慮參數(shù)不確定的魯棒軌跡優(yōu)化仿真分析

為驗證所設(shè)計算法的可行性,本節(jié)開展參數(shù)不確定條件下的再入滑翔魯棒軌跡優(yōu)化仿真,通過調(diào)整權(quán)重系數(shù)kC,kJ以滿足不同的可靠性和魯棒性需求[13]。采用“DO”、“RO1”與“RO2”分別代表確定性、可靠及魯棒優(yōu)化結(jié)果,設(shè)置RO1參數(shù)kC=3,kJ=0;RO2參數(shù)kC=3,kJ=1,仿真結(jié)果如圖5所示。

圖5 DO、RO1與RO2的優(yōu)化結(jié)果對比

圖5對比了3種條件下優(yōu)化所得H-V與傾側(cè)角指令曲線。由圖5a)可知,RO1與RO2對應(yīng)的高度曲線相較于DO更遠離過程約束邊界,即優(yōu)化軌跡的可靠性得到明顯提升;與RO1相比,RO2的優(yōu)化軌跡在中段出現(xiàn)2次拉起,表明這種軌跡形式有利于降低軌跡對不確定因素的敏感度,提升軌跡魯棒性。由圖5b)可知,相較于DO,RO1與RO2初始下降段對應(yīng)的傾側(cè)角幅值基本維持在0°附近,大小基本一致,表明RO1與RO2對熱流峰值的可靠性相當;RO1與RO2繼續(xù)保持較小的傾側(cè)角幅值以保持更高的飛行高度,提升過程可靠性。

利用圖5b)中的傾側(cè)角優(yōu)化結(jié)果進行5 000次蒙特卡羅打靶仿真,以驗證魯棒優(yōu)化算法的有效性,所得部分參數(shù)的統(tǒng)計結(jié)果如表2與圖6所示。其中,Δμ(·)為均值偏差,Pe為超出過程約束邊界的概率,即約束違反率。

圖6 RO2優(yōu)化控制量打靶結(jié)果

表2 3組仿真打靶數(shù)據(jù)對比

由表2可知,RO1與RO2的終端經(jīng)緯度統(tǒng)計精度基本一致,且相較于DO有明顯改善,而由于考慮了終端性能的魯棒性,RO2對應(yīng)的終端經(jīng)緯度在三者中最小。由圖6a)～6b)及表2可知,相較于DO,RO1與RO2的熱流、動壓、過載均遠離約束邊界,約束超出的概率由78.5%分別降低到0.32%與0.56%,可見兩者的可靠性均有所提升。與RO1結(jié)果相比,RO2在保證軌跡可靠性的同時,對不確定性的敏感度最低,具有更好的魯棒性。此外,相較于DO結(jié)果,RO1與RO2對應(yīng)的飛行時間均值較大,這是由于不確定性的引入與可靠性約束的施加會造成優(yōu)化求解的可行域減小,迫使軌跡高度提升導致飛行時間增大;而相較于RO1,RO2的飛行時間均值更大,這是由于RO2的目標函數(shù)中綜合考慮了飛行時間的均值與終端經(jīng)緯度方差,而RO1僅考慮了飛行時間均值。可見,軌跡的可靠性與魯棒性的提升是以犧牲一定的最優(yōu)性為代價的,需要根據(jù)設(shè)計需求合理調(diào)整權(quán)重系數(shù),以權(quán)衡軌跡的可靠性、魯棒性與最優(yōu)性。

4.4 2種魯棒軌跡優(yōu)化算法仿真對比

為驗證所設(shè)計算法的計算效率優(yōu)勢,將本文算法(convex optimization-polynomial chaos,CO-PC)與現(xiàn)有文獻[13]中基于偽譜法與PC的典型魯棒軌跡優(yōu)化方法(guess pseudospectral-polynomial chaos,GP-PC)進行對比。設(shè)置權(quán)重系數(shù)kC=kJ=3,其余仿真條件同上,所得優(yōu)化結(jié)果如圖7所示。

圖7 GP-PC與CO-PC的優(yōu)化結(jié)果對比

隨后,利用圖7b)中的傾側(cè)角優(yōu)化結(jié)果進行5 000次蒙特卡羅打靶仿真,以進一步對比2種魯棒優(yōu)化算法的性能,所得部分參數(shù)的統(tǒng)計結(jié)果如表3所示。其中“DO”為確定性優(yōu)化結(jié)果。

表3 2種方法優(yōu)化結(jié)果的打靶數(shù)據(jù)對比

由圖7a)～7b)可知,2種方法的魯棒軌跡優(yōu)化結(jié)果基本一致,且對應(yīng)的高度曲線相較于DO更遠離過程約束邊界;相較于GP-PC,CO-PC的高度曲線前段更高,且飛行全程存在2次高度拉起過程。由表3可知,兩者的約束違反概率及終端經(jīng)緯度均值精度基本一致,而相較于GP-PC,CO-PC的終端經(jīng)緯度標準差以及時間均值則更小,表明CO-PC的優(yōu)化結(jié)果魯棒性更強,目標函數(shù)更優(yōu)。

表4則進一步對比了2種優(yōu)化算法的計算耗時及待優(yōu)化變量數(shù)?？梢钥闯?CO-PC方法在優(yōu)化變量數(shù)及迭代步數(shù)更多的情況下,優(yōu)化耗時僅為GP-PC方法的10%,表明本文所提出的CO-PC方法能夠有效處理狀態(tài)維數(shù)擴展的魯棒軌跡優(yōu)化問題,顯著提升該問題求解的效率。

表4 GP-PC與CO-PC的數(shù)值計算性能對比

5 結(jié) 論

針對存在參數(shù)不確定的再入滑翔軌跡設(shè)計問題,本文研究了一種基于非嵌入式混沌多項式與凸優(yōu)化相結(jié)合的魯棒軌跡優(yōu)化算法,理論分析與仿真結(jié)果表明:

1) 與確定性軌跡優(yōu)化和現(xiàn)有典型魯棒軌跡優(yōu)化算法相比,本文方法能夠顯著提升再入滑翔軌跡的魯棒性與可靠性,降低軌跡對參數(shù)不確定性的敏感度,且具有相當?shù)木群透叩挠嬎阈蕛?yōu)勢。

2) 軌跡魯棒性與可靠性的提升以犧牲一定的最優(yōu)性為代價,通過調(diào)整權(quán)重系數(shù),能夠權(quán)衡三者關(guān)系,以綜合滿足軌跡的可靠性、魯棒性和最優(yōu)性需求。

3) 所構(gòu)建的基于高斯求積與非嵌入式混沌多項式的不確定量化傳播算法,可實現(xiàn)隨機問題的有效轉(zhuǎn)化,且計算精度與MC方法相當情況下,具有顯著的計算效率優(yōu)勢。

4) 本文僅考慮了三維隨機變量的不確定問題,若進一步增加隨機變量維數(shù),會導致擴展狀態(tài)維數(shù)呈指數(shù)型增加,難于求解。因此,在高維不確定條件下,如何實現(xiàn)魯棒軌跡的快速優(yōu)化求解可作為后續(xù)的研究方向。