基于多目标优化算法的多弹道轨迹规划及仿真 *

张耀华 ，徐亚宁 ，高少杰 ，张鹏 ，刘仁体 ，石晓龙 ，杨曦

（1.上海机电工程研究所，上海 201109；2.中国人民解放军93128 部队，北京 100854）

0 引言

随着兵器科学技术的不断发展，各军事强国空天防御体系不断完善变化，未来战争的形态将颠覆人们对常规战争的认知观感。其中，协同作战［1-4］将成为未来战争的主流，可产生1 + 1 ≫2 的效果。导弹作为现代武器的重要一员，在战争中发挥着重要作用，研究如何进行多导弹协同作战，是未来导弹武器发展的重要趋势。对于多弹协同作战，合理规划多个导弹的弹道，具有重要意义。

弹道规划是指根据既定战术指标，建立飞行动力学方程，在满足状态量约束和控制量约束的前提下，求解最优控制问题的过程，属于飞行器轨迹优化［5-7］的范畴。

文献［8］基于多枚弹道导弹的空中动态组网技术，提出了多弹协同规划作战应用方法，为多弹协同作战规划技术的应用提供了思路；文献［9］提出了一种面向突防的多导弹攻击时间/攻击角度协同的弹道规划方法，得到了满足攻击时间/攻击角度协同的弹道；文献［10］提出了一种在满足最小末速约束前提下能够实现指定攻击时间的弹道导弹协同飞行策略，实现了多弹道导弹飞行时间的协同；文献［11］提出了一种多飞行器再入段时间协同弹道规划方法，可规划出满足到达时间与终端约束的协同再入轨迹。目前针对多弹协同中弹道规划领域的研究多是基于时间协同对2 个或3 个导弹进行弹道规划，而对于规划更多数目导弹的飞行弹道的研究相对较少，本文以此为切入点展开介绍。

对于实际的弹道优化设计问题，根据不同的战术指标要求往往会有不同的优化目标，比如弹道设计要求射程尽可能大，同时希望飞行过程中的气动加热量尽量小，末端攻击速度尽量大等。上述优化目标之间可能存在较大冲突，甚至可能自相矛盾。对于这种具有多个优化目标的弹道优化问题，很难找到在所有目标中均为最优的解，求解此类多目标优化问题，要权衡各个目标，需用多目标优化算法求解。多目标优化问题的解通常被称为非支配解或 Pareto 最优解（Pareto optimal solution）［12］，其特点是无法在不削弱其中一个目标的同时去改善其他目标。对于弹道优化问题，目前常用的多目标优化算法包括带精英策略的非支配排序遗传算法（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］等，多目标优化算法在飞行器轨迹优化方面应用前景广阔。

本文着眼于规划多个导弹的飞行弹道，将MOEA/D 算法应用于助推滑翔导弹的滑翔段弹道规划，在平衡多个冲突目标（最大射程与最大末速度）的同时，获得多条优质的导弹飞行弹道。此外，考虑到初始种群对多目标优化问题的影响，利用伪谱法进行弹道优化获得的控制量作为多目标优化算法的初始种群，提高初始种群的质量。

1 MOEA/D 算法应用

1.1 多目标优化问题概述

1.1.1 MOP 的数学定义

一般，以最小化研究问题为例，多目标优化问题（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 的最优解及最优解集

一般而言，多目标优化问题很难获得一个唯一的全局最优解，而是得到多个最优解的集合，通常将MOP 的最优解称为Pareto 最优解，得到的解集称为Pareto 最优解集（Pareto set，PS）［12］。

1.2 MOEA/D 算法的基本原理

MOEA/D 是基于分解的多目标进化算法，其基本思想是为各个冲突的目标设定权重向量，借助聚合方法（分解方法）将多目标优化问题分解为单目标优化问题，进而利用进化算法求解获得的多个单目标优化问题［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 最优解也是原多目标优化问题的Pareto 最优解。换而言之，对于任意一个Pareto 最优点x*，总可以找到一个与之对应的权重向量λ*，使得x*为式（2）的Pareto 最优解。

1.2.2 基本MOEA/D 算法流程

