楊光宇，陳國(guó)光，王 捷，朱豪坤

(中北大學(xué) 機(jī)電工程學(xué)院， 太原 030051)

空地導(dǎo)彈是軍事戰(zhàn)爭(zhēng)中使用范圍廣、消耗量大、用途多的武器類(lèi)型[1]。隨著作戰(zhàn)任務(wù)及打擊目標(biāo)的多樣化，對(duì)空地導(dǎo)彈彈道軌跡也提出了更高的要求。彈道優(yōu)化設(shè)計(jì)是現(xiàn)代導(dǎo)彈設(shè)計(jì)中的重要內(nèi)容，有利于提高導(dǎo)彈飛行品質(zhì)以滿(mǎn)足既定任務(wù)要求。對(duì)導(dǎo)彈設(shè)計(jì)有著十分重要的意義及實(shí)際工程價(jià)值[2-3]。

彈道優(yōu)化設(shè)計(jì)，實(shí)質(zhì)就是求解在滿(mǎn)足各種約束條件下的最優(yōu)控制問(wèn)題，同時(shí)也是一個(gè)動(dòng)態(tài)優(yōu)化問(wèn)題。彈道優(yōu)化的方法主要分為直接法和間接法兩大類(lèi)[4-5]。直接法是將連續(xù)的優(yōu)化問(wèn)題直接離散并進(jìn)行參數(shù)化，不需要求解出最優(yōu)值的條件。間接法是基于極大值原理將最優(yōu)控制問(wèn)題轉(zhuǎn)換為哈密頓邊值問(wèn)題[6]。偽譜法是數(shù)值法的典型代表方法，大多數(shù)配點(diǎn)是基于高斯積分規(guī)則[7]。高斯偽譜法作為一種基于全局插值多項(xiàng)式的直接配點(diǎn)法，它相對(duì)于一般直接配點(diǎn)法的優(yōu)勢(shì)是可以用較少的節(jié)點(diǎn)獲得較高的精度[8]。

本文利用高斯偽譜法將最優(yōu)控制問(wèn)題(動(dòng)態(tài)優(yōu)化問(wèn)題)轉(zhuǎn)化為非線(xiàn)性規(guī)劃問(wèn)題，SQP算法因其穩(wěn)定性和精度優(yōu)勢(shì)在求解此類(lèi)問(wèn)題上應(yīng)用最為廣泛。但SQP算法對(duì)初值依賴(lài)度很高，如若隨意選取初值，很難求得最優(yōu)解[9]。粒子群算法對(duì)目標(biāo)函數(shù)和應(yīng)用條件的限制較少，而且收斂速度快，具有更大的靈活性和適應(yīng)性[10]，因此本文利用全局尋優(yōu)能力較強(qiáng)的粒子群算法快速求解一條優(yōu)化彈道作為初值，而后代入SQP算法求解該問(wèn)題。

1 問(wèn)題描述

空地導(dǎo)彈彈道優(yōu)化設(shè)計(jì)問(wèn)題是一類(lèi)復(fù)雜的多約束動(dòng)態(tài)優(yōu)化問(wèn)題，問(wèn)題描述包括目標(biāo)函數(shù)、約束條件和優(yōu)化變量等。一般最優(yōu)控制問(wèn)題的Bolza形式描述如下式：

(1)

其中，u(t)∈Rm為控制變量；x(t)∈Rn為狀態(tài)變量；t0和tf為時(shí)間的初值和末值，可固定也可自由。

系統(tǒng)滿(mǎn)足的動(dòng)力學(xué)約束、邊界條件約束和路徑約束分別為

(2)

期望得到一條用時(shí)最短且滿(mǎn)足約束的彈道，故目標(biāo)函數(shù)選取終端時(shí)刻，即

minJ=tf

(3)

方案彈道初步設(shè)計(jì)過(guò)程中簡(jiǎn)化了運(yùn)動(dòng)模型，將導(dǎo)彈描述為理想的質(zhì)點(diǎn)，暫時(shí)不考慮彈體繞質(zhì)心轉(zhuǎn)動(dòng)的影響，以縱向平面內(nèi)的質(zhì)心運(yùn)動(dòng)方程為動(dòng)力學(xué)約束[11]。

動(dòng)力學(xué)約束方程如式(4)所示，式中:m為導(dǎo)彈質(zhì)量；P為發(fā)動(dòng)機(jī)推力；ISP為發(fā)動(dòng)機(jī)比沖；mc為質(zhì)量流量；X=qsCx、Y=qsCy分別為阻力和升力，q為動(dòng)壓，s為參考面積，Cx、Cy分別為阻力和升力系數(shù)，兩者均是馬赫數(shù)Ma和攻角α的函數(shù)，可由吹風(fēng)試驗(yàn)和數(shù)值計(jì)算獲得；?為導(dǎo)彈俯仰角；θ為彈道傾角；x、y分別為射程和飛行高度。

(4)

邊界約束條件主要包括終端時(shí)刻的射程x和飛行高度y。

(5)

路徑約束條件主要包括飛行過(guò)程中的動(dòng)壓q、法向過(guò)載ny和控制變量(攻角)α的約束。

