運(yùn)載火箭推力故障不確定性下軌跡可達(dá)包絡(luò)分析

李超兵，包為民，李忠奎，程曉明，魏才盛

(1.北京大學(xué)工學(xué)院，北京 100871；2.中國航天科技集團(tuán)有限公司，北京 100048；3.北京航天自動控制研究所，北京 100854；4.中南大學(xué)自動化學(xué)院，長沙 410083)

0 引 言

運(yùn)載火箭是執(zhí)行航天任務(wù)的重要基礎(chǔ)和保障。隨著我國航天強(qiáng)國建設(shè)的步伐逐漸加快，對運(yùn)載火箭的可靠性和智能化提出了更高的要求[1-2]。不同類型、不同級別的故障對任務(wù)的影響程度不同。動力系統(tǒng)故障是運(yùn)載火箭最容易遭受到的故障，輕則導(dǎo)致火箭入軌精度降低，重則直接導(dǎo)致飛行任務(wù)的失敗[3]。隨著運(yùn)載火箭在線計算能力的提高和火箭智能化相關(guān)技術(shù)的發(fā)展，當(dāng)動力系統(tǒng)出現(xiàn)故障時，有望通過在線故障診斷和軌跡重規(guī)劃，在線調(diào)整飛行程序，完成或部分完成既定任務(wù)。

近幾年，許多學(xué)者面向推力故障下軌跡重規(guī)劃問題開展研究。譚述君等[4]通過分別求解故障下面向目標(biāo)軌道和救援軌道的離線重規(guī)劃問題，構(gòu)建了離線樣本集，在線通過調(diào)用樣本集完成實(shí)時軌跡重規(guī)劃。馬昊磊等[6]提出一種基于自適應(yīng)偽譜法的故障火箭軌跡重規(guī)劃方法。Li等[7]針對運(yùn)載火箭突發(fā)電力故障問題，考慮規(guī)劃方法的實(shí)時性，提出了一種基于凸優(yōu)化的軌跡重規(guī)劃方法。Song等[8]通過構(gòu)造狀態(tài)相關(guān)的評價指標(biāo)，使得推力下降故障發(fā)生后，運(yùn)載火箭能夠自主決策和規(guī)劃到目標(biāo)軌道的軌跡。李遠(yuǎn)東等[9]針對可重復(fù)利用運(yùn)載火箭的推力下降問題，提出了基于粒子群算法和優(yōu)化ZEM/ZEV制導(dǎo)律的軌跡優(yōu)化方案。Ma等[10]針對運(yùn)載火箭突發(fā)的推力下降問題，提出了一種改進(jìn)的并行牛頓制導(dǎo)算法，能夠有效解決多約束下的制導(dǎo)問題。He等[11]基于凸優(yōu)化方法和深度神經(jīng)網(wǎng)絡(luò)離線計算和構(gòu)建故障狀態(tài)到最優(yōu)救援軌道的數(shù)據(jù)庫，在線突發(fā)故障時可快速獲得最優(yōu)重規(guī)劃救援軌道。Diao等[12]提出了推力故障運(yùn)載火箭的在軌自主決策方法和基于目標(biāo)自適應(yīng)更新的迭代制導(dǎo)方法。值得注意的是，上述研究雖然能夠針對特定故障，在線獲得重規(guī)劃的目標(biāo)軌道或救援軌道，但并未考慮故障存在不確定性的實(shí)際工況。當(dāng)運(yùn)載火箭突發(fā)推力故障時，故障信息由故障診斷子系統(tǒng)估計獲得，存在較強(qiáng)的不確定性。定量評估運(yùn)載火箭推力故障不確定性的影響，分析不確定性下軌跡可達(dá)包絡(luò)，對于推力故障下運(yùn)載火箭的能力分析、在線自主決策和軌跡重規(guī)劃具有重要的參考價值。

