劉 旭,李 響,張后軍,郭宇恒,王曉鵬

(1.北京理工大學(xué)宇航學(xué)院,北京 100081；2.北京機(jī)電工程研究所,北京 100083；3.上海機(jī)電工程研究所,上海 200233)

0 引 言

傳統(tǒng)的再入軌跡優(yōu)化不考慮不確定性因素,即默認(rèn)優(yōu)化模型中所有參數(shù)都是準(zhǔn)確的[1-3],而實(shí)際飛行過程中不確定性因素客觀存在,如飛行器氣動(dòng)參數(shù)誤差、大氣參數(shù)誤差等,這些誤差的影響一般通過后續(xù)的控制系統(tǒng)設(shè)計(jì)消除。近年來,有學(xué)者在飛行軌跡優(yōu)化階段就將這些不確定性因素考慮進(jìn)去,提出了飛行軌跡的魯棒優(yōu)化(Robust design optimization,RDO)[4]與可靠性優(yōu)化(Reliability-based design optimization,RBDO)[5]。RDO在最小化目標(biāo)函數(shù)期望的同時(shí),盡量使目標(biāo)函數(shù)對(duì)不確定性因素不敏感；RBDO著重于在不確定性因素影響下將失效概率限制在某個(gè)值以內(nèi)。通過RDO或RBDO優(yōu)化后得到的飛行軌跡對(duì)不確定性因素的敏感度降低,自身就有一定抗干擾的能力。

文獻(xiàn)[4]指出在考慮不確定性的情況下,相較于飛行器靜態(tài)魯棒優(yōu)化(如氣動(dòng)優(yōu)化、結(jié)構(gòu)優(yōu)化等),飛行器動(dòng)態(tài)魯棒優(yōu)化(如飛行軌跡優(yōu)化)的研究比較有限,因此提出一種不確定條件下的飛行軌跡魯棒優(yōu)化方法；文獻(xiàn)[6]指出飛機(jī)軌跡規(guī)劃和預(yù)測(cè)中的不確定性給未來空中交通管理系統(tǒng)帶來了重大挑戰(zhàn),因此提出一種考慮氣象不確定性的基于最優(yōu)控制的飛機(jī)軌跡規(guī)劃方法；文獻(xiàn)[7]注意到不確定性(尤其是風(fēng)場(chǎng)不確定性)對(duì)于無人機(jī)動(dòng)態(tài)滑翔軌跡的影響,在量化不確定性后提出一種動(dòng)態(tài)滑翔軌跡魯棒優(yōu)化方法；火星大氣密度具有的強(qiáng)不確定性是火星進(jìn)入段飛行任務(wù)的主要干擾源,對(duì)此文獻(xiàn)[8]提出一種火星大氣進(jìn)入段的軌跡魯棒優(yōu)化,得到的最優(yōu)軌跡在不確定條件下的魯棒性明顯優(yōu)于確定性優(yōu)化,而文獻(xiàn)[5]則是提出用序列可靠性優(yōu)化(Reliability-based sequential optimization,RBSO)來降低參數(shù)不確定性對(duì)火星進(jìn)入段軌跡的影響,同樣得到優(yōu)于確定性優(yōu)化的結(jié)果；文獻(xiàn)[9]應(yīng)用RBDO優(yōu)化飛行器在不規(guī)則小行星軟著陸的軌跡,將不確定動(dòng)態(tài)模型的軌跡優(yōu)化轉(zhuǎn)化為帶有可靠性約束的參數(shù)優(yōu)化,并應(yīng)用序列優(yōu)化和可靠性評(píng)估(Sequential optimization and reliability assessment,SORA)流程求解；文獻(xiàn)[10]利用脫敏最優(yōu)控制(Desensitized optimal control,DOC),將標(biāo)稱軌跡設(shè)計(jì)和標(biāo)稱軌跡跟蹤結(jié)合在一起,得到的火星下降段脫敏軌跡對(duì)不確定性和擾動(dòng)的敏感性大幅降低。

