基于隨機響應面與混沌多項式的魯棒軌跡優化

王培臣, 閆循良,*, 王 寬, 鄭 雄

(1. 西北工業大學航天學院, 陜西 西安 710072; 2. 中國運載火箭技術研究院, 北京 100076)

0 引 言

作為飛行器設計的關鍵技術之一,軌跡優化逐漸成為高超聲速等航空航天領域研究的熱點問題[1-2]。然而,現有的軌跡優化研究通常僅注重標稱情況下的軌跡性能優化提升,沒有考慮真實情況下軌跡的魯棒性和可靠性[3-6]。在諸如高超聲速再入等實際飛行任務中,初始狀態、動力學及環境參數等均存在諸多不確定性,勢必使得飛行器沿標稱軌跡飛行時產生較大的過程及終端偏差,進而增大約束違反的風險,降低飛行可靠性和終端任務精度。因此,為了提高優化軌跡的可靠性和抗干擾能力,部分學者逐步將研究重心聚焦于考慮不確定性的魯棒軌跡優化[7-13],即在軌跡設計階段事先考慮不確定性的影響,從而降低由軌跡規劃迭代耗時和軌跡偏差修正帶來的制導控制系統負擔。

求解魯棒軌跡優化問題的典型策略為:先將含有隨機微分方程的隨機最優控制問題轉化為包含擴展確定性微分方程的確定最優控制問題,進而采用高效的軌跡優化算法對該高維問題進行數值求解[8],即包含不確定性量化傳播和確定性軌跡優化兩部分關鍵技術。典型的不確定性量化傳播理論和方法包括蒙特卡羅(Monte Carlo,MC)方法[14]、線性協方差分析法[15-16]、無跡變換法[17]及混沌多項式展開(polynomial chaos expansion,PCE)方法[18-19]等。MC方法因簡單易行而被廣泛應用于不確定性量化及傳播,但大量的數值模擬會導致統計結果收斂速度較慢[14];線性協方差分析法可對線性高斯過程模型中的不確定性進行量化傳播,且具有較高的計算效率,但其不適用于高度非線性的軌跡優化問題[15];相比之下,PCE方法將隨機變量表示為隨機正交多項式的加權求和,可對具有任意分布類型的隨機變量實現較為精確的逼近,同時能夠以較低的計算代價達到與MC方法相當的精度,近年來在火星進入軌跡不確定性量化[18-21]、流體力學不確定性分析[22-25]中得到了廣泛應用。PCE方法與系統模型結合的方式有兩種,即嵌入式[17]與非嵌入式[19,23-25],嵌入式混沌多項式(intrusive polynomial chaos, IPC)將隨機變量的混沌多項式展開代入動力學微分方程中,并利用伽遼金投影法將隨機微分方程轉換為一組以混沌多項式系數為狀態量的高維確定性微分方程組,最終通過求解該微分方程組完成不確定性傳播[21-22],具有較高的計算效率;非IPC(non-IPC, NIPC)則將原系統模型視為“黑箱”,進而根據模型的輸入輸出逼近混沌多項式系數,無需針對具體問題進行推導建模,因此能夠處理一般非線性系統的不確定性量化傳播問題,具有較好的通用性及工程應用價值。

而對于另一關鍵技術,即確定性軌跡優化而言,常見的方法主要包括間接法和直接法。間接法因對解的初值猜測較敏感、需要對最優性必要條件進行推導且無法處理復雜約束等原因,難以在復雜約束軌跡優化問題中得到廣泛應用。近年來,以偽譜法[2]、凸優化[4-6]等為代表的直接法已被廣泛應用于高超聲速等航空航天領域的軌跡優化任務,并展現出了良好的數值計算優勢和普適性等特點。其中,凸優化方法因具有多項式復雜度、收斂速度快,且其局部最優解理論上即為全局最優解等特點,已逐步被用于求解復雜約束下的軌跡優化問題以及考慮不確定性的高維狀態擴展優化問題。

