大氣層內固體火箭多約束魯棒三維能量管理制導

劉 飛，王松艷，楊 明，晁 濤

(哈爾濱工業大學航天學院控制與仿真中心，哈爾濱 150001)

0 引 言

近些年來，固體運載火箭(Solid launch vehicle，SLV)上升段制導受到了國內外學者的廣泛關注[1]。基于優化理論，文獻[2]利用間接法對固體火箭的上升段軌跡進行了優化，將軌跡優化問題轉化為Hamilton兩點邊值問題，并用數值方法求解。然而，間接法存在兩個明顯的缺點：由最優條件組成的約束方程非常復雜；數值方法對初值猜想的敏感度較高。相比于間接法，直接法(如高斯偽譜法[3])對于具有快時變和非線性特性的SLV上升段具有更大的工程實踐意義。自適應制導算法是近年來的研究熱點[4]，但在動態參數變化較大的極端情況下，制導算法的可靠性無法得到保證。為了提高軌跡生成效率，許多學者對無損凸優化技術進行了研究，已成功應用到運載火箭[5]、航天器[6]、高超聲速滑翔飛行器[7]的軌跡優化中。然而，大多數過程約束和非線性動態都不能無損凸化。近些年來，計算制導引起了學者們的廣泛關注。預測校正方法是一種典型的計算制導方法[8]，但其對算法的實時性要求很高。

SLV的一個顯著特點是發動機推力和關機時間不可控。對此，文獻[9]采用了耗盡關機技術，通過改變飛行姿態角來調節所需剩余速度。在過去的幾十年里，學者們致力于SLV能量管理算法的研究[10]。Yao等[11]提出了一種單向姿態能量管理方法。文獻[12]研究了基于零射程線的耗盡關機能量管理算法，其優點是降低了算法對關機點的靈敏度，但并沒有考慮過程約束。Patha等[13]研究了基于交替俯仰角變化的姿態能量管理(Attitude energy management，AEM)方法，但姿態角變化曲線需要離線計算。Wang等[14]提出了一種改進的姿態能量管理，減小由于忽略角速度而產生的跟蹤誤差。上述方法均未考慮不確定性，在有參數攝動的情況下魯棒性較差。Zarchan[15]提出的通用能量管理(General energy management，GEM)方法是一種閉環制導算法，其速度控制精度較高，但起始點和結束點的攻角較大，可能會違反過程約束。文獻[16]采用樣條能量管理方法構建約束方程，利用數值方法在線求解。文獻[17]提出了一種基于動態逆的能量管理方法。這兩種方法都能實現對終端速度的精確控制，但不能調整終端高度。文獻[18]提出了基于推力矢量控制的能量管理的完整方案。上述算法在縱向平面進行能量管理，避免了與側向運動的耦合，但減弱了速度調節能力。文獻[19]引入側向能量管理消耗剩余能量，從而精確控制終端速度。然而由于縱、側向機動的耦合，側向制導指令可能會對終端高度產生影響。文獻[20]針對滿足多約束條件的標稱軌跡，提出了一種離線軌跡優化方法。雖然該算法的精度較高，但不具備在線自適應能力。現有的能量管理算法僅在大氣層外進行[21]，這導致終端速度調節能力較差。因此，需要開發一種大氣內能量管理算法以延長能量管理時間，提高速度調節能力。

結合上述國內外研究現狀分析，本文研究的大氣層內三維能量管理存在如下難點：

1)由于氣動力的存在，難以準確計算剩余速度能力。

2)由于發動機工作環境的差異，導致發動機的比沖和秒流量與地面試驗存在差異。理論計算的局限性以及風洞實驗與實際飛行條件不同，使得計算時用到的氣動系數與真實值存在偏差。這些參數不確定性將影響制導精度。

3)與大氣層外能量管理相比，大氣層內能量管理需要考慮動壓約束。對于在真空中飛行的SLV，過載約束只與當前速度有關。而動壓約束不僅與速度有關，還與高度有關。動壓約束下的控制量可行域的計算比由過載約束構成的可行域[22]更為復雜。

