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

(1. 西北工業(yè)大學(xué)航天學(xué)院, 陜西 西安 710072; 2. 中國(guó)運(yùn)載火箭技術(shù)研究院, 北京 100076)

0 引 言

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

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

而對(duì)于另一關(guān)鍵技術(shù),即確定性軌跡優(yōu)化而言,常見的方法主要包括間接法和直接法。間接法因?qū)獾某踔挡聹y(cè)較敏感、需要對(duì)最優(yōu)性必要條件進(jìn)行推導(dǎo)且無(wú)法處理復(fù)雜約束等原因,難以在復(fù)雜約束軌跡優(yōu)化問(wèn)題中得到廣泛應(yīng)用。近年來(lái),以偽譜法[2]、凸優(yōu)化[4-6]等為代表的直接法已被廣泛應(yīng)用于高超聲速等航空航天領(lǐng)域的軌跡優(yōu)化任務(wù),并展現(xiàn)出了良好的數(shù)值計(jì)算優(yōu)勢(shì)和普適性等特點(diǎn)。其中,凸優(yōu)化方法因具有多項(xiàng)式復(fù)雜度、收斂速度快,且其局部最優(yōu)解理論上即為全局最優(yōu)解等特點(diǎn),已逐步被用于求解復(fù)雜約束下的軌跡優(yōu)化問(wèn)題以及考慮不確定性的高維狀態(tài)擴(kuò)展優(yōu)化問(wèn)題。

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

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

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

因此,本文針對(duì)航空航天領(lǐng)域魯棒軌跡優(yōu)化算法存在的問(wèn)題,提出一種高效的魯棒軌跡優(yōu)化設(shè)計(jì)方法。通過(guò)構(gòu)造基于NIPC與SRSM的不確定性量化傳播模型,實(shí)現(xiàn)了不確定性隨機(jī)問(wèn)題的有效轉(zhuǎn)化;隨后,將該模型與凸優(yōu)化方法相結(jié)合,設(shè)計(jì)了一種基于序列凸優(yōu)化算法的求解策略,以實(shí)現(xiàn)對(duì)該高維確定性問(wèn)題的快速求解。最終,以高超聲速滑翔再入為例進(jìn)行仿真,驗(yàn)證了本文方法的有效性。

1 魯棒軌跡優(yōu)化問(wèn)題描述

一般的軌跡優(yōu)化問(wèn)題可以歸結(jié)為如下確定性最優(yōu)控制問(wèn)題PD:

(1)

當(dāng)系統(tǒng)參數(shù)向量w及狀態(tài)初值x(t0)存在不確定性時(shí),目標(biāo)函數(shù)J、動(dòng)力學(xué)函數(shù)f及約束函數(shù)均變?yōu)殡S機(jī)函數(shù)。因此,原確定性優(yōu)化問(wèn)題PD轉(zhuǎn)變?yōu)殡S機(jī)優(yōu)化問(wèn)題。為了降低優(yōu)化軌跡對(duì)不確定性的敏感度,即當(dāng)存在不確定性時(shí),預(yù)期的軌跡離散度在統(tǒng)計(jì)意義上盡可能小,通常需要在目標(biāo)函數(shù)、約束函數(shù)中引入隨機(jī)量的統(tǒng)計(jì)矩。因此,一般魯棒軌跡優(yōu)化問(wèn)題PR可描述為

(2)

式中:下標(biāo)μ和σ分別代表均值和標(biāo)準(zhǔn)差;kJ、kg為待設(shè)計(jì)權(quán)重系數(shù);x0μ與x0σ分別為初始狀態(tài)的均值和標(biāo)準(zhǔn)差;εf為狀態(tài)終值標(biāo)準(zhǔn)差的上界。當(dāng)采用上述模型進(jìn)行軌跡優(yōu)化時(shí),需要引入不確定性量化傳播技術(shù)計(jì)算相關(guān)統(tǒng)計(jì)量并處理隨機(jī)動(dòng)力學(xué)方程。

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

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

2 不確定量化傳播建模與問(wèn)題轉(zhuǎn)化

本節(jié)將NIPC與SRSM相結(jié)合,同時(shí)引入拉丁超立方采樣策略,構(gòu)建了不確定性量化傳播模型,從而將魯棒軌跡優(yōu)化問(wèn)題轉(zhuǎn)化為高維狀態(tài)空間中的等價(jià)確定性優(yōu)化問(wèn)題。以下給出算法原理及實(shí)現(xiàn)步驟。

2.1 隨機(jī)狀態(tài)方程轉(zhuǎn)化

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

(3)

(4)

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

