運(yùn)載火箭推力故障下軌跡在線重規(guī)劃方法*

高天域，何睿智，溫長(zhǎng)新，湯國(guó)建

(1.國(guó)防科技大學(xué) 空天科學(xué)學(xué)院空天工程系·長(zhǎng)沙·410073；2.中國(guó)人民解放軍95899部隊(duì)·北京·100000)

0 引 言

近年來，隨著各國(guó)運(yùn)載火箭技術(shù)不斷提升，運(yùn)載能力不斷增強(qiáng)，對(duì)外層空間的開發(fā)與利用也在不斷深化，運(yùn)載火箭需要更強(qiáng)的自主控制能力從而妥善處理突發(fā)故障[1]。發(fā)動(dòng)機(jī)是運(yùn)載火箭研制過程中的重中之重，其可靠性關(guān)乎整個(gè)飛行任務(wù)的成敗，然而新技術(shù)、新材料、新理論的應(yīng)用不可避免地造成其可靠性降低，據(jù)統(tǒng)計(jì)，運(yùn)載火箭主動(dòng)段動(dòng)力系統(tǒng)故障占所發(fā)故障的60%[2]。例如日本H-2運(yùn)載火箭由于一級(jí)發(fā)動(dòng)機(jī)提前關(guān)機(jī)，升空后嚴(yán)重偏離預(yù)定軌道，導(dǎo)致發(fā)射任務(wù)失敗，造成巨大經(jīng)濟(jì)損失；印度GSLV-F06運(yùn)載火箭發(fā)射后，2臺(tái)游動(dòng)發(fā)動(dòng)機(jī)點(diǎn)火失敗，導(dǎo)致任務(wù)失??；同樣，我國(guó)也曾發(fā)生此類事故，長(zhǎng)征五號(hào)遙二火箭因芯一級(jí)2臺(tái)發(fā)動(dòng)機(jī)其中一臺(tái)故障導(dǎo)致發(fā)射任務(wù)失敗，對(duì)我國(guó)后續(xù)航天任務(wù)造成影響。隨著自適應(yīng)制導(dǎo)和軌跡在線規(guī)劃技術(shù)的不斷成熟，發(fā)動(dòng)機(jī)非致命故障的情況下，仍能通過調(diào)整火箭飛行程序完成任務(wù)。SpaceX公司發(fā)射的Falcon 9運(yùn)載火箭曾在一臺(tái)發(fā)動(dòng)機(jī)因故障提前關(guān)機(jī)的情況下，通過增加一二級(jí)飛行時(shí)間，將主要載荷送入預(yù)定軌道，次要載荷進(jìn)入較低軌道，基本完成任務(wù)；美國(guó)德爾塔4運(yùn)載火箭曾利用制導(dǎo)控制系統(tǒng)重新生成故障后的飛行軌跡，對(duì)推力下降及時(shí)進(jìn)行了補(bǔ)償，充分利用剩余燃料完成入軌。

針對(duì)火箭故障條件下的軌跡在線重規(guī)劃和制導(dǎo)問題，國(guó)內(nèi)外開展了一定的研究。韓業(yè)鵬[3]針對(duì)大氣層外動(dòng)力故障，研究了入軌點(diǎn)更新迭代制導(dǎo)方法，提高了控制系統(tǒng)適應(yīng)故障的能力。王志祥等[4]以某型號(hào)捆綁火箭為研究對(duì)象，分析了單臺(tái)發(fā)動(dòng)機(jī)推力下降對(duì)火箭飛行軌跡和姿態(tài)的影響，并利用Gauss偽譜法展開仿真研究。宋征宇等[5]提出了一種在線軌跡規(guī)劃策略，判斷火箭不可達(dá)時(shí)自主尋找最優(yōu)救援軌道。李源[6]針對(duì)二級(jí)運(yùn)載火箭在不同飛行段發(fā)生不同程度的故障的軌跡規(guī)劃方法以及制導(dǎo)方法進(jìn)行了研究，為運(yùn)載火箭軌跡在線規(guī)劃和自適應(yīng)制導(dǎo)技術(shù)應(yīng)用提供了一定參考。H.Lee等[7]通過牛頓迭代法提升了火箭應(yīng)對(duì)機(jī)構(gòu)故障的能力。H.Sun[8]提出了一種基于有限微分的能夠適應(yīng)發(fā)動(dòng)機(jī)推力故障的上升段制導(dǎo)方法。B.Beneditkter等[9]采用無損凸化技術(shù)和連續(xù)凸化技術(shù)相結(jié)合的方式，并在優(yōu)化過程中明確考慮自由終端時(shí)間的要求，提出了一種多級(jí)運(yùn)載火箭主動(dòng)段凸規(guī)劃方法。G.A.Dukema[10]和A.J.Calise等[11]提出了一種基于預(yù)測(cè)-校正的初值生成方法，能夠處理故障后返回地面的問題。除此之外，離線設(shè)計(jì)應(yīng)急軌道的方法也是提高運(yùn)載火箭應(yīng)對(duì)故障能力的一種選擇，“挑戰(zhàn)者”號(hào)曾應(yīng)用離線設(shè)計(jì)好的應(yīng)急軌道，成功在一臺(tái)發(fā)動(dòng)機(jī)故障后進(jìn)入較低的安全軌道。A.L.Cowling[12]研究了不同時(shí)刻發(fā)動(dòng)機(jī)停止工作后的中止軌跡。隨著箭載計(jì)算機(jī)的計(jì)算能力不斷提高，未來基于自主制導(dǎo)方法的新型閉環(huán)控制系統(tǒng)架構(gòu)以及多級(jí)分層軌跡優(yōu)化將成為主要的研究方向[13-15]。

