基于多目標(biāo)優(yōu)化算法的多彈道軌跡規(guī)劃及仿真 *

張耀華 ，徐亞寧 ，高少杰 ，張鵬 ，劉仁體 ，石曉龍 ，楊曦

（1.上海機(jī)電工程研究所，上海 201109；2.中國(guó)人民解放軍93128 部隊(duì)，北京 100854）

0 引言

隨著兵器科學(xué)技術(shù)的不斷發(fā)展，各軍事強(qiáng)國(guó)空天防御體系不斷完善變化，未來戰(zhàn)爭(zhēng)的形態(tài)將顛覆人們對(duì)常規(guī)戰(zhàn)爭(zhēng)的認(rèn)知觀感。其中，協(xié)同作戰(zhàn)［1-4］將成為未來戰(zhàn)爭(zhēng)的主流，可產(chǎn)生1 + 1 ?2 的效果。導(dǎo)彈作為現(xiàn)代武器的重要一員，在戰(zhàn)爭(zhēng)中發(fā)揮著重要作用，研究如何進(jìn)行多導(dǎo)彈協(xié)同作戰(zhàn)，是未來導(dǎo)彈武器發(fā)展的重要趨勢(shì)。對(duì)于多彈協(xié)同作戰(zhàn)，合理規(guī)劃多個(gè)導(dǎo)彈的彈道，具有重要意義。

彈道規(guī)劃是指根據(jù)既定戰(zhàn)術(shù)指標(biāo)，建立飛行動(dòng)力學(xué)方程，在滿足狀態(tài)量約束和控制量約束的前提下，求解最優(yōu)控制問題的過程，屬于飛行器軌跡優(yōu)化［5-7］的范疇。

文獻(xiàn)［8］基于多枚彈道導(dǎo)彈的空中動(dòng)態(tài)組網(wǎng)技術(shù)，提出了多彈協(xié)同規(guī)劃作戰(zhàn)應(yīng)用方法，為多彈協(xié)同作戰(zhàn)規(guī)劃技術(shù)的應(yīng)用提供了思路；文獻(xiàn)［9］提出了一種面向突防的多導(dǎo)彈攻擊時(shí)間/攻擊角度協(xié)同的彈道規(guī)劃方法，得到了滿足攻擊時(shí)間/攻擊角度協(xié)同的彈道；文獻(xiàn)［10］提出了一種在滿足最小末速約束前提下能夠?qū)崿F(xiàn)指定攻擊時(shí)間的彈道導(dǎo)彈協(xié)同飛行策略，實(shí)現(xiàn)了多彈道導(dǎo)彈飛行時(shí)間的協(xié)同；文獻(xiàn)［11］提出了一種多飛行器再入段時(shí)間協(xié)同彈道規(guī)劃方法，可規(guī)劃出滿足到達(dá)時(shí)間與終端約束的協(xié)同再入軌跡。目前針對(duì)多彈協(xié)同中彈道規(guī)劃領(lǐng)域的研究多是基于時(shí)間協(xié)同對(duì)2 個(gè)或3 個(gè)導(dǎo)彈進(jìn)行彈道規(guī)劃，而對(duì)于規(guī)劃更多數(shù)目導(dǎo)彈的飛行彈道的研究相對(duì)較少，本文以此為切入點(diǎn)展開介紹。

對(duì)于實(shí)際的彈道優(yōu)化設(shè)計(jì)問題，根據(jù)不同的戰(zhàn)術(shù)指標(biāo)要求往往會(huì)有不同的優(yōu)化目標(biāo)，比如彈道設(shè)計(jì)要求射程盡可能大，同時(shí)希望飛行過程中的氣動(dòng)加熱量盡量小，末端攻擊速度盡量大等。上述優(yōu)化目標(biāo)之間可能存在較大沖突，甚至可能自相矛盾。對(duì)于這種具有多個(gè)優(yōu)化目標(biāo)的彈道優(yōu)化問題，很難找到在所有目標(biāo)中均為最優(yōu)的解，求解此類多目標(biāo)優(yōu)化問題，要權(quán)衡各個(gè)目標(biāo)，需用多目標(biāo)優(yōu)化算法求解。多目標(biāo)優(yōu)化問題的解通常被稱為非支配解或 Pareto 最優(yōu)解（Pareto optimal solution）［12］，其特點(diǎn)是無(wú)法在不削弱其中一個(gè)目標(biāo)的同時(shí)去改善其他目標(biāo)。對(duì)于彈道優(yōu)化問題，目前常用的多目標(biāo)優(yōu)化算法包括帶精英策略的非支配排序遺傳算法（nondominated sorting genetic algorithm II ，NSGA-II）［13］、多目標(biāo)粒子群算法（multi-objective particle swarm optimization，MOPSO）［14］和 基 于 分 解 的 多 目 標(biāo) 進(jìn) 化算法（multi-objective evolutionary algorithm based on decomposition，MOEA/D）［15］等，多目標(biāo)優(yōu)化算法在飛行器軌跡優(yōu)化方面應(yīng)用前景廣闊。