在飛行軌跡魯棒優(yōu)化或可靠性優(yōu)化的實(shí)際計(jì)算中,不確定性傳播(Uncertainty propagation,UP),即如何根據(jù)輸入不確定量計(jì)算得到輸出不確定量的統(tǒng)計(jì)特征(一般指均值和方差)是關(guān)鍵步驟。傳統(tǒng)蒙特卡洛(Monte Carlo,MC)仿真法易于操作,對(duì)各種復(fù)雜問題均有很強(qiáng)的適應(yīng)性,得到了較廣泛的應(yīng)用,也常用于檢驗(yàn)其他方法的準(zhǔn)確性。但MC仿真法的精度依賴于大量樣本點(diǎn)上的仿真次數(shù),為得到可靠結(jié)果需要的高昂計(jì)算代價(jià)限制了它在軌跡魯棒優(yōu)化中的應(yīng)用?；煦缍囗?xiàng)式展開(Polynomial chaos expansion,PCE)是一種基于Askey正交多項(xiàng)式隨機(jī)展開的不確定性傳播分析法,能高效地計(jì)算輸出不確定性變量的統(tǒng)計(jì)特征[11],近年來在飛行軌跡魯棒優(yōu)化中得到了較廣泛的應(yīng)用[4,6-8]。基于可靠性的飛行軌跡優(yōu)化[5,9]多采用最可能失效點(diǎn)(Most probable point,MPP)法進(jìn)行不確定性分析,包括一次可靠度法(First-order reliability method,FORM)[12]和二次可靠度法(Second-order reliability method,SORM)[13],其基本思路是在MPP點(diǎn)進(jìn)行局部展開來近似極限狀態(tài)函數(shù),并以此計(jì)算失效概率。文獻(xiàn)[14]引入人工神經(jīng)網(wǎng)絡(luò)模型,用訓(xùn)練后的神經(jīng)網(wǎng)絡(luò)作為復(fù)雜系統(tǒng)的替代模型,改善單獨(dú)使用MC法的不足。文獻(xiàn)[15] 通過Kriging代理模型來構(gòu)建目標(biāo)變量和不確定性輸入?yún)?shù)的響應(yīng)關(guān)系,以此來分析目標(biāo)變量的不確定性信息。

除了前述PCE方法與基于MPP的方法,協(xié)方差分析描述函數(shù)技術(shù)(Covariance analysis describing function technique,CADET)是另外一種高效的隨機(jī)系統(tǒng)統(tǒng)計(jì)特征計(jì)算方法。該方法最早由Gelb提出[16],先采用統(tǒng)計(jì)線性化方法將非線性系統(tǒng)轉(zhuǎn)化為近似線性系統(tǒng),再用協(xié)方差分析法導(dǎo)出系統(tǒng)隨機(jī)狀態(tài)量均值和協(xié)方差的傳播方程,只需一次求解就能確定特定隨機(jī)系統(tǒng)的統(tǒng)計(jì)特征,在飛行器制導(dǎo)系統(tǒng)性能統(tǒng)計(jì)分析中得到了廣泛應(yīng)用。文獻(xiàn)[17]應(yīng)用CADET進(jìn)行單發(fā)武器的射擊效能分析；文獻(xiàn)[18] 基于CADET建立了制導(dǎo)炸彈半實(shí)物仿真系統(tǒng)的誤差傳播理論,分析了仿真姿態(tài)角和姿態(tài)角速度的誤差對(duì)仿真結(jié)果的影響；文獻(xiàn)[19]基于閉環(huán)控制系統(tǒng)下線性化的相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)模型,采用CADET對(duì)定義航天器誤差橢球的協(xié)方差矩陣進(jìn)行分析；文獻(xiàn)[20] 針對(duì)火星探測(cè)捕獲的誤差傳播,利用CADET分析單一誤差和綜合誤差對(duì)不同捕獲策略入軌精度的影響?，F(xiàn)有的文獻(xiàn)中,CADET主要用于制導(dǎo)精度分析,本研究將CADET引入再入軌跡魯棒優(yōu)化中,實(shí)現(xiàn)了再入軌跡魯棒優(yōu)化的高效計(jì)算。

1 飛行軌跡魯棒優(yōu)化模型

不失一般性,先以傳統(tǒng)確定性條件下Bolza型飛行軌跡優(yōu)化問題構(gòu)建飛行軌跡魯棒優(yōu)化模型,可以描述為：將u(t)作為控制量,使得

(1a)

(1b)

C[x(t),u(t),t]≤0

(1c)

H[x(t0),u(t0),x(tf),u(tf),tf]=0

(1d)

式中：x(t)是狀態(tài)向量；u(t)是控制向量；f是動(dòng)態(tài)約束向量函數(shù)；C是過程約束向量函數(shù)；H是端點(diǎn)約束向量函數(shù)。

現(xiàn)在考慮不確定性因素的影響,將隨機(jī)量s引入式(1),得到：將u(t)作為控制量,使得

u(t),t]dt

(2a)

(2b)

C[x(t,s),u(t),t]≤0

(2c)

H[x(t0,s),u(t0),x(tf,s),u(tf),tf]=0(2d)