表1 正交基函數(shù)與隨機(jī)變量分布類型對(duì)應(yīng)關(guān)系Table 1 Correspondence relationship between distribution type of orthogonal basis function and random variables

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

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

(5)

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

(6)

式中:X與W分別表示隨機(jī)狀態(tài)x與參數(shù)w的PC系數(shù)向量,且有

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)

設(shè)Ω(Δ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,則可通過(guò)最小二乘回歸法求解式的微分方程組,即有

(11)

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

2.2 過(guò)程約束及目標(biāo)函數(shù)的PC展開

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

(12)

(13)

式中:gj(x(t),u(t),t)與Jj(x(t),u(t),t)為對(duì)應(yīng)函數(shù)的PC多項(xiàng)式系數(shù)。將隨機(jī)樣本點(diǎn){Δ1,Δ2,…,Δns}代入方程(12)、方程(13)可得

(14)

(15)

(16)

(17)

(18)

(19)

類似地,有

(20)

(21)

(22)

2.3 統(tǒng)計(jì)量計(jì)算

根據(jù)文獻(xiàn)[15],可將式(2)中隨機(jī)狀態(tài)x、過(guò)程約束函數(shù)g(x(t),u(t),t,w)及目標(biāo)函數(shù)J(x(t),u(t),t,w)的均值和方差表示為PC多項(xiàng)式系數(shù)的函數(shù),即

(23)

(24)

(25)

式中:〈·〉表示內(nèi)積運(yùn)算;gj,Jj(j=1,2,…,P)可通過(guò)式(21)、式(22)表示為X、u的非線性函數(shù),進(jìn)而可通過(guò)以上公式計(jì)算狀態(tài)、過(guò)程約束與目標(biāo)函數(shù)的統(tǒng)計(jì)量。

至此,已將魯棒軌跡優(yōu)化問(wèn)題PR轉(zhuǎn)化為關(guān)于擴(kuò)展?fàn)顟B(tài)X及控制量u的高維確定性優(yōu)化問(wèn)題PED,進(jìn)而可利用凸優(yōu)化算法進(jìn)行處理求解,限于篇幅,問(wèn)題PED在此不再贅述。

3 基于凸優(yōu)化的魯棒優(yōu)化求解策略

考慮到復(fù)雜軌跡優(yōu)化問(wèn)題大多是非凸的,如何將非凸問(wèn)題進(jìn)行凸化處理是應(yīng)用凸優(yōu)化技術(shù)求解最優(yōu)軌跡的關(guān)鍵,并吸引了廣大學(xué)者對(duì)此開展了深入研究,相關(guān)成果亦在航空航天領(lǐng)域最優(yōu)控制及軌跡優(yōu)化問(wèn)題上得到了較為廣泛的應(yīng)用[4-6]。

本節(jié)引入逐次線性化[27]等技術(shù)對(duì)第2節(jié)中的高維確定性優(yōu)化問(wèn)題PED進(jìn)行凸化處理,并將該問(wèn)題離散處理為參數(shù)凸優(yōu)化問(wèn)題,進(jìn)而建立相應(yīng)的序列凸優(yōu)化迭代算法,對(duì)該問(wèn)題進(jìn)行求解。序列凸優(yōu)化方法通過(guò)外環(huán)迭代更新參考軌跡,內(nèi)環(huán)求解凸問(wèn)題逐步逼近原問(wèn)題的最優(yōu)解[27]。以下給出凸化、離散處理過(guò)程以及迭代算法的原理和實(shí)現(xiàn)步驟。

3.1 凸化處理

3.1.1 狀態(tài)方程凸化

分析式(11)可知,H為常值矩陣,F(X,u,t,W)為關(guān)于擴(kuò)展?fàn)顟B(tài)X的非線性函數(shù),因此該式為非線性微分方程。下面采用逐次線性化方法對(duì)該式進(jìn)行凸化處理。

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

(26)

定義:

(27)

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

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

(28)

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

3.1.2 過(guò)程約束及目標(biāo)函數(shù)凸化

類似地,考慮式(2)中的過(guò)程約束及目標(biāo)函數(shù)均為關(guān)于擴(kuò)展?fàn)顟B(tài)X及控制量u的非線性函數(shù),仍采用逐次線性化方法對(duì)其進(jìn)行凸化處理。

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

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 狀態(tài)不等式約束凸化

考慮式(2)中的狀態(tài)不等式約束為關(guān)于擴(kuò)展?fàn)顟B(tài)X的非線性函數(shù),其凸化處理方法[11]如下。

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

(33)

Q=LTL

(34)

