一種改進的運載火箭迭代制導方法

馬宗占，許志,*，唐碩，張遷

1. 西北工業大學 航天學院，西安 710072 2. 陜西省空天飛行器設計技術重點實驗室，西安 710072

提高可靠性并降低成本已成為當前運載火箭發展的兩大趨勢。SpaceX的“獵鷹9”火箭在主動段飛行時，當9臺發動機中的2臺出現故障時依然能夠準確將有效載荷送入預定軌道[1]，這要求其制導算法必須具備高精度、強魯棒性和對非災難故障狀況的適應能力[2-3]。

從20世紀60年代開始，基于最優控制理論的迭代制導算法便已經應用于運載火箭大氣層外制導，該算法由于不依賴于標準軌跡，只基于當前飛行狀態和目標軌道參數迭代計算出最優入軌指令，因此該算法對參數不確定性都具有較強的適應能力。早期土星火箭采用基于線性正切形式的迭代制導算法(IGM)[4-5]，該算法將引力場假設為常值，通過把制導指令分解成俯仰和偏航兩個平面(假設偏航指令為小量)來簡化迭代流程。20世紀 70年代初，為了適應航天飛機多任務制導要求，Brand[6]和Long[7]等在IGM算法的基礎上，開發了動力顯式制導(PEG)算法[6-7]，該算法將三維推力矢量作為控制量，以正交假設代替橫截條件，對推力積分做線性化展開，以剩余速度作為迭代變量，并引入推力可調模式等處理手段，這使得該算法相比于IGM在最優性和任務適應性上具有更大的優勢。Jaggers[8-9]、Von der Porten[10]等對PEG算法進行了改進，采用線性引力場假設，控制量限幅和修改雅克比矩陣等措施進一步增強了算法的最優性和收斂性。而陳新民和余夢倫[11]在推導最優解時用3個速度分量和2個位置分量約束代替軌道根數約束從而簡化了橫截條件，茹家欣[12]在此基礎上增加了迭代入軌點緯度幅角，進一步保障了算法的收斂性。上述迭代制導算法都是降維迭代的方式求解最優控制量從而保障算法收斂，但由于采用了簡化假設代替準確的橫截條件而難以保證算法最優性，并且在進行推力積分時假設與位置約束相關的控制量為小量進行線性化處理導致降低了積分精度，這兩點不足使傳統迭代制導難以在推力大幅下降情況下保障其最優性和故障適應性。

與迭代制導方法不同，Brown等在1967年提出的最優制導(OPGUID)方法[13-14]采用線性引力場假設并嚴格滿足所有最優性條件，以協態變量和剩余時間作為迭代變量，采用多維Newton-Raphson法進行求解。Lu和Dukeman等則對最優制導算法在計算和求解上進行了改進[15-17]，并將其應用拓展到了大氣層內和助推-滑翔-助推模式。而文獻[18-20]則通過推導了牛頓法中Jacobi矩陣的解析表達式對求解過程進行了改進，并推導了基于軌道要素形式終端約束的邊界條件。與迭代制導相比，最優制導方法盡管保留了完整的最優性條件，但對初值的敏感性和多維迭代的低收斂特性也給算法在線應用帶來了較大的風險[21]，因此嚴重制約了其工程應用。

針對火箭推力故障情況下迭代制導在最優性和故障適應性上存在的局限性，本文提出了一種改進迭代制導算法。該算法重新推導了基于火箭最優入軌問題的具有5個終端軌道根數約束的橫截條件，提高了算法最優性；采用高斯-勒讓德方法提高推力積分精度；最后該算法保留了迭代制導的降維求解模式，并結合控制量限幅算法，進一步保障算法實時性和收斂性。本文算法通過能夠迭代出新的入軌點真近點角進而保障火箭入軌，對推力故障具有更強適應性。

1 火箭入軌最優控制問題

1.1 火箭動力學模型

本文在入軌點坐標系O-XYZ下建立的火箭飛行動力學模型，入軌點坐標系定義如圖1所示：坐標原點取為地心，OY軸指向參考入軌點，在慣性空間內固定，OX軸在軌道平面內，垂直于OY軸，指向衛星運行方向，OZ軸、OX軸和OY軸組成右手坐標系。圖中O-XIYIZI為地心赤道慣性系；i*為目標軌道的軌道傾角；w*為近地點幅角；Ω*為升交點赤經；f*為參考入軌點真近點角。

火箭在大氣層外飛行動力學方程為

(1)

式中：u為單位推力矢量；T為常值推力大小；Vex為發動機比沖。采用線性引力場假設

g=-ω2r

(2)

式中：ω為平均角速度，取作當前位置矢徑r0與入軌點矢徑rf模的平均值處的圓軌道角速度，即

(3)

式中：μ為地球引力常量。

1.2 最優控制問題描述