由于引入了隨機(jī)量s,式(1)中的狀態(tài)量x(t)由確定量轉(zhuǎn)變?yōu)殡S機(jī)量x(t,s)。相應(yīng)地,式(1)中確定性的目標(biāo)函數(shù)式(1a)、動(dòng)力學(xué)方程式(1b)、過程約束式(1c)和端點(diǎn)約束式(1d)分別轉(zhuǎn)變?yōu)殡S機(jī)條件下的式(2a)、(2b)、(2c)和(2d)。

由于含有隨機(jī)量s,式(2)中的目標(biāo)函數(shù)式(2a)和約束條件式(2c-2d)與傳統(tǒng)軌跡優(yōu)化中的目標(biāo)函數(shù)式(1a)和約束條件式(1c-1d)有本質(zhì)的不同?，F(xiàn)有的軌跡優(yōu)化方法都是針對(duì)如式(1)所示的確定性問題發(fā)展起來的,不能直接處理式(2)所示的含有隨機(jī)量的優(yōu)化問題。為此,引入飛行器氣動(dòng)和結(jié)構(gòu)優(yōu)化中魯棒優(yōu)化的理念[21-22],將傳統(tǒng)確定條件下飛行軌跡優(yōu)化模型轉(zhuǎn)化為如式(3)所示不確定條件下的飛行軌跡魯棒優(yōu)化模型：將u(t)作為控制量,使得

min μ(J)+kJ·σ(J)

(3a)

(3b)

μ(C[x(t,s),u(t),t])+kC·σ(C[x(t,s),

u(t),t])≤0

(3c)

μ(H[x(t0,s),u(t0),x(tf,s),u(tf),tf])=0

(3d)

σ(H[x(t0,s),u(t0),x(tf,s),u(tf),tf])≤εH

(3e)

式中：μ(·)和σ(·)分別表示隨機(jī)量的均值和標(biāo)準(zhǔn)差；kJ和kC是權(quán)重系數(shù)；εH是隨機(jī)端點(diǎn)約束標(biāo)準(zhǔn)差的容許值。

式(3)所表達(dá)的具體含義如下：

式(3a)與目標(biāo)函數(shù)式(1a)對(duì)應(yīng),是隨機(jī)目標(biāo)函數(shù)的均值μ(J)和標(biāo)準(zhǔn)差σ(J)的加權(quán)和,通過最小化隨機(jī)目標(biāo)函數(shù)的均值來表達(dá)確定性目標(biāo)函數(shù)式(1a)的設(shè)計(jì)意圖,同時(shí)通過最小化隨機(jī)目標(biāo)函數(shù)的標(biāo)準(zhǔn)差來限制它的離散程度。目標(biāo)函數(shù)的最優(yōu)性和魯棒性可以通過選取權(quán)重系數(shù)達(dá)到設(shè)計(jì)要求,理論上kJ值越大,目標(biāo)函數(shù)的離散程度越小。

式(3c)與過程約束式(1c)對(duì)應(yīng),式(1c)是不等式約束,在魯棒優(yōu)化中,不等式約束通過隨機(jī)約束的均值μ(C)和標(biāo)準(zhǔn)差σ(C)的加權(quán)和來表達(dá),kC越大,所得最優(yōu)軌跡的魯棒性越高,滿足過程約束的概率越大,但求解難度也會(huì)相應(yīng)增加。

式(3d)、(3e)與端點(diǎn)約束式(1d)對(duì)應(yīng),式(1d)是等式約束,在魯棒優(yōu)化中,等式約束通過隨機(jī)約束的均值式(3d)和標(biāo)準(zhǔn)差式(3e)兩組約束來表達(dá),這與不等式約束是不同的。均值約束要求等于原問題中的值,而標(biāo)準(zhǔn)差約束要求小于εH,εH越小,魯棒性越高。

總體來看,魯棒優(yōu)化模型式(3)通過要求隨機(jī)優(yōu)化模型式(2)中隨機(jī)量的均值滿足確定性優(yōu)化模型式(1)中相應(yīng)的約束,來表達(dá)原有的設(shè)計(jì)意圖；同時(shí)通過約束隨機(jī)優(yōu)化式(2)中隨機(jī)量的標(biāo)準(zhǔn)差,限制其離散程度,提高最優(yōu)軌跡的魯棒性。這樣從統(tǒng)計(jì)角度來看,式(3)得到的最優(yōu)軌跡在隨機(jī)擾動(dòng)的影響下偏離設(shè)計(jì)值的程度要小于式(1),從而達(dá)到魯棒優(yōu)化設(shè)計(jì)的目標(biāo)。

值得注意的是,魯棒優(yōu)化模型中式(3b)目前仍是含有隨機(jī)量s的動(dòng)力學(xué)方程,現(xiàn)有軌跡優(yōu)化方法(如偽譜法等)難以直接處理這類問題,在下節(jié)中將討論如何利用協(xié)方差分析描述函數(shù)技術(shù)(CADET)將含有隨機(jī)量的魯棒優(yōu)化問題式(3)轉(zhuǎn)化為等效的確定性優(yōu)化問題,從而可以采用偽譜法等現(xiàn)有軌跡優(yōu)化方法求解。