本文著眼于規(guī)劃多個(gè)導(dǎo)彈的飛行彈道，將MOEA/D 算法應(yīng)用于助推滑翔導(dǎo)彈的滑翔段彈道規(guī)劃，在平衡多個(gè)沖突目標(biāo)（最大射程與最大末速度）的同時(shí)，獲得多條優(yōu)質(zhì)的導(dǎo)彈飛行彈道。此外，考慮到初始種群對(duì)多目標(biāo)優(yōu)化問題的影響，利用偽譜法進(jìn)行彈道優(yōu)化獲得的控制量作為多目標(biāo)優(yōu)化算法的初始種群，提高初始種群的質(zhì)量。

1 MOEA/D 算法應(yīng)用

1.1 多目標(biāo)優(yōu)化問題概述

1.1.1 MOP 的數(shù)學(xué)定義

一般，以最小化研究問題為例，多目標(biāo)優(yōu)化問題（multi-objective optimization problem，MOP）的數(shù)學(xué)描述可以表示為如下形式：

式中：f1(x)，f2(x)，…，fn(x)為n個(gè)目標(biāo)函數(shù)，它們共同構(gòu)成目標(biāo)函數(shù)向量f(x)；gi(x) ≤0 和hi(x) = 0 分別表示MOP 的k個(gè)不等式約束和l個(gè)等式約束；x=(x1，x2，…，xm)表示決策向量；若存在f：Ω→Φ，則Ω∈Rm，Φ∈Rn，稱Ω和Φ分別為MOP 的決策空間和目標(biāo)空間。

1.1.2 MOP 的最優(yōu)解及最優(yōu)解集

一般而言，多目標(biāo)優(yōu)化問題很難獲得一個(gè)唯一的全局最優(yōu)解，而是得到多個(gè)最優(yōu)解的集合，通常將MOP 的最優(yōu)解稱為Pareto 最優(yōu)解，得到的解集稱為Pareto 最優(yōu)解集（Pareto set，PS）［12］。

1.2 MOEA/D 算法的基本原理

MOEA/D 是基于分解的多目標(biāo)進(jìn)化算法，其基本思想是為各個(gè)沖突的目標(biāo)設(shè)定權(quán)重向量，借助聚合方法（分解方法）將多目標(biāo)優(yōu)化問題分解為單目標(biāo)優(yōu)化問題，進(jìn)而利用進(jìn)化算法求解獲得的多個(gè)單目標(biāo)優(yōu)化問題［15］。

1.2.1 聚合方法

常用的聚合方法主要有加權(quán)和方法、切比雪夫（Tchebycheff）方法以及邊界交叉法等，本文主要針對(duì)Tchebycheff 方法進(jìn)行介紹，Tchebycheff 聚合方法的數(shù)學(xué)描述如下所示：

式中：m為MOP 的目標(biāo)函數(shù)個(gè)數(shù)；λ為權(quán)重向量；z*=(z1，z2，…，zm) 被 稱 為 參 考 點(diǎn) ，滿 足z*i=min {fi(x)|x∈Ω}，i= 1，2，…，m。Tchebycheff 聚合方法的相關(guān)理論研究表明，式（2）的Pareto 最優(yōu)解也是原多目標(biāo)優(yōu)化問題的Pareto 最優(yōu)解。換而言之，對(duì)于任意一個(gè)Pareto 最優(yōu)點(diǎn)x*，總可以找到一個(gè)與之對(duì)應(yīng)的權(quán)重向量λ*，使得x*為式（2）的Pareto 最優(yōu)解。

