基于速度增量序列凸化的固體運(yùn)載火箭制導(dǎo)算法

陳之潤 鄭 新 趙鵬雷 于煜斌

北京宇航系統(tǒng)工程研究所，北京 100076

0 引言

近年來，隨著微小衛(wèi)星市場的蓬勃發(fā)展，固體運(yùn)載火箭憑借高可靠、快履約、快發(fā)射、低成本等優(yōu)勢迅速占領(lǐng)市場，成為世界航天大國運(yùn)載技術(shù)發(fā)展的重點(diǎn)領(lǐng)域[1]。但不同于液體運(yùn)載火箭推力可調(diào)、隨時(shí)關(guān)機(jī)的特性，固體運(yùn)載火箭為提高可靠性、降低成本、增大干質(zhì)比，一般取消推力終止裝置而采用耗盡關(guān)機(jī)制導(dǎo)方式，一方面需進(jìn)行固體發(fā)動(dòng)機(jī)能量管理，另一方面需考慮發(fā)動(dòng)機(jī)推力、秒耗量、比沖、關(guān)機(jī)時(shí)間等偏差大的問題，這給固體運(yùn)載火箭制導(dǎo)精度帶來了新的挑戰(zhàn)[2]。

目前，針對(duì)耗盡關(guān)機(jī)制導(dǎo)算法的研究，以閉路制導(dǎo)算法最為經(jīng)典，該方法以視速度模值為變量，通過閉路制導(dǎo)段、能量管理段和常姿態(tài)導(dǎo)引段實(shí)現(xiàn)關(guān)機(jī)點(diǎn)制導(dǎo)精度。文獻(xiàn)[3-5,6]在經(jīng)典閉路制導(dǎo)算法上進(jìn)行改進(jìn)，考慮了終端速度管控、引力攝動(dòng)J2項(xiàng)、需用速度增益曲面以及閉路制導(dǎo)與零射程線結(jié)合等問題。文獻(xiàn)[7-9]針對(duì)能量管理算法分別提出了姿態(tài)調(diào)制制導(dǎo)(AEM)、通用能量管理(GEM)、以及樣條能量管理(SEM)等方式，有效降低了終端速度偏差。但是，上述算法對(duì)終端多約束問題適應(yīng)性較差，且無法處理過程約束，難以滿足主動(dòng)段多過程及終端約束制導(dǎo)需求。

另外，有學(xué)者將固體火箭耗盡關(guān)機(jī)制導(dǎo)問題轉(zhuǎn)化為最優(yōu)控制問題[10]進(jìn)行求解，可處理過程及終端約束問題，但多約束的最優(yōu)控制問題難以采用解析法求解，常采用數(shù)值方法求解，受限于求解速度，較少應(yīng)用于彈上制導(dǎo)。近年來，凸優(yōu)化算法由于其求解快速性等優(yōu)點(diǎn)，逐漸應(yīng)用于固體火箭耗盡關(guān)機(jī)最優(yōu)控制問題的求解。文獻(xiàn)[11]提出了一種基于改進(jìn)序列凸優(yōu)化的固體火箭制導(dǎo)算法，將滑行時(shí)間擴(kuò)展為控制變量，實(shí)現(xiàn)了高精度入軌需求，但未考慮高低溫下發(fā)動(dòng)機(jī)內(nèi)彈道變化；文獻(xiàn)[12-13]利用二階錐規(guī)劃求解運(yùn)載火箭上升段入軌問題，文獻(xiàn)[14]提出了在推力下降時(shí)，基于序列凸優(yōu)化的在線重規(guī)劃入軌算法，但文獻(xiàn)[12-14]均為液體火箭發(fā)動(dòng)機(jī)，推力可調(diào)節(jié)，無法直接應(yīng)用于采用耗盡關(guān)機(jī)方式的固體運(yùn)載火箭制導(dǎo)。

本文針對(duì)上述耗盡關(guān)機(jī)制導(dǎo)方法的不足，提出一種以速度模值增量為變量的固體火箭序列凸化制導(dǎo)算法。主要?jiǎng)?chuàng)新點(diǎn)如下：