2 采用CADET重構(gòu)軌跡魯棒優(yōu)化模型

對(duì)于式(3)所示的魯棒優(yōu)化模型,如何高效計(jì)算式(3a,3c-3e)中的統(tǒng)計(jì)特征(均值μ(·)和標(biāo)準(zhǔn)差σ(·))是關(guān)鍵步驟。蒙特卡洛法程序結(jié)構(gòu)簡(jiǎn)單,易于實(shí)現(xiàn),得到了廣泛應(yīng)用。根據(jù)大數(shù)定律,樣本數(shù)N→∞時(shí),蒙特卡洛法的結(jié)果以概率1收斂于真實(shí)值。一般認(rèn)為,樣本點(diǎn)數(shù)量達(dá)到104時(shí),蒙特卡洛法的結(jié)果才是穩(wěn)定收斂和可靠的[23],工程上通常要成千上萬樣本點(diǎn)才能得到滿足精度要求的結(jié)果。對(duì)單次隨機(jī)過程的統(tǒng)計(jì)特征分析而言,蒙特卡洛法是可以接受的；而軌跡優(yōu)化計(jì)算是一個(gè)迭代過程,在每一迭代步都需要計(jì)算均值μ(·)和標(biāo)準(zhǔn)差σ(·)的值,在這種情況下,蒙特卡洛法將導(dǎo)致巨大的計(jì)算量,在實(shí)際計(jì)算中是不可行的。根據(jù)上述情況,本文采用CADET計(jì)算隨機(jī)系統(tǒng)的統(tǒng)計(jì)特征,并將含有隨機(jī)量的式(3)轉(zhuǎn)化為等效的確定性優(yōu)化模型。

2.1 CADET方法

一般連續(xù)時(shí)域非線性隨機(jī)系統(tǒng)可以表示為：

(4)

式中：x(t,s)是n維隨機(jī)狀態(tài)向量；w(t)是m維強(qiáng)迫函數(shù)向量,包含作用于系統(tǒng)的隨機(jī)擾動(dòng)s和控制輸入u；f(x,t)是x(t,s)的非線性矢量函數(shù)；G(t)是n×m的連續(xù)函數(shù)矩陣。

協(xié)方差分析描述函數(shù)法,就是先利用描述函數(shù)理論對(duì)非線性系統(tǒng)進(jìn)行統(tǒng)計(jì)線性化,然后對(duì)線性化后的系統(tǒng)進(jìn)行協(xié)方差分析。所謂統(tǒng)計(jì)線性化,是根據(jù)狀態(tài)量x(t,s)的概率密度函數(shù)形式,在較大的x(t,s)變化范圍內(nèi),用一個(gè)線性函數(shù)來逼近非線性函數(shù)f(x,t)。其本質(zhì)是對(duì)系統(tǒng)的每一個(gè)非線性用一個(gè)或幾個(gè)近似的、對(duì)輸入幅度靈敏的線性增益來代替[24]。這種方法的優(yōu)點(diǎn)是不要求f(x,t)連續(xù)可微,使許多不能用真線性化進(jìn)行線性化處理的實(shí)際函數(shù)有了線性化的可能[25]。

(5)

(6)

(7)

式中：b(t)是w(t)的確定性分量,Q(t)是w(t)隨機(jī)分量的譜密度矩陣。只需確定初始條件m(0)和p(0),利用式(7)進(jìn)行一次計(jì)算就能獲得任意時(shí)刻系統(tǒng)狀態(tài)量的統(tǒng)計(jì)特征,而且精度與幾千次的蒙特卡洛仿真相當(dāng),極大程度上節(jié)省了計(jì)算時(shí)間與成本。

由式(6)可知,描述函數(shù)的計(jì)算依賴于狀態(tài)量的概率密度函數(shù),一般來說十分復(fù)雜,甚至不能求得解析式。工程上常假定x(t,s)服從聯(lián)合正態(tài)分布,從而描述函數(shù)的計(jì)算可簡(jiǎn)化為

(8)

