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

0 引 言

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

1 主動段彈道模型

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

1.1 彈道計算模型

飛行器主動段運動過程是一個變質(zhì)量飛行過程，其三自由度質(zhì)心運動方程為[1]

式中 O為發(fā)射坐標(biāo)系原點與發(fā)射點的連接點；x，y，z為飛行器在發(fā)射坐標(biāo)系中的位置三分量；Ox軸為在發(fā)射點水平面內(nèi)指向發(fā)射瞄準(zhǔn)方向；Oy軸為垂直于發(fā)射點水平面指向上方；Oz軸與xOy面相垂直并構(gòu)成右手直角坐標(biāo)系； Vx， Vy， Vz為飛行器在發(fā)射坐標(biāo)系中的速度三分量； ROx， ROy， ROz為發(fā)射點地心矢徑在發(fā)射坐標(biāo)系中的三分量； gx， gy， gz為引力加速度在發(fā)射坐標(biāo)系中的三分量；P為發(fā)動機推力；X，Y，Z分別為氣動阻力、升力和側(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è)計

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

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

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

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

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

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

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

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

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

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

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

2.1 性能指標(biāo)

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

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

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

2.3 約束條件

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

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

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

同時落點航程偏差LΔ和落點橫程偏差HΔ滿足：

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

彈道優(yōu)化分為靜態(tài)優(yōu)化和動態(tài)優(yōu)化問題。嚴(yán)格來講，彈道優(yōu)化是一個動態(tài)優(yōu)化問題，但由于彈道優(yōu)化自身的復(fù)雜性，采用動態(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)化為兩點邊值問題，其求解的精度較高，且其解滿足一階最優(yōu)性必要條件。但由于飛行器數(shù)學(xué)模型復(fù)雜，協(xié)態(tài)方程與橫截條件推導(dǎo)過程較為復(fù)雜和繁瑣；同時協(xié)態(tài)量初值無實際物理意義，其初值難以估計。直接法采用離散化方法將連續(xù)的彈道優(yōu)化問題轉(zhuǎn)化為非線性規(guī)劃問題，然后采用數(shù)值方法對非線性規(guī)劃問題進行求解。直接法不需求解最優(yōu)解的必要條件，操作簡單；同時具有較好的通用性和魯棒性[7～10]。

求解非線性規(guī)劃問題的數(shù)值方法包括以可行方向法、梯度下降法、罰函數(shù)方法、動態(tài)規(guī)劃法、SQP法等為代表的精確算法和以遺傳算法、進化策略、模擬退火算法、神經(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）更新算法完成計算[15～18]。非線性規(guī)劃問題可表示為

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

式中 d為搜索方向； ? J (uk)， ? c (uk)， ? p (uk)分別為函數(shù) J (u)， c (u)和 p (u)在 uk處的梯度； Bk矩陣為Lagrange函數(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 仿真分析

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

圖3 主動段攻角隨時間變化曲線Fig.3 Curve of Attack Angle for Boost-phase Trajectory

圖4 主動段動壓隨時間變化曲線Fig.4 Curve of Dynamic Pressure for Boost-phase Trajectory

圖5 主動段動壓攻角乘積隨時間變化曲線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算法對基于動壓變化的二次下壓彈道進行優(yōu)化設(shè)計，在滿足最大飛行高度和落點等飛行約束的同時，第2次下壓攻角絕對值最大值減小，最大動壓減小5%，動壓與攻角乘積絕對值的最大值減小8%，驗證了所提彈道設(shè)計方法的有效性。

5 結(jié) 論

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

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

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