MOEA/D 算法利用预先生成的一组均匀权重向量λ将原多目标优化问题分解为N个单目标优化子问题（以标量形式呈现），每个单目标优化子问题均视为一个个体，可形成由N个个体组成的种群，N即为MOEA/D 算法设置的种群大小。MOEA/D 采用进化算法对整个种群进行处理，即相当于同时优化N个单目标优化子问题，对种群每一代的处理过程中，挑选每个单目标优化子问题的最优解组成种群。值得注意的是，gte对权重向量λ是连续的，对于取值很接近的两个权重向量，其最优解也是相近的，单个子问题的优化主要依靠其邻居子问题的当前解进行处理。

以Tchebycheff 聚合方法进行分解，将原MOP 分解为N个单目标优化子问题，取均匀分布的N个权重向量λ1，λ2，…，λN，则第l个子问题的目标函数可表述为

式中：z*一般取全局最优点（又称理想点）；λl为第l个子问题对应的权重向量。

MOEA/D 算法主要操作步骤包括初始化操作、个体及邻域的迭代更新以及算法终止的判定等流程。工作流程如下：

算法输入：

MOP；

算法停止准则；

tgen：算法设置的最大进化代数

N：算法定义的种群规模

对应种群规模的N个权重向量：λ1，λ2，…，λN

T：每个权重向量对应的邻居权重向量个数算 法 输 出：记 为EP，由x1，x2，…，xN和FV1，FV2，…，FVN组成

（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′，若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′，FVj=F(y′)；

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

（3） 判断算法是否停止

如若算法满足停止准则，则输出EP；否则，转至2.1）继续进行迭代更新。

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 种算法在实际运行时耗费时间区别很大，以测试函数ZDT1 为例，采用2 种算法分别进行10 次重复仿真，MOEA/D 算法平均耗时约为8.23 s，NSGA-II 算 法 平 均 耗 时 约 为23.18 s，可 见MOEA/D 算法平均计算时间相比NSGA-II 算法小很多，这与文献［15］中对2 种算法每一代计算复杂度的分析是一致的。对于弹道优化问题的求解，需要进行大量计算，MOEA/D 算法计算耗时少是算法的一大优势所在。

图2 测试函数ZDT3 仿真结果Fig.2 Simulation results of test function ZDT3

2 弹道模型及气动参数

2.1 弹道模型

导弹弹道模型［16］可描述为

式中：m为导弹质量；v为导弹速度；P为发动机推力；α为攻角；β为侧滑角；θ为弹道倾角；γv为速度滚转角；ψv为弹道偏角；X为阻力；Y为升力；Z为侧向力；x为导弹射程；y为导弹飞行高度；z为导弹横程；g为重力加速度。

为了便于弹道设计，常对式（4）进行简化，得到纵向平面弹道方程为

以某型助推滑翔导弹为研究对象，上升段采用分段常值攻角控制，针对滑翔段进行弹道规划，控制变量α并不在运动模型中显含，而是隐含在空气动力当中，即

式中：Cx，Cy分别为升力系数和阻力系数；q为导弹所受空气动力与来流的动压；S为导弹参考面积（特征面积）；ρ为空气密度。

导弹总质量为2 000 kg，在滑翔段无发动机推力，质量为400 kg，参考面积S为0.204 8 m2。导弹飞行过程中，控制变量α隐含在升力系数Cx和阻力系数Cy中，控制α的取值可改导弹所受的空气动力。

2.2 大气参数模型

导弹全程在大气层内飞行，其飞行过程受大气参数变化的影响，研究导弹的飞行弹道需要获取基本的大气变化规律。

本文在进行弹道设计时直接利用美国1976 标准大气参数［17］进行仿真计算，现给出不同参数（大气压强、密度、声速值）随海拔高度的变化趋势，如图3 所示。其中，大气压强的变化与发动机推力计算相关，大气密度的变化则影响动压的计算，进而对空气动力的计算产生影响，声速值的变化则直接影响马赫数的计算。

3 基于MOEA/D 算法的多弹道轨迹仿真实现

本文研究对象为某型助推滑翔导弹，对无动力滑翔段［18］（从弹道顶点到20 km 高度这一阶段）进行仿真研究，规划多条具备不同目标值的仿真弹道。

3.1 弹道优化问题的离散化