針對(duì)運(yùn)載火箭主動(dòng)段非入軌飛行段出現(xiàn)推力下降的問題，本文提出了一種基于能量最優(yōu)的凸優(yōu)化軌跡規(guī)劃方法，將原非凸問題通過線性化估算和常值假設(shè)轉(zhuǎn)為凸優(yōu)化問題，并通過迭代得到原問題的最優(yōu)解，實(shí)現(xiàn)該問題下軌跡快速求解。

1 數(shù)學(xué)模型描述

1.1 坐標(biāo)系定義

首先定義以下2個(gè)坐標(biāo)系：1)發(fā)射慣性坐標(biāo)系O-xyz，坐標(biāo)原點(diǎn)為發(fā)射點(diǎn)，Ox軸在火箭發(fā)射點(diǎn)水平面內(nèi)指向發(fā)射瞄準(zhǔn)方向，Oy軸垂直于發(fā)射點(diǎn)水平指向上方，Oz軸滿足右手坐標(biāo)系；2)近焦點(diǎn)坐標(biāo)系Op-xpypzp，坐標(biāo)原點(diǎn)為地心，xp軸為地心和目標(biāo)軌道近地點(diǎn)的連線，指向近地點(diǎn)，yp為xp在目標(biāo)軌道平面內(nèi)沿軌道方向旋轉(zhuǎn)90°，zp滿足右手坐標(biāo)系。近焦點(diǎn)坐標(biāo)系如圖1所示。

圖1 近焦點(diǎn)坐標(biāo)系Fig.1 Perifocal coordinate system

1.2 動(dòng)力學(xué)模型

本文主要針對(duì)三級(jí)火箭的二級(jí)飛行段進(jìn)行研究。假設(shè)故障發(fā)生于一二級(jí)分離時(shí)刻，此時(shí)火箭高度已達(dá)100km左右，大氣稀薄，可忽略氣動(dòng)力的影響。假設(shè)地球?yàn)闄E球，此時(shí)在發(fā)射慣性坐標(biāo)系下建立運(yùn)載火箭運(yùn)動(dòng)方程

(1)

式中，向量(r,v,a)分別為火箭在慣性系中的位置、速度和加速度；T為火箭在慣性系中3個(gè)方向的推力分量；Treal為運(yùn)載火箭實(shí)際推力。一般情況下，火箭推力為常值，m為火箭質(zhì)量，隨時(shí)間變化，Isp為發(fā)動(dòng)機(jī)比沖，g為引力加速度，g0為海平面重力加速度。

1.3 發(fā)動(dòng)機(jī)推力故障模型

本文考慮的發(fā)動(dòng)機(jī)推力下降模式為：由于秒耗量下降發(fā)動(dòng)機(jī)推力迅速下降至某一固定值，且不考慮二次故障。表征故障的主要參數(shù)有故障發(fā)生時(shí)刻和發(fā)動(dòng)機(jī)故障后輸出推力的大小，以某型號(hào)芯二級(jí)4臺(tái)發(fā)動(dòng)機(jī)為例，各發(fā)動(dòng)機(jī)發(fā)生推力下降后的數(shù)學(xué)模型為

Ti=kT(i=1,2,3,4)

(2)

其中，k為推力下降系數(shù)；Ti為故障后的推力；T為標(biāo)準(zhǔn)推力。

