沖壓發(fā)動(dòng)機(jī)導(dǎo)彈爬升軌跡與推力調(diào)節(jié)規(guī)律優(yōu)化

尚騰，谷良賢，趙吉松，龔春林

(西北工業(yè)大學(xué)航天學(xué)院，陜西西安 710072)

引言

沖壓發(fā)動(dòng)機(jī)導(dǎo)彈的彈道設(shè)計(jì)不同于固體火箭發(fā)動(dòng)機(jī)導(dǎo)彈的彈道設(shè)計(jì)，主要體現(xiàn)在沖壓發(fā)動(dòng)機(jī)推力與彈道的強(qiáng)耦合性。一方面，沖壓發(fā)動(dòng)機(jī)的內(nèi)部參數(shù)和性能指標(biāo)隨導(dǎo)彈的飛行速度、飛行高度、攻角及實(shí)際進(jìn)入發(fā)動(dòng)機(jī)的空氣流量而變化;另一方面，沖壓發(fā)動(dòng)機(jī)的比沖及推力直接影響導(dǎo)彈的質(zhì)量和飛行性能。所以，在設(shè)計(jì)以沖壓發(fā)動(dòng)機(jī)為動(dòng)力的導(dǎo)彈彈道時(shí)必須考慮推力調(diào)節(jié)規(guī)律的設(shè)計(jì)［1］。沖壓發(fā)動(dòng)機(jī)的推力調(diào)節(jié)多采用控制余氣系數(shù)的方法，其特點(diǎn)是實(shí)現(xiàn)方便且有利于保證沖壓發(fā)動(dòng)機(jī)燃燒的穩(wěn)定性［2］。

以沖壓發(fā)動(dòng)機(jī)為動(dòng)力的超聲速、中遠(yuǎn)程戰(zhàn)術(shù)導(dǎo)彈，通常用助推器將導(dǎo)彈加速到達(dá)預(yù)定的高度和速度，沖壓發(fā)動(dòng)機(jī)再接力工作，將導(dǎo)彈加速到規(guī)定的巡航速度和巡航高度進(jìn)行巡航飛行。爬升彈道的設(shè)計(jì)優(yōu)化需要解決如何最優(yōu)地使導(dǎo)彈從接力點(diǎn)(沖壓發(fā)動(dòng)機(jī)開(kāi)始工作)爬升到巡航狀態(tài)。文獻(xiàn)［3-4］在進(jìn)行沖壓發(fā)動(dòng)機(jī)導(dǎo)彈爬升彈道優(yōu)化的研究中，余氣系數(shù)只根據(jù)導(dǎo)彈飛行馬赫數(shù)進(jìn)行了適當(dāng)調(diào)節(jié)，并沒(méi)有參與優(yōu)化。而沖壓發(fā)動(dòng)機(jī)導(dǎo)彈爬升過(guò)程中飛行速度、高度、導(dǎo)彈姿態(tài)等均在較大范圍內(nèi)變化，只有把沖壓發(fā)動(dòng)機(jī)的推力調(diào)節(jié)規(guī)律和導(dǎo)彈彈道設(shè)計(jì)協(xié)調(diào)起來(lái)，進(jìn)行一體化優(yōu)化設(shè)計(jì)，才能獲得合理的結(jié)果。

直接法對(duì)目標(biāo)函數(shù)的解析性質(zhì)不作苛刻要求，且適合各種精度的彈道模型，本文采用直接法求解沖壓發(fā)動(dòng)機(jī)導(dǎo)彈爬升這一最優(yōu)控制問(wèn)題，建立了爬升軌跡與推力調(diào)節(jié)規(guī)律共同參與優(yōu)化的一體化優(yōu)化設(shè)計(jì)模型。為了克服直接法對(duì)初始點(diǎn)選取敏感，容易收斂到局部最優(yōu)解的缺點(diǎn)，將粒子群算法(PSO)與變尺度算法(BFGS)串聯(lián)結(jié)合，尋優(yōu)前期利用PSO較強(qiáng)的全局搜索能力，搜索至全局最優(yōu)點(diǎn)附近;尋優(yōu)后半階段利用BFGS較強(qiáng)的局部尋優(yōu)能力，以提高搜索精度。設(shè)計(jì)結(jié)果表明，本文提出的爬升彈道一體化設(shè)計(jì)方法能夠成功地求解沖壓發(fā)動(dòng)機(jī)爬升彈道問(wèn)題，具有較好的工程應(yīng)用價(jià)值。