1.2.2 基本MOEA/D 算法流程

MOEA/D 算法利用預(yù)先生成的一組均勻權(quán)重向量λ將原多目標(biāo)優(yōu)化問題分解為N個(gè)單目標(biāo)優(yōu)化子問題（以標(biāo)量形式呈現(xiàn)），每個(gè)單目標(biāo)優(yōu)化子問題均視為一個(gè)個(gè)體，可形成由N個(gè)個(gè)體組成的種群，N即為MOEA/D 算法設(shè)置的種群大小。MOEA/D 采用進(jìn)化算法對(duì)整個(gè)種群進(jìn)行處理，即相當(dāng)于同時(shí)優(yōu)化N個(gè)單目標(biāo)優(yōu)化子問題，對(duì)種群每一代的處理過程中，挑選每個(gè)單目標(biāo)優(yōu)化子問題的最優(yōu)解組成種群。值得注意的是，gte對(duì)權(quán)重向量λ是連續(xù)的，對(duì)于取值很接近的兩個(gè)權(quán)重向量，其最優(yōu)解也是相近的，單個(gè)子問題的優(yōu)化主要依靠其鄰居子問題的當(dāng)前解進(jìn)行處理。

以Tchebycheff 聚合方法進(jìn)行分解，將原MOP 分解為N個(gè)單目標(biāo)優(yōu)化子問題，取均勻分布的N個(gè)權(quán)重向量λ1，λ2，…，λN，則第l個(gè)子問題的目標(biāo)函數(shù)可表述為

式中：z*一般取全局最優(yōu)點(diǎn)（又稱理想點(diǎn)）；λl為第l個(gè)子問題對(duì)應(yīng)的權(quán)重向量。

MOEA/D 算法主要操作步驟包括初始化操作、個(gè)體及鄰域的迭代更新以及算法終止的判定等流程。工作流程如下：

算法輸入：

MOP；

算法停止準(zhǔn)則；

tgen：算法設(shè)置的最大進(jìn)化代數(shù)

N：算法定義的種群規(guī)模

對(duì)應(yīng)種群規(guī)模的N個(gè)權(quán)重向量：λ1，λ2，…，λN

T：每個(gè)權(quán)重向量對(duì)應(yīng)的鄰居權(quán)重向量個(gè)數(shù)算 法 輸 出：記 為EP，由x1，x2，…，xN和FV1，F(xiàn)V2，…，F(xiàn)VN組成

（1） 初始化操作

1.1） 設(shè)置算法輸出EP為空集；

1.2） 計(jì)算任意兩個(gè)權(quán)重向量λi和λj之間的歐氏距離，找出權(quán)重向量λi的T個(gè)鄰居權(quán)重向量記為{λi1，λi2，…，λiT}，將鄰居權(quán)重向量的索引值記錄為B(i) ={i1，i2，…，iT}(i= 1，2，…，N)；

1.3） 采用隨機(jī)取樣方式獲得初始化種群x1，x2，…，xN∈Ω，計(jì)算對(duì)應(yīng)xi的目標(biāo)函數(shù)值F(xi)；

1.4） 初始化z=(z1，z2，…，zm)T，取zi= min {fi(x)|x∈Ω}，i= 1，2，…，m。

（2） 迭代計(jì)算并更新對(duì)于i= 1，2，…，N，進(jìn)行如下操作：

2.1） 重新組合：從記錄Bi中隨機(jī)挑選2 個(gè)索引k和l，與 之 對(duì) 應(yīng) 的2 個(gè) 個(gè) 體 為xk和xl，利 用 遺 傳 算 子由二者產(chǎn)生新個(gè)體y，對(duì)y進(jìn)行修復(fù)和改進(jìn)得到y(tǒng)′，若y′超出Ω范圍，則用邊界內(nèi)的隨機(jī)數(shù)代替；

2.2） 更新z：對(duì)于j= 1，2，…，m，計(jì)算得到的新個(gè)體y′對(duì)應(yīng)的函數(shù)值fj(y′)，若zj＞fj(y′)，則需更新z，令zj=fj(y′)；