這里假設(shè)隨機(jī)狀態(tài)量x(t,s)服從正態(tài)分布,只有對(duì)高斯輸入的線性系統(tǒng)才是嚴(yán)格成立的,而對(duì)非高斯輸入的線性和非線性系統(tǒng)往往是近似的。對(duì)于非線性系統(tǒng),盡管輸入是高斯的,輸出一般來說是非高斯的。由中心極限定理可知：當(dāng)經(jīng)過低通濾波器后,隨機(jī)過程趨于高斯過程。因此,當(dāng)信號(hào)通過系統(tǒng)傳遞時(shí),可以借助系統(tǒng)的線性部分來保證非高斯的非線性輸出結(jié)果近似是正態(tài)的[24]。前述高斯假設(shè)的有效性已經(jīng)得到了廣泛的研究和檢驗(yàn)。

在附錄A中,給出了采用CADET方法得到的再入軌跡非線性動(dòng)力學(xué)模型的統(tǒng)計(jì)線性化模型。

2.2 軌跡魯棒優(yōu)化的等效確定性優(yōu)化模型

式(3)所示的軌跡魯棒優(yōu)化模型盡管表達(dá)了均值、標(biāo)準(zhǔn)差等引入隨機(jī)量后的性能信息,但式(3b)仍是關(guān)于x(t,s)的隨機(jī)微分方程,含有隨機(jī)量s,且其余各式中關(guān)于均值和標(biāo)準(zhǔn)差的求解方法也沒有給出,式(3)是無法用于實(shí)際計(jì)算的?，F(xiàn)在應(yīng)用CADET,將式(3)中含有的隨機(jī)變量s消去,得到等效的確定性優(yōu)化模型,具體步驟為將式(3b)轉(zhuǎn)化為關(guān)于x(t,s)均值和協(xié)方差的常微分方程,同時(shí)將式(3)中目標(biāo)函數(shù)、過程約束和端點(diǎn)約束表示為關(guān)于x(t,s)均值和協(xié)方差的函數(shù),從而得到等效的確定性軌跡魯棒優(yōu)化模型：將u(t)作為控制量,使得

minJμ(m(t),p(t),u(t),t)+kJ·Jσ(m(t),

p(t),u(t),t)

(9a)

(9b)

(9c)

Cμ(m(t),p(t),u(t),t)+kC·Cσ(m(t),

p(t),u(t),t)≤0

(9d)

Hμ(m(t),p(t),u(t),t)=0

(9e)

Hσ(m(t),p(t),u(t),t)≤εH

(9f)

式中：m(t)和p(t)分別是隨機(jī)狀態(tài)量x(t,s)的均值和協(xié)方差,根據(jù)CADET中的式(7)計(jì)算；(·)μ和(·)σ分別表示用于計(jì)算(·)的均值和標(biāo)準(zhǔn)差的函數(shù)。

對(duì)比式(9)和式(3),可以發(fā)現(xiàn),采用CADET后,含有隨機(jī)量s的魯棒優(yōu)化模型式(3)轉(zhuǎn)化為不含隨機(jī)變量的等效確定性優(yōu)化模型式(9),應(yīng)用現(xiàn)有的軌跡優(yōu)化算法完全可以求解。換言之,通過求解確定性優(yōu)化問題式(9),就能得到不確定條件下,魯棒優(yōu)化意義下的最優(yōu)控制量和最優(yōu)軌跡。

3 算例與分析

3.1 確定性再入軌跡優(yōu)化模型

文獻(xiàn)[26]中的再入軌跡優(yōu)化問題是一例經(jīng)典的軌跡優(yōu)化問題,希望在滿足終端和過程約束的條件下,最大化終點(diǎn)的橫向距離(即緯度),在確定性條件下可以描述為：將α(t)和σ(t)作為控制量,使得

maxJ=φ(tf)

(10a)

(10b)

(10c)

(10d)

式中：r是到地球球心的距離；θ是經(jīng)度；φ是緯度；v是速度；γ是航跡角；ψ是航向角；α是攻角；σ是傾側(cè)角；m是質(zhì)量；μ是地心引力常數(shù)；Re是地球半徑；D是阻力加速度；L是升力加速度；q是熱流密度。式(10c)表示端點(diǎn)約束,式(10d)表示過程約束。關(guān)于該模型中各參數(shù)的詳細(xì)說明可參見文獻(xiàn)[26],這里不再詳述。

進(jìn)一步,阻力加速度和升力加速度表示為：

(11)

式中：CD和CL分別是阻力系數(shù)和升力系數(shù)；ρ是大氣密度S是特征面積；h是高度；H是密度尺度高；ρ0是海平面大氣密度。

在傳統(tǒng)飛行軌跡優(yōu)化中,氣動(dòng)參數(shù)CD·和CL·視為確定值,根據(jù)實(shí)驗(yàn)得到的氣動(dòng)數(shù)據(jù)求得。

3.2 再入軌跡魯棒優(yōu)化模型

由于制造工藝誤差、測(cè)量誤差等各方面因素,將氣動(dòng)參數(shù)CD·和CL·視為不確定性量更為合理。本文中假設(shè)：