鑒于PCE方法與直接法的優勢,諸多研究將二者進行有機結合,設計了多種基于PCE方法與直接法的魯棒軌跡優化算法和求解策略。Fisher等[8]首次提出了數值求解隨機最優控制問題的框架,該框架利用IPC將隨機最優控制問題轉化為高維確定性最優控制問題,并采用配點法對該問題進行求解。Li等[9]將NIPC與高斯偽譜法相結合,用于解決氣動系數不確定的高超聲速飛機的爬升魯棒軌跡優化問題。Jiang等[10]將NIPC與偽譜法相結合,通過調用OpenCossan工具箱進行不確定量化傳播,并采用并行計算技術提升高維確定性優化問題的求解效率,解決了初始狀態及模型參數均不確定的火星再入魯棒軌跡優化設計問題。Wang等[11]提出一種將IPC與凸優化結合的魯棒軌跡優化方法,并與基于IPC與高斯偽譜法的方法進行對比,證明了所提方法在計算精度相當的情況下,計算效率顯著提高。楊奔等[12]基于最小二乘法建立了NIPC與凸優化相結合的魯棒軌跡優化模型,提升了氣動參數不確定條件下軌跡的可靠性。

基于PCE和直接法的隨機問題求解技術雖然已被逐步用于魯棒軌跡優化,但仍存在計算效率低、通用性較差以及難于處理多維不確定性等問題,這些問題進而限制了其在復雜約束及強非線性魯棒優化問題(如航空航天領域魯棒軌跡優化)中的應用。首先,以高斯偽譜法為代表的傳統直接法雖然解決了部分魯棒軌跡優化問題,但其仍面臨計算效率低、難于快速求解多維不確定問題等局限,而凸優化技術與PCE的結合仍然不夠完善。其次,基于IPC的優化方法雖然具有較高的計算效率,但需要針對不同問題進行推導建模和代碼改寫,由于推導過程存在模型簡化近似,只適用于處理多項式非線性問題[9];基于NIPC的優化方法雖然能夠避免上述問題,但PCE系數求解大都基于求積、采樣的方法,例如高斯求積(Gauss quadrature,GQ)方法[9-10]、稀疏網格求積法[13],而隨著系統隨機變量維數及非線性程度增大,這類方法的狀態擴展維數會呈指數增長,使得狀態離散節點數量顯著增加,從而導致解決大規模NLP問題時非常困難,計算成本高昂且效率低下,尤其是當不確定因素維數n>3時,大量的內存消耗和計算成本會嚴重阻礙優化過程,甚至導致優化失敗。

值得注意的是,除去上述基于求積的方法,NIPC系數也可采用隨機響應面法(stochastic response surface method, SRSM)[23]進行計算。該方法基于代理模型思想,通過擬合系統輸出與輸入之間的關系獲得混沌多項式模型,并結合典型抽樣方法進行少量采樣即可實現多項式系數計算,可以顯著降低模型轉化維數和規模,因此具有較高的精度和計算效率。目前,SRSM已被廣泛應用于流體力學不確定性分析中[24-26],但尚未在魯棒軌跡優化設計中廣泛應用。

因此,本文針對航空航天領域魯棒軌跡優化算法存在的問題,提出一種高效的魯棒軌跡優化設計方法。通過構造基于NIPC與SRSM的不確定性量化傳播模型,實現了不確定性隨機問題的有效轉化;隨后,將該模型與凸優化方法相結合,設計了一種基于序列凸優化算法的求解策略,以實現對該高維確定性問題的快速求解。最終,以高超聲速滑翔再入為例進行仿真,驗證了本文方法的有效性。

1 魯棒軌跡優化問題描述

一般的軌跡優化問題可以歸結為如下確定性最優控制問題PD:

(1)

當系統參數向量w及狀態初值x(t0)存在不確定性時,目標函數J、動力學函數f及約束函數均變為隨機函數。因此,原確定性優化問題PD轉變為隨機優化問題。為了降低優化軌跡對不確定性的敏感度,即當存在不確定性時,預期的軌跡離散度在統計意義上盡可能小,通常需要在目標函數、約束函數中引入隨機量的統計矩。因此,一般魯棒軌跡優化問題PR可描述為

(2)

式中:下標μ和σ分別代表均值和標準差;kJ、kg為待設計權重系數;x0μ與x0σ分別為初始狀態的均值和標準差;εf為狀態終值標準差的上界。當采用上述模型進行軌跡優化時,需要引入不確定性量化傳播技術計算相關統計量并處理隨機動力學方程。