2 凸化算法

2.1 約束條件

2.1.1 初始狀態(tài)約束

假設(shè)發(fā)動(dòng)機(jī)出現(xiàn)故障時(shí)刻為t0，后續(xù)軌跡優(yōu)化初始狀態(tài)即為t0時(shí)刻火箭運(yùn)動(dòng)狀態(tài)，X為火箭狀態(tài)變量，下標(biāo)0表示優(yōu)化起始時(shí)刻，起始點(diǎn)等式約束即為

X0=[r0v0m0]

(3)

2.1.2 終端約束

由于研究對(duì)象為三級(jí)運(yùn)載火箭，二級(jí)飛行段的主要任務(wù)為提升高度和速度，對(duì)終端約束并沒有嚴(yán)格的要求，為了保證后續(xù)飛行段的穩(wěn)定，增加兩個(gè)終端約束。

(1)軌道面約束

為防止后續(xù)飛行出現(xiàn)偏離目標(biāo)軌道的情況，需要將運(yùn)載火箭約束在目標(biāo)軌道平面內(nèi)飛行。由于在發(fā)射坐標(biāo)系下軌道平面相關(guān)參數(shù)計(jì)算復(fù)雜，因此將終端軌道面約束轉(zhuǎn)到近焦點(diǎn)坐標(biāo)系下，可以有效提高計(jì)算效率，軌道面約束表現(xiàn)形式如下

rpfz=0
vpfz=0

(4)

式中，vpfz、rpfz分別為近焦點(diǎn)坐標(biāo)系下終端速度和高度。

近焦點(diǎn)坐標(biāo)系和發(fā)射慣性坐標(biāo)系間的轉(zhuǎn)換矩陣為

(5)

式中，Ω為升交點(diǎn)經(jīng)度；ω為近地點(diǎn)幅角；i為軌道傾角；λ0為發(fā)射點(diǎn)經(jīng)度；B0為地理緯度；A0為射擊方位角。

(2)終端高度約束

由于發(fā)動(dòng)機(jī)秒耗量下降，火箭在燃料冗余的情況下可以適當(dāng)延長(zhǎng)二級(jí)飛行時(shí)間，火箭在三級(jí)一次飛行段過后會(huì)進(jìn)入一個(gè)過渡軌道，飛行時(shí)間的延長(zhǎng)可能導(dǎo)致飛行高度超過或者接近這一過渡軌道，不利于后續(xù)三級(jí)一次飛行段的軌跡規(guī)劃和制導(dǎo)，因此在二級(jí)飛行段軌跡規(guī)劃時(shí)增加飛行高度約束。以標(biāo)準(zhǔn)軌跡二級(jí)飛行結(jié)束時(shí)刻地心距作為約束邊界條件。具體表現(xiàn)形式如下

(6)

其中，rf為故障條件下二級(jí)飛行結(jié)束時(shí)火箭的地心距；Rmax為標(biāo)準(zhǔn)軌跡二級(jí)飛行結(jié)束時(shí)火箭的地心距。

2.1.3 過程約束

當(dāng)出現(xiàn)推力下降情況時(shí)，由于火箭推力下降，同一時(shí)刻的火箭速度、加速度和標(biāo)準(zhǔn)彈道相比也出現(xiàn)了下降的情況，且二級(jí)飛行段已飛出了稠密大氣層，大氣密度可以忽略不計(jì)。因此，與這三項(xiàng)相關(guān)的動(dòng)壓、熱流、過載、彎矩等約束在同一時(shí)刻均能滿足約束要求，在優(yōu)化過程中不予考慮，主要考慮飛行過程中的推力約束與質(zhì)量約束。

(1)推力約束

該液體火箭發(fā)動(dòng)機(jī)不具備推力調(diào)節(jié)能力，其推力固定，因此視為等式約束。

(7)

(2)質(zhì)量約束

為防止二級(jí)飛行段結(jié)束時(shí)出現(xiàn)計(jì)算錯(cuò)誤，保證仿真符合物理規(guī)律，運(yùn)載火箭在二級(jí)飛行時(shí)質(zhì)量必須大于其結(jié)構(gòu)質(zhì)量，具體形式如下

mdry≤m(t)

(8)

其中，mdry為火箭結(jié)構(gòu)質(zhì)量。

2.1.4 最優(yōu)化性能指標(biāo)

本文所設(shè)計(jì)的最優(yōu)化性能指標(biāo)為終端能量最大，即飛行過程中的能量損耗最小，具體指標(biāo)形式如下

(9)