1)針對(duì)固體火箭發(fā)動(dòng)機(jī)在高低溫下的內(nèi)彈道(推力、秒耗量、耗盡時(shí)間)散布大的特性，提出一種以速度模值增量為變量的序列凸化制導(dǎo)算法，解決了以飛行時(shí)間為自變量的序列凸化制導(dǎo)算法精度差且魯棒性低的問題；

2)提出一種以速度模值增量為變量的高低溫內(nèi)彈道修正模型；通過對(duì)比速度模值增量在線辨識(shí)高低溫工況，以此修正內(nèi)彈道模型，提高了耗盡關(guān)機(jī)的制導(dǎo)精度。

1 問題建模

1.1 運(yùn)動(dòng)學(xué)與動(dòng)力學(xué)模型

為保證稠密大氣層內(nèi)飛行時(shí)的載荷等約束，固體運(yùn)載火箭助推段第一級(jí)或兩級(jí)一般采用攝動(dòng)制導(dǎo)，在助推段最后一級(jí)及末修級(jí)采用閉路制導(dǎo)，實(shí)現(xiàn)最終的入軌精度要求。由于火箭實(shí)施閉路制導(dǎo)的飛行時(shí)間較短，在不影響制導(dǎo)精度的前提下，可忽略地球自轉(zhuǎn)的影響。因此，基于非旋轉(zhuǎn)圓球地球模型，建立固體運(yùn)載火箭助推段運(yùn)動(dòng)學(xué)與動(dòng)力學(xué)方程如下：

(1)

其中：r為地心距，λ和φ分別為經(jīng)緯度，V為速度，θ為當(dāng)?shù)貜椀纼A角，σ為彈道偏角；α和β為攻角和側(cè)滑角；L、D和Z分別為升力、阻力和側(cè)力；m為質(zhì)量，F(xiàn)為推力，g為引力加速度。上式(1)中動(dòng)力學(xué)模型可簡記為：

(2)

其中：X(t)=[r,λ,f,V,θ,σ]T，U=[α,β]T。

由于固體運(yùn)載火箭實(shí)施閉路制導(dǎo)時(shí)已飛出稠密大氣層，氣動(dòng)力對(duì)彈道影響較小，且高超聲速飛行段的氣動(dòng)特性相對(duì)穩(wěn)定，在不影響制導(dǎo)精度的前提下，可對(duì)氣動(dòng)力系數(shù)進(jìn)行多項(xiàng)式擬合，以提升求解速度，插值公式選取如下：

L=L(α,Ma)，D=D(α,Ma)，Z=Z(β,Ma)

(3)

1.2 最優(yōu)控制問題

固體運(yùn)載火箭耗盡關(guān)機(jī)制導(dǎo)問題主要為滿足地心距、當(dāng)?shù)貜椀纼A角和彈道偏角的約束，以及飛行過程中姿控能力約束，其可轉(zhuǎn)換為控制變量攻角和側(cè)滑角的約束。該問題可轉(zhuǎn)換為如下最優(yōu)控制問題：

(4)

其中:Umax為控制量U(U=[α,β]T)的二范數(shù)約束，θf為耗盡點(diǎn)處目標(biāo)當(dāng)?shù)貜椀纼A角，rf為耗盡點(diǎn)處目標(biāo)地心距，σf為耗盡點(diǎn)處目標(biāo)彈道偏角。

1.3 內(nèi)彈道模型

固體發(fā)動(dòng)機(jī)內(nèi)彈道模型主要包含發(fā)動(dòng)機(jī)的推力、室壓和流量等參數(shù)。發(fā)動(dòng)機(jī)的室壓Pc、推力F和流量qm的理論計(jì)算公式如下：

(5)