2.3） 更 新 鄰 域 解：對(duì) 于 每 個(gè)j∈B(i)，如 果gte(y′|λj，z) ≤gte(xj|λj，z)，則說明新個(gè)體y′相 比xj更好，應(yīng)進(jìn)行更新，令xj=y′，F(xiàn)Vj=F(y′)；

2.4） 更 新 算 法 輸 出EP：對(duì) 于 輸 出EP中 受F(y′)支配的，將其移除出去，若F(y′)不受EP中向量的支配，將F(y′)存入EP。

（3） 判斷算法是否停止

如若算法滿足停止準(zhǔn)則，則輸出EP；否則，轉(zhuǎn)至2.1）繼續(xù)進(jìn)行迭代更新。

1.2.3 測(cè)試函數(shù)仿真實(shí)驗(yàn)

為檢驗(yàn)MOEA/D 算法的性能，選取文獻(xiàn)［7］中的兩目標(biāo)測(cè)試函數(shù)ZDT1 和ZDT3 進(jìn)行仿真測(cè)試，分別采用MOEA/D 算法和NSGA-II 算法進(jìn)行對(duì)比仿真實(shí)驗(yàn)。仿真時(shí)2 種算法均設(shè)置種群數(shù)目為200，迭代次數(shù)為300，決策變量個(gè)數(shù)取為30，交叉概率取為1，變異概率取為1/30，結(jié)果如圖1，2 所示。

圖1 測(cè)試函數(shù)ZDT1 仿真結(jié)果Fig.1 Simulation results of test function ZDT1

分析圖1 和圖2 可知，在上述算法設(shè)置條件下，MOEA/D 算法和NSGA-II 算法都可獲得較好的仿真結(jié)果，與測(cè)試函數(shù)ZDT1 和ZDT3 的理想Pareto 前沿很接近。然而，2 種算法在實(shí)際運(yùn)行時(shí)耗費(fèi)時(shí)間區(qū)別很大，以測(cè)試函數(shù)ZDT1 為例，采用2 種算法分別進(jìn)行10 次重復(fù)仿真，MOEA/D 算法平均耗時(shí)約為8.23 s，NSGA-II 算 法 平 均 耗 時(shí) 約 為23.18 s，可 見MOEA/D 算法平均計(jì)算時(shí)間相比NSGA-II 算法小很多，這與文獻(xiàn)［15］中對(duì)2 種算法每一代計(jì)算復(fù)雜度的分析是一致的。對(duì)于彈道優(yōu)化問題的求解，需要進(jìn)行大量計(jì)算，MOEA/D 算法計(jì)算耗時(shí)少是算法的一大優(yōu)勢(shì)所在。

圖2 測(cè)試函數(shù)ZDT3 仿真結(jié)果Fig.2 Simulation results of test function ZDT3

2 彈道模型及氣動(dòng)參數(shù)

2.1 彈道模型

導(dǎo)彈彈道模型［16］可描述為

式中：m為導(dǎo)彈質(zhì)量；v為導(dǎo)彈速度；P為發(fā)動(dòng)機(jī)推力；α為攻角；β為側(cè)滑角；θ為彈道傾角；γv為速度滾轉(zhuǎn)角；ψv為彈道偏角；X為阻力；Y為升力；Z為側(cè)向力；x為導(dǎo)彈射程；y為導(dǎo)彈飛行高度；z為導(dǎo)彈橫程；g為重力加速度。

為了便于彈道設(shè)計(jì)，常對(duì)式（4）進(jìn)行簡(jiǎn)化，得到縱向平面彈道方程為

以某型助推滑翔導(dǎo)彈為研究對(duì)象，上升段采用分段常值攻角控制，針對(duì)滑翔段進(jìn)行彈道規(guī)劃，控制變量α并不在運(yùn)動(dòng)模型中顯含，而是隱含在空氣動(dòng)力當(dāng)中，即

式中：Cx，Cy分別為升力系數(shù)和阻力系數(shù)；q為導(dǎo)彈所受空氣動(dòng)力與來流的動(dòng)壓；S為導(dǎo)彈參考面積（特征面積）；ρ為空氣密度。

