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

（1.上海機電工程研究所，上海 201109；2.中國人民解放軍93128 部隊，北京 100854）

0 引言

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

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

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

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

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

1 MOEA/D 算法應用

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

1.1.1 MOP 的數學定義

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

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

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

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

1.2 MOEA/D 算法的基本原理

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

1.2.1 聚合方法

常用的聚合方法主要有加權和方法、切比雪夫（Tchebycheff）方法以及邊界交叉法等，本文主要針對Tchebycheff 方法進行介紹，Tchebycheff 聚合方法的數學描述如下所示：

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

1.2.2 基本MOEA/D 算法流程

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

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

式中：z*一般取全局最優(yōu)點（又稱理想點）；λl為第l個子問題對應的權重向量。

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

算法輸入：

MOP；

算法停止準則；

tgen：算法設置的最大進化代數

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

對應種群規(guī)模的N個權重向量：λ1，λ2，…，λN

T：每個權重向量對應的鄰居權重向量個數算 法 輸 出：記 為EP，由x1，x2，…，xN和FV1，F(xiàn)V2，…，F(xiàn)VN組成

（1） 初始化操作

1.1） 設置算法輸出EP為空集；

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

1.3） 采用隨機取樣方式獲得初始化種群x1，x2，…，xN∈Ω，計算對應xi的目標函數值F(xi)；

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

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

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

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

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

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

（3） 判斷算法是否停止

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

1.2.3 測試函數仿真實驗

為檢驗MOEA/D 算法的性能，選取文獻［7］中的兩目標測試函數ZDT1 和ZDT3 進行仿真測試，分別采用MOEA/D 算法和NSGA-II 算法進行對比仿真實驗。仿真時2 種算法均設置種群數目為200，迭代次數為300，決策變量個數取為30，交叉概率取為1，變異概率取為1/30，結果如圖1，2 所示。

圖1 測試函數ZDT1 仿真結果Fig.1 Simulation results of test function ZDT1

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

圖2 測試函數ZDT3 仿真結果Fig.2 Simulation results of test function ZDT3

2 彈道模型及氣動參數

2.1 彈道模型

導彈彈道模型［16］可描述為

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

為了便于彈道設計，常對式（4）進行簡化，得到縱向平面彈道方程為

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

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

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

2.2 大氣參數模型

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

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

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

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

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

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

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

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

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

3.2.1 優(yōu)化設置

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

式中：J1和J2為目標函數；xtf為射程；vtf為末端速度。

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

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

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

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

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

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

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

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

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

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

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

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

3.2.3 仿真結果分析

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

4 結束語

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