其中:a為燃速系數(shù)，在發(fā)動(dòng)機(jī)裝藥配方確定下，其主要與高低溫相關(guān)；ρ為推進(jìn)劑密度；C*為推進(jìn)劑特征速度；Ab為裝藥燃面，裝藥方案確定后為常量；At為喉部面積；n為壓強(qiáng)指數(shù)，通常為常數(shù)。推力系數(shù)CF=ηCFth，η為比沖效率，發(fā)動(dòng)機(jī)方案確定后為常數(shù)；CFth為理論推力系數(shù)，與喉部面積相關(guān)。

通過式(5)關(guān)系可知，在發(fā)動(dòng)機(jī)方案、裝藥配方、裝藥量和飛行彈道等確定后，內(nèi)彈道的室壓、推力和秒耗量模型主要受推進(jìn)劑溫度和喉部燒蝕率影響，將兩種偏差作極限組合，可得高溫、標(biāo)溫和低溫3種內(nèi)彈道模型，其推力模型如圖1所示：

從圖1可知，在實(shí)際飛行狀態(tài)下，內(nèi)彈道模型在3種工況下推力模值和耗盡時(shí)間偏差很大，由于該偏差無法提前準(zhǔn)確辨識(shí)，所以難以保證固體發(fā)動(dòng)機(jī)耗盡關(guān)機(jī)的制導(dǎo)精度。

針對(duì)上述問題，考慮到圖1中三種工況下推力曲線與X軸圍成面積偏差較小，即不同工況下推力速度模值增量基本相等，故本文提出以推力速度模值增量替代傳統(tǒng)的飛行時(shí)間為自變量，建立新的動(dòng)力學(xué)模型，并在實(shí)際飛行中，通過對(duì)比速度增量模值與高低溫三種工況下的預(yù)示模值，在線修正當(dāng)前內(nèi)彈道模型，最后采用序列凸優(yōu)化算法求解該最優(yōu)控制問題。

2 模型轉(zhuǎn)換

2.1 動(dòng)力學(xué)轉(zhuǎn)換

根據(jù)固體發(fā)動(dòng)機(jī)預(yù)示的內(nèi)彈道模型，建立速度模值增量E(即推力產(chǎn)生的加速度模值積分)的動(dòng)力學(xué)模型如下：

(6)

其中:m為質(zhì)量，F(xiàn)(t)為預(yù)示的內(nèi)彈道推力，qm(t)為預(yù)示的內(nèi)彈道秒耗量。由于內(nèi)彈道推力F(t)大于0，可知速度模值增量E為單調(diào)遞增。

(7)

為減小耗盡時(shí)間偏差對(duì)制導(dǎo)精度的影響，本節(jié)提出將方程(1)中的自變量時(shí)間t替換為速度模值增量E，修正后的模型如下：

(8)

(9)

其中：X(E)=[r,λ,f,V,θ,σ]T，U=[α,β]T。

飛行過程中真實(shí)推力速度模值增量計(jì)算值為：

(10)

2.2 內(nèi)彈道模型修正

考慮工程實(shí)踐，將理論上發(fā)動(dòng)機(jī)溫度的連續(xù)性辨識(shí)問題轉(zhuǎn)換為三分類(高溫、標(biāo)溫、低溫)的離散辨識(shí)問題。本節(jié)提出在飛行前裝訂式(7)的內(nèi)彈道預(yù)示模型，通過式(7)中飛行時(shí)間與速度模值增量的函數(shù)，可得當(dāng)前速度增量下高溫、標(biāo)溫、低溫內(nèi)彈道的飛行時(shí)間，通過對(duì)比飛行時(shí)間修正發(fā)動(dòng)機(jī)內(nèi)彈道模型。內(nèi)彈道修正模型如下所示：

F(E)=

(11)

(12)

qm(E)=

(13)

2.3 最優(yōu)控制問題轉(zhuǎn)換

通過上述轉(zhuǎn)換，可將式(4)中最優(yōu)控制問題轉(zhuǎn)換為如下最優(yōu)控制問題：

(14)