其中，μ為地球引力常數(shù)；vf為終端時(shí)刻的速度。

由于終端時(shí)刻在近焦點(diǎn)坐標(biāo)系下z方向的速度、位置分量均為0，為了簡(jiǎn)便計(jì)算，提高計(jì)算效率，將優(yōu)化命題在近焦點(diǎn)坐標(biāo)系下表示，即為

(10)

其中，rpxf、rpyf、vpxf、vpyf分別為終端時(shí)刻運(yùn)載火箭在近焦點(diǎn)坐標(biāo)系下位置和速度在x軸和y軸的分量。

2.2 最優(yōu)控制問題的建立與求解

根據(jù)上述的性能指標(biāo)及初始狀態(tài)、終端和過程約束，以及飛行過程中的動(dòng)力學(xué)約束，同時(shí)定義狀態(tài)變量X=[rvm]和控制變量U=[TxTyTzTreal]，其中Tx、Ty、Tz分別為推力在近焦點(diǎn)坐標(biāo)系下3個(gè)方向的分量。

建立如下最優(yōu)控制問題Problem 1：

(11)

式中的非凸項(xiàng)主要來源有：1)性能指標(biāo)為非凸函數(shù)；2)動(dòng)力學(xué)方程中的引力相關(guān)項(xiàng)為非線性。

2.2.1 非凸項(xiàng)凸化處理

由于性能指標(biāo)函數(shù)為非凸，因此通過線性化手段進(jìn)行凸化處理，將性能指標(biāo)進(jìn)行一階泰勒展開并做歸一化處理，省略掉無意義的常數(shù)項(xiàng)和高階項(xiàng)后，性能指標(biāo)轉(zhuǎn)化為如下形式

(12)

其中，上標(biāo)0表示標(biāo)準(zhǔn)軌跡下各個(gè)相關(guān)參數(shù)。

除此之外將控制量等式約束松弛為二階錐約束。即為

(13)

對(duì)于引力相關(guān)項(xiàng)，由于引力所引起的加速度與推力加速度相比為小量，同時(shí)軌跡規(guī)劃的高度變化一般在幾十千米以內(nèi)，遠(yuǎn)小于地心距，所以對(duì)軌跡規(guī)劃的影響程度并不大，因此可以基于初始軌跡進(jìn)行引力加速度的快速計(jì)算。根據(jù)初始軌跡將引力相關(guān)加速度項(xiàng)轉(zhuǎn)化為與地心距相關(guān)的線性函數(shù)，在軌跡規(guī)劃過程中根據(jù)各離散點(diǎn)運(yùn)載火箭地心距插值得到引力加速度。

gk=f(r)

(14)

其中，gk表示在離散點(diǎn)處的引力加速度；f(·)為引力加速度和地心距之間的關(guān)系。

由于推進(jìn)劑秒耗量為常值，因此各離散點(diǎn)的質(zhì)量mk可以直接計(jì)算

(15)

針對(duì)動(dòng)力學(xué)方程，采用梯形離散法進(jìn)行動(dòng)力學(xué)方程的凸化，將飛行時(shí)間[t0tf]平均分為N段，離散間隔為Δt=(tf-t0)/N。具體表現(xiàn)形式如下

(16)

2.2.2 最優(yōu)控制模型

結(jié)合上述分析，建立如下二級(jí)飛行段軌跡規(guī)劃問題Problem 2：

(17)

2.2.3 求解流程

通過求解上述最優(yōu)控制問題，可以得到推力故障下的最優(yōu)軌跡，具體求解流程如圖2所示。

1)給定離散區(qū)間N以及終端精度要求εr、εv；

2)根據(jù)運(yùn)載火箭故障時(shí)間及推力故障系數(shù)計(jì)算故障后飛行時(shí)間；

(18)

其中，twork為二級(jí)工作時(shí)長(zhǎng)；t0為發(fā)生故障時(shí)刻；tf為二級(jí)飛行段結(jié)束時(shí)刻；k為推力下降系數(shù)。

3)根據(jù)所給參數(shù)積分求解初始軌跡；

4)根據(jù)初始軌跡進(jìn)行引力項(xiàng)的凸化；

5)應(yīng)用原對(duì)偶內(nèi)點(diǎn)法，求解問題Problem 2，得到最優(yōu)解X*=[r*v*m*]，U*；

3 仿真驗(yàn)證與結(jié)果分析

3.1 仿真條件設(shè)置