1 數(shù)學(xué)模型

1.1 爬升彈道運(yùn)動(dòng)模型

將導(dǎo)彈看作可控質(zhì)點(diǎn)，導(dǎo)彈在鉛垂平面內(nèi)的運(yùn)動(dòng)方程組為:

式中，V為導(dǎo)彈速度;P為推力;m為導(dǎo)彈質(zhì)量;X為阻力;Y為升力;α，θ分別為攻角和彈道傾角;x為縱向射程;y為高度;mf為發(fā)動(dòng)機(jī)燃料消耗率。

設(shè)沖壓發(fā)動(dòng)機(jī)導(dǎo)彈的爬升彈道光滑且二階連續(xù)，則可用三次樣條曲線插值擬合其爬升軌跡。下面推導(dǎo)導(dǎo)彈按照預(yù)設(shè)飛行軌跡y*(x)飛行時(shí)的需用攻角。

由式(3)、式(4)可得:

給定導(dǎo)彈爬升飛行的預(yù)設(shè)飛行軌跡y*(x)，則可以利用式(6)～式(9)求出按此軌跡飛行時(shí)的需用攻角。

對(duì)式(1)～式(5)進(jìn)行積分求解，飛行攻角取需用攻角，即可得到導(dǎo)彈的實(shí)際飛行彈道y(x)。其中，沖壓發(fā)動(dòng)機(jī)的性能隨導(dǎo)彈的飛行速度、飛行高度、迎角及實(shí)際進(jìn)入發(fā)動(dòng)機(jī)的空氣流量而變化。

1.2 沖壓發(fā)動(dòng)機(jī)推力模型及推力調(diào)節(jié)規(guī)律

沖壓發(fā)動(dòng)機(jī)的推力計(jì)算公式為:

式中，ma為空氣流量;ve為發(fā)動(dòng)機(jī)噴口氣流速度;v為飛行器速度;pe為噴口壓強(qiáng);pa為環(huán)境壓強(qiáng);αf為余氣系數(shù)，αf=ma/(Lmf);L為單位燃料消耗的理論空氣量。

以某頭部進(jìn)氣沖壓發(fā)動(dòng)機(jī)為例，采用一元流理論，得到?jīng)_壓發(fā)動(dòng)機(jī)的推力、燃料秒流量隨飛行高度、馬赫數(shù)、攻角及余氣系數(shù)的變化關(guān)系［2］，即:

式中，f1，f2以插值表的形式給出。發(fā)動(dòng)機(jī)的推力調(diào)節(jié)就是通過(guò)調(diào)節(jié)供油量，控制余氣系數(shù)按照所需的規(guī)律進(jìn)行變化。

2 優(yōu)化問(wèn)題表達(dá)與優(yōu)化算法

2.1 優(yōu)化問(wèn)題描述

以沖壓發(fā)動(dòng)機(jī)為動(dòng)力的導(dǎo)彈爬升段設(shè)計(jì)優(yōu)化的目標(biāo)是使導(dǎo)彈達(dá)到巡航高度并以巡航馬赫數(shù)飛行，且要求燃料消耗量少，飛行時(shí)間短。時(shí)間最短和燃料最省是兩個(gè)相互矛盾的彈道優(yōu)化目標(biāo)［4］，本文研究燃料最省的情況。