其中:Umax為控制量U的二范數(shù)約束，θf為耗盡點(diǎn)處目標(biāo)當(dāng)?shù)貜椀纼A角，rf為耗盡點(diǎn)處目標(biāo)地心距，σf為耗盡點(diǎn)處目標(biāo)彈道偏角。

3 序列凸化算法

采用序列凸優(yōu)化算法對(duì)式(14)中最優(yōu)控制問題進(jìn)行求解。針對(duì)動(dòng)力學(xué)模型的非凸等式約束，對(duì)動(dòng)力學(xué)模型進(jìn)行線性化處理，可得如下的線性系統(tǒng)方程：

(15)

為避免初值迭代時(shí)不可行問題，引入虛擬控制量w，可將線性系統(tǒng)方程轉(zhuǎn)化為：

(16)

對(duì)式(16)采用梯形積分離散化，可得如下離散動(dòng)力學(xué)方程：

(17)

其中:上角標(biāo)k為第k次迭代，下角標(biāo)i為速度增量模值離散后第i個(gè)值。

優(yōu)化目標(biāo)函數(shù)轉(zhuǎn)化為：

(18)

為了提升系統(tǒng)動(dòng)力學(xué)方程序列凸化的有效性，對(duì)迭代狀態(tài)施加信賴域約束，約束形式如下：

(19)

由于式(14)中不等式約束均為凸約束，故最優(yōu)控制問題已轉(zhuǎn)換為凸優(yōu)化問題。

4 仿真驗(yàn)證

為驗(yàn)證本文提出制導(dǎo)算法的精度與可靠性，選擇某型固體運(yùn)載火箭的三級(jí)飛行段作為本文仿真對(duì)象。本節(jié)分別對(duì)基于速度增量和基于時(shí)間的序列凸優(yōu)化方法進(jìn)行仿真，并對(duì)內(nèi)彈道修正前后的制導(dǎo)精度進(jìn)行了對(duì)比仿真。

4.1 仿真參數(shù)

以離線規(guī)劃的三級(jí)彈道作為初始參考軌跡，氣動(dòng)系數(shù)按風(fēng)洞試驗(yàn)數(shù)據(jù)的插值表給出，質(zhì)心位置和轉(zhuǎn)動(dòng)慣量按質(zhì)量特性插值表給出。

仿真初始參數(shù)分布和關(guān)機(jī)點(diǎn)約束條件如下表所示：

表1 初始值參數(shù)及偏差

表2 關(guān)機(jī)點(diǎn)約束條件

仿真軟件采用Matlab R2017a，凸優(yōu)化求解器采用SDPT3，制導(dǎo)采樣周期為50ms。

4.2 自變量轉(zhuǎn)換對(duì)比

本節(jié)仿真中，對(duì)比以速度模值增量和以飛行時(shí)間為自變量的序列凸優(yōu)化方法的耗盡關(guān)機(jī)處制導(dǎo)精度。仿真中，裝訂標(biāo)溫內(nèi)彈道預(yù)示模型，未采用內(nèi)彈道修正。以速度增量和時(shí)間為變量的推力模型如下圖所示：

圖3 推力與飛行時(shí)間關(guān)系

從圖2～3可知，耗盡點(diǎn)速度模值增量偏差更小，而耗盡時(shí)間偏差更大，且以飛行時(shí)間為自變量時(shí)，推力與溫度的非線性關(guān)系更強(qiáng)，難以修正。

圖2 推力與速度模值增量關(guān)系

在仿真中，由于耗盡點(diǎn)速度模值增量存在偏差，故在高、低溫工況下耗盡點(diǎn)Ef均取標(biāo)溫下的值。

采用基于速度增量的序列凸優(yōu)化制導(dǎo)算法，高度和當(dāng)?shù)貜椀纼A角仿真結(jié)果如下：

圖5 高低溫下當(dāng)?shù)貜椀纼A角仿真曲線

采用基于時(shí)間的序列凸優(yōu)化制導(dǎo)算法，高度和當(dāng)?shù)貜椀纼A角仿真結(jié)果如下：

圖6 高低溫下高度