4)當發動機關機時，SLV需要以零側向速度回到發射平面。三維制導帶來了終端側向位移和側向速度約束，這使得終端約束的數量增多，從而增加了制導的復雜性。

基于上述問題，本文提出一種改進的上升段制導算法——基于能量管理的魯棒三維制導算法(Robust three-dimensional energy management, R3DEM)，用以解決考慮模型參數不確定的上升段在線軌跡規劃制導問題。引入側向能量管理對發動機剩余速度能力進行耗散，提高速度調節能力。將閉環制導轉化為求解關于攻角和側向速度能力曲線參數的非線性方程組問題。此外，還考慮了控制量變化率約束、過載約束和動壓約束。針對氣動系數、發動機比沖和秒流量存在的攝動，利用容積卡爾曼濾波進行模型參數辨識提高了制導算法的魯棒性。

1 問題描述

1.1 上升段運動模型

在發射坐標下建立SLV上升段運動模型，將三維空間運動分解為縱向平面運動和側向平面運動，如圖1所示。v為SLV的速度；vxy和vxz分別為縱向平面和側向平面上的速度投影；Pe為發動機推力，其方向沿SLV體軸；Pey和Pez分別為推力在縱向平面和側向平面上的投影；Θ為當地彈道傾角；Φ為SLV彈體軸線與當地水平面的夾角，稱為當地俯仰角；α為攻角；σ為當地彈道偏角；Ψ為體軸與當地縱向平面的夾角，稱為當地偏航角；β為側滑角；D,L和R分別為阻力、升力和側力；G為重力。

圖1 上升段受力分析示意圖Fig.1 Diagram of force analysis in ascending phase

為了簡化計算，對模型進行了若干假設。第一，將SLV視為質點，忽略旋轉。第二，不考慮擺動噴管的影響。第三，地球被視為一個球體。第四，忽略科氏力和離心力。第五，重力加速度為常值。基于上述假設，大氣層內SLV上升段三維運動方程為

(1)

式中：h為高度；z為側向位移；Pe和m分別表示發動機推力大小和SLV質量；g為重力加速度。

氣動力的計算方式如下：

(2)

式中：cD,cL和cR分別為升力系數、阻力系數和側向力系數；ρ為大氣密度；Sm為氣動有效面積。

1.2 約束條件

初始條件：根據上升段三維動力學模型(1)，初始條件可表示為

(3)

過程約束：對于大氣層內的上升段運動，SLV受過載、攻角和側滑角變化率約束，以及動壓約束。

其中，縱向平面和側向平面的法向過載約束分別為

(4)

式中：nv和nh分別為縱向平面和側向平面的法向過載；nvmax和nhmax分別為nv和nh允許最大值絕對值。

攻角側滑角的變化率約束分別為

(5)

動壓約束為

q=0.5ρv2≤qmax

(6)

式中：qmax為允許的最大動壓值。

終端約束：包含終端高度、速度、終端當地彈道傾角、側向位移約束以及當地彈道偏角約束

(7)

不確定性：SLV飛行過程中面臨的主要不確定性來自于發動機參數，即發動機比沖和質量秒流量

(8)

對于大氣層內的飛行，除了式(8)中的不確定性外，還考慮氣動系數不確定性。為了簡化問題，本文只考慮cL和cD的不確定性：

(9)

2 制導算法設計

R3DEM算法的制導方案示意圖如圖2所示，圖中物理量的意義將在后文進行說明。R3DEM各部分功能如下：

在線參數辨識：根據加速度計的測量值，將ΔcD, ΔcL和Pe/m視作擴展狀態量。通過對模型的合理簡化，建立上升段參數辨識的近似模型。利用容積卡爾曼濾波對模型中的不確定性參數進行辨識，并實時修正模型參數。詳細過程見3.1。

在線軌跡規劃：實現軌跡的在線更新。構建攻角和剩余速度能力曲線的數學形式，將終端高度、速度、當彈道傾角、側向位移和當地彈道偏角表示為關于曲線參數的代數方程。提出一種改進的側向能量管理，在傳統真空能量管理的基礎上，考慮氣動力的影響，分別給出發動機和氣動力產生的剩余速度能力的計算方法。詳細過程見3.3。