(12)

式中：(·)N表示(·)的參考值；s1,s2表示隨機(jī)變量。本文仿真實(shí)驗(yàn)中,s1～N(0.01,0.012),s2～N(-0.01,0.012)且s1,s2相互獨(dú)立。需要說明的是,這里s1,s2的均值不為0,目的是使干擾條件更苛刻,以顯示方法的有效性。關(guān)于s1,s2的概率分布規(guī)律等因素本身也是值得研究的問題,本研究著眼于再入軌跡魯棒優(yōu)化方法,因此在這里對(duì)此不作過多討論。

引入隨機(jī)變量s1,s2后,確定性軌跡優(yōu)化問題式(10)轉(zhuǎn)化為如式(2)所示的隨機(jī)軌跡優(yōu)化問題,在此基礎(chǔ)上進(jìn)一步轉(zhuǎn)化為式(3)所示的軌跡魯棒優(yōu)化問題,具體如下：

式(10)中的目標(biāo)函數(shù)式(10a)變?yōu)椋?/p>

maxJ=μ(φ(tf,s))-kJ·σ(φ(tf,s))

(13)

動(dòng)力學(xué)方程式(10b)轉(zhuǎn)化為含隨機(jī)變量的動(dòng)力學(xué)方程：

(14)

式中：x(t,s)表示隨機(jī)狀態(tài)向量；f表示式(10b)的右端項(xiàng)函數(shù)。

相應(yīng)地,式(10)中的確定性端點(diǎn)約束式(10c)和過程約束式(10d)轉(zhuǎn)化為用均值和標(biāo)準(zhǔn)差描述的統(tǒng)計(jì)意義上的約束。其中,等式約束式(10c)轉(zhuǎn)化為由均值和標(biāo)準(zhǔn)差表示的兩組約束式(15a)和式(15b),不等式約束(10d)轉(zhuǎn)化為式(15c)：

(15a)

(15b)

(15c)

需要說明：式(15b)僅對(duì)原確定性問題中有終端約束的狀態(tài)變量施加了標(biāo)準(zhǔn)差約束；所考慮的氣動(dòng)不確定性不會(huì)影響到初始狀態(tài),因此式(15b)中沒有設(shè)置對(duì)初值標(biāo)準(zhǔn)差的約束；此外,用來限制離散程度的標(biāo)準(zhǔn)差容許值和權(quán)重系數(shù)由設(shè)計(jì)者根據(jù)實(shí)際情況設(shè)定,本研究中設(shè)置如下：εhf=4500 m,εvf=250 m/s,εγf=3°,kJ=3,kCr=3,kCθ=3,kCv=3,kCγ=3,kCq=3。

綜上,再入軌跡魯棒優(yōu)化模型(RO)表示為：將α(t)和σ(t)作為控制量,使得

(16)

至此,已構(gòu)建起2個(gè)再入軌跡優(yōu)化問題,即傳統(tǒng)確定性優(yōu)化問題DO(式(10))以及同時(shí)考慮均值和標(biāo)準(zhǔn)差約束的魯棒優(yōu)化問題RO(式(16))。在下節(jié)中,將分別求解這兩個(gè)優(yōu)化問題,并對(duì)結(jié)果進(jìn)行對(duì)比分析。

3.3 計(jì)算結(jié)果及討論

魯棒優(yōu)化模型式(16)含有式(14),而式(14)是含有隨機(jī)變量s的隨機(jī)飛行動(dòng)力學(xué)方程,在實(shí)際計(jì)算求解式(16)時(shí),需要采用CADET方法將式(14)轉(zhuǎn)化為形如式(7)的關(guān)于均值和協(xié)方差的確定性動(dòng)力學(xué)方程,具體的轉(zhuǎn)化結(jié)果可參考附錄A。

3.3.1最優(yōu)控制量的對(duì)比分析

采用偽譜法求解這兩個(gè)軌跡優(yōu)化問題(關(guān)于偽譜法的相關(guān)論述這里不再展開,可參看相關(guān)文獻(xiàn)[27-29]),得到傳統(tǒng)優(yōu)化DO的最優(yōu)控制量：αDO(t)和σDO(t),魯棒優(yōu)化RO的最優(yōu)控制量：αRO(t)和σRO(t),如圖1所示,為了更清楚地顯示,兩個(gè)小方框給出了相應(yīng)時(shí)間段的放大圖。

圖1 最優(yōu)控制量(方框內(nèi)為放大圖)Fig.1 Optimal control

觀察圖1中控制量的時(shí)間歷程曲線,可得如下信息：