盡管對故障不確定性影響的研究較少，但已有部分研究開展了環(huán)境、測量、執(zhí)行等環(huán)節(jié)不確定性的影響分析，對本文的研究具有重要參考價值。Luo等[13]綜述了當(dāng)前主流的不確定性影響分析方法，包含蒙特卡洛仿真法、局部線性化法、混沌多項(xiàng)式法、狀態(tài)轉(zhuǎn)移張量法、微分代數(shù)法等。Jones等[14]針對軌道不確定性的傳播過程，基于混沌多項(xiàng)式法進(jìn)行了求解。Valli等[15]基于微分代數(shù)法進(jìn)行了不確定性的傳播分析，以較小的計算代價實(shí)現(xiàn)了較高的估計精度。Yang等[16]考慮了航天器推力器的執(zhí)行不確定性，基于狀態(tài)轉(zhuǎn)移張量法分析了不確定性對軌跡可達(dá)包絡(luò)的影響。Jin等[17]面向航天器的交會問題，考慮了環(huán)境不確定性、測量不確定性、執(zhí)行不確定性的影響，基于狀態(tài)轉(zhuǎn)移張量法構(gòu)建了不確定性的傳播和更新方程。在不確定性影響分析方法中，狀態(tài)轉(zhuǎn)移張量法由于不需要對不確定性進(jìn)行隨機(jī)采樣[18]，因此非常適合用于推力故障運(yùn)載火箭這類復(fù)雜動力學(xué)系統(tǒng)的不確定性分析。

本文針對突發(fā)推力故障的運(yùn)載火箭軌跡重規(guī)劃問題，考慮故障和環(huán)境存在不確定性，分析了運(yùn)載火箭的軌跡可達(dá)包絡(luò)。首先，建立了推力故障火箭軌跡重規(guī)劃問題模型，并利用序列二次規(guī)劃法進(jìn)行求解。然后，通過對不確定性傳播的建模，提出了基于狀態(tài)轉(zhuǎn)移張量法的運(yùn)載火箭軌跡可達(dá)包絡(luò)分析方法。最后，通過與傳統(tǒng)蒙特卡洛仿真方法的對比，驗(yàn)證了本文方法的有效性和優(yōu)勢，有望為推力故障運(yùn)載火箭的能力分析、自主決策和軌跡重規(guī)劃提供技術(shù)支撐。

1 推力故障火箭運(yùn)動模型與軌跡優(yōu)化

運(yùn)載火箭在地心慣性系下關(guān)于時間t的運(yùn)動方程為：

(1)

式中：r為火箭位置矢量，v為火箭速度矢量，μ為萬有引力常數(shù)，T為推力大小，u為推力單位矢量，m為運(yùn)載火箭的當(dāng)前質(zhì)量，D為氣動力矢量，g0為引力加速度，Isp為火箭發(fā)動機(jī)的比沖。

運(yùn)載火箭上升過程中可能遭遇多類型故障。本文選擇發(fā)生頻率較高的推力下降故障進(jìn)行分析?；鸺l(fā)生推力下降故障后，真實(shí)推力大小為：

Treal=k0T

(2)

式中：k0為推力下降比例系數(shù)。

當(dāng)運(yùn)載火箭突發(fā)故障后，執(zhí)行原定推力序列無法使火箭到達(dá)設(shè)計軌道的期望位置r(tf)和期望速度v(tf)，其中tf是入軌時間。為使推力故障對飛行任務(wù)的影響最小，需進(jìn)行運(yùn)載火箭軌跡的在線重規(guī)劃。本文假設(shè)運(yùn)載火箭在t0時刻突發(fā)故障，以燃料最省(到達(dá)設(shè)計軌道后火箭剩余質(zhì)量最大)為優(yōu)化目標(biāo)，考慮動力學(xué)約束、起始與終端條件約束、推力約束等，建立如下推力故障火箭軌跡重規(guī)劃問題模型：

(3)

式中：J為燃料最省性能指標(biāo)，Re表示地球半徑。

推力故障運(yùn)載火箭的軌跡重規(guī)劃問題可通過求解最優(yōu)問題的一般方法，如序列二次規(guī)劃[19]、凸優(yōu)化[20]等方法進(jìn)行求解，本文采用序列二次規(guī)劃法對問題進(jìn)行求解。

值得注意的是，軌跡重規(guī)劃問題假設(shè)系統(tǒng)是完全已知的，忽略了環(huán)境不確定性(如氣動不確定性)和故障不確定性對運(yùn)動軌跡的影響。特別地，當(dāng)發(fā)生推力下降故障后，火箭控制系統(tǒng)依賴故障診斷模塊獲得估計的推力下降比例系數(shù)k，當(dāng)k存在估計不確定性時，小范圍的推力大小變化將引起較大范圍的軌跡波動，使得重規(guī)劃軌跡失效。本文將基于狀態(tài)轉(zhuǎn)移張量法分析火箭軌跡優(yōu)化受不確定性的影響，獲得火箭軌跡的可達(dá)包絡(luò)。

2 系統(tǒng)不確定性建模

包含不確定性因素下運(yùn)載火箭的開環(huán)動態(tài)系統(tǒng)方程可以建模為：

(4)

E[w(t)wT(t′)]=Sw(t)δ(t-t′)

