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

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

0 引 言

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

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

結(jié)合上述國內(nèi)外研究現(xiàn)狀分析，本文研究的大氣層內(nèi)三維能量管理存在如下難點(diǎn)：

1)由于氣動力的存在，難以準(zhǔn)確計(jì)算剩余速度能力。

2)由于發(fā)動機(jī)工作環(huán)境的差異，導(dǎo)致發(fā)動機(jī)的比沖和秒流量與地面試驗(yàn)存在差異。理論計(jì)算的局限性以及風(fēng)洞實(shí)驗(yàn)與實(shí)際飛行條件不同，使得計(jì)算時(shí)用到的氣動系數(shù)與真實(shí)值存在偏差。這些參數(shù)不確定性將影響制導(dǎo)精度。

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

4)當(dāng)發(fā)動機(jī)關(guān)機(jī)時(shí)，SLV需要以零側(cè)向速度回到發(fā)射平面。三維制導(dǎo)帶來了終端側(cè)向位移和側(cè)向速度約束，這使得終端約束的數(shù)量增多，從而增加了制導(dǎo)的復(fù)雜性。

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

1 問題描述

1.1 上升段運(yùn)動模型

在發(fā)射坐標(biāo)下建立SLV上升段運(yùn)動模型，將三維空間運(yùn)動分解為縱向平面運(yùn)動和側(cè)向平面運(yùn)動，如圖1所示。v為SLV的速度；vxy和vxz分別為縱向平面和側(cè)向平面上的速度投影；Pe為發(fā)動機(jī)推力，其方向沿SLV體軸；Pey和Pez分別為推力在縱向平面和側(cè)向平面上的投影；Θ為當(dāng)?shù)貜椀纼A角；Φ為SLV彈體軸線與當(dāng)?shù)厮矫娴膴A角，稱為當(dāng)?shù)馗┭鼋?；α為攻角；σ為?dāng)?shù)貜椀榔?；Ψ為體軸與當(dāng)?shù)乜v向平面的夾角，稱為當(dāng)?shù)仄浇?；β為?cè)滑角；D,L和R分別為阻力、升力和側(cè)力；G為重力。

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

為了簡化計(jì)算，對模型進(jìn)行了若干假設(shè)。第一，將SLV視為質(zhì)點(diǎn)，忽略旋轉(zhuǎn)。第二，不考慮擺動噴管的影響。第三，地球被視為一個(gè)球體。第四，忽略科氏力和離心力。第五，重力加速度為常值?；谏鲜黾僭O(shè)，大氣層內(nèi)SLV上升段三維運(yùn)動方程為

(1)

式中：h為高度；z為側(cè)向位移；Pe和m分別表示發(fā)動機(jī)推力大小和SLV質(zhì)量；g為重力加速度。

氣動力的計(jì)算方式如下：

(2)

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

1.2 約束條件

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

(3)

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

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

(4)

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

攻角側(cè)滑角的變化率約束分別為

(5)

動壓約束為

q=0.5ρv2≤qmax

(6)

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

終端約束：包含終端高度、速度、終端當(dāng)?shù)貜椀纼A角、側(cè)向位移約束以及當(dāng)?shù)貜椀榔羌s束

(7)

不確定性：SLV飛行過程中面臨的主要不確定性來自于發(fā)動機(jī)參數(shù)，即發(fā)動機(jī)比沖和質(zhì)量秒流量

(8)

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

(9)

2 制導(dǎo)算法設(shè)計(jì)

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

在線參數(shù)辨識：根據(jù)加速度計(jì)的測量值，將ΔcD, ΔcL和Pe/m視作擴(kuò)展?fàn)顟B(tài)量。通過對模型的合理簡化，建立上升段參數(shù)辨識的近似模型。利用容積卡爾曼濾波對模型中的不確定性參數(shù)進(jìn)行辨識，并實(shí)時(shí)修正模型參數(shù)。詳細(xì)過程見3.1。

在線軌跡規(guī)劃：實(shí)現(xiàn)軌跡的在線更新。構(gòu)建攻角和剩余速度能力曲線的數(shù)學(xué)形式，將終端高度、速度、當(dāng)彈道傾角、側(cè)向位移和當(dāng)?shù)貜椀榔潜硎緸殛P(guān)于曲線參數(shù)的代數(shù)方程。提出一種改進(jìn)的側(cè)向能量管理，在傳統(tǒng)真空能量管理的基礎(chǔ)上，考慮氣動力的影響，分別給出發(fā)動機(jī)和氣動力產(chǎn)生的剩余速度能力的計(jì)算方法。詳細(xì)過程見3.3。

