基于混沌多項式的再入滑翔魯棒軌跡凸優化

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

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

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

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

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

1 軌跡優化問題描述

1.1 確定性再入軌跡優化建模

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

(1)

式中,自變量為無量綱能量e=1/r-V2/2,無量綱狀態量X=[r,θ,φ,γ,ψ]T;r為無量綱地心距,V為無量綱速度,θ為地心經度,φ為地心緯度,γ為當地速度傾角,ψ為航跡偏航角,υ為傾側角,控制量取為攻角α和傾側角υ,即u=[α,υ]T。

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

(2)

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

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

υmin≤|υ|≤υmax

(3)

軌跡優化的邊界條件包括初始狀態約束和終端狀態約束

(4)

式中:X0與Xf由具體任務給定;e0與ef分別為初始狀態和終端狀態對應的能量。

(5)

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

(6)

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

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

r(e)≥L(e)

(7)

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

(8)

式中,τ為無量綱時間。

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

1.2 不確定性因素建模

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

(9)

且有

(10)

(11)

1.3 魯棒再入軌跡優化問題描述

在問題P0中引入參數不確定性,則其狀態量、目標函數、約束條件均為關于隨機變量S的隨機函數,為盡可能降低隨機干擾對尋優結果的影響,即當存在不確定性時,預期軌跡離散度在統計意義上盡可能小,需要在目標函數、約束函數中引入隨機量的統計矩,故可將問題P0轉化為P1進行描述

(12)

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

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

2 再入滑翔隨機最優控制問題轉化

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

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

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

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

(13)

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

(14)

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

(15)

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

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

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

2.2 再入滑翔隨機問題轉化

(20)

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

結合2.1節的非嵌入式混沌多項式算法和上述統計量計算公式,可將問題P1轉化為狀態擴展的確定性軌跡優化問題P2,即

(21)

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

3 基于序列凸優化的軌跡優化求解

3.1 凸化處理

為避免抖振問題[4]出現,此處定義新的控制變量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)為與地球自轉相關的小量,且f0(Xm,sm)關于Xm是非線性的。以下采用逐次線性化方法,對相關非線性方程進行凸化處理。

1) 動力學方程凸化

(27)

且有

(28)

(29)

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

2) 目標函數凸化

為使班會的內容更為飽滿，形式更為生動，有效實現主題，我嘗試著將寓理式微型班會與體驗式微型班會的結構模式相結合。下面我將從主題的導入、展開及深化三個方面來談談這堂課的環節設計：

(30)

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

(31)

可得凸函數

(32)

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

3) 控制約束凸化

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

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

(33)

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

(34)

4) 過程約束凸化

(35)

且有

(36)

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

3.2 離散化處理

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

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

(37)

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

(38)

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

3.3 序列凸優化求解策略

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

(39)

以相鄰2次迭代的凸優化解對應的狀態量最大偏差作為收斂準則,即

(40)

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

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

1) 根據隨機向量S服從的分布確定配點sm和對應積分權重τm,其中m=1,2,…,Q;

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

4 仿真分析

4.1 仿真條件設置

表1 邊界條件取值

(41)

所有仿真均在搭載Intel Core i7-8700 3.20 GHz Intel處理器的臺式機完成,仿真環境為MATLAB 2016b平臺。基于自行開發的代碼進行不確定性量化傳播,同時基于CVX工具包進行軌跡優化算法開發,并調用SDPT3求解器求解凸優化子問題P4。

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

(42)

4.2 確定性軌跡優化及不確定性量化傳播仿真分析

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

圖1 標稱情況下滑翔軌跡優化結果

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

圖2 標稱情況下優化控制量打靶結果

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

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

圖3 不確定性量化傳播的經度均值對比

圖4 不確定性量化傳播的經度標準差對比

4.3 考慮參數不確定的魯棒軌跡優化仿真分析

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

圖5 DO、RO1與RO2的優化結果對比

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

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

圖6 RO2優化控制量打靶結果

表2 3組仿真打靶數據對比

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

4.4 2種魯棒軌跡優化算法仿真對比

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

圖7 GP-PC與CO-PC的優化結果對比

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

表3 2種方法優化結果的打靶數據對比

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

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

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

5 結 論

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

1) 與確定性軌跡優化和現有典型魯棒軌跡優化算法相比,本文方法能夠顯著提升再入滑翔軌跡的魯棒性與可靠性,降低軌跡對參數不確定性的敏感度,且具有相當的精度和更高的計算效率優勢。

2) 軌跡魯棒性與可靠性的提升以犧牲一定的最優性為代價,通過調整權重系數,能夠權衡三者關系,以綜合滿足軌跡的可靠性、魯棒性和最優性需求。

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

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