需要注意的是:①εf需要基于先驗統計信息確定,若取值不當可能造成問題無解,為了提升優化算法的收斂性,可將xσ(tf)以罰函數的形式引入目標函數中;② 上述魯棒軌跡優化模型同時考慮了軌跡的魯棒性及可靠性,其中,在目標函數及端點約束中引入均值、標準差項可提升軌跡的魯棒性,而通過引入考慮統計矩的過程約束模型,則可以保證以一定概率滿足過程約束從而提升軌跡可靠性;③ 通過設計權重系數kJ、kg可實現軌跡的魯棒性及可靠性之間的權衡,具體討論見后文仿真部分。

下文將詳細給出考慮初始狀態x(t0)及系統參數向量w不確定性的魯棒軌跡優化算法。首先,采用基于響應面法的PCE技術將原始魯棒軌跡優化問題PR轉化為狀態擴展的高維確定性優化問題;隨后,對該問題進行凸化、離散處理,并采用序列凸優化算法進行求解。

2 不確定量化傳播建模與問題轉化

本節將NIPC與SRSM相結合,同時引入拉丁超立方采樣策略,構建了不確定性量化傳播模型,從而將魯棒軌跡優化問題轉化為高維狀態空間中的等價確定性優化問題。以下給出算法原理及實現步驟。

2.1 隨機狀態方程轉化

假設隨機狀態向量x為n維,系統參數向量w為s維,則每個狀態變量xi(t)及參數wk可展開為PCE多項式[20]:

(3)

(4)

式中:Δ=[Δ1,Δ2,…,Δd]為d維隨機向量;Φj(Δ)為正交多項式函數,可通過與Δ分布類型相對應的一元正交基函數φj(Δl)(l=1,2,…,d)求張量積得到;xij(t)和wkj為PCE系數,展開項數P+1=[(p+d)!]/(p!d!),即由正交基函數φj(Δl)的階數p和隨機變量維數d決定;根據Askey法則[20],不同隨機變量Δl的分布類型對應不同的正交基函數{φj(Δl)},具體對應形式如表1所示。

表1 正交基函數與隨機變量分布類型對應關系Table 1 Correspondence relationship between distribution type of orthogonal basis function and random variables

PC(polynomial chaos,PC)系數wkj及xij(t0)可根據wk及xi(t0)的分布確定,而隨機狀態xi(t)的PC系數xij(t)未知,下面采用SRSM導出求解系數xij(t)的表達式。

將式(3)、式(4)代入式(2)的微分方程中,可得

(5)

采用拉丁超立方[24]進行ns次隨機采樣,并將隨機樣本點{Δ1,Δ2,…,Δns}代入方程(5),可得微分方程組:

(6)

式中:X與W分別表示隨機狀態x與參數w的PC系數向量,且有

X=[x10(t),x11(t),…,x1P(t),x20(t),x21(t),…,x2P(t),…,
xn0(t),xn1(t),…,xnP(t)]T

(7)

W=[w10,w11,…,w1P,w20,w21,…,w2P,…,ws0,ws1,…,wsP]T

(8)

(9)

設Ω(Δm)=[Φ0(Δm),Φ1(Δm),…,ΦP(Δm)],(m=1, 2,…,ns),定義

Θ(Δm)=In?Ω(Δm)

(10)

則Ψ=[Θ(Δ1)T,Θ(Δ2)T,…,Θ(Δns)T]T∈Rn·ns×n(P+1)。

若ns≥P+1,則可通過最小二乘回歸法求解式的微分方程組,即有

(11)

式中:矩陣H=(ΨTΨ)-1ΨT∈Rn(P+1)×n·ns。為保證該模型的逼近精度,本文統一取采樣點ns=2(P+1)[22],后文不再贅述。式(11)即為關于擴展系數X與控制量u的高維狀態方程。

2.2 過程約束及目標函數的PC展開

由于非線性約束函數g(x(t),u(t),t,w)與目標函數J(x(t),u(t),t,w)均為隨機狀態變量x(t)的函數,故也可采用類似方法進行混沌多項式展開,具體如下:

(12)

(13)