軌跡解算：采用牛頓法實時求解攻角參數和速度能力曲線，實現在線軌跡規劃。詳細過程見3.3。

過程約束：將動壓、過載和控制量幅值及變化率約束轉化為攻角曲線和側向速度能力曲線的可行域。若在線軌跡規劃得到的攻角曲線和側向速度能力曲線在可行界內，則按照規劃的路徑飛行；否則，沿著可行域邊界進行飛行，保證軌跡實時滿足路徑約束。詳細過程見3.4。

圖2 制導方案示意圖Fig.2 Diagram of the guidance scheme

2.1 基于容積卡爾曼濾波的參數在線辨識

由于大氣層內上升段面臨發動機參數和氣動系數等多種不確定性，常用的擴張狀態觀測器[22]對所有不確定性的總體影響進行觀測，無法將每種不確定性分開觀測。因此，本文采用濾波方法對各參數不確定性進行辨識。參數辨識中常用的非線性濾波方法，如擴展卡爾曼濾波(Extended Kalman filter，EKF)，是將非線性方程線性化處理。EKF算法需要計算擴展狀態方程的雅可比矩陣，不但增加了計算量，且容易引起濾波發散[23]。相比之下，容積卡爾曼濾波(Cubature Kalman filter，CKF)通過三階容積積分規則，利用容積點的加權對高斯積分進行逼近，可達到三階Taylor展開的精度，具有較高的濾波穩定性[24]。

(10)

由于SLV在第二級時質量較大，有如下近似關系：

(11)

在實際飛行中，上升段動力學參數的攝動具有隨機性，難以用精確的模型描述。為了簡化問題，本節將參數攝動模型近似如下：

(12)

式中：δD,δL>0，為常數。

(13)

將待辨識的模型參數與原系統狀態量合并，擴展成新的狀態量。將式(1)和式(12)組合后，采用四階龍格-庫塔法則離散化，得到擴展后的狀態方程：

xk=F(xk-1,uk-1)

(14)

式中：xk=[v,Θ,σ,h,z,ΔcD,ΔcL,ΔPe/m]T|t=k是系統在k時刻的狀態向量；uk=[α,β]T|t=k是系統在k時刻的輸入向量。

采用四階龍格-庫塔法則將式(13)離散化，得到觀測方程：

zk=H(xk,uk)+wk

(15)

利用文獻[24]中的容積卡爾曼濾波，對由式(14)和(15)構成的離散系統進行辨識，可實現對ΔPe/m, ΔcD和ΔcL的在線估計。

2.2 大氣層內側向能量管理的制導

2.2.1大氣層內側向能量管理的基本原理

大氣層內能量管理的基本原理如圖3所示。橫軸和縱軸分別代表當地水平和當地側向速度。vz和vf分別為當前時刻側向速度矢量和終端時刻期望側向速度矢量；Δvf為期望速度增量；DP和RP分別是阻力和側向力在側向平面推力方向上的分量；對于耗盡式關機的SLV，側向平面內能量管理將發動機推力和氣動力在側向平面所提供的剩余速度能力進行耗散。VCAPz是側向平面剩余速度能力，可以分解為VTCAPE(由發動機推力提供的剩余速度能力)、VACAPP和VACAPL(由氣動力提供的剩余速度能力)。根據能量管理基本原理[22]，有如下關系式：

(16)

式中:VTCAPz為側向平面內,由發動機推力和側向力共同產生的負荷剩余速度能力曲線。

圖3 側向能量管理原理Fig.3 The principle of lateral energy management

注1.根據SLV在上升過程中的飛行特點，v在推力作用下增大。SLV在發射的過程中，姿態由垂直到水平，當地彈道傾角從90°降至0。因此，vx=vcosΘ是單調遞增的，可以用來作為飛行過程的自變量。

2.2.2氣動力的分解