狀態(tài)變量xi的標(biāo)準(zhǔn)差可表示為

(35)

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

(36)

式中:ximax與ximin分別為第i個(gè)狀態(tài)量對(duì)應(yīng)的上下邊界。

至此,問(wèn)題PED轉(zhuǎn)化為以下等價(jià)的高維確定性軌跡凸優(yōu)化問(wèn)題PED-Convex:

(37)

需要說(shuō)明的是,問(wèn)題PED-Convex的待優(yōu)化狀態(tài)為X,擴(kuò)展為問(wèn)題PR的(P+1)倍。

3.2 離散化處理

為了應(yīng)用數(shù)值方法對(duì)問(wèn)題PED-Convex進(jìn)行求解,需要對(duì)狀態(tài)量和控制量進(jìn)行離散化處理?；谔菪畏▌t,將自變量變化域[t0,tf]等間距離散為N個(gè)間隔,自變量離散為{t0,t1,…,tN},狀態(tài)量離散為{X0,X1,…,XN},控制量離散為{u0,u1,…,uN}。

對(duì)式(37)的線性動(dòng)力學(xué)方程進(jìn)行離散化,可得:

(38)

類似地,對(duì)問(wèn)題PED-Convex的其他約束及目標(biāo)函數(shù)進(jìn)行離散化處理,最終可得參數(shù)凸問(wèn)題PED-disConvex:

(39)

該問(wèn)題的待優(yōu)化狀態(tài)擴(kuò)展為問(wèn)題PR的(P+1)(N+1)倍。

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

綜上所述,非線性隨機(jī)最優(yōu)控制問(wèn)題PR的解可通過(guò)迭代求解問(wèn)題PED-disConvex進(jìn)行逼近,并以當(dāng)前獲得軌跡不斷更新參考軌跡,以提高解的收斂性和精度,即令

(40)

直到滿足約束

(41)

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

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

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

4 仿真分析

本節(jié)首先通過(guò)定量分析對(duì)兩種不確定性量化傳播算法的狀態(tài)擴(kuò)展維數(shù)及計(jì)算效率進(jìn)行比較,之后以典型的高超聲速滑翔問(wèn)題為例,開展多不確定性條件下的魯棒軌跡快速優(yōu)化仿真,以驗(yàn)證本文所設(shè)計(jì)算法的有效性。

4.1 不確定性量化傳播算法定量對(duì)比分析

如前所述,當(dāng)采用非嵌入式混沌多項(xiàng)式進(jìn)行參數(shù)不確定性量化傳播時(shí),時(shí)常用到兩種典型模型轉(zhuǎn)化方法,即SRSM方法和GQ方法,以下通過(guò)定量分析對(duì)這兩種方法的計(jì)算效率進(jìn)行比較。

一般而言,對(duì)于具有p維隨機(jī)變量的不確定問(wèn)題,若PC多項(xiàng)式展開的階數(shù)為p,根據(jù)文獻(xiàn)[22]可知,GQ方法的狀態(tài)擴(kuò)展倍數(shù)為p,而SRSM的狀態(tài)擴(kuò)展倍數(shù)為p。表2給出了兩種方法在不同p、d取值情況下的狀態(tài)擴(kuò)展情況。

表2 不同條件下GQ與SRSM的狀態(tài)擴(kuò)展性能對(duì)比Table 2 Comparison of state expansion between GQ and SRSM under different conditions

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

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

4.2 高超聲速滑翔魯棒軌跡優(yōu)化仿真

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

4.2.1 問(wèn)題描述

假設(shè)地球?yàn)樾D(zhuǎn)圓球,標(biāo)稱條件下的高超聲速滑翔質(zhì)心運(yùn)動(dòng)模型采用文獻(xiàn)[5]的形式。針對(duì)目標(biāo)點(diǎn)固定情況下的飛行時(shí)間最短滑翔問(wèn)題,若結(jié)合任務(wù)需求預(yù)先選定攻角剖面,而僅將傾側(cè)角作為待優(yōu)化控制量,可將該軌跡優(yōu)化對(duì)應(yīng)的確定性最優(yōu)控制問(wèn)題描述為

(42)

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

(43)

且有

(44)

式中:c為動(dòng)力學(xué)參數(shù)向量;ρ為大氣密度;CL,CD為升力、阻力系數(shù);δ為服從特定分布且相互獨(dú)立的隨機(jī)變量;上劃線“-”代表標(biāo)稱情況。

為了盡可能降低終端位置散布,將終端經(jīng)緯度的標(biāo)準(zhǔn)差引入魯棒軌跡優(yōu)化問(wèn)題的目標(biāo)函數(shù),則考慮初始狀態(tài)及動(dòng)力學(xué)參數(shù)不確定的滑翔魯棒軌跡優(yōu)化問(wèn)題為