圖7 當(dāng)?shù)貜椀纼A角偏差

從圖4～7中可知，本文基于速度模值增量的序列凸優(yōu)化制導(dǎo)算法在發(fā)動(dòng)機(jī)耗盡點(diǎn)處的高度和當(dāng)?shù)貜椀纼A角偏差均更小，制導(dǎo)精度更高。

圖4 高低溫下高度仿真曲線

4.3 內(nèi)彈道修正對(duì)比

本節(jié)仿真中，對(duì)比內(nèi)彈道修正前后的耗盡關(guān)機(jī)點(diǎn)的制導(dǎo)精度。仿真中，裝訂高溫、標(biāo)溫、低溫內(nèi)彈道預(yù)示模型，考慮到高低溫內(nèi)彈道在發(fā)動(dòng)機(jī)工作下降段難以修正，且下降段時(shí)間短，對(duì)制導(dǎo)精度影響不大，故在下降段均采用標(biāo)溫內(nèi)彈道模型。修正前后的內(nèi)彈道模型如下圖所示：

圖9 修正后推力與速度增量關(guān)系

圖10 修正前時(shí)間對(duì)速度增量導(dǎo)數(shù)與速度增量關(guān)系

圖11 修正后時(shí)間對(duì)速度增量導(dǎo)數(shù)與速度增量關(guān)系

圖12 修正前質(zhì)量與速度增量關(guān)系

從圖8～13可知，內(nèi)彈道預(yù)示模型修正后，可有效降低相對(duì)實(shí)際內(nèi)彈道的偏差，且易于實(shí)際工程應(yīng)用。

圖8 修正前推力與速度增量關(guān)系

考慮內(nèi)彈道高低溫偏差，對(duì)三級(jí)飛行段初始點(diǎn)質(zhì)量、高度、速度、當(dāng)?shù)貜椀纼A角等進(jìn)行極限拉偏仿真。采用基于視速度增量的序列凸優(yōu)化制導(dǎo)算法，內(nèi)彈道修正前耗盡關(guān)機(jī)點(diǎn)制導(dǎo)偏差如下圖所示：

圖15 修正前關(guān)機(jī)點(diǎn)當(dāng)?shù)貜椀纼A角偏差

采用基于視速度增量的序列凸優(yōu)化制導(dǎo)算法，內(nèi)彈道修正后關(guān)機(jī)點(diǎn)制導(dǎo)偏差如下圖所示：

圖16 修正后關(guān)機(jī)點(diǎn)高度偏差

圖17 修正后關(guān)機(jī)點(diǎn)當(dāng)?shù)貜椀纼A角偏差

從圖14～17中可知，內(nèi)彈道修正前耗盡關(guān)機(jī)點(diǎn)高度最大偏差為550m，當(dāng)?shù)貜椀纼A角最大偏差為0.11°；內(nèi)彈道修正后高度耗盡關(guān)機(jī)點(diǎn)最大偏差為280m，當(dāng)?shù)貜椀纼A角最大偏差為0.1°。可知，修正后耗盡關(guān)機(jī)點(diǎn)高度和當(dāng)?shù)貜椀纼A角偏差顯著降低，故內(nèi)彈道修正方式對(duì)制導(dǎo)精度有明顯提高。

圖14 修正前關(guān)機(jī)點(diǎn)高度偏差

5 結(jié)論

針對(duì)固體發(fā)動(dòng)機(jī)在高低溫下的耗盡時(shí)間、推力模型和秒耗量散布大的問題，提出了一種以速度模值增量為變量的序列凸優(yōu)化算法，并通過對(duì)比速度增量模值，對(duì)當(dāng)前內(nèi)彈道模型進(jìn)行在線修正。通過仿真驗(yàn)證，該算法比以時(shí)間為變量的序列凸優(yōu)化算法的制導(dǎo)精度更高，且以速度模值增量為判據(jù)的內(nèi)彈道修正算法能顯著提高制導(dǎo)精度，具有較高的工程應(yīng)用價(jià)值。