傳統的縱向平面大氣層外能量管理分別計算推力和重力產生的剩余速度能力[16]。基于該思想，分別將D和R沿著側向平面內發動機推力方向和當地側向方向分解，如圖4所示。DP,DL,RP和RL分別是D和R沿著這兩個方向的分量。Dxoz表示阻力在當地水平面的分量。

圖4 側向平面氣動力分解Fig.4 Aerodynamic force decompositionin lateral plane

通過上述分解，可以得到R和D在側向平面內，沿推力方向和當地側向方向上的分量為

(17)

2.2.3剩余速度能力分量的計算

將積分變量由t替換為vx，則的表達式如下：

(18)

式中：M(vx)為質量關于vx的函數。

由于氣動力遠小于飛行第二階段的推力且固體SLV采用耗盡式關機，可以得到以下近似關系：

(19)

由式(18)可知，在發動機剩余速度能力計算時，將Pe/M(vx)作為一個整體。可通過2.1節對Pe/m進行更新，獲得較為精確的發動機剩余速度能力。

參考VTCAPE的計算方法，有：

(20)

式中：

(21)

2.2.4VTCAPz的起始點與終點

設VTCAPz的表達式為S(vx)，則S(vx)當前時刻下的起始點為(vtcosΘtcosσt,vtcosΘtsinσt)

S(vtcosΘtcosσt)=vtcosΘtsinσt

(22)

(23)

為了保證終端當地偏航角Ψ=0,S(vx)在終點(vf,0)的導數為：

S′(vf)=0

(24)

當速度到達vf時，為了使推力和氣動力在側向平面的分量產生的剩余速度能力完全耗盡，有：

(25)

2.2.5側向位移增量Δz計算

設A(vxA,S(vxA))和B(vxB,S(vxB))分別為VTCAPz上的兩點。當速度沿著VTCAPz從A至B時，側向位移ΔzA→B可以按照如下計算：

(26)

由于直接求取?ΔzA→B/?vx較為困難，可以將其寫成導數乘積的形式：

(27)

在圖1中，有如下近似關系：

β=Ψ-σ

(28)

(29)

將式(27)和(28)代入到式(1)，有

(30)

因此，當速度向量沿著VTCAPz由A點運動到B點時，產生的側向位移ΔzA→B為

ΔzA→B=

(31)

則從當前時刻到終端時刻，由VTCAPz和VACAPP共同引起的側向位移Δz為

(32)

2.2.6高度增量和當地彈道傾角增量的計算

(33)

(34)

因此，Θ和h均可以表示成關于vx的方程：

(35)

終端高度和當地彈道傾角約束可表示為

(36)

由于終端高度和當地彈道傾角通過攻角進行控制，為了保證攻角解的唯一性，構造了一條包含兩個待定參數的攻角曲線。假設SLV將在如下攻角的數學形式下進行剩余飛行：

α(vx)=γ1vx+γ2

(37)

式中：γ1和γ2是制導周期內需要求解的參數。

2.3 軌跡的解算

2.3.1約束代數方程的建立

由上述設計過程可知，VTCAPz=VTCAPE+VACAPP。根據式(36),(22)-(25)和(32)，三維軌跡應滿足：

(38)

由于SLV在終端時刻需要返回發射面，因此終端側向位移為零，即Δzf=0。

注3.為了保證式(38)有唯一解，α(vx)和S(vx)中未知數的總數應該等于式(37)中的方程個數。由于α(vx)包含了2個未知數[γ1,γ2]T,S(vx)應包含5個未知數。有無窮多數學形式的S(vx)可以滿足式(38)。對于該問題，多項式曲線有兩個優點：首先，多項式函數的特點是它可以任意增加未知數的數目來匹配方程組的數目。其次，它的導數形式比較簡單。因此，可用5階多項式曲線來擬合S(vx)：

S(vx)=A(vx-vxf)4+B(vx-vxf)3+C(vx-vxf)2+

D(vx-vxf)+E

(39)

通過將攻角和剩余速度能力曲線設計為特定的數學形式，將上升段能量管理制導問題轉化為代數方程求解問題，進而在線獲得制導指令。