(5)

式中:δ(·)為Kroneckerδ函數(shù)，Sw(t)為相應(yīng)的協(xié)方差矩陣。

系統(tǒng)的控制量輸入依據(jù)導(dǎo)航狀態(tài)計算得到，在該運(yùn)載火箭的動力學(xué)模型中的物理意義為火箭推力，廣義表達(dá)式為：

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

E[η(t)ηT(t′)]=Sη(t)δ(t-t′)

(16)

式中:Sη(t)為相應(yīng)的協(xié)方差矩陣。

在突發(fā)推力故障工況下，雖然可以通過建立和求解在線重規(guī)劃問題獲得優(yōu)化軌跡，但軌跡的重規(guī)劃過程未考慮不確定性因素的影響。當(dāng)發(fā)生推力下降故障后，火箭控制系統(tǒng)依賴故障診斷模塊獲得估計的推力下降比例系數(shù)k0,k0的估計值往往存在較強(qiáng)的不確定性。此外，大氣密度ρ往往也存在不確定性，進(jìn)而影響為氣動力矢量D。最后，在故障發(fā)生時刻運(yùn)載火箭所處位置狀態(tài)也會存在一定不確定偏差，該偏差在運(yùn)動過程中會隨著時間累積而對運(yùn)載火箭整體狀態(tài)的計算產(chǎn)生一定影響。

為表示出系統(tǒng)動態(tài)方程(4)在考慮此類不確定性下的具體形式，需對系統(tǒng)的狀態(tài)增廣處理，也即定義運(yùn)載火箭運(yùn)動系統(tǒng)的真實(shí)狀態(tài)為：

x=[rT,vT,cT,eT,dT,p]T

(17)

結(jié)合式(1)和式(4)，在考慮系統(tǒng)不確定因素作用下真實(shí)狀態(tài)的動態(tài)方程具體可表示為：

(18)

(19)

vrel=(v+c)-ω×(r+e)

(20)

式中：Cd為阻力系數(shù),Aref為參考面積，c∈R3,e∈R3,d∈R3和p∈R分別表示作用在遙測速度、遙測位置、推力加速度及大氣密度上的不確定因素；ω×為地球自轉(zhuǎn)角速度矢量ω的反對稱矩陣。這些不確定性變量的動態(tài)過程滿足一階馬爾克夫過程，即某時刻的狀態(tài)只與上一時刻的狀態(tài)有關(guān)。τc,τe,τd及τp分別對應(yīng)偏差在整個系統(tǒng)中的時間常數(shù)。并且系統(tǒng)加入的不確定偏差滿足高斯噪聲的性質(zhì)：

(21)

(22)

(23)

(24)

式中:σj(j=c,e,d,p)為對應(yīng)不確定因素的穩(wěn)定方差。

此外定義火箭運(yùn)動系統(tǒng)的導(dǎo)航狀態(tài)中包含火箭的位置、速度以及大氣密度，具體為：

(25)

(26)

(27)

系統(tǒng)選擇的測量狀態(tài)為地心慣性系下的位置信息。也即：

(28)

至此建立了推力故障后運(yùn)載火箭在不確定性因素作用下的系統(tǒng)狀態(tài)動態(tài)方程。系統(tǒng)狀態(tài)偏差的求解特別是運(yùn)載火箭軌跡可達(dá)包絡(luò)的求解還需對本章建立的不確定性模型作進(jìn)一步處理。

3 運(yùn)載火箭軌跡可達(dá)包絡(luò)分析

(29)

(30)

(31)

(32)

式中：

(33)

(34)

(35)

同樣地，對狀態(tài)更新中的狀態(tài)偏差量，即式(8)進(jìn)行線性化可以得到：

(36)

(37)

(38)

可以得到增廣狀態(tài)的狀態(tài)方程為：

(39)

(40)

式中：

(41)

基于式(39)和(40)中給出的線性化后的系統(tǒng)狀態(tài)方程，由狀態(tài)轉(zhuǎn)移張量法計算系統(tǒng)的線性協(xié)方差及其更新方程為：

CA=E[X(t)XT(t)]

(42)

(43)

(44)

考慮到系統(tǒng)狀態(tài)在模型不確定因素作用下的誤差傳遞狀況可以通過系統(tǒng)真實(shí)狀態(tài)的協(xié)方差矩陣來分析，因此提取協(xié)方差矩陣為：

(45)

基于系統(tǒng)狀態(tài)偏差的協(xié)方差矩陣Dtrue，可以定量評估狀態(tài)偏離期望狀態(tài)的軌跡包絡(luò)范圍。

