一種基于動(dòng)壓變化的火箭助推段彈道設(shè)計(jì)方法

王 智，馬保海，張 婕，熊 偉，宋劍爽

0 引 言

為滿足特定的技術(shù)驗(yàn)證目的，諸如新型熱防護(hù)材料的考核驗(yàn)證，需要飛行器全程在大氣層內(nèi)高速飛行來創(chuàng)造需要的飛行環(huán)境。飛行器全程在大氣層內(nèi)高速飛行時(shí)，一方面主動(dòng)段動(dòng)壓數(shù)值大，變化幅度大；另一方面，傳統(tǒng)主動(dòng)段一級(jí)彈道一次負(fù)攻角轉(zhuǎn)彎下壓的形式已不能滿足飛行器全程大氣層內(nèi)飛行的要求。針對(duì)以上問題，根據(jù)飛行器的飛行特點(diǎn)，本文在彈道設(shè)計(jì)時(shí)采用基于動(dòng)壓變化的二次下壓特殊低彈道形式，同時(shí)采用序列二次規(guī)劃算法（Sequential Quadratic Programming，SQP），對(duì)基于動(dòng)壓變化的二次下壓彈道進(jìn)行優(yōu)化設(shè)計(jì)，可以避免在大動(dòng)壓同時(shí)出現(xiàn)大攻角，為飛行器提供良好的飛行環(huán)境。通過與傳統(tǒng)彈道設(shè)計(jì)方法的對(duì)比分析，驗(yàn)證了所提彈道設(shè)計(jì)方法在滿足給定約束的情況下，減小了主動(dòng)段最大動(dòng)壓，縮小了動(dòng)壓變化范圍，具有較好的工程參考價(jià)值。

1 主動(dòng)段彈道模型

考慮地球?yàn)樾D(zhuǎn)的均質(zhì)橢球體，忽略風(fēng)的干擾，采用標(biāo)準(zhǔn)大氣模型，在發(fā)射坐標(biāo)系內(nèi)建立飛行器三自由度質(zhì)心運(yùn)動(dòng)模型。

1.1 彈道計(jì)算模型

飛行器主動(dòng)段運(yùn)動(dòng)過程是一個(gè)變質(zhì)量飛行過程，其三自由度質(zhì)心運(yùn)動(dòng)方程為[1]

式中 O為發(fā)射坐標(biāo)系原點(diǎn)與發(fā)射點(diǎn)的連接點(diǎn)；x，y，z為飛行器在發(fā)射坐標(biāo)系中的位置三分量；Ox軸為在發(fā)射點(diǎn)水平面內(nèi)指向發(fā)射瞄準(zhǔn)方向；Oy軸為垂直于發(fā)射點(diǎn)水平面指向上方；Oz軸與xOy面相垂直并構(gòu)成右手直角坐標(biāo)系； Vx， Vy， Vz為飛行器在發(fā)射坐標(biāo)系中的速度三分量； ROx， ROy， ROz為發(fā)射點(diǎn)地心矢徑在發(fā)射坐標(biāo)系中的三分量； gx， gy， gz為引力加速度在發(fā)射坐標(biāo)系中的三分量；P為發(fā)動(dòng)機(jī)推力；X，Y，Z分別為氣動(dòng)阻力、升力和側(cè)向力；m為飛行器質(zhì)量； m0為飛行器初始質(zhì)量；m˙為質(zhì)量秒消耗量； GB為箭體坐標(biāo)系至發(fā)射坐標(biāo)系的轉(zhuǎn)換矩陣； GV為速度坐標(biāo)系至發(fā)射坐標(biāo)系的轉(zhuǎn)換矩陣；A為離心慣性力矩陣；B為哥氏慣性力矩陣。

1.2 飛行程序角設(shè)計(jì)

飛行程序角設(shè)計(jì)是彈道優(yōu)化設(shè)計(jì)的關(guān)鍵，也是彈道參數(shù)優(yōu)化問題的主要對(duì)象[2～4]。

飛行器全程在大氣層中高速飛行時(shí)，主動(dòng)段動(dòng)壓數(shù)值大，變化幅度大，為便于程序角設(shè)計(jì)，采用跟蹤程序攻角的方式。同時(shí)在進(jìn)行程序攻角設(shè)計(jì)時(shí)，為使飛行程序能反映動(dòng)壓及其變化情況，主要根據(jù)動(dòng)壓的變化情況進(jìn)行設(shè)計(jì)。主動(dòng)段程序攻角設(shè)計(jì)曲線見圖1。

圖1 主動(dòng)段程序攻角與動(dòng)壓曲線Fig.1 Curve of Attack Angle and Dynamic Pressure for Boost-phase Trajectory