利用MOEA/D 算法进行滑翔段弹道设计时，需要进行大量的计算，从减小计算量角度出发，将纵向平面运动模型作为弹道模型进行仿真，滑翔段无推力纵向平面运动模型可由式（5）去掉推力得到。

在借助MOEA/D 算法进行优化之前，需要将弹道优化问题离散化处理，本文采用直接打靶法［5］进行处理，将其转化为NLP（非线性规划）问题。直接打靶法对控制量的离散过程描述为

式中：飞行时间分为M个时间段；t0与tf分别为该段飞行的起始时刻和终端时刻；u(t)表示控制量。完成控制变量的离散化处理后，借助MOEA/D 算法对控制变量u(t)进行优化，即MOEA/D 算法中种群的每个个体对应控制变量u(t)。

3.2 多弹道轨迹仿真实现

3.2.1 优化设置

进行滑翔段的多目标优化设计时，以导弹射程最大和末速最大为目标函数，表述为式（9）所示，优化过程中以攻角作为控制量，对于不等式约束采用惩罚函数的形式在目标函数中体现。

式中：J1和J2为目标函数；xtf为射程；vtf为末端速度。

由于MOEA/D 算法初始化对决策变量的选取具有很大随机性，会增大优化难度，本文借助Radau 伪谱法［19］对滑翔段弹道进行优化，利用得到的优化结果作为初始种群。伪谱法优化时以射程最大为目标函数，以末速度作为约束条件，末速度约束设为500 m/s，取攻角作为控制量，范围为［-18°，+18°］，对弹道顶点（约50 km）至20 km 高度这一滑翔段弹道进行优化（弹道顶点前为上升段），仿真结果如图4～7 所示。

图4 中橙色线条为采用伪谱法优化的滑翔段弹道，弹道顶点为49.8 km，末端设定高度为20 km，优化所得最大射程为790 km；图5 给出了速度随时间变化曲线，末端速度为500 m/s；图7 给出了控制量攻角随时间变化曲线，为后续多目标优化提供算法所需初始种群。

图5 速度-时间变化曲线Fig.5 Velocity versus time

图6 弹道倾角-时间变化曲线Fig.6 Trajectory inclination angle versus time

图7 伪谱法优化攻角曲线Fig.7 Optimized angle of attack curve by pseudo-spectral method

3.2.2 多弹道轨迹仿真实现

采用MOEA/D 算法进行多目标优化时，初始种群在3.2.1 节伪谱法优化弹道获取的攻角取值基础上进行选择，振荡幅度选为±5°。仿真时弹道参数的初始条件如表1 所示，多目标优化设置子问题（种群）数目N=400，迭代次数It=20，进行弹道仿真，获得多组不同控制量的取值，进而获得多条弹道轨迹。MOEA/D 算法优化结果如图8 所示，图8 中a）为仿真所得初始Pareto 前沿，b）为筛选末速度小于900 m/s 后所得Pareto 前沿。在优化所得的Pareto 最优解中，选择末速度指标小于600 m/s 的对应解，以此作为控制量反求滑翔段飞行弹道，得到8 条飞行弹道，弹道仿真结果如图9～11 所示。

表1 初始状态参数Table 1 Parameters of initial state

图8 MOEA/D 算法优化结果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 伪谱法所得弹道对比可知，在单个目标最优性方面，MOEA/D 算法稍差于伪谱法，但MOEA/D 算法的优势在于可以同时优化多个目标（如本文中的射程与速度），并且可以得到多个Pareto 最优解，获取多种兼具不同优化目标的飞行弹道。

4 结束语

本文以射程最优和速度最优为目标函数，采用MOEA/D 算法进行滑翔段弹道规划，并借助伪谱法为多目标优化算法提供初始种群，得到了多条具备不同射程指标和速度指标的飞行弹道，可为多弹协同飞行提供参考弹道，文中给出的仿真实例说明了该方法的有效性。多目标优化算法种类很多，未来进行下一步研究还可采用其他算法（如MOPSO 算法）用于弹道规划；对于优化目标的选取，也不局限于射程、速度，还可选择横程最大、总加热量最小等目标，甚至增加至3 个目标函数进行优化。利用多目标优化算法获取多条兼具不同目标效果飞行弹道的方法，为飞行器弹道规划领域的研究提供了新思路。