為了充分驗(yàn)證本文推力故障下運(yùn)載火箭非入軌段軌跡規(guī)劃算法的正確性和收斂性，重點(diǎn)針對(duì)運(yùn)載火箭二級(jí)飛行段發(fā)生故障的情況進(jìn)行仿真驗(yàn)證。假設(shè)發(fā)生故障時(shí)刻的火箭位置速度參數(shù)均可根據(jù)導(dǎo)航系統(tǒng)或地面基站獲得。

圖2 求解流程Fig.2 Solving process

本文使用的仿真實(shí)驗(yàn)軟硬件條件為Inter Corel i5-2300 CPU，3.10GHz，操作系統(tǒng)為Windows 7，求解軌跡規(guī)劃過程中的二階錐規(guī)劃問題。

仿真過程中的無故障下運(yùn)載火箭參數(shù)設(shè)置如表1所示。

表1 初始任務(wù)參數(shù)

由于研究目標(biāo)為三級(jí)運(yùn)載火箭，在三級(jí)一次飛行段過后運(yùn)載火箭會(huì)進(jìn)入一個(gè)軌道滑行，三級(jí)一次飛行段進(jìn)入的目標(biāo)軌道參數(shù)如表2所示。

表2 目標(biāo)軌道參數(shù)

發(fā)射場(chǎng)選為海南文昌衛(wèi)星發(fā)射中心，其相關(guān)參數(shù)如表3所示。

表3 發(fā)射場(chǎng)參數(shù)

運(yùn)載火箭在整個(gè)三級(jí)飛行過程中引力和地心距的關(guān)系在發(fā)射坐標(biāo)系下如圖3所示，在整個(gè)優(yōu)化過程中可以根據(jù)圖3的關(guān)系插值得到引力加速度，從而實(shí)現(xiàn)引力項(xiàng)的凸化。

圖3 引力變化曲線Fig.3 Gravitational force change curve

3.2 軌跡重規(guī)劃仿真驗(yàn)證

針對(duì)二級(jí)點(diǎn)火時(shí)刻立即發(fā)生不同程度故障進(jìn)行模擬，驗(yàn)證基于能量最優(yōu)的軌跡在線重規(guī)劃算法。由于二次飛行段為運(yùn)載火箭的非入軌段，其主要任務(wù)為爬高和加速，因此并無嚴(yán)格的終端精度要求，只需速度和高度達(dá)到要求即可。能量最優(yōu)性能指標(biāo)可以使運(yùn)載火箭充分利用所攜帶的燃料，以最小的能量損耗滿足運(yùn)載火箭終端速度和高度的約束。取離散區(qū)間N=100，εr=10-5，εv=10-5。

針對(duì)二級(jí)點(diǎn)火時(shí)發(fā)動(dòng)機(jī)發(fā)生故障，對(duì)推力、秒耗量同時(shí)下降10%的情況進(jìn)行分析，二級(jí)飛行段在線軌跡重規(guī)劃仿真實(shí)驗(yàn)結(jié)果如圖4所示。

(a)高度變化曲線

具體終端誤差如表4所示。

表4 終端誤差對(duì)比

為驗(yàn)證推力經(jīng)松弛后其幅值不隨時(shí)間改變，推力幅值隨時(shí)間的變化曲線如圖5所示。

圖5 推力幅值變化曲線Fig.5 Thrust amplitude variation curve

仿真結(jié)果表明，應(yīng)用故障條件下基于能量最優(yōu)的軌跡在線重規(guī)劃算法得到的軌跡終端速度和高度誤差均控制在0.2%以內(nèi)，具有良好的準(zhǔn)確度；且僅需迭代5次即可達(dá)到收斂條件，單次迭代時(shí)間在0.3s左右，具有良好的魯棒性和穩(wěn)定性。推力幅值經(jīng)松弛后不隨時(shí)間改變，驗(yàn)證了其有效性。

4 結(jié) 論

本文研究了基于能量最優(yōu)的運(yùn)載火箭非入軌段在發(fā)生推力故障后的軌跡重規(guī)劃問題。通過無損凸化、線性化等方法，將最優(yōu)控制問題轉(zhuǎn)化為二階錐規(guī)劃問題后，應(yīng)用原對(duì)偶內(nèi)點(diǎn)法精確求解。仿真結(jié)果表明，該方法能夠滿足終端速度和高度的要求，能夠完成運(yùn)載火箭在該類故障下軌跡重規(guī)劃，計(jì)算效率高，能夠?qū)崿F(xiàn)軌跡的快速收斂，提高了運(yùn)載火箭適應(yīng)故障的能力，具有在線應(yīng)用的潛力。