導(dǎo)彈總質(zhì)量為2 000 kg，在滑翔段無(wú)發(fā)動(dòng)機(jī)推力，質(zhì)量為400 kg，參考面積S為0.204 8 m2。導(dǎo)彈飛行過程中，控制變量α隱含在升力系數(shù)Cx和阻力系數(shù)Cy中，控制α的取值可改導(dǎo)彈所受的空氣動(dòng)力。

2.2 大氣參數(shù)模型

導(dǎo)彈全程在大氣層內(nèi)飛行，其飛行過程受大氣參數(shù)變化的影響，研究導(dǎo)彈的飛行彈道需要獲取基本的大氣變化規(guī)律。

本文在進(jìn)行彈道設(shè)計(jì)時(shí)直接利用美國(guó)1976 標(biāo)準(zhǔn)大氣參數(shù)［17］進(jìn)行仿真計(jì)算，現(xiàn)給出不同參數(shù)（大氣壓強(qiáng)、密度、聲速值）隨海拔高度的變化趨勢(shì)，如圖3 所示。其中，大氣壓強(qiáng)的變化與發(fā)動(dòng)機(jī)推力計(jì)算相關(guān)，大氣密度的變化則影響動(dòng)壓的計(jì)算，進(jìn)而對(duì)空氣動(dòng)力的計(jì)算產(chǎn)生影響，聲速值的變化則直接影響馬赫數(shù)的計(jì)算。

3 基于MOEA/D 算法的多彈道軌跡仿真實(shí)現(xiàn)

本文研究對(duì)象為某型助推滑翔導(dǎo)彈，對(duì)無(wú)動(dòng)力滑翔段［18］（從彈道頂點(diǎn)到20 km 高度這一階段）進(jìn)行仿真研究，規(guī)劃多條具備不同目標(biāo)值的仿真彈道。

3.1 彈道優(yōu)化問題的離散化

利用MOEA/D 算法進(jìn)行滑翔段彈道設(shè)計(jì)時(shí)，需要進(jìn)行大量的計(jì)算，從減小計(jì)算量角度出發(fā)，將縱向平面運(yùn)動(dòng)模型作為彈道模型進(jìn)行仿真，滑翔段無(wú)推力縱向平面運(yùn)動(dòng)模型可由式（5）去掉推力得到。

在借助MOEA/D 算法進(jìn)行優(yōu)化之前，需要將彈道優(yōu)化問題離散化處理，本文采用直接打靶法［5］進(jìn)行處理，將其轉(zhuǎn)化為NLP（非線性規(guī)劃）問題。直接打靶法對(duì)控制量的離散過程描述為

式中：飛行時(shí)間分為M個(gè)時(shí)間段；t0與tf分別為該段飛行的起始時(shí)刻和終端時(shí)刻；u(t)表示控制量。完成控制變量的離散化處理后，借助MOEA/D 算法對(duì)控制變量u(t)進(jìn)行優(yōu)化，即MOEA/D 算法中種群的每個(gè)個(gè)體對(duì)應(yīng)控制變量u(t)。

3.2 多彈道軌跡仿真實(shí)現(xiàn)

3.2.1 優(yōu)化設(shè)置

進(jìn)行滑翔段的多目標(biāo)優(yōu)化設(shè)計(jì)時(shí)，以導(dǎo)彈射程最大和末速最大為目標(biāo)函數(shù)，表述為式（9）所示，優(yōu)化過程中以攻角作為控制量，對(duì)于不等式約束采用懲罰函數(shù)的形式在目標(biāo)函數(shù)中體現(xiàn)。

式中：J1和J2為目標(biāo)函數(shù)；xtf為射程；vtf為末端速度。

由于MOEA/D 算法初始化對(duì)決策變量的選取具有很大隨機(jī)性，會(huì)增大優(yōu)化難度，本文借助Radau 偽譜法［19］對(duì)滑翔段彈道進(jìn)行優(yōu)化，利用得到的優(yōu)化結(jié)果作為初始種群。偽譜法優(yōu)化時(shí)以射程最大為目標(biāo)函數(shù)，以末速度作為約束條件，末速度約束設(shè)為500 m/s，取攻角作為控制量，范圍為［-18°，+18°］，對(duì)彈道頂點(diǎn)（約50 km）至20 km 高度這一滑翔段彈道進(jìn)行優(yōu)化（彈道頂點(diǎn)前為上升段），仿真結(jié)果如圖4～7 所示。