式中:gj(x(t),u(t),t)與Jj(x(t),u(t),t)為對應函數的PC多項式系數。將隨機樣本點{Δ1,Δ2,…,Δns}代入方程(12)、方程(13)可得

(14)

(15)

(16)

(17)

(18)

(19)

類似地,有

(20)

(21)

(22)

2.3 統計量計算

根據文獻[15],可將式(2)中隨機狀態x、過程約束函數g(x(t),u(t),t,w)及目標函數J(x(t),u(t),t,w)的均值和方差表示為PC多項式系數的函數,即

(23)

(24)

(25)

式中:〈·〉表示內積運算;gj,Jj(j=1,2,…,P)可通過式(21)、式(22)表示為X、u的非線性函數,進而可通過以上公式計算狀態、過程約束與目標函數的統計量。

至此,已將魯棒軌跡優化問題PR轉化為關于擴展狀態X及控制量u的高維確定性優化問題PED,進而可利用凸優化算法進行處理求解,限于篇幅,問題PED在此不再贅述。

3 基于凸優化的魯棒優化求解策略

考慮到復雜軌跡優化問題大多是非凸的,如何將非凸問題進行凸化處理是應用凸優化技術求解最優軌跡的關鍵,并吸引了廣大學者對此開展了深入研究,相關成果亦在航空航天領域最優控制及軌跡優化問題上得到了較為廣泛的應用[4-6]。

本節引入逐次線性化[27]等技術對第2節中的高維確定性優化問題PED進行凸化處理,并將該問題離散處理為參數凸優化問題,進而建立相應的序列凸優化迭代算法,對該問題進行求解。序列凸優化方法通過外環迭代更新參考軌跡,內環求解凸問題逐步逼近原問題的最優解[27]。以下給出凸化、離散處理過程以及迭代算法的原理和實現步驟。

3.1 凸化處理

3.1.1 狀態方程凸化

分析式(11)可知,H為常值矩陣,F(X,u,t,W)為關于擴展狀態X的非線性函數,因此該式為非線性微分方程。下面采用逐次線性化方法對該式進行凸化處理。

將式(11)在參考軌跡{Xk,uk}附近進行一階泰勒展開,可得

(26)

定義:

(27)

式中:?F/?X|(Xk,uk,t,W)與?F/?u|(Xk,uk,t,W)分別為參考軌跡上函數F(X,u,t,W)相對于狀態X和控制量u的偏導數。其中,參考軌跡{Xk,uk}可通過迭代算法上一步優化結果得到,第一步優化采用的參考軌跡初值{X0,u0}可以X(t0)作為初始狀態,以標稱條件下的確定性優化結果作為控制量u(t)積分式(11)得到。

為保證線性化的精度,增加信賴域約束:

(28)

式中:δX和δu表示信賴域半徑。

3.1.2 過程約束及目標函數凸化

類似地,考慮式(2)中的過程約束及目標函數均為關于擴展狀態X及控制量u的非線性函數,仍采用逐次線性化方法對其進行凸化處理。

在參考軌跡{Xk,uk}附近進行一階泰勒展開,可得:

gμ(X,u,t,W)+kggσ(X,u,t,W)=
GAX+GBu+GC

(29)

Jμ(X,u,t,W)+kJJσ(X,u,t,W)=
JAX+JBu+JC

(30)

定義:

(31)

(32)

3.1.3 狀態不等式約束凸化

考慮式(2)中的狀態不等式約束為關于擴展狀態X的非線性函數,其凸化處理方法[11]如下。

由式(23)可知,隨機狀態變量xi(t)的方差可表示為

(33)

Q=LTL

(34)

狀態變量xi的標準差可表示為

(35)

因此,式(2)中的狀態不等式約束可無損凸化為如下二階錐約束:

(36)

式中:ximax與ximin分別為第i個狀態量對應的上下邊界。

至此,問題PED轉化為以下等價的高維確定性軌跡凸優化問題PED-Convex:

(37)

需要說明的是,問題PED-Convex的待優化狀態為X,擴展為問題PR的(P+1)倍。

3.2 離散化處理

為了應用數值方法對問題PED-Convex進行求解,需要對狀態量和控制量進行離散化處理。基于梯形法則,將自變量變化域[t0,tf]等間距離散為N個間隔,自變量離散為{t0,t1,…,tN},狀態量離散為{X0,X1,…,XN},控制量離散為{u0,u1,…,uN}。