基于狀態(tài)轉(zhuǎn)移張量法的推力故障運(yùn)載火箭軌跡可達(dá)包絡(luò)分析算法如圖1所示。從輸入推力故障后重構(gòu)的優(yōu)化軌跡開始，在依據(jù)系統(tǒng)初始條件下確定的不確定性因素后建立對應(yīng)狀態(tài)偏差的協(xié)方差微分方程，最后利用數(shù)值計算的方法即可求解各個時刻下的狀態(tài)偏差量。通過運(yùn)載火箭狀態(tài)與軌道位置的對應(yīng)關(guān)系即可求出推力故障下不確定性作用的可達(dá)域包絡(luò)。

圖1 推力故障運(yùn)載火箭軌跡可達(dá)包絡(luò)分析算法Fig.1 Reachable envelope analysis algorithm for thrust fault launch vehicle trajectory

4 仿真校驗(yàn)

為驗(yàn)證本文提出的推力故障、氣動以及位置信息等不確定下運(yùn)載火箭軌跡可達(dá)包絡(luò)分析方法的有效性,本節(jié)組織了兩種不同推力故障及不確定因素作用下運(yùn)載火箭軌跡包絡(luò)優(yōu)化計算算例，其中推力故障包括小推力故障損失以及較大的推力故障損失。運(yùn)載火箭仿真相關(guān)參數(shù)如表1所示。

表1 運(yùn)載火箭數(shù)據(jù)Table 1 Launch vehicle data

4.1 故障類型1計算仿真

故障類型1的仿真場景假定運(yùn)載火箭在時間t0=200 s時突發(fā)推力下降故障，推力下降比例系數(shù)k0=70%，并且主要考慮推力及大氣密度的不確定性因素?；谛蛄卸我?guī)劃方法求解式(3)中推力故障火箭軌跡重規(guī)劃問題，可獲得故障后的重規(guī)劃軌跡，如圖2所示。重規(guī)劃軌跡與標(biāo)稱軌跡的入軌軌道根數(shù)對比如表2所示。從圖2和表2的仿真結(jié)果可以得出，入軌在t0=200 s發(fā)生推力損失30%，經(jīng)過軌跡重規(guī)劃，重規(guī)劃軌跡終點(diǎn)的軌道根數(shù)盡管與標(biāo)稱軌跡的軌道根數(shù)存在一定的誤差，但基本能夠滿足飛行任務(wù)的軌道要求。

圖2 故障類型1重規(guī)劃軌跡與標(biāo)稱軌跡對比圖Fig.2 Comparison between the replanning and nominal trajectories after fault type 1

表2 重規(guī)劃軌跡與標(biāo)稱軌跡入軌軌道根數(shù)對比Table 2 Terminal orbital elements comparison between the replanned trajectory and the nominal trajectory

為進(jìn)一步分析故障發(fā)生后不確定性的影響，仿真中引入推力不確定性和氣動不確定性，相關(guān)參數(shù)如表3所示。

表3 故障類型1不確定性參數(shù)設(shè)置Table 3 Uncertainty parameter setting of fault type 1

基于本文提出的推力故障運(yùn)載火箭軌跡可達(dá)包絡(luò)算法，可以獲得故障發(fā)生后運(yùn)載火箭的狀態(tài)偏移包絡(luò)(3σ包絡(luò))，如圖3到圖5所示。其中圖3給出了運(yùn)載火箭位置矢量的包絡(luò)分量，圖4給出了運(yùn)載火箭速度矢量的包絡(luò)分量，圖5給出了推力故障火箭重規(guī)劃軌跡的局部軌跡偏差范圍。

圖3 運(yùn)載火箭位置偏差隨時間變化曲線Fig.3 Deviation curve of launch vehicle position with time

圖4 運(yùn)載火箭速度偏差隨時間變化曲線Fig.4 Deviation curve of launch vehicle velocity with time

圖5 推力故障火箭重規(guī)劃軌跡的局部軌跡偏差示意圖Fig.5 Schematic diagram of local trajectory deviation of rocket replanning trajectory with thrust fault

為驗(yàn)證所提方法的優(yōu)勢，本文基于表2中不確定性參數(shù)設(shè)置，隨機(jī)進(jìn)行了100次獨(dú)立的蒙特卡洛仿真，仿真結(jié)果同樣繪制于圖3到圖4中以方便進(jìn)行對比分析。