彈道優(yōu)化的初始條件是給定的接力點(diǎn)位置、接力馬赫數(shù);終端約束為給定的巡航高度、彈道傾角(0°)、最小巡航馬赫數(shù);在導(dǎo)彈爬升過(guò)程中有最大攻角的限制。彈道優(yōu)化的目標(biāo)函數(shù)和約束條件的數(shù)學(xué)描述為:

式中，ycruise，Mamin為給定的巡航高度和最小巡航馬赫數(shù)。

2.2 彈道樣條離散

在三次樣條曲線中，兩個(gè)節(jié)點(diǎn)之間的曲線由三次多項(xiàng)式擬合生成，整個(gè)樣條曲線由各段三次多項(xiàng)式曲線組成，相鄰兩段在節(jié)點(diǎn)處二階連續(xù)。本文構(gòu)造三次樣條時(shí)采用自然邊界條件，即兩端點(diǎn)處的二階導(dǎo)數(shù)為零。采用三次樣條插值方法可以對(duì)彈道優(yōu)化這一最優(yōu)控制問(wèn)題進(jìn)行離散并參數(shù)化［5-6］。

將彈道沿射程x分為n段，形成n+1個(gè)節(jié)點(diǎn)［x0，x1，…，xn］，優(yōu)化變量為各節(jié)點(diǎn)處的飛行軌跡縱坐標(biāo)yi和余氣系數(shù)αfi(i=0，1，…，n)。通過(guò)樣條插值，可構(gòu)建出導(dǎo)彈在豎直平面內(nèi)的軌跡y*(x)和余氣系數(shù)變化規(guī)律α*f(x)。

2.3 PSO-BFGS串聯(lián)優(yōu)化方法

PSO算法有著較強(qiáng)的全局搜索能力，但局部搜索能力較差，搜索精度不夠高［7］。BFGS變尺度法是求解無(wú)約束優(yōu)化問(wèn)題最有效的算法之一，搜索精度較高，數(shù)值穩(wěn)定性好。結(jié)合兩者優(yōu)點(diǎn)，將PSO算法與BFGS算法串聯(lián)形成PSO-BFGS優(yōu)化算法。

PSO算法首先在可行解空間和速度空間隨機(jī)初始化粒子群(Particle Swarm)，即確定粒子(Particle)的初始位置和初始速度。其中，位置用于表征問(wèn)題解。設(shè)d維搜索空間中的第i個(gè)粒子的位置和速度分別表示為 xi= ［xi，1，xi，2，…，xi，d］和 vi= ［vi，1，vi，2，…，vi，d］。通過(guò)評(píng)價(jià)各粒子的適應(yīng)度值，確定每個(gè)粒子所經(jīng) 過(guò) 的最佳位置 pi= ［pi，1，pi，2，…，pi，d］，以及群體所發(fā)現(xiàn)的最佳位置 pg= ［pg，1，pg，2，…，pg，d］，再按下式更新各粒子的速度和位置:

式中，w為慣性權(quán)重因子;c1和c2為正的加速常數(shù);r1和r2為在0～1之間均勻分布的隨機(jī)數(shù)。

設(shè)粒子的速度區(qū)間為［vmin，j，vmax，j］，位置范圍為［xmin，j，xmax，j］，當(dāng)粒子速度超出限制時(shí)取邊界值。位置超出范圍時(shí)，按下式確定［8］:

式中，f為彈性因子(0≤ f≤1)，表示粒子在超出邊界時(shí)，被彈回到可行解空間內(nèi)的強(qiáng)弱程度。采用此式有利于求解最優(yōu)解在邊界附近的問(wèn)題。

BFGS的基本思想為［8］:迭代過(guò)程中，在 x(k+1)處按下式生成對(duì)稱正定矩陣Hk+1，以代替海森矩陣的逆陣:

式中，g(k)=▽ f(x(k))，由中心差分格式給出;Δg(k)=g(k+1)－ g(k);Δx(k)=x(k+1)－ x(k)。然后，以p(k+1)= －Hk+1g(k+1)為x(k+1)處的搜索方向進(jìn)行搜索。BFGS對(duì)一維搜索算法精度要求較高，文中采用三次曲線擬合法［8］。