魯棒優(yōu)化最優(yōu)控制量αRO(t)、σRO(t),分別與對(duì)應(yīng)的確定性優(yōu)化的最優(yōu)控制量αDO(t)、σDO(t)的時(shí)間歷程趨勢(shì)基本一致,但αDO(t)和σDO(t)相對(duì)平緩,而αRO(t)和σRO(t)的波動(dòng)較大,且越接近終端時(shí)刻波動(dòng)越明顯。

αRO(t)和σRO(t)所需的時(shí)間比αDO(t)和σDO(t)的更長(zhǎng)(RO為2289.4 s,DO為2200.1 s)。

上述現(xiàn)象表明,魯棒優(yōu)化為了降低不確定性,比傳統(tǒng)確定性優(yōu)化付出了更多的控制成本,這是符合客觀實(shí)際的。此外,需要說明的是,圖1顯示攻角在末端變化較快,這是圖形比例原因?qū)е碌?總的飛行時(shí)間較長(zhǎng)),攻角實(shí)際最大變化率絕對(duì)值為3.43°/s,工程上是完全可以實(shí)現(xiàn)的。

3.3.2目標(biāo)和約束的魯棒性對(duì)比分析

為驗(yàn)證兩種優(yōu)化模型得到的最優(yōu)控制量在不確定因素影響下的效果,將αD0(t)和σD0(t)、αRO(t)和σRO(t)分別代入式(14)所示的含隨機(jī)變量的動(dòng)力學(xué)方程中,進(jìn)行10000次的MC仿真,并對(duì)仿真結(jié)果進(jìn)行對(duì)比分析,主要從三個(gè)方面進(jìn)行：目標(biāo)函數(shù)的魯棒性對(duì)比；終端約束的魯棒性對(duì)比；過程約束的魯棒性對(duì)比。

首先給出相應(yīng)的統(tǒng)計(jì)結(jié)果：

目標(biāo)函數(shù)終端時(shí)刻緯度的統(tǒng)計(jì)特性及散布情況如表1和圖2所示；終端約束終端時(shí)刻高度、速度和航跡角的統(tǒng)計(jì)特性及散布情況如表2、表3和圖3所示；過程約束熱流密度約束的情況如圖4和表4所示。圖2和圖3中橫坐標(biāo)表示MC仿真次數(shù),對(duì)圖表的解釋及分析見后。

表1 終端時(shí)刻緯度的統(tǒng)計(jì)特性Table 1 Statistical property of latitude at terminal time

圖2 終端時(shí)刻緯度的散布情況Fig.2 Dispersion of latitude at terminal time

表2 終端時(shí)刻高度、速度和航跡角的均值和相對(duì)誤差Table 2 Mean and relative error of height,speed and flight path angle at terminal time

表3 終端時(shí)刻高度、速度和航跡角的標(biāo)準(zhǔn)差Table 3 Standard deviation of height,speed and flight path angle at terminal time

圖3 終端時(shí)刻高度、速度和航跡角的散布情況Fig.3 Dispersion of height,speed and flight path angle at terminal time

圖4 熱流密度的散布情況Fig.4 Dispersion of heat flux density

通過分析上述圖2到圖4中的結(jié)果及表1到表4中的數(shù)據(jù),可以得到以下結(jié)論：

表4 熱流密度約束的滿足情況Table 4 Satisfaction of heat flux density constraint

1)目標(biāo)函數(shù)結(jié)果對(duì)比

表1中第一、二行數(shù)據(jù)顯示,對(duì)于目標(biāo)函數(shù)φ(tf)的均值,DO大于RO,而對(duì)于φ(tf)的標(biāo)準(zhǔn)差,RO小于DO。

φ(tf)的均值越大而標(biāo)準(zhǔn)差越小越好,因?yàn)闃O大化φ(tf)是設(shè)計(jì)目標(biāo),而較小的標(biāo)準(zhǔn)差則意味著抗干擾能力較強(qiáng)。因此上述現(xiàn)象表明,RO為提高系統(tǒng)對(duì)不確定性因素的魯棒性,犧牲掉一部分φ(tf)的均值,這是符合設(shè)計(jì)預(yù)期的。

為了顯示出均值和標(biāo)準(zhǔn)差的綜合效果,表1中第三行數(shù)據(jù)顯示了均值和標(biāo)準(zhǔn)差的加權(quán)和,kJ取3時(shí),RO表現(xiàn)優(yōu)于DO。

圖2顯示了氣動(dòng)不確定條件下,終端時(shí)刻緯度的散布情況,可以直觀地看出RO中的樣本點(diǎn)分布更為緊密,表明RO的目標(biāo)函數(shù)對(duì)不確定性因素的敏感程度更低。

2)終端約束結(jié)果對(duì)比