火箭入軌最優控制問題描述為：尋找u(t)，使火箭沿著該推力方向可以在滿足一定約束條件下以最小燃料消耗量進入指定的軌道。由于推力大小一定時，性能指標最小燃料消耗可以等效為最短入軌時間，因此該最優控制問題的狀態變量微分方程為方程組(1)中前兩式，性能指標為

(4)

式中：tgo為剩余飛行時間。終端約束(通常以軌道根數形式給出)用狀態量表示為

Ci(vf,rf)=0i=1,2,…,n(n<6)

(5)

式中：下標f表示終端值(下文同)。控制量約束為

(6)

對于該最優控制問題，可以利用變分法相關理論，將其轉化為終端時間自由的兩點邊值問題。

1.3 最優控制問題求解

本文應用極大值原理[22]求解最優控制量，選擇的哈密爾頓函數為

(7)

式中：λr、λv為協態變量。根據極大值原理解得控制量為

(8)

協態變量為

(9)

(10)

1.4 5個軌道根數約束下的橫截條件推導

火箭入軌一般需要滿足5個軌道根數形式的終端約束，即半長軸a、偏心率e、軌道傾角i、升交點赤經Ω以及近地點幅角w。最優控制解不僅需要滿足5個軌道根數的終端約束，還需要滿足相應的橫截條件。對于終端時間自由的兩點邊值問題，其橫截條件的表達式為

(11)

Hf=0

(12)

式中：ηi為拉格朗日乘子，終端約束Ci的表達式見附錄A。其中約束(12)可以由式(13)代替[14]

(13)

而式(11)是一組6維方程組，無法直接用于迭代制導，因此下文將對其進行降維處理。將Ci(i=1,2,…,5)的表達式代入式(11)中的前兩式得到

λf=-η1vf-η2(rf×vf)×rf-η4e3-

η3[2vf(rfN)-rf(vfN)-(rfvf)N]

(14)

(15)

(16)

在式(15)兩端同時點乘vf，得到

(17)

用式(16)減去式(17)即可消去拉格朗日乘子，得到

(18)

式(18)即為第一個橫截條件的表達式。傳統迭代制導方法通常用簡化假設代替這一約束條件，如PEG采用正交性假設[7]

(19)

而文獻[12]用3個速度分量和2個位置分量約束代替5個軌道根數約束，即令

(20)

式(20)和下文中標量加下標X、Y、Z均分別表示對應矢量在入軌點坐標系下的X、Y、Z坐標。由式(11)和式(20)得到

(21)

從而代替式(18)。本文則保留了式給出的準確橫截條件，消除了PEG及文獻[12]等對橫截條件的假設，從理論層面上提升了算法最優性。

2 迭代過程推力積分計算

定義推力積分為

(22)

(23)

由式(8)和式(10)得到最優推力單位矢量的表達式為

(24)

考慮到約束(13)，式(24)的分母項(定義為den)可以轉化為

(25)

其中：

(26)

定義推力積分系數為

(27)

(28)

有

(29)

(30)

對于積分

(31)

做變換

(32)

于是積分可以轉化為

(33)

其中：

(34)

根據GLQF可得積分式(33)的值為

(35)

式中：si為高斯-勒讓德求積公式的節點；Ai為相應的系數。同理令

(36)

可以得到

(37)

對于二重積分

(38)

做變換

(39)

Jacobi行列式為

(40)

于是積分式(38)可以轉化為

(41)

其中

gS(α,β)=

(42)

積分式(41)的值為

(43)

同理令

gQ(α,β)=

(44)

可以得到

(45)

GLQF方法精度取決于節點個數。

圖2 推力積分對比圖Fig.2 Comparison chart of thrust integral

3 迭代過程中引力積分計算

文獻[12]基于平均引力場假設進行引力積分運算，在發動機故障導致的小推力長燃弧情況下計算精度較差。本文采用球極坐標系下泰勒多項式逼近的方式給出了一種適用于小推力長燃弧飛行模式的引力積分近似解[8]：

定義引力積分為

(46)

(47)

定義平均位置rv和rr，滿足：

(48)

(49)

在線性引力場假設下，得到rv和rr的近似表達式

(50)

(51)

在球極坐標系下估算平均位置更適用于小推力長燃弧模式。令vd為滿足終端約束的目標入軌點速度，rd為目標入軌點位置，建立如圖3所示的球極坐標系。圖3中erdy為沿rd的單位矢量，er0y為沿r0的單位矢量，erdz為沿入軌點坐標系OZ軸的單位矢量，有

erdx=erdy×erdz

(52)

圖3 球極坐標系示意圖Fig.3 Diagram of polar coordinate system

er0x=er0y×er0z

(53)

er0z=er0x×er0y

(54)

令r0的極坐標形式為

(55)

其中旋轉角為