點(diǎn)火起飛后，零攻角垂直飛行至1t時(shí)刻。1t時(shí)刻后，動(dòng)壓Q呈現(xiàn)增大趨勢(shì)，為避免在大動(dòng)壓時(shí)出現(xiàn)大攻角，在動(dòng)壓增大一定幅度后，即2t時(shí)刻，攻角恢復(fù)到零，從而實(shí)現(xiàn)飛行器在出現(xiàn)動(dòng)壓最大值及其附近時(shí)的攻角數(shù)值比較小，為飛行器提供良好的飛行環(huán)境。在1t～2t時(shí)間段，進(jìn)行第1次程序負(fù)攻角轉(zhuǎn)彎飛行，程序攻角設(shè)計(jì)為如下形式：

式中mα為第1次轉(zhuǎn)彎攻角絕對(duì)值的最大值；mt為攻角達(dá)到極值mα的時(shí)間。

在23～tt時(shí)間段，飛行器保持零攻角重力轉(zhuǎn)彎飛行，3t取在接近動(dòng)壓最大值處。

在35～tt時(shí)間段，飛行器進(jìn)行第2次程序負(fù)攻角轉(zhuǎn)彎飛行。為有效調(diào)整飛行動(dòng)壓的變化情況，通過時(shí)間變量4t將第2次程序負(fù)攻角轉(zhuǎn)彎飛行分成兩段，即第2次程序負(fù)攻角轉(zhuǎn)彎的第1段和第2段。5t取在動(dòng)壓下降段變平緩處。程序角設(shè)計(jì)為如下形式：

式中 ?cx0，?cx1分別為t3時(shí)刻和t4時(shí)刻的程序角；?˙cx1，?˙cx2分別為二次下壓第一段和第二段俯仰程序角變化率；ωz為地球自轉(zhuǎn)角速度在發(fā)射坐標(biāo)系z(mì)軸的分量；θ為彈道傾角。

在56～tt時(shí)間段，飛行器進(jìn)行攻角歸零飛行，以使在分離時(shí)刻攻角為零，為分離創(chuàng)造良好的條件。6t取在發(fā)動(dòng)機(jī)關(guān)機(jī)前4～5 s之間。

綜上所述，時(shí)間點(diǎn)1t～6t、mα、二次下壓第1段和第2段俯仰程序角變化率cx1?˙和cx2?˙決定了主動(dòng)段程序角的變化規(guī)律。依據(jù)各種設(shè)計(jì)指標(biāo)要求，mα、4t、1cx?˙和cx2?˙由彈道優(yōu)化得到。由于第一次下壓段動(dòng)壓較小，第二次下壓段動(dòng)壓較大，為避免第二次負(fù)攻角轉(zhuǎn)彎出現(xiàn)較大的攻角，需要 ?˙cx1＜?˙cx2。

2 優(yōu)化問題建模

2.1 性能指標(biāo)

為給全程大氣層內(nèi)高速飛行器創(chuàng)造較好的飛行環(huán)境，減小姿態(tài)控制系統(tǒng)和載荷設(shè)計(jì)難度，需要避免在大動(dòng)壓時(shí)出現(xiàn)大攻角。因此，在進(jìn)行二次下壓彈道優(yōu)化設(shè)計(jì)時(shí)，目標(biāo)函數(shù)選擇為

2.2 優(yōu)化控制變量

選取變量mα，4t，cx1?˙，cx2?˙和射向角0A作為優(yōu)化參數(shù)，因此，優(yōu)化控制變量為

2.3 約束條件

考慮可實(shí)現(xiàn)性，控制變量應(yīng)滿足如下的約束條件：

此外，由于彈道采用低彈道形式，要求飛行高度滿足：

式中ch為給定的最大飛行高度。

同時(shí)落點(diǎn)航程偏差LΔ和落點(diǎn)橫程偏差HΔ滿足：

3 主動(dòng)段彈道優(yōu)化

彈道優(yōu)化分為靜態(tài)優(yōu)化和動(dòng)態(tài)優(yōu)化問題。嚴(yán)格來講，彈道優(yōu)化是一個(gè)動(dòng)態(tài)優(yōu)化問題，但由于彈道優(yōu)化自身的復(fù)雜性，采用動(dòng)態(tài)優(yōu)化求解面臨較大難度。因此通常將彈道優(yōu)化轉(zhuǎn)化為靜態(tài)優(yōu)化問題，再通過數(shù)值方法求解[5,6]。