結(jié)合該故障類型下的重規(guī)劃結(jié)果以及計算出的火箭狀態(tài)偏差量，為了能夠更加直觀地觀察到系統(tǒng)狀態(tài)偏離值對火箭終端入軌能力的影響程度，將狀態(tài)偏差結(jié)果轉(zhuǎn)換為火箭入軌時的軌道偏差，其偏差結(jié)果如表4所示。

表4 重規(guī)劃軌跡入軌軌道根數(shù)偏差Table 4 Deviation of the terminal orbit elements in replanning trajectory

4.2 故障類型2計算仿真

在故障類型2的仿真場景中,假定運(yùn)載火箭在時間t0=500 s時突發(fā)較大的推力下降故障，推力下降比例系數(shù)k0=50%，并且在考慮了推力及大氣密度的不確定性外還加入了運(yùn)載火箭故障位置不確定性因素。故障后的重規(guī)劃軌跡如圖6所示,重規(guī)劃軌跡與標(biāo)稱軌跡的入軌軌道根數(shù)對比如表5所示。

圖6 故障類型2重規(guī)劃軌跡與標(biāo)稱軌跡對比圖Fig.6 Comparison between the replanning and nominal trajectories after fault type 2

表5 重規(guī)劃軌跡與標(biāo)稱軌跡入軌軌道根數(shù)對比Table 5 Terminal orbital elements comparison between the replanned trajectory and the nominal trajectory

在該故障類型中添加了額外的運(yùn)載火箭火箭位置不確定偏差，各不確定性相關(guān)參數(shù)如表6所示。

表6 故障類型2不確定性參數(shù)設(shè)置Table 6 Uncertainty parameter setting of fault type 2

在此故障類型下同樣仿真計算出了運(yùn)載火箭在位置和速度上各分量的偏差分布狀況，計算仿真結(jié)果如圖7-9所示。

圖7 運(yùn)載火箭位置偏差隨時間變化曲線Fig.7 Deviation curve of launch vehicle position with time

圖8 運(yùn)載火箭速度偏差隨時間變化曲線Fig.8 Deviation curve of launch vehicle velocity with time

圖9 推力故障火箭重規(guī)劃軌跡的局部軌跡偏差示意圖Fig.9 Schematic diagram of local trajectory deviation of rocket replanning trajectory with thrust fault

結(jié)合該故障類型下的重規(guī)劃結(jié)果以及計算出的火箭狀態(tài)偏差量，將狀態(tài)偏差結(jié)果轉(zhuǎn)換為火箭入軌時的軌道根數(shù)偏差，其偏差結(jié)果如表7所示。

表7 重規(guī)劃軌跡入軌軌道根數(shù)偏差Table 7 Deviation of the terminal orbital elements in replanning trajectory

從以上2種故障類型的仿真計算結(jié)果中可以得出：本文提出的推力故障不確定性下運(yùn)載火箭軌跡可達(dá)包絡(luò)分析方法給出的位置矢量包絡(luò)邊界和速度矢量包絡(luò)邊界與100次蒙特卡洛仿真得到的偏差上界結(jié)果一致，驗(yàn)證了所提方法的有效性。區(qū)別于多次蒙特卡洛方法，本文方法能夠一次性給出運(yùn)載火箭軌跡可達(dá)包絡(luò)，對于后續(xù)的運(yùn)載火箭能力分析、故障發(fā)生后的決策等工作具有明確、直接、及時的參考價值。

5 結(jié) 論

針對突發(fā)推力故障的運(yùn)載火箭的軌跡優(yōu)化與設(shè)計問題，提出了一種基于狀態(tài)轉(zhuǎn)移張量法的故障火箭軌跡可達(dá)包絡(luò)分析方法。通過建立推力下降故障的運(yùn)載火箭運(yùn)動模型，在不考慮不確定因素影響下，采用序列二次規(guī)劃方法優(yōu)化節(jié)求解出推力故障下的重規(guī)劃軌跡?？紤]故障和環(huán)境不確定性，基于線性協(xié)方差對非線性系統(tǒng)的傳播過程進(jìn)行建模和模型線性化，并采用狀態(tài)轉(zhuǎn)移張量法實(shí)現(xiàn)了對運(yùn)載火箭軌跡可達(dá)包絡(luò)分析。仿真結(jié)果表面，所提方法能夠一次性給出運(yùn)載火箭軌跡可達(dá)包絡(luò)邊界，與100次獨(dú)立蒙特卡洛仿真得到的包絡(luò)邊界一致。論文所研方法有望為推力故障運(yùn)載火箭的能力分析、自主決策和軌跡重規(guī)劃提供技術(shù)支撐。