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

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

引言

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

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

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

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

1.1 爬升彈道運動模型

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

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

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

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

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

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

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

沖壓發(fā)動機的推力計算公式為:

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

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

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

2 優(yōu)化問題表達與優(yōu)化算法

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

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

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

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

2.2 彈道樣條離散

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

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

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

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

PSO算法首先在可行解空間和速度空間隨機初始化粒子群(Particle Swarm)，即確定粒子(Particle)的初始位置和初始速度。其中，位置用于表征問題解。設(shè)d維搜索空間中的第i個粒子的位置和速度分別表示為 xi= ［xi，1，xi，2，…，xi，d］和 vi= ［vi，1，vi，2，…，vi，d］。通過評價各粒子的適應(yīng)度值，確定每個粒子所經(jīng) 過 的最佳位置 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之間均勻分布的隨機數(shù)。

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

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

BFGS的基本思想為［8］:迭代過程中，在 x(k+1)處按下式生成對稱正定矩陣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)處的搜索方向進行搜索。BFGS對一維搜索算法精度要求較高，文中采用三次曲線擬合法［8］。

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

圖1 PSO-BFGS算法流程圖

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

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

3 設(shè)計實例

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

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

此方案類似文獻［3-4］中采用的方案。

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

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

圖2 方案2的求解收斂過程

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

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

圖5 爬升中攻角變化

圖6 爬升中推力變化

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

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

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

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

4 結(jié)束語

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

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

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

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

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

［5］ 陳勝琪.基于樣條曲線的機動彈頭再入軌跡優(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)化計算原理與算法程序設(shè)計［M］.長沙:國防科技大學(xué)出版社，2001.