對式(37)的線性動力學方程進行離散化,可得:

(38)

類似地,對問題PED-Convex的其他約束及目標函數進行離散化處理,最終可得參數凸問題PED-disConvex:

(39)

該問題的待優化狀態擴展為問題PR的(P+1)(N+1)倍。

3.3 序列凸優化求解策略

綜上所述,非線性隨機最優控制問題PR的解可通過迭代求解問題PED-disConvex進行逼近,并以當前獲得軌跡不斷更新參考軌跡,以提高解的收斂性和精度,即令

(40)

直到滿足約束

(41)

式中:εX和εu為收斂誤差限。

圖1給出了基于序列凸優化的魯棒軌跡優化算法原理框圖。

圖1 魯棒軌跡優化算法原理Fig.1 Principle of robust trajectory optimization algorithm

4 仿真分析

本節首先通過定量分析對兩種不確定性量化傳播算法的狀態擴展維數及計算效率進行比較,之后以典型的高超聲速滑翔問題為例,開展多不確定性條件下的魯棒軌跡快速優化仿真,以驗證本文所設計算法的有效性。

4.1 不確定性量化傳播算法定量對比分析

如前所述,當采用非嵌入式混沌多項式進行參數不確定性量化傳播時,時常用到兩種典型模型轉化方法,即SRSM方法和GQ方法,以下通過定量分析對這兩種方法的計算效率進行比較。

一般而言,對于具有p維隨機變量的不確定問題,若PC多項式展開的階數為p,根據文獻[22]可知,GQ方法的狀態擴展倍數為p,而SRSM的狀態擴展倍數為p。表2給出了兩種方法在不同p、d取值情況下的狀態擴展情況。

表2 不同條件下GQ與SRSM的狀態擴展性能對比Table 2 Comparison of state expansion between GQ and SRSM under different conditions

由表2可知,在p、d較小的情況下,兩種方法的狀態擴展倍數比較接近;隨著p的增大,兩種方法的狀態擴展倍數均有所增大。相較于GQ方法,SRSM的增長速度較為緩慢;隨著隨機變量維數d的不斷增大,GQ方法的擴展倍數呈指數增長,而SRSM增長較為緩慢。

可見,相較于GQ方法,本文采用的SRSM具有更高的計算效率,尤其在處理高維不確定性問題時具有顯著優勢,而GQ方法一般只能適用于d≤3的情況。此外,下文仿真部分亦給出了兩種方法在精度與計算效率方面的進一步對比分析。

4.2 高超聲速滑翔魯棒軌跡優化仿真

高超聲速滑翔飛行具有航程遠、速度高、飛行環境復雜等典型特性,且面臨多種復雜約束以及參數不確定性和干擾的影響,這使得其軌跡優化問題成為一個典型的快時變、非線性、強耦合、非凸的隨機最優控制問題。文獻[4-7]研究表明,標稱條件下的再入滑翔軌跡優化具有較大的挑戰性和應用價值,而考慮真實情況的多參數不確定性的魯棒軌跡優化問題則具有更大的難度和挑戰性。因此,本節以高超聲速滑翔問題為例,開展多不確定性條件下的魯棒軌跡快速優化仿真,以驗證本文所設計算法的有效性。此外,由于文獻[11]已將基于PC的高斯偽譜法與凸優化方法進行了對比,證明了基于PC的凸優化方法在計算精度相當的情況下,計算效率顯著提高。因此,本文不再將以上兩種方法進行比較,而是將所提基于隨機響應面和混沌多項式的凸優化方法(簡稱為SRSM-PC-CO)與基于GQ及混沌多項式的凸優化方法(簡稱為GQ-PC-CO)進行比較,以體現本文算法在計算效率方面的優勢。

4.2.1 問題描述

假設地球為旋轉圓球,標稱條件下的高超聲速滑翔質心運動模型采用文獻[5]的形式。針對目標點固定情況下的飛行時間最短滑翔問題,若結合任務需求預先選定攻角剖面,而僅將傾側角作為待優化控制量,可將該軌跡優化對應的確定性最優控制問題描述為

(42)

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

(43)

且有

(44)