ψ0=arcsin(er0y·er0z)

(56)

由于ψ0通常是一小量，所以可以近似得到

φ0=arcsin(er0y·er0x)

(57)

同理可以得到rd的極坐標形式

(58)

以及它們對時間的一階導數

(59)

(60)

由泰勒展開式可知任意時刻的位置矢量r(t)的極坐標形式為

(61)

其對時間一階導數為

(62)

令平均位置rv和rr的極坐標形式分別為pv和pr，其近似值取做

(63)

(64)

結合式(61)計算得到

(65)

把pv和pr轉換到入軌點坐標系下，得到

rv= [cospv2(erdycospv3+erdxcospv3)+

erdzsinpv2]pv1

(66)

rr= [cospr2(erdycospr3+erdxcospr3)+

erdzsinpr2]pr1

(67)

式中：pr1表示pr的第1個坐標，其余同理。最終由式(48)和式(49)計算出引力積分值。

4 迭代制導計算流程

圖4 制導算法流程圖Fig.4 Flow chart of guidance algorithm

4.1 計算剩余時間tgo和K

定義剩余速度vgn和剩余位置rgn

vgn≡vd-v0-vgrav

(68)

rgn≡rd-r0-rgrav-vtgo

(69)

對式(1)前兩式進行兩次積分，可以看出

(70)

即

(71)

(72)

對推力積分進行簡化，認為ω(t-K)為小量，并假設[8]

(73)

對式(27)進行泰勒展開并忽略二階小量，得到

(74)

其中：LP和JP是近似的推力積分系數，僅應用于本環節，其表達式為

(75)

將式(74)代入到式(71)中，得到

(76)

令

(77)

從而得到

JT≈0

(78)

因此有

(79)

所以剩余時間為

(80)

根據式(68)和式(80)，可以迭代求出剩余時間tgo的值，而K的值依然由式(77)決定。

4.2 計算控制變量λK和

由4.1節求出tgo和K后，根據第2節中式(35)、 式(37)、式(43)和式(45)可以計算出推力積分系數LT、JT、ST、QT，由式(71)～式(72)可以得到

(81)

(82)

其中：

A=ωsin[ω(tgo-K)]vd-

(83)

B=cos[ω(tgo-K)]vd+

(84)

在式(72)兩端同時點乘B，結合式(82)可以得到修正后的ST表達式

(85)

最終得到

(86)

在式(10)中令ω趨于0，得到控制矢量的近似表達式為

(87)

圖5 控制量限幅和入軌點校正示意圖Fig.5 Diagram of control limits and injection point correction

首先計算出當前時刻單位推力矢量

(88)

當erdx·u0≥0時，不需要做任何調整，否則令

erdx·u0=0

(89)

(90)

4.4 入軌點預測和校正

根據式(29)和式(30)計算得到推力積分vthrust和rthrust，則預測入軌速度vp和位置rp為

(91)

令rd與入軌點坐標系的Y軸的夾角為θ，其表達式為

(92)

參考入軌點的真近點角為f*，則目標入軌點的真近點角為

fd=f*+θ

(93)

可以得到rd的模值為

(94)

于是

(95)

(96)

于是

(97)

當校正后入軌點位置收斂時，認為本周期迭代制導算法收斂，當前時刻制導指令由式(88)給出。

5 仿真驗證

為了驗證本文所改進迭代制導方法的最優性、發動機故障狀態下算法適應能力以及各種干擾條件下算法的魯棒性，本文針對某型火箭大氣層外末級入軌飛行段進行迭代制導算法評估。

5.1 仿真初始和終端條件

火箭初始質量為184 334.0 kg，發動機參數軸向推力為999 036 N，秒流量為56.029 kg/s。初始點和標準入軌點軌道要素見表1。

表1 仿真初始點和標準入軌點軌道根數

5.2 干擾條件下制導算法魯棒性驗證

為了驗證本文所改進的迭代制導算法的魯棒性和入軌精度，針對火箭末級給出的各種干擾和偏差(如表2所示)，進行了1 000次蒙特卡羅打靶仿真，打靶結果見圖6、表3和表4。

從蒙特卡羅打靶統計結果中可以看出，在表3 所示的偏差條件下，采用本文方法可以使火箭滿足5個終端約束精確入軌且迭代步數不大于5，說明算法可以同時適應多種偏差，具有較強的收斂性和魯棒性。

表2 蒙特卡羅打靶偏差項設置Table 2 Setting of error terms in Monte Carlo simulation

5.3 發動機故障情況下制導算法適應性驗證