(6)

2 高斯偽譜法

高斯偽譜法以L(fǎng)egendre多項(xiàng)式的根為離散點(diǎn)，將連續(xù)最優(yōu)控制問(wèn)題的狀態(tài)變量和控制變量離散化，并以離散點(diǎn)為節(jié)點(diǎn)采用全區(qū)間Lagrange插值多項(xiàng)式來(lái)近似狀態(tài)變量和控制變量，從而將動(dòng)力學(xué)方程等約束轉(zhuǎn)化為代數(shù)約束，形成非線(xiàn)性規(guī)劃問(wèn)題[12-13]。

高斯偽譜法的基本算法過(guò)程如下：

1) 對(duì)時(shí)間變量t作變換到τ∈[-1，1]上。

(7)

動(dòng)力學(xué)方程約束式(2)轉(zhuǎn)化為

(8)

2) 取N階Legendre多項(xiàng)式N的個(gè)根τk∈(-1，1)，k=1，2，…，N；再令τ0=-1，τf=1，則τ0，τk，τf構(gòu)成彈道的N+2個(gè)離散點(diǎn)。

3) 將狀態(tài)變量與控制變量在離散點(diǎn)上離散化，從而將動(dòng)力學(xué)方程約束、邊界條件約束和路徑約束轉(zhuǎn)化為代數(shù)約束。

由高斯偽譜法的算法過(guò)程可得：

(9)

對(duì)式(9)求導(dǎo)：

(10)

將式(10)代入式(8)，得到狀態(tài)向量在配點(diǎn)上應(yīng)滿(mǎn)足的代數(shù)約束：

(11)

由式(8)可得：

(12)

采用高斯積分公式近似式(12)中的積分項(xiàng)，可得離散化終端邊界條件：

(13)

把邊界條件和過(guò)程約束也離散化，將彈道優(yōu)化問(wèn)題轉(zhuǎn)化為非線(xiàn)性規(guī)劃問(wèn)題。其一般形式描述如下式：

(14)

其中，c為非線(xiàn)性不等式約束;ceq為非線(xiàn)性等式約束。

3 粒子群算法

粒子群算法是計(jì)算智能領(lǐng)域，除了蟻群算法、魚(yú)群算法之外的一種群體智能的優(yōu)化算法。該算法最早由Kennedy和Eberhart在1995年提出的。粒子群算法源于對(duì)鳥(niǎo)類(lèi)捕食行為的研究，鳥(niǎo)類(lèi)捕食時(shí)，找到食物最簡(jiǎn)單有效的策略就是搜尋當(dāng)前距離食物最近的鳥(niǎo)的周?chē)鷧^(qū)域。

粒子群尋優(yōu)算法中每個(gè)粒子追隨個(gè)體極值和全局極值的過(guò)程進(jìn)行信息更新[14]，同時(shí)也代表優(yōu)化問(wèn)題的一個(gè)潛在解，用位置、速度和適應(yīng)度值三項(xiàng)指標(biāo)表示該粒子特征。

粒子的適應(yīng)度值由適應(yīng)度函數(shù)決定，其值表示了粒子的優(yōu)劣。粒子的速度決定了粒子移動(dòng)的方向和距離，速度隨自身及其他粒子的移動(dòng)經(jīng)驗(yàn)進(jìn)行動(dòng)態(tài)調(diào)整。粒子的位置隨著自身位置及速度進(jìn)行動(dòng)態(tài)調(diào)整。從而實(shí)現(xiàn)個(gè)體在可行解空間中的尋優(yōu)。

粒子的速度和位置更新公式分別為[15]：

(15)

(16)

設(shè)向量X={α1，α2，…，αn，tf}為優(yōu)化問(wèn)題的一個(gè)潛在解，即為粒子群算法中的一個(gè)粒子。其中前n維為攻角的變化值，最后一維tf為終端飛行時(shí)間(粒子的維數(shù)為n+1)，換言之每個(gè)粒子都代表了攻角隨時(shí)間變化的一種規(guī)律。

初始化粒子群，設(shè)定種群規(guī)模(Sizepop)、最大進(jìn)化代數(shù)(Maxgen)、慣性權(quán)重、加速度常數(shù)、位置和速度的上下限等。通過(guò)適度函數(shù)計(jì)算出n維的所有狀態(tài)變量值，進(jìn)而求出每個(gè)粒子的適度值，保存?zhèn)€體和群體經(jīng)歷的最優(yōu)位置，更新粒子的速度和位置，重復(fù)上述過(guò)程，直到滿(mǎn)足終止條件。將得到的最優(yōu)值擬合成一條連續(xù)彈道，然后在得到的彈道曲線(xiàn)上用高斯偽譜法的取點(diǎn)規(guī)則取點(diǎn)，得到SQP算法的初值X0，代入求解非線(xiàn)性規(guī)劃問(wèn)題，得到最終優(yōu)化設(shè)計(jì)彈道，算法流程如圖1所示。

圖1 算法流程框圖

4 仿真結(jié)果與分析