彈道優(yōu)化轉(zhuǎn)化的數(shù)值方法主要有間接法和直接法兩種。間接法根據(jù) Pontryagin極小值原理將最優(yōu)控制問題轉(zhuǎn)化為兩點(diǎn)邊值問題，其求解的精度較高，且其解滿足一階最優(yōu)性必要條件。但由于飛行器數(shù)學(xué)模型復(fù)雜，協(xié)態(tài)方程與橫截條件推導(dǎo)過程較為復(fù)雜和繁瑣；同時(shí)協(xié)態(tài)量初值無實(shí)際物理意義，其初值難以估計(jì)。直接法采用離散化方法將連續(xù)的彈道優(yōu)化問題轉(zhuǎn)化為非線性規(guī)劃問題，然后采用數(shù)值方法對(duì)非線性規(guī)劃問題進(jìn)行求解。直接法不需求解最優(yōu)解的必要條件，操作簡(jiǎn)單；同時(shí)具有較好的通用性和魯棒性[7～10]。

求解非線性規(guī)劃問題的數(shù)值方法包括以可行方向法、梯度下降法、罰函數(shù)方法、動(dòng)態(tài)規(guī)劃法、SQP法等為代表的精確算法和以遺傳算法、進(jìn)化策略、模擬退火算法、神經(jīng)網(wǎng)絡(luò)算法等為代表的現(xiàn)代啟發(fā)式算法[11～14]。其中，SQP算法是一種有效的求解非線性規(guī)劃問題的工具，是飛行器彈道優(yōu)化的主流方法。

SQP算法將有約束的優(yōu)化問題通過拉格朗日函數(shù)轉(zhuǎn)化為無約束的優(yōu)化問題，并將原優(yōu)化問題通過求解Lagrange極小值方程轉(zhuǎn)化為二次規(guī)劃問題，最后通過擬牛頓迭代算法與 BFGS（Broyden Fletcher Goldfarb Shanno）更新算法完成計(jì)算[15～18]。非線性規(guī)劃問題可表示為

非線性規(guī)劃問題對(duì)應(yīng)的二次規(guī)劃子問題可表示為

式中 d為搜索方向； ? J (uk)， ? c (uk)， ? p (uk)分別為函數(shù) J (u)， c (u)和 p (u)在 uk處的梯度； Bk矩陣為L(zhǎng)agrange函數(shù)的Hessian陣的良好近似。

本文通過選擇程序角中的特征參數(shù)作為優(yōu)化變量，將彈道優(yōu)化轉(zhuǎn)化為參數(shù)優(yōu)化問題，利用SQP算法求解。彈道優(yōu)化流程如圖2所示。

圖2 彈道優(yōu)化流程Fig.2 Flow Chart of Trajectory Optimization

4 仿真分析

分別對(duì)基于動(dòng)壓變化的二次下壓優(yōu)化彈道和傳統(tǒng)彈道設(shè)計(jì)結(jié)果進(jìn)行對(duì)比分析。主要特征參數(shù)比較見圖 3～6。

圖3 主動(dòng)段攻角隨時(shí)間變化曲線Fig.3 Curve of Attack Angle for Boost-phase Trajectory

圖4 主動(dòng)段動(dòng)壓隨時(shí)間變化曲線Fig.4 Curve of Dynamic Pressure for Boost-phase Trajectory

圖5 主動(dòng)段動(dòng)壓攻角乘積隨時(shí)間變化曲線Fig.5 Curve of Product of Dynamic Pressure and Attack Angle for Boost-phase Trajectory

圖6 高度隨航程變化曲線Fig.6 Curve of Height and Range

由圖3～6可以看出，采用SQP算法對(duì)基于動(dòng)壓變化的二次下壓彈道進(jìn)行優(yōu)化設(shè)計(jì)，在滿足最大飛行高度和落點(diǎn)等飛行約束的同時(shí)，第2次下壓攻角絕對(duì)值最大值減小，最大動(dòng)壓減小5%，動(dòng)壓與攻角乘積絕對(duì)值的最大值減小8%，驗(yàn)證了所提彈道設(shè)計(jì)方法的有效性。

5 結(jié) 論

本文首先闡述了基于動(dòng)壓變化的二次下壓特殊低彈道優(yōu)化設(shè)計(jì)方法，然后通過仿真分析，獲得以下結(jié)論：

a）采用基于動(dòng)壓變化的二次下壓程序攻角設(shè)計(jì)，可使飛行程序反映動(dòng)壓及其變化情況；

b）采用基于動(dòng)壓變化的二次下壓彈道優(yōu)化設(shè)計(jì)時(shí)，在飛行約束滿足的同時(shí)，第二次下壓攻角絕對(duì)值最大值、最大動(dòng)壓和動(dòng)壓與攻角乘積絕對(duì)值最大值均有所減小，驗(yàn)證了所提彈道設(shè)計(jì)方法的有效性，具有較好的工程應(yīng)用價(jià)值。