由于本文只對發動機出現非災難故障狀態時(包括燃料堵塞、燃料泄漏、燃燒室壓強異常、伺服機構卡死等狀況，但總沖依然滿足入軌要求)驗證制導算法適應能力，此時發動機故障對制導的影響可歸結為秒流量大幅下降和等效比沖下降，因此設計了表5所示的3種典型故障模式。其中，當等效比沖與預裝訂值一致時可以通過實時測量視加速度獲取秒流量變化，因此故障模式1主要驗證當故障精確已知模式下算法的最優性；而某些特殊發動機故障模式會導致秒流量和等效比沖同時出現大幅下降，這對算法魯棒性和適應性同時提出挑戰，表5中的故障模式2和故障模式3對該模式下算法性能進行驗證。

表3 蒙特卡羅仿真結果Table 3 Results of Monte Carlo simulation

表4 最大迭代步數統計Table 4 Statistics of maximum iteration steps

表5 不同故障模式下入軌時間統計Table 5 Entry time in different fault modes

圖7 推力故障情況下制導指令對比Fig.7 Comparison of guidance instructions with thrust fault

采用本文算法與文獻[12]中的迭代制導算法進行仿真對比，仿真結果見表4～表5和圖7～圖8。 仿真結果顯示，本文改進的制導算法在模式1條件下相對于文獻[12]方法關機時間提前了約13.36 s，說明本文采用的包括改進推力積分和引力積分等措施可以提高算法最優性；在模式2和模式3條件下，隨著等效比沖的下降幅度增大，文獻[12]方法會出現迭代不收斂的狀況，而本文方法依然可以保障火箭精確入軌。此外，從表4中可以看出故障模式2下本文方法相對于文獻[12]方法將關機時間提前了約25.86 s，說明在更嚴重的故障狀況下本文方法最優性也更明顯。

為進一步獲得2種方法對故障的適應邊界，本文將秒流量下降幅度設置為額定的35%，文獻[12]方法最多可以適應16%的等效比沖下降，而本文改進的制導方法可以適應25%的等效比沖下降，進一步驗證了本文所提方法對故障具有更強的魯棒性和適應性。

圖8 推力大偏差情況下軌道根數曲線對比Fig.8 Comparison of orbital elements curves with thrust fault

由于算法采用預裝訂的迭代初值，在故障發生之后的第一個制導周期內收斂時間最長，因此在故障模式2的第一個制導周期將本文方法與文獻[12] 方法就算法收斂過程進行對比，結果如圖9 所示。其中niter表示迭代步數，xmiss表示本步與上一步迭代的目標入軌位置X坐標之差的絕對值。結果說明在故障模式下本文方法收斂速度更快。

表6顯示在不同故障模式下對于任意給定的真近點角初值，算法均能在有限步數內收斂(表6中邊界模式指秒流量下降35%，比沖下降25%的算法適應邊界狀況)，說明算法對于初值不敏感并且可以保障實時性。

圖9 故障模式2下收斂過程曲線對比Fig.9 Comparison of convergence curves in fault mode 2

表6 不同故障模式下最大迭代步數統計Table 6 Maximum iteration steps in different fault modes

6 結 論

1) 本文基于文獻[12]與PEG方法，重新推導了滿足5個軌道根數約束的準確的橫截條件，采用了精度更高的推力和引力積分方法，因此進一步完善了算法流程，具有較強的最優性和適應性。

2) 本文方法相對于文獻[12]方法在相同推力故障情況下具有更強的最優性與收斂性，且對于更嚴重的推力故障魯棒性和適應性更強。

3) 本文對多變量參數最優化問題進行降維迭代求解，并結合控制量限幅算法，保障了故障條件下算法的收斂性與實時性，且能適應多種偏差，具有一定的工程實用價值。

[21] AHMAD N, HAWKINS M, VON DERPORTEN P, et al. Closed loop guidance trade study for Space Launch System Block-1B vehicle[C]∥41st AAS GNC Conference, 2018.

附錄A:

由于開普勒軌道滿足:

(A1)

式中：vf和rf分別為vf和rf的模值，上標*表示軌道根數的目標值。所以半長軸約束為

(A2)

由于軌道角動量模值h滿足:

h2=μ·p=μ·a*(1-e*2)

(A3)

式中：p為軌道半通徑，所以式(A2)和式(A4)共同約束了目標軌道偏心率大小。

C2=(rf×vf)·(rf×vf)-μ·a*(1-e*2)=0

(A4)

令N表示軌道平面與赤道平面交線(由地心指向升交點)單位矢量，則近地點幅角的表達式為

(A5)

式中:e*為偏心率矢量，其表達式為

(A6)

所以式(A2)、式(A4)和式(A7)共同約束了目標軌道的近地點幅角。

C3=Ne-cos(w*)e*=

cos(w*)e*=0

(A7)

在入軌點坐標系下，令入軌點Z向速度和位置為0即可約束目標軌道的軌道傾角和升交點赤經[12]，因此軌道傾角約束和升交點赤經約束為

C4=vfe3=0

(A8)

C5=rfe3=0

(A9)

其中：