軌跡解算：采用牛頓法實(shí)時(shí)求解攻角參數(shù)和速度能力曲線，實(shí)現(xiàn)在線軌跡規(guī)劃。詳細(xì)過程見3.3。

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

圖2 制導(dǎo)方案示意圖Fig.2 Diagram of the guidance scheme

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

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

(10)

由于SLV在第二級時(shí)質(zhì)量較大，有如下近似關(guān)系：

(11)

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

(12)

式中：δD,δL>0，為常數(shù)。

(13)

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

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

(14)

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

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

zk=H(xk,uk)+wk

(15)

利用文獻(xiàn)[24]中的容積卡爾曼濾波，對由式(14)和(15)構(gòu)成的離散系統(tǒng)進(jìn)行辨識，可實(shí)現(xiàn)對ΔPe/m, ΔcD和ΔcL的在線估計(jì)。

2.2 大氣層內(nèi)側(cè)向能量管理的制導(dǎo)

2.2.1大氣層內(nèi)側(cè)向能量管理的基本原理

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

(16)

式中:VTCAPz為側(cè)向平面內(nèi),由發(fā)動機(jī)推力和側(cè)向力共同產(chǎn)生的負(fù)荷剩余速度能力曲線。

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

注1.根據(jù)SLV在上升過程中的飛行特點(diǎn)，v在推力作用下增大。SLV在發(fā)射的過程中，姿態(tài)由垂直到水平，當(dāng)?shù)貜椀纼A角從90°降至0。因此，vx=vcosΘ是單調(diào)遞增的，可以用來作為飛行過程的自變量。

2.2.2氣動力的分解

傳統(tǒng)的縱向平面大氣層外能量管理分別計(jì)算推力和重力產(chǎn)生的剩余速度能力[16]?；谠撍枷耄謩e將D和R沿著側(cè)向平面內(nèi)發(fā)動機(jī)推力方向和當(dāng)?shù)貍?cè)向方向分解，如圖4所示。DP,DL,RP和RL分別是D和R沿著這兩個(gè)方向的分量。Dxoz表示阻力在當(dāng)?shù)厮矫娴姆至俊?/p>

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

通過上述分解，可以得到R和D在側(cè)向平面內(nèi)，沿推力方向和當(dāng)?shù)貍?cè)向方向上的分量為

(17)

2.2.3剩余速度能力分量的計(jì)算

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

(18)

式中：M(vx)為質(zhì)量關(guān)于vx的函數(shù)。

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

(19)

由式(18)可知，在發(fā)動機(jī)剩余速度能力計(jì)算時(shí)，將Pe/M(vx)作為一個(gè)整體?？赏ㄟ^2.1節(jié)對Pe/m進(jìn)行更新，獲得較為精確的發(fā)動機(jī)剩余速度能力。

參考VTCAPE的計(jì)算方法，有：

(20)

式中：

(21)

2.2.4VTCAPz的起始點(diǎn)與終點(diǎn)

設(shè)VTCAPz的表達(dá)式為S(vx)，則S(vx)當(dāng)前時(shí)刻下的起始點(diǎn)為(vtcosΘtcosσt,vtcosΘtsinσt)

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

(22)

(23)

為了保證終端當(dāng)?shù)仄浇铅?0,S(vx)在終點(diǎn)(vf,0)的導(dǎo)數(shù)為：

S′(vf)=0

(24)

當(dāng)速度到達(dá)vf時(shí)，為了使推力和氣動力在側(cè)向平面的分量產(chǎn)生的剩余速度能力完全耗盡，有：

(25)

2.2.5側(cè)向位移增量Δz計(jì)算

設(shè)A(vxA,S(vxA))和B(vxB,S(vxB))分別為VTCAPz上的兩點(diǎn)。當(dāng)速度沿著VTCAPz從A至B時(shí)，側(cè)向位移ΔzA→B可以按照如下計(jì)算：

(26)

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

(27)

在圖1中，有如下近似關(guān)系：

β=Ψ-σ

(28)

(29)

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

(30)

因此，當(dāng)速度向量沿著VTCAPz由A點(diǎn)運(yùn)動到B點(diǎn)時(shí)，產(chǎn)生的側(cè)向位移ΔzA→B為

ΔzA→B=

(31)

則從當(dāng)前時(shí)刻到終端時(shí)刻，由VTCAPz和VACAPP共同引起的側(cè)向位移Δz為

(32)

2.2.6高度增量和當(dāng)?shù)貜椀纼A角增量的計(jì)算

(33)

(34)

因此，Θ和h均可以表示成關(guān)于vx的方程：

(35)