2.3.2S(vx)和α(vx)參數計算

假設所有狀態量都是可測的，則式(38)是關于[γ1,γ2,A,B,C,D,E]T的方程組。在式(38)中：

S(vf)=VACAPL,S′(vf)=0

(40)

顯然有：

(41)

因此，式(38)可以化簡為

(42)

設K=[K1,K2,K3,K4,K5]T,Π=[γ1,γ2,A,B,C]T。為了便于計算，采用Simpson法對積分進行近似。

注4.一些改進的數值迭代算法同樣可以有效地解決本問題，兼顧計算效率和收斂性能，本文選擇牛頓法進行迭代計算：

(43)

2.4 過程約束

本文考慮過載、和控制量變化率、控制量幅值以及動壓約束約束。

注6.側向平面法向過載約束和側滑角的變化率約束可行域求解方法與文獻[22]的相同，這里不再贅述。

2.4.1攻角變化率約束

根據導數分解原則，有

(44)

則

(45)

式中：

(46)

2.4.2縱向平面法向過載約束

縱向平面內的法向過載約束如下：

(47)

其邊界情況分為：(a)α>0和(b)α<0。考慮到篇幅的限制，這里僅以α>0的情況為例進行說明。

根據式(47)，當觸發法向過載約束時滿足：

cosΨmax1cosΦmax1sinΘt+D

(48)

設υ=sinΦmax，則式(48)可以寫作：

S2(vxt))cLSm

(49)

對式(49)進行牛頓迭代，可以獲得υ1：

υ(k-1)

(50)

式中：υ(k-1)和υk分別表示相鄰兩次牛頓迭代的υ值。

根據式(44)，有：

|Φ|=|γ1vx+γ2+Θt|<|γ1vx+Θt|+|γ2|=

γ2

(51)

因此，有：

|γ2|

(52)

由此可以得出，縱向平面內的法向過載約束對|γ2|施加了限制。

2.4.3動壓約束

大氣層內能量管理不但要考慮過載約束、攻角和側滑角的幅值、變化率約束，還要滿足動壓約束。由于側向位移相對高度而言較小，近似認為動壓只與高度和速度有關。假設大氣密度隨高度變化如下：

ρ(h)=ρ0e-ζh

(53)

式中：ρ0是零海拔處的大氣密度；ζ>0為常數。

根據式(53)，需要實時滿足：

(54)

定義動壓約束觸發判斷方程(55)。第一項表示在某一速度vxt+τ時，觸發最大動壓約束qmax的最低臨界高度。第二項表示速度為vxt+τ時的預測高度。這里，τ=(vfcos Θf-vxt)/100是關于vx的積分步長。當預測高度高于臨界高度時，不違反動壓約束。

h(vxt+τ)

(55)

若Jq>0，則說明當vx=vxt+τ時，動壓超出約束值qmax。此時，對式(42)中的K5(vx)進行替換，使得vx=vxt+τ時滿足動壓約束：

(56)

式中：

(57)

另一方面，若當前時刻Jq<0，則K5(vx)保持不變。

3 仿真校驗

本節通過仿真校驗來驗證R3DEM算法的有效性。首先，給出仿真的初始條件。然后，驗證該算法在不同任務中的適用性，以及在含有模型參數不確定性情況下的魯棒性。最后，通過與模型預測靜態規劃算法和改進粒子群優化算法的比較，驗證算法的優越性。

3.1 仿真條件

本文的研究背景是三級SLV的第二級，表1給出了制導任務滿足的部分約束條件。

表1 初始條件、終端約束和路徑約束Table 1 Initial, terminal and path constraints

3.2 適應性驗證

為了驗證R3DEM算法的適應性，對終端高度63000～72000 m，終端速度2550～3000 m/s下的四種制導任務進行測試。終端當地彈道傾角設置為6.5°。