式中:c為動力學參數向量;ρ為大氣密度;CL,CD為升力、阻力系數;δ為服從特定分布且相互獨立的隨機變量;上劃線“-”代表標稱情況。

為了盡可能降低終端位置散布,將終端經緯度的標準差引入魯棒軌跡優化問題的目標函數,則考慮初始狀態及動力學參數不確定的滑翔魯棒軌跡優化問題為

(45)

式中:rminμ(e),rminσ(e)為地心距下邊界的均值與方差,在算法迭代求解過程中,可基于上一步參考軌跡并結合式(24)計算得到。上述問題符合一般魯棒軌跡優化問題PR的形式,可通過本文所設計的魯棒軌跡優化算法進行求解。

同時,對于上述考慮不確定性的滑翔軌跡魯棒優化問題,初始狀態不確定性為5,過程參數不確定性為3,若同時考慮這些不確定性,則隨機變量維數d=8。考慮該問題具有較強非線性,設置混沌多項式階數p=2,若采用GQ-PC-CO方法進行魯棒軌跡優化,則需將待優化狀態量擴展為標稱情況的38(6 561倍),原始隨機常微分方程隨之擴展為6 561×5維耦合的確定性微分方程組。如考慮凸優化過程中的離散化處理,待優化參數量將會更大,以上高維擴展會導致計算效率下降甚至無法收斂等計算困難;而采用本文所提出SRSM-PC-CO方法進行求解時,待優化狀態只需擴展(8+2)!/(8!2!)(45倍),相較而言,顯著降低了待優化變量個數,因而具有良好的優化計算效率和收斂性。

4.2.2 仿真條件設置

表3 邊界條件取值Table 3 Boundary condition values

表4 參數不確定性設置Table 4 Parameters uncertainty setting

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

4.2.3 不確定性量化傳播仿真分析

首先,基于標稱條件進行確定性滑翔軌跡優化仿真,所得結果如圖2所示。由圖2(a)、圖2(b)可以看出,所得速度-高度曲線連續光滑并且滿足過程約束要求,熱流密度和動壓曲線在一段時間內靠近約束邊界,具有超出約束邊界的潛在風險;由圖2(a)、圖2(c)可以看出,優化結果與積分結果一致,且控制指令反轉次數較少,優化結果具有較好的可行性。

圖2 標稱情況下滑翔軌跡優化結果Fig.2 Optimization results of glide trajectory under nominal conditions

僅考慮表4中的大氣密度、升力和阻力系數3個動力學參數不確定性,基于上述優化結果開展GQ-PC、SRSM-PC方法的不確定性量化傳播仿真,并與MC打靶結果進行對比,以分析不確定性對優化軌跡的影響,同時驗證所給出不確定性量化算法的有效性。設置PC多項式階數p=2,MC打靶次數為3 000次,部分仿真結果如圖3～圖5及表5所示。

圖3 標稱情況下優化控制量MC仿真結果Fig.3 MC simulation result of optimized control quantity under nominal condition

圖4 不確定性量化傳播的狀態均值對比Fig.4 Comparison of state mean using uncertainty propagation

圖5 不確定性量化傳播的狀態標準差對比Fig.5 Comparison of state standard deviation using uncertainty propagation

表5 不確定量化傳播的終端狀態統計量對比Table 5 Comparison of terminal state statistics using uncertainty propagation algorithm

由圖3可知,標稱情況下的軌跡優化結果易受到參數不確定性的影響,從而出現過程約束超出限制、終端精度下降等問題,即軌跡可靠性及魯棒性均面臨不同程度的降低。因此,需要進一步開展不確定條件下的魯棒軌跡優化,以提升軌跡的抗干擾能力及可靠性,進而降低制導控制系統負擔和任務失敗風險。

圖4和圖5分別給出了不同方法對應的部分狀態均值及標準差比較曲線,表5則進一步對比了終端狀態的統計精度。

可以看出,3種方法所得統計結果較為一致,相較于MC方法,GQ-PC與SRSM-PC方法所計算的終端高度均值與標準差誤差小于0.1 km,經緯度誤差小于0.01°,說明GQ-PC與SRSM-PC方法均可逼近打靶結果,即可替代打靶方法進行不確定性量化傳播。在計算精度相當的單次量化傳播過程中,SRSM-PC方法僅需進行20次動力學方程積分,而GQ-PC方法需要進行27次積分,其計算量為SRSM-PC的1.35倍;若進一步增加不確定性維數并考慮優化算法的迭代計算過程,SRSM-PC方法在計算效率方面的優勢將會更加明顯。