(45)

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

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

4.2.2 仿真條件設(shè)置

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

表4 參數(shù)不確定性設(shè)置Table 4 Parameters uncertainty setting

所有仿真均在搭載Intel Core i7-8700 3.20 GHz Intel處理器的臺(tái)式機(jī)完成,仿真環(huán)境為Matlab 2016b平臺(tái)。基于自行開發(fā)的代碼進(jìn)行不確定性量化傳播,同時(shí)基于CVX工具包[29]進(jìn)行軌跡優(yōu)化算法開發(fā),并調(diào)用SDPT3求解器求解凸優(yōu)化子問(wèn)題PED-disConvex。

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

首先,基于標(biāo)稱條件進(jìn)行確定性滑翔軌跡優(yōu)化仿真,所得結(jié)果如圖2所示。由圖2(a)、圖2(b)可以看出,所得速度-高度曲線連續(xù)光滑并且滿足過(guò)程約束要求,熱流密度和動(dòng)壓曲線在一段時(shí)間內(nèi)靠近約束邊界,具有超出約束邊界的潛在風(fēng)險(xiǎn);由圖2(a)、圖2(c)可以看出,優(yōu)化結(jié)果與積分結(jié)果一致,且控制指令反轉(zhuǎn)次數(shù)較少,優(yōu)化結(jié)果具有較好的可行性。

圖2 標(biāo)稱情況下滑翔軌跡優(yōu)化結(jié)果Fig.2 Optimization results of glide trajectory under nominal conditions

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

圖3 標(biāo)稱情況下優(yōu)化控制量MC仿真結(jié)果Fig.3 MC simulation result of optimized control quantity under nominal condition

圖4 不確定性量化傳播的狀態(tài)均值對(duì)比Fig.4 Comparison of state mean using uncertainty propagation

圖5 不確定性量化傳播的狀態(tài)標(biāo)準(zhǔn)差對(duì)比Fig.5 Comparison of state standard deviation using uncertainty propagation

表5 不確定量化傳播的終端狀態(tài)統(tǒng)計(jì)量對(duì)比Table 5 Comparison of terminal state statistics using uncertainty propagation algorithm

由圖3可知,標(biāo)稱情況下的軌跡優(yōu)化結(jié)果易受到參數(shù)不確定性的影響,從而出現(xiàn)過(guò)程約束超出限制、終端精度下降等問(wèn)題,即軌跡可靠性及魯棒性均面臨不同程度的降低。因此,需要進(jìn)一步開展不確定條件下的魯棒軌跡優(yōu)化,以提升軌跡的抗干擾能力及可靠性,進(jìn)而降低制導(dǎo)控制系統(tǒng)負(fù)擔(dān)和任務(wù)失敗風(fēng)險(xiǎn)。

圖4和圖5分別給出了不同方法對(duì)應(yīng)的部分狀態(tài)均值及標(biāo)準(zhǔn)差比較曲線,表5則進(jìn)一步對(duì)比了終端狀態(tài)的統(tǒng)計(jì)精度。

可以看出,3種方法所得統(tǒng)計(jì)結(jié)果較為一致,相較于MC方法,GQ-PC與SRSM-PC方法所計(jì)算的終端高度均值與標(biāo)準(zhǔn)差誤差小于0.1 km,經(jīng)緯度誤差小于0.01°,說(shuō)明GQ-PC與SRSM-PC方法均可逼近打靶結(jié)果,即可替代打靶方法進(jìn)行不確定性量化傳播。在計(jì)算精度相當(dāng)?shù)膯未瘟炕瘋鞑ミ^(guò)程中,SRSM-PC方法僅需進(jìn)行20次動(dòng)力學(xué)方程積分,而GQ-PC方法需要進(jìn)行27次積分,其計(jì)算量為SRSM-PC的1.35倍;若進(jìn)一步增加不確定性維數(shù)并考慮優(yōu)化算法的迭代計(jì)算過(guò)程,SRSM-PC方法在計(jì)算效率方面的優(yōu)勢(shì)將會(huì)更加明顯。

4.2.4 兩種算法的魯棒軌跡優(yōu)化仿真對(duì)比

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

圖6 GQ-PC-CO與SRSM-PC-CO的魯棒軌跡優(yōu)化結(jié)果Fig.6 Robust trajectory optimization results of GQ-PC-CO and SRSM-PC-CO

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