結(jié)合某型空地導(dǎo)彈彈道仿真結(jié)果和算法流程圖，粒子群算法各個(gè)參數(shù)的初值如下：Maxgen=50，Sizepop=20，c1=c2=1.45，ω=0.8，t0=0?？盏貜椀婪桨竷?yōu)化設(shè)計(jì)條件為：x(τ0)=0，y(τ0)=7 000，m(τ0)=225，x(τf)=14 000，y(τf)=0，nymax=15，αmax=0.2 rad。適應(yīng)度函數(shù)的功能是通過(guò)粒子X(jué)={α1，α2，…，αn，tf}計(jì)算出狀態(tài)變量的變化規(guī)律，最后適度值的大小表明了初始彈道最終與目標(biāo)點(diǎn)距離的遠(yuǎn)近。

粒子的維數(shù)為n+1，n的值可任意選取，但不宜過(guò)大或過(guò)小。

以n取為20、30和40為例，算法進(jìn)化曲線(xiàn)如圖2所示，仿真結(jié)果如圖3和圖4所示。

圖2 粒子群算法進(jìn)化曲線(xiàn)

圖3 n=20、n=30、n=40速度-時(shí)間仿真曲線(xiàn)

圖4 n=20、n=30、n=40高度-射程仿真曲線(xiàn)

圖2表明：通過(guò)粒子群算法可以較快的算出一條初始彈道，計(jì)算用時(shí)與n的大小有關(guān)。計(jì)算用時(shí)可以通過(guò)程序內(nèi)的相應(yīng)命令得到。當(dāng)n分別取20、30和40時(shí)，平均用時(shí)為7.76 s、11.42 s和15.67 s。

圖3和圖4表明：初始彈道的可行性、精度和計(jì)算量也與n的大小有關(guān)。n的值太小會(huì)導(dǎo)致彈道不可行、精度較低，太大又會(huì)導(dǎo)致計(jì)算量增加。通過(guò)多次試驗(yàn)發(fā)現(xiàn)，當(dāng)n≤20時(shí)，彈道的可行性和精度都很差。n≥30發(fā)現(xiàn)再增加n的值所計(jì)算出的結(jié)果差別不大，而且繼續(xù)增加n的值，彈道的計(jì)算量會(huì)顯著的增加。這表明在滿(mǎn)足方案設(shè)計(jì)要求下，選取n=30為宜，可以做到初始彈道的可行性、精度和計(jì)算量三者兼顧。此外，由于粒子群算法的隨機(jī)性，每次解得的初始彈道也不完全相同。

以攻角為控制變量，選取n=30，得到初始彈道，利用高斯偽譜法的取點(diǎn)規(guī)則取點(diǎn)，得到SQP算法的初值，代入求解得到最終優(yōu)化設(shè)計(jì)彈道。通過(guò)反復(fù)計(jì)算發(fā)現(xiàn)，取點(diǎn)數(shù)N取10到20之間為宜。為了更好的在圖上表現(xiàn)出高斯偽譜法的取點(diǎn)規(guī)則，這里選取N=12(加上τ0，τf，構(gòu)成圖中彈道的14個(gè)離散點(diǎn))。同等條件下，隨意選取SQP算法的初值(本文直接選用了一個(gè)初始化粒子)，仿真結(jié)果如圖5、圖6、圖7所示。

圖5 攻角-時(shí)間仿真曲線(xiàn)

圖6 速度-時(shí)間仿真曲線(xiàn)

圖7 高度-射程仿真曲線(xiàn)

圖5、圖6、圖7表明：在滿(mǎn)足約束條件的前提下，使用PSO-SQP算法完全優(yōu)于僅使用SQP算法計(jì)算出的彈道。此外通過(guò)比較圖3(n=30的速度曲線(xiàn))和圖6(PSO-SQP的速度曲線(xiàn))，可以發(fā)現(xiàn)最終彈道相較于初始彈道，飛行時(shí)間從32.20 s縮短到了31.63 s。

仿真結(jié)果表明：僅用SQP算法求解彈道優(yōu)化設(shè)計(jì)問(wèn)題，由于選取問(wèn)題，很難求得最優(yōu)解，利用粒子群算法可以很好地解決SQP算法的初值依賴(lài)的問(wèn)題，增加了得到最優(yōu)解的可能性。并且計(jì)算出的最終設(shè)計(jì)彈道相較于初始彈道，也更好地符合了方案設(shè)計(jì)的預(yù)期效果。

5 結(jié)論

針對(duì)空地導(dǎo)彈快速打擊地面目標(biāo)的彈道優(yōu)化設(shè)計(jì)問(wèn)題，提出了一種基于PSO-SQP算法的優(yōu)化方法，通過(guò)研究分析和仿真結(jié)果得出：利用高斯偽譜法離散約束條件，將最優(yōu)控制問(wèn)題轉(zhuǎn)化為非線(xiàn)性規(guī)劃問(wèn)題，可以用較少的節(jié)點(diǎn)獲得較高的精度。利用粒子群算法可快速、有效的選取初始彈道。本文提出的算法流程可以更快更有效地解決彈道優(yōu)化設(shè)計(jì)問(wèn)題，可為工程設(shè)計(jì)提供參考。