PSO-BFGS串聯(lián)混合算法的優(yōu)化流程如圖1所示。

圖1 PSO-BFGS算法流程圖

針對(duì)爬升彈道的末端約束，通過(guò)在適應(yīng)度函數(shù)中增加懲罰函數(shù)來(lái)處理。適應(yīng)度函數(shù)選為:

式中，cy，cθ，cMa分別為終端高度、彈道傾角和速度的懲罰系數(shù)。則爬升彈道優(yōu)化問(wèn)題轉(zhuǎn)換為無(wú)約束優(yōu)化問(wèn)題，可由PSO-BFGS混合算法尋優(yōu)求解。

3 設(shè)計(jì)實(shí)例

采用上述彈道優(yōu)化方法對(duì)沖壓發(fā)動(dòng)機(jī)導(dǎo)彈的爬升彈道進(jìn)行優(yōu)化。算例設(shè)計(jì)參考了文獻(xiàn)［1］和文獻(xiàn)［3］，主要計(jì)算條件為:接力點(diǎn)馬赫數(shù)2.2，接力點(diǎn)高度2 km，巡航高度15 km，終點(diǎn)彈道傾角0°，終點(diǎn)馬赫數(shù)大于2.4，爬升終點(diǎn)距接力點(diǎn)的水平距離為25 km。為了觀察軌跡與推力調(diào)節(jié)一體化設(shè)計(jì)在節(jié)省燃料方面的貢獻(xiàn)，本文對(duì)以下兩種方案進(jìn)行了優(yōu)化。

方案1:僅爬升軌跡參與彈道優(yōu)化，余氣系數(shù)變化規(guī)律預(yù)先設(shè)定為:

此方案類似文獻(xiàn)［3-4］中采用的方案。

方案2:爬升軌跡和余氣系數(shù)都參與彈道優(yōu)化，即前文提出的一體化設(shè)計(jì)方案。

方案1中，彈道采用6個(gè)節(jié)點(diǎn)作為控制點(diǎn)，沿x向均勻分布，y1取 2 000.0，y2～ y5∈［2 000.0，15 000.0］，y6∈［14 500.0，15 500.0］。方案 2 中，彈道離散同方案1，余氣系數(shù)離散控制點(diǎn)也取6個(gè)，沿 x向均勻分布，αf1取0.9，αf2～ αf6∈［0.9，2.0］。方案2的優(yōu)化問(wèn)題規(guī)模是方案1的2倍，圖2給出了方案2的求解收斂過(guò)程。可以看出，PSO尋優(yōu)前期，適應(yīng)度值(J1)迅速下降，更新n=40代以后，適應(yīng)度值就無(wú)明顯改進(jìn)，收斂至37.536 0;BFGS方法經(jīng)過(guò)23次迭代，收斂至35.149 7，較PSO算法的優(yōu)化結(jié)果提高6.36%。優(yōu)化求解得到的兩種方案的爬升軌跡和飛行參數(shù)見(jiàn)圖3～圖6。

圖2 方案2的求解收斂過(guò)程

圖3 兩種方案最優(yōu)爬升彈道

圖4 爬升中馬赫數(shù)變化

圖5 爬升中攻角變化

圖6 爬升中推力變化

圖6方案1曲線中點(diǎn)劃線表示:推力值在導(dǎo)彈飛行馬赫數(shù)接近2.41時(shí)出現(xiàn)震蕩。兩種方案爬升飛行的末端參數(shù)及飛行過(guò)程中的攻角絕對(duì)值的最大值見(jiàn)表1。

表1 兩種方案的優(yōu)化結(jié)果比較