表2中數(shù)據(jù)顯示,在氣動(dòng)不確定因素影響下,DO的終端高度、速度和航跡角的均值與約束要求的值有較大偏差,最大達(dá)到51.40%。而RO的終端狀態(tài)均值與終端約束要求的值基本吻合。這說明,在氣動(dòng)不確定因素影響下RO能使終端狀態(tài)以約束要求值為中心散布,而DO的終端狀態(tài)散布偏離了約束要求值。

表3中數(shù)據(jù)顯示,RO中的高度和航跡角的標(biāo)準(zhǔn)差小于DO中對(duì)應(yīng)量的標(biāo)準(zhǔn)差,這說明,在氣動(dòng)不確定影響下,RO中的高度和航跡角的分布更緊密。此外,盡管DO中速度的標(biāo)準(zhǔn)差小于RO中速度的標(biāo)準(zhǔn)差,但由于DO中速度的均值偏離了要求的均值,較小的標(biāo)準(zhǔn)差并未帶來實(shí)際的收益,這在圖3中得到了體現(xiàn),具體如下。

圖3表明了氣動(dòng)不確定條件下,高度、速度和航跡角終端狀態(tài)的散布情況,紅實(shí)線是設(shè)計(jì)要求,表示均值,上、下紅虛線表示均值加減標(biāo)準(zhǔn)差,則落入上、下虛線之間部分的樣本點(diǎn)越多,說明抗干擾能力越強(qiáng),魯棒性越好。DO的10000個(gè)樣本點(diǎn)中分別有5937(高度)、6599(速度)和6310(航跡角)個(gè)在這個(gè)范圍內(nèi),占比為59.37%、65.99%和63.10%；而RO的10000個(gè)樣本點(diǎn)中,分別有6883(高度)、6766(速度)和8267(航跡角)個(gè)在這個(gè)范圍內(nèi),占比為68.83%、67.66%和82.67%。顯然,RO所得到的結(jié)果比DO的魯棒性更好。

3)過程約束結(jié)果對(duì)比

本算例中共有7個(gè)過程約束,其中熱流密度的過程約束是主動(dòng)約束,因此這里以熱流密度過程約束來討論不確定因素的影響。

圖4a和圖4b分別顯示了DO和RO中q(t)時(shí)間歷程樣本的分布情況,紅線代表飛行過程中熱流密度的約束上限。對(duì)比圖4a和圖4b,直觀地顯示出DO中大量熱流密度時(shí)間歷程樣本超出了約束限制,而在RO中,絕大多數(shù)樣本都是滿足約束的。表4中給出了具體值,在DO的10000條q(t)時(shí)間歷程樣本中,有17條完全滿足約束,占比為0.17%；而在RO的10000條q(t)時(shí)間歷程樣本中,有9933條完全滿足約束,占比為99.33%。上述結(jié)果表明,在氣動(dòng)不確定條件下,RO中的熱流密度過程約束有極強(qiáng)的抗干擾能力,而DO中的熱流密度過程約束的抗干擾能力很弱。

總的來說,在不確定性因素的影響下,再入軌跡魯棒優(yōu)化的最優(yōu)軌跡相較于確定性優(yōu)化的最優(yōu)軌跡表現(xiàn)出更強(qiáng)的抗干擾能力,而對(duì)于不同的物理性能指標(biāo),抗干擾的能力有所不同,在本算例中,對(duì)熱流密度的抗干擾能力最強(qiáng)。

4 總結(jié)與展望

傳統(tǒng)的再入軌跡優(yōu)化是在確定性條件下展開的,本研究將不確定性參數(shù)引入再入軌跡優(yōu)化中,構(gòu)建了再入軌跡魯棒優(yōu)化模型,使得到的飛行軌跡具有一定的抗干擾能力。將CADET作為不確定性分析手段,引入再入軌跡魯棒優(yōu)化的實(shí)際計(jì)算中,實(shí)現(xiàn)了再入軌跡魯棒優(yōu)化的高效計(jì)算。仿真結(jié)果表明,相比于確定性軌跡優(yōu)化,本文方法得到的最優(yōu)軌跡在不確定性因素的干擾下具有更強(qiáng)的魯棒性。

本研究的著眼點(diǎn)是再入軌跡魯棒優(yōu)化方法,而軌跡優(yōu)化中的不確定性來源眾多,將來還應(yīng)對(duì)輸入不確定性的合理建模展開深入研究。此外,協(xié)方差分析法的精度會(huì)隨著系統(tǒng)非線性程度的增加而降低,對(duì)于非高斯輸入的有效性也有待考察,因此準(zhǔn)確高效的不確定性傳播方法依然是未來關(guān)注的重點(diǎn)。