4.2.4 兩種算法的魯棒軌跡優化仿真對比

采用本文所設計的SRSM-PC-CO算法進行魯棒軌跡優化仿真,同時與現有的GQ-PC-CO優化方法及確定性優化(deterministic optimization,DO)方法進行對比,以驗證本文算法的有效性及計算效率優勢。設置權重系數kJ=kg=3,其余仿真條件同上,結果對比如圖6所示。隨后,利用圖6(d)中的傾側角指令分別進行MC打靶仿真,得到部分參數的統計結果,如表6所示。其中,Δθμf和Δφμf分別表示終端經度、緯度均值偏差,θσf和φσf分別表示終端經度、緯度標準差;tμf和tσf分別表示飛行時間的均值及標準差,Pe表示軌跡超出過程約束邊界的概率,即事故發生率。

圖6 GQ-PC-CO與SRSM-PC-CO的魯棒軌跡優化結果Fig.6 Robust trajectory optimization results of GQ-PC-CO and SRSM-PC-CO

由圖6(a)和圖6(b)可以看出,兩組魯棒軌跡優化結果基本一致,且高度、熱流密度和動壓曲線相較于DO更遠離過程約束邊界,因而約束超出概率降低,優化軌跡的可靠性得到了明顯提升。

由表6可知,兩組魯棒優化控制量的打靶統計結果較為一致,相較于DO,其終端經、緯度分布更加集中,均值偏差及標準差均有所降低,過程約束超出概率顯著降低。以上結果表明,本文方法與GQ-PC-CO的優化結果較為一致,其軌跡的魯棒性和可靠性均顯著提升,驗證了本文方法的有效性。

表6 兩種方法優化結果的MC仿真結果對比Table 6 Comparison of MC simulation results of optimization results for two methods

此外,由圖6(d)可知,相較于DO,魯棒優化使得傾側角指令在初始下降段基本維持在0°附近,隨后繼續以相對較小的幅值飛行,以增大縱向平面內的法向爬升過載,進而提升前段軌跡高度,最終降低熱流密度以遵守過程約束,而為了完成航向調整,在飛行中后期則需要采用相對更大的傾側角幅值實現快速轉彎,如圖6(c)所示。

表7給出了兩種算法的數值處理性能對比。可見,在相同條件下,本文算法對應的待優化變量數目顯著降低,在迭代步數相當情況下的計算耗時降低了約62%,進一步說明了SRSM-PC-CO魯棒軌跡優化算法在計算效率方面具有顯著優勢。

表7 兩種優化算法的數值計算性能對比Table 7 Numerical computation performance comparison of two optimization algorithms

4.2.5 基于權重參數調整的軌跡優化可靠性與魯棒性能分析

為了進一步探討所提出軌跡優化算法在可靠性和魯棒性方面的權衡能力,以更好地滿足不同任務的需求側重,下文開展基于權重參數調整的魯棒優化仿真與分析。由前文魯棒優化建模可知,調整權重參數kg、kJ可以滿足不同的可靠性和魯棒性需求。若令kg=3,kJ=0,則優化模型演化為基于可靠性的優化[10](reliability-based optimization, RBO);若令kg=kJ=3,則問題演化為魯棒優化(robust optimization, RO)。

以表4中的3個動力學參數及5個初始狀態作為不確定參數,采用SRSM-PC-CO算法進行仿真,同時將其與DO對比,以探究參數調整所帶來的軌跡特性改變。表8給出了不同條件下的權重系數取值,其他仿真條件與前文一致。仿真結果如圖7所示。

表8 權重系數設置Table 8 Weight factors setting

圖7 DO、RBO與RO的結果參數對比Fig.7 Parameter comparison of results of DO, RBO, and RO