仿真結(jié)果顯示，兩種設(shè)計(jì)方案的設(shè)計(jì)結(jié)果均達(dá)到了設(shè)計(jì)要求。方案2中，爬升過(guò)程總消耗燃料量為35.149 3 kg，比方案 1 減少 5.102 4 kg，減少了12.68%，爬升時(shí)間增加5.0 s。燃料消耗的減少對(duì)于減少導(dǎo)彈發(fā)射質(zhì)量或增加導(dǎo)彈的射程是十分有意義的，且接力爬升段在沖壓發(fā)動(dòng)機(jī)的整個(gè)飛行過(guò)程中所占比例較小，爬升段時(shí)間的少量增加是可以接受的。

對(duì)比兩種方案下的爬升彈道:在爬升前半段，方案1中導(dǎo)彈作加速爬升飛行，速度已達(dá)到末端速度要求;方案2通過(guò)優(yōu)化余氣系數(shù)，使導(dǎo)彈作減速爬升飛行。在爬升后半段，方案1中，導(dǎo)彈速度略有增加;而方案2中，導(dǎo)彈呈加速飛行狀態(tài)，由于導(dǎo)彈已經(jīng)接近平飛高度，加速效率高?？梢?jiàn)，采用方案2，導(dǎo)彈以較小的平均飛行速度完成爬升過(guò)程，可以實(shí)現(xiàn)節(jié)省燃料，進(jìn)而達(dá)到減小導(dǎo)彈發(fā)射質(zhì)量或增加射程的目的。

4 結(jié)束語(yǔ)

本文采用直接法建立了飛行軌跡形狀和推力調(diào)節(jié)規(guī)律共同參與優(yōu)化的沖壓發(fā)動(dòng)機(jī)導(dǎo)彈的爬升彈道一體化優(yōu)化模型，并采用PSO-BFGS混合算法進(jìn)行求解。優(yōu)化結(jié)果表明:采用樣條插值直接離散方法和PSO-BFGS混合算法能夠較好地求解沖壓發(fā)動(dòng)機(jī)導(dǎo)彈爬升彈道優(yōu)化問(wèn)題。在沖壓發(fā)動(dòng)機(jī)導(dǎo)彈爬升彈道設(shè)計(jì)優(yōu)化中，推力調(diào)節(jié)規(guī)律與爬升軌跡一體化優(yōu)化設(shè)計(jì)，可獲得先減速再加速的爬升方案，采用此方案，導(dǎo)彈能以較小的平均速度完成爬升過(guò)程，可以顯著節(jié)省燃料，提高導(dǎo)彈性能。

［1］ 劉恒軍，沙建科，王華.沖壓發(fā)動(dòng)機(jī)導(dǎo)彈彈道多目標(biāo)優(yōu)化［J］.系統(tǒng)仿真學(xué)報(bào)，2009，21(9):2764-2766.

［2］ 劉興洲，于守志.飛航導(dǎo)彈動(dòng)力裝置(上)［M］.北京:宇航出版社，1992.

［3］ 王華，楊存富，劉恒軍.沖壓發(fā)動(dòng)機(jī)為動(dòng)力導(dǎo)彈爬升彈道優(yōu)化［J］.彈箭與制導(dǎo)學(xué)報(bào)，2008，28(3):185-188.

［4］ 王華，楊存富，劉恒軍.以沖壓發(fā)動(dòng)機(jī)為動(dòng)力的導(dǎo)彈爬升彈道研究［J］.現(xiàn)代防御技術(shù)，2008，36(4):27-30.

［5］ 陳勝琪.基于樣條曲線的機(jī)動(dòng)彈頭再入軌跡優(yōu)化研究［J］.飛行力學(xué)，2009，27(1):63-65.

［6］ Ken Cote.Complex 3D flight trajectory generation and tracking using cubic splines［C］//AIAA’s 1st Technical Conference and Workshop on Unmanned Aerospace Vehicles.Portsmouth，USA，2002.

［7］ Rhonald M J，Roy J H.Hybrid particle swarm-pattern search optimizer for aerospace propulsion applications//46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit.Nashville，USA，2010.

［8］ 粟塔山.最優(yōu)化計(jì)算原理與算法程序設(shè)計(jì)［M］.長(zhǎng)沙:國(guó)防科技大學(xué)出版社，2001.