終端高度和當(dāng)?shù)貜椀纼A角約束可表示為

(36)

由于終端高度和當(dāng)?shù)貜椀纼A角通過攻角進(jìn)行控制，為了保證攻角解的唯一性，構(gòu)造了一條包含兩個(gè)待定參數(shù)的攻角曲線。假設(shè)SLV將在如下攻角的數(shù)學(xué)形式下進(jìn)行剩余飛行：

α(vx)=γ1vx+γ2

(37)

式中：γ1和γ2是制導(dǎo)周期內(nèi)需要求解的參數(shù)。

2.3 軌跡的解算

2.3.1約束代數(shù)方程的建立

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

(38)

由于SLV在終端時(shí)刻需要返回發(fā)射面，因此終端側(cè)向位移為零，即Δzf=0。

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

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

D(vx-vxf)+E

(39)

通過將攻角和剩余速度能力曲線設(shè)計(jì)為特定的數(shù)學(xué)形式，將上升段能量管理制導(dǎo)問題轉(zhuǎn)化為代數(shù)方程求解問題，進(jìn)而在線獲得制導(dǎo)指令。

2.3.2S(vx)和α(vx)參數(shù)計(jì)算

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

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

(40)

顯然有：

(41)

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

(42)

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

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

(43)

2.4 過程約束

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

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

2.4.1攻角變化率約束

根據(jù)導(dǎo)數(shù)分解原則，有

(44)

則

(45)

式中：

(46)

2.4.2縱向平面法向過載約束

縱向平面內(nèi)的法向過載約束如下：

(47)

其邊界情況分為：(a)α>0和(b)α<0?？紤]到篇幅的限制，這里僅以α>0的情況為例進(jìn)行說明。

根據(jù)式(47)，當(dāng)觸發(fā)法向過載約束時(shí)滿足：

cosΨmax1cosΦmax1sinΘt+D

(48)

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

S2(vxt))cLSm

(49)

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

υ(k-1)

(50)

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

根據(jù)式(44)，有：

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

γ2

(51)

因此，有：

|γ2|

(52)

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

2.4.3動壓約束

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

ρ(h)=ρ0e-ζh

(53)

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

根據(jù)式(53)，需要實(shí)時(shí)滿足：

(54)

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

h(vxt+τ)

(55)

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

(56)

式中：

(57)

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

3 仿真校驗(yàn)

本節(jié)通過仿真校驗(yàn)來驗(yàn)證R3DEM算法的有效性。首先，給出仿真的初始條件。然后，驗(yàn)證該算法在不同任務(wù)中的適用性，以及在含有模型參數(shù)不確定性情況下的魯棒性。最后，通過與模型預(yù)測靜態(tài)規(guī)劃算法和改進(jìn)粒子群優(yōu)化算法的比較，驗(yàn)證算法的優(yōu)越性。

3.1 仿真條件

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

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

3.2 適應(yīng)性驗(yàn)證

為了驗(yàn)證R3DEM算法的適應(yīng)性，對終端高度63000～72000 m，終端速度2550～3000 m/s下的四種制導(dǎo)任務(wù)進(jìn)行測試。終端當(dāng)?shù)貜椀纼A角設(shè)置為6.5°。

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

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

3.3 魯棒性驗(yàn)證

圖6(a)和(b)分別為蒙特卡洛仿真試驗(yàn)下，h-v和z-v的包絡(luò)線。所有軌跡均滿足終端高度、速度和當(dāng)?shù)貜椀纼A角約束，且攝動軌跡均在標(biāo)稱情況附近。終端高度、速度、當(dāng)?shù)貜椀纼A角最大誤差值分別為13.3 m, 5.44 m/s, 0.013°，滿足精度指標(biāo)要求。

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

3.4 與現(xiàn)有方法比較

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

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

4 結(jié) 論

本文提出了一種基于能量管理的三維在線軌跡規(guī)劃制導(dǎo)算法，仿真結(jié)果表明，該算法適用于耗盡關(guān)機(jī)的SLV大氣層內(nèi)上升段。本文的主要結(jié)論如下：

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

2) 本文提出的制導(dǎo)算法將傳統(tǒng)的縱向平面制導(dǎo)擴(kuò)展到三維空間，擴(kuò)大了終端速度的可調(diào)范圍。當(dāng)終端在縱向平面上達(dá)到規(guī)定的高度和速度時(shí)，其側(cè)向位移和速度收斂至零。與模型預(yù)測靜態(tài)規(guī)劃方法相比，拓寬了終端速度調(diào)整范圍。與改進(jìn)的粒子群算法相比，計(jì)算效率得到了提高。

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