由表6可知,兩組魯棒優(yōu)化控制量的打靶統(tǒng)計(jì)結(jié)果較為一致,相較于DO,其終端經(jīng)、緯度分布更加集中,均值偏差及標(biāo)準(zhǔn)差均有所降低,過(guò)程約束超出概率顯著降低。以上結(jié)果表明,本文方法與GQ-PC-CO的優(yōu)化結(jié)果較為一致,其軌跡的魯棒性和可靠性均顯著提升,驗(yàn)證了本文方法的有效性。

表6 兩種方法優(yōu)化結(jié)果的MC仿真結(jié)果對(duì)比Table 6 Comparison of MC simulation results of optimization results for two methods

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

表7給出了兩種算法的數(shù)值處理性能對(duì)比?？梢?在相同條件下,本文算法對(duì)應(yīng)的待優(yōu)化變量數(shù)目顯著降低,在迭代步數(shù)相當(dāng)情況下的計(jì)算耗時(shí)降低了約62%,進(jìn)一步說(shuō)明了SRSM-PC-CO魯棒軌跡優(yōu)化算法在計(jì)算效率方面具有顯著優(yōu)勢(shì)。

表7 兩種優(yōu)化算法的數(shù)值計(jì)算性能對(duì)比Table 7 Numerical computation performance comparison of two optimization algorithms

4.2.5 基于權(quán)重參數(shù)調(diào)整的軌跡優(yōu)化可靠性與魯棒性能分析

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

以表4中的3個(gè)動(dòng)力學(xué)參數(shù)及5個(gè)初始狀態(tài)作為不確定參數(shù),采用SRSM-PC-CO算法進(jìn)行仿真,同時(shí)將其與DO對(duì)比,以探究參數(shù)調(diào)整所帶來(lái)的軌跡特性改變。表8給出了不同條件下的權(quán)重系數(shù)取值,其他仿真條件與前文一致。仿真結(jié)果如圖7所示。

表8 權(quán)重系數(shù)設(shè)置Table 8 Weight factors setting

圖7 DO、RBO與RO的結(jié)果參數(shù)對(duì)比Fig.7 Parameter comparison of results of DO, RBO, and RO

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

利用圖7(c)中的傾側(cè)角優(yōu)化結(jié)果進(jìn)行MC打靶仿真,得到部分參數(shù)的統(tǒng)計(jì)結(jié)果,如圖8及表9所示,其中,圖8給出了打靶結(jié)果中部分參數(shù)的近似概率密度函數(shù)。

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

表9 DO、RBO與RO的MC仿真數(shù)據(jù)對(duì)比Table 9 Data comparison of MC simulation of DO, RBO and RO

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

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

表10 DO、RBO與RO的數(shù)值計(jì)算性能對(duì)比Table 10 Numerical computation performance comparison of DO, RBO and RO

5 結(jié) 論

針對(duì)存在多參數(shù)不確定性的軌跡優(yōu)化問(wèn)題,本文研究了一種通用、高效的魯棒軌跡優(yōu)化方法,并以高超聲速再入滑翔軌跡優(yōu)化為例進(jìn)行了仿真分析與對(duì)比驗(yàn)證,相關(guān)研究結(jié)論如下:

(1) 與確定性軌跡優(yōu)化和傳統(tǒng)魯棒軌跡優(yōu)化算法相比,所設(shè)計(jì)魯棒軌跡優(yōu)化方法可以獲得兼具可靠性和魯棒性的優(yōu)化軌跡,且具有相當(dāng)?shù)木群惋@著的計(jì)算效率優(yōu)勢(shì)。

(2) 軌跡魯棒性和可靠性的提升是以犧牲一定最優(yōu)性為代價(jià)的,在實(shí)際應(yīng)用時(shí),可根據(jù)設(shè)計(jì)需求調(diào)整權(quán)重系數(shù),以綜合滿足軌跡的可靠性、魯棒性和最優(yōu)性需求。

(3) 所構(gòu)建的基于隨機(jī)響應(yīng)面法及NIPC的不確定性量化傳播模型,可實(shí)現(xiàn)初始狀態(tài)和過(guò)程動(dòng)力學(xué)等多參數(shù)不確定隨機(jī)問(wèn)題的有效轉(zhuǎn)化,且具有狀態(tài)擴(kuò)展維數(shù)少、計(jì)算效率高的優(yōu)勢(shì)。

(4) 與基于IPC的優(yōu)化方法相比,本文方法將原系統(tǒng)模型視為“黑箱”,無(wú)需針對(duì)具體問(wèn)題進(jìn)行推導(dǎo)建模,因而具有較高的工程應(yīng)用價(jià)值。