由圖7(a)、圖7(b)可知,相較于圖3(b)的確定性優化結果,RO與RBO對應的高度、熱流密度和動壓曲線更遠離約束邊界,即優化軌跡的可靠性得到了明顯提升;與RBO相比,RO對應的高度在第一個波峰處和末段更高,說明這種軌跡形式有利于降低隨機不確定性的影響,提升終端位置和飛行時間指標的魯棒性。由圖7(c)可知,與前文分析類似,為了降低高度曲線在第一個波谷處的熱流峰值,即提升軌跡可靠性,RBO與RO對應初始下降段的傾側角幅值均較小,且大小基本一致,這表明RBO與RO對熱流峰值的可靠性相當;在波谷之后,RO繼續保持更小的傾側角幅值,以保持更大的飛行高度;而在飛行后段,RO則采用了更大的傾側角幅值,以實現航向快速調整,這一現象在圖7(d)的地面航跡曲線結果中亦得到了進一步印證。事實上,軌跡可靠性和魯棒性的提升對飛行器在不同飛行階段的控制能力提出了更高要求,且沒有統一的規律可循,需要針對不同的問題進行專門分析。

利用圖7(c)中的傾側角優化結果進行MC打靶仿真,得到部分參數的統計結果,如圖8及表9所示,其中,圖8給出了打靶結果中部分參數的近似概率密度函數。

圖8 DO、RBO與RO的MC仿真結果概率密度函數Fig.8 Probability density functions of DO, RBO, and RO’s MC simulation result

表9 DO、RBO與RO的MC仿真數據對比Table 9 Data comparison of MC simulation of DO, RBO and RO

由圖8(a)、圖8(b)及表9可知,RBO與RO的終端經、緯度統計精度基本一致,且相較于DO有明顯改善,而由于考慮了終端性能和優化指標的魯棒性,RO對應的終端經、緯度和飛行時間散布在三者中最小。由圖8(c)～圖8(e)及表9可知,相較于DO結果,RBO與RO的熱流、動壓、過載峰值分布均遠離約束邊界,約束超出的概率由DO的67.4%分別降低到0.90%和0.88%,可見兩者的可靠性均有所提升。以上結果表明,與RBO結果相比,RO在保證軌跡可靠性的同時,對不確定性的敏感度最低,具有更好的魯棒性。由圖8(f)及表9可知,相較于DO結果,RBO與RO對應的飛行時間均值較大,這是由于不確定性的引入和可靠性約束的施加會造成優化求解的可行域減小,迫使軌跡高度提升,進而導致飛行時間增大。而相較于RBO,RO的飛行時間均值更大,這是由于RO的目標函數綜合考慮了飛行時間的均值及標準差,而RBO的目標函數僅考慮了飛行時間的均值。因此,軌跡魯棒性與可靠性的提升是以犧牲一定最優性為代價的。在實際應用時,可根據設計需求調整權重系數,以滿足軌跡的可靠性、魯棒性和最優性綜合需求。

表10對比了3種仿真條件下的計算耗時及待優化變量數。由表10可以看出,即使同時考慮初始狀態及參數不確定性,RBO及RO的優化耗時均在可接受范圍內,若采用并行計算,優化耗時可進一步縮短。

表10 DO、RBO與RO的數值計算性能對比Table 10 Numerical computation performance comparison of DO, RBO and RO

5 結 論

針對存在多參數不確定性的軌跡優化問題,本文研究了一種通用、高效的魯棒軌跡優化方法,并以高超聲速再入滑翔軌跡優化為例進行了仿真分析與對比驗證,相關研究結論如下:

(1) 與確定性軌跡優化和傳統魯棒軌跡優化算法相比,所設計魯棒軌跡優化方法可以獲得兼具可靠性和魯棒性的優化軌跡,且具有相當的精度和顯著的計算效率優勢。

(2) 軌跡魯棒性和可靠性的提升是以犧牲一定最優性為代價的,在實際應用時,可根據設計需求調整權重系數,以綜合滿足軌跡的可靠性、魯棒性和最優性需求。

(3) 所構建的基于隨機響應面法及NIPC的不確定性量化傳播模型,可實現初始狀態和過程動力學等多參數不確定隨機問題的有效轉化,且具有狀態擴展維數少、計算效率高的優勢。

(4) 與基于IPC的優化方法相比,本文方法將原系統模型視為“黑箱”,無需針對具體問題進行推導建模,因而具有較高的工程應用價值。