4種終端高度和終端速度下的仿真結果如圖5所示。圖5(a)和(c)為h-v曲線和當地彈道傾角隨時間變化曲線。可以得出，R3DEM算法可以滿足4種不同任務下的終端狀態約束。圖5(e)為攻角隨時間變化曲線。由圖5(d)中當地彈道偏角隨時間變化曲線可知，當速度達到終端約束值時，當地彈道偏角收斂至0。由圖5(b)中的z-v曲線可知，側向位移z在終端時刻收斂至0，這表明SLV最終回到了發射平面。此外，圖5(b)中側向位移變化率收斂至0，圖5(d)中終端當地彈道偏角收斂至0，二者均表明終端時刻側向速度收斂至0。終端側滑角收斂至0，表示速度方向與彈體體軸方向重合，體軸調整到發射平面，滿足了第三級助推的起始條件。圖5(f)表明，為了達到精確控制速度的目的，在側向平面內機動消耗了多余的能量，且終端高度和速度越小，需要消耗的發動機能量越多。圖5(g)和(h)分別為縱向平面法向過載絕對值和側向平面法向過載絕對值隨時間變化曲線。圖5(i)為動壓隨時間變化曲線。第二級飛行過程始終滿足所有的過程約束，且觸發了過程約束的邊界條件。盡管如此，由于速度能力曲線始終在其可行域內，因此約束始終沒有超出邊界值。

圖5 不同高度、速度任務下的仿真結果Fig.5 Simulation results for various combinations of altitude and velocity

3.3 魯棒性驗證

圖6(a)和(b)分別為蒙特卡洛仿真試驗下，h-v和z-v的包絡線。所有軌跡均滿足終端高度、速度和當地彈道傾角約束，且攝動軌跡均在標稱情況附近。終端高度、速度、當地彈道傾角最大誤差值分別為13.3 m, 5.44 m/s, 0.013°，滿足精度指標要求。

圖6 含有不確定性的蒙特卡洛仿真結果Fig.6 Monte Carlo simulation results with uncertainty

3.4 與現有方法比較

為了檢驗R3DEM的速度調節能力，分別給出了R3DEM和MPSP算法的終端可達區域。由圖7(a)可知，R3DEM的可達區域完全覆蓋了MPSP。因此，與MPSP相比，R3DEM具有更強的終端適應能力。IPSO適應度隨迭代代數的變化曲線如圖7(b)所示。IPSO經過120次迭代后收斂，得到最優解，平均收斂時間超過5 min，這顯然不適用于在線制導。相比之下，R3DEM僅需10 s左右。圖7(c)和(d)分別為高度和速度隨時間變化曲線，三種算法都可以達到期望的終端速度和高度。由圖7(e)可知，MPSP產生的攻角大于R3DEM和IPSO。R3DEM和IPSO算法在縱、側向兩個平面對速度進行耗散，因此產生的攻角和側滑角較小。由于MPSP的速度耗散完全依賴于縱向平面，導致了大范圍的攻角變化，進而產生圖7(f)中較大的法向過載。過載約束是制約速度調節能力的重要因素，MPSP算法觸發過載約束的時間較長，導致其速度調節能力較低。R3DEM和IPSO算法的過載在整個過程中都處于約束值之下。

圖7 R3DEM，MPSP，IPSO的仿真結果對比Fig.7 Simulation results comparison among R3DEM, MPSP and IPSO

4 結 論

本文提出了一種基于能量管理的三維在線軌跡規劃制導算法，仿真結果表明，該算法適用于耗盡關機的SLV大氣層內上升段。本文的主要結論如下：

1) 提出了一種SLV上升段制導算法，將求解空間軌跡曲線的問題轉化為求解攻角和速度能力曲線參數的問題，能以較高的精度滿足終端約束。

2) 本文提出的制導算法將傳統的縱向平面制導擴展到三維空間，擴大了終端速度的可調范圍。當終端在縱向平面上達到規定的高度和速度時，其側向位移和速度收斂至零。與模型預測靜態規劃方法相比，拓寬了終端速度調整范圍。與改進的粒子群算法相比，計算效率得到了提高。

3) 可在過載、動壓和控制量變化率約束，以及模型參數不確定性情況下，能以較高精度完成制導任務。