圖4 中橙色線條為采用偽譜法優(yōu)化的滑翔段彈道，彈道頂點(diǎn)為49.8 km，末端設(shè)定高度為20 km，優(yōu)化所得最大射程為790 km；圖5 給出了速度隨時(shí)間變化曲線，末端速度為500 m/s；圖7 給出了控制量攻角隨時(shí)間變化曲線，為后續(xù)多目標(biāo)優(yōu)化提供算法所需初始種群。

圖5 速度-時(shí)間變化曲線Fig.5 Velocity versus time

圖6 彈道傾角-時(shí)間變化曲線Fig.6 Trajectory inclination angle versus time

圖7 偽譜法優(yōu)化攻角曲線Fig.7 Optimized angle of attack curve by pseudo-spectral method

3.2.2 多彈道軌跡仿真實(shí)現(xiàn)

采用MOEA/D 算法進(jìn)行多目標(biāo)優(yōu)化時(shí)，初始種群在3.2.1 節(jié)偽譜法優(yōu)化彈道獲取的攻角取值基礎(chǔ)上進(jìn)行選擇，振蕩幅度選為±5°。仿真時(shí)彈道參數(shù)的初始條件如表1 所示，多目標(biāo)優(yōu)化設(shè)置子問題（種群）數(shù)目N=400，迭代次數(shù)It=20，進(jìn)行彈道仿真，獲得多組不同控制量的取值，進(jìn)而獲得多條彈道軌跡。MOEA/D 算法優(yōu)化結(jié)果如圖8 所示，圖8 中a）為仿真所得初始Pareto 前沿，b）為篩選末速度小于900 m/s 后所得Pareto 前沿。在優(yōu)化所得的Pareto 最優(yōu)解中，選擇末速度指標(biāo)小于600 m/s 的對(duì)應(yīng)解，以此作為控制量反求滑翔段飛行彈道，得到8 條飛行彈道，彈道仿真結(jié)果如圖9～11 所示。

表1 初始狀態(tài)參數(shù)Table 1 Parameters of initial state

圖8 MOEA/D 算法優(yōu)化結(jié)果Fig.8 Optimization results of MOEA/D algorithm

圖9 高度-射程曲線Fig.9 Altitude versus range

圖10 速度-時(shí)間變化曲線Fig.10 Velocity versus time

圖11 彈道傾角-時(shí)間變化曲線Fig.11 Trajectory inclination angle versus time

3.2.3 仿真結(jié)果分析

在仿真得到的8 條飛行彈道中，射程最大為761.4 km，射程最小為740.3 km，末速最大為566.4 m/s，末速最小為484.6 m/s。與Radau 偽譜法所得彈道對(duì)比可知，在單個(gè)目標(biāo)最優(yōu)性方面，MOEA/D 算法稍差于偽譜法，但MOEA/D 算法的優(yōu)勢(shì)在于可以同時(shí)優(yōu)化多個(gè)目標(biāo)（如本文中的射程與速度），并且可以得到多個(gè)Pareto 最優(yōu)解，獲取多種兼具不同優(yōu)化目標(biāo)的飛行彈道。

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

本文以射程最優(yōu)和速度最優(yōu)為目標(biāo)函數(shù)，采用MOEA/D 算法進(jìn)行滑翔段彈道規(guī)劃，并借助偽譜法為多目標(biāo)優(yōu)化算法提供初始種群，得到了多條具備不同射程指標(biāo)和速度指標(biāo)的飛行彈道，可為多彈協(xié)同飛行提供參考彈道，文中給出的仿真實(shí)例說明了該方法的有效性。多目標(biāo)優(yōu)化算法種類很多，未來進(jìn)行下一步研究還可采用其他算法（如MOPSO 算法）用于彈道規(guī)劃；對(duì)于優(yōu)化目標(biāo)的選取，也不局限于射程、速度，還可選擇橫程最大、總加熱量最小等目標(biāo)，甚至增加至3 個(gè)目標(biāo)函數(shù)進(jìn)行優(yōu)化。利用多目標(biāo)優(yōu)化算法獲取多條兼具不同目標(biāo)效果飛行彈道的方法，為飛行器彈道規(guī)劃領(lǐng)域的研究提供了新思路。