一种基于系统状态约束的低能耗弹道优化方法与流程

1.本发明涉及一种基于系统状态约束的低能耗弹道优化方法，属于制导火箭弹弹道优化技术领域。


背景技术：

2.制导火箭弹在普通无控火箭弹的基础上增加探测制导装置，具有更高的命中精度和较低的成本，可有效填补无控火箭弹和导弹之间的空缺。如图1所示，制导火箭弹的飞行分为无控段、中制导段、末制导段三部分，弹体在中制导段按照方案弹道进行方案飞行，末制导按照制导律飞行实现对目标的追击。因此，末制导是实现精准打击的关键，为保证其顺利进行，中制导和末制导交接的状态是至关重要的。
3.针对目标位置区域，设定终端海拔和速度约束并选取耗能小为目标适应度函数，保证末制导阶段有充足的能量对目标进行追击，实现精准制导；设定终端弹道倾角约束有利于末制导阶段角度的变化以及打击时终端攻击角约束的实现；设定过程约束保证弹体不会失稳，顺利进入末制导阶段。通过优化算法对弹道进行优化，优化后的弹道满足终端约束、过程约束以及优化目标，为顺利进入末制导并实现精准打击打下基础。


技术实现要素：

4.为保证制导火箭弹顺利进入末制导阶段，并有良好的末制导初始状态，本发明提供一种基于系统状态约束的低能耗弹道优化方法，可降低能耗并满足过程约束、终端约束，为末制导提供良好的初始状态，可提高制导精度并保证终端攻击角约束的顺利进行。本发明的步骤包括建立弹体的运动学方程组、采用龙格
‑
库塔法求解微分方程组、设定约束条件及目标适应度函数、通过自适应遗传算法对弹道进行优化，具体如下：
5.步骤一、建立弹道运动学方程组
6.优选地，方案飞行一般选取一个固定铅垂平面飞行，因此在纵平面建立弹道运动学方程组，经过推导化简得到如下方程组：
[0007][0008]
式(1)中式中m为弹体质量、v为弹体速度、α为攻角、θ为弹道倾角、g为重力加速度、t为飞行时刻；x、y分别为弹体的水平位置和海拔高度；c
x0
为零升阻力、c
x1
为诱导阻力；为升力系数导数、s为特征面积、q为动压q＝ρv2，ρ为大气密度。
[0009]
步骤二、求解微分方程组
[0010]
在本发明中选取攻角α为控制变量，对攻角离散化处理，将弹道优化问题转化为非线性规划问题，对时间进行r等分，得到u＝[α1,α2,...α
r
,α
r 1
]作为待优化的控制变量。为保证控制的连续性，采用五次多项式插值法对控制量平滑处理。通过α(t)以及初始状态采用四阶龙格
‑
库塔法求解微分方程组，得到弹道的位置、速度、弹道倾角等参数信息。
[0011]
步骤三、设定约束条件及目标适应度函数
[0012]
选取约束如下：
[0013]
(1)攻角控制容许约束：
[0014]
为防止弹体失稳，不能以过大攻角飞行，选取控制容许范围约束为：
[0015]
|α|≤α
max
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0016]
式(2)中α
max
为允许控制攻角绝对值的最大值。
[0017]
(2)终端速度约束：
[0018]
为保证制导火箭弹在末制导阶段对目标顺利追击，弹体速度不能过低，设定约束为：
[0019]
v(t
f
)≥v
min
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0020]
式(3)中v(t
f
)为终端时刻的速度值，v
min
为终端速度允许的最小值。
[0021]
(3)终端位置约束：
[0022]
为保证制导火箭弹在末制导阶段对目标顺利追击，终端位置要接近目标区域，为便于追击目标，终端海拔不能过低，选取约束为：
[0023]
x
min
≤x(t
f
)≤x
max
,y
min
≤y(t
f
)≤y
max
ꢀꢀꢀꢀꢀ
(4)
[0024]
式(4)中x(t
f
)、x
min
、x
max
分别为终端水平飞行位置及其约束的最小值和最大值，y(t
f
)、y
min
、y
max
分别为终端海拔及其约束的最小值和最大值。
[0025]
(4)终端弹道倾角约束：
[0026]
为保证制导火箭弹在末制导阶段实现终端攻击角约束，设定较小的弹道倾角约束，有利于末制导终端约束角度的实现，选取约束为：
[0027]
θ
min
≤θ(t
f
)≤θ
max
ꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0028]
式(5)中θ(t
f
)、θ
min
、θ
max
分别为终端弹道倾角及其约束的最小值和最大值。
[0029]
选取目标适应度函数为：
[0030]
为保证弹道飞行能量消耗最低，本发明选取优化目标为在满足约束条件下，弹体保留能量最大即消耗能量最小，选取目标适应度函数为：
[0031][0032]
式(6)为能量消耗的相反数，j为适应度函数，为在α为控制量时使j为最大y(t0)、v(t0)分别为海拔和速度的初始值。
[0033]
步骤四、使用自适应遗传算法对弹道进行优化
[0034]
算法流程由弹道参数与适应度计算和遗传算法两部分组成，其中弹道参数与适应度计算具体步骤为：
[0035]
(1)采用上述五次多项式插值法对初代种群进行连续处理，即可将离散的控制变量近似为连续的控制函数；
[0036]
(2)将控制函数代入到系统状态微分方程组即式(1)中，使用龙格
‑
库塔法计算弹体每一时刻的运动参数；
[0037]
(3)依据求解的运动参数计算目标适应度函数，对于不满足步骤三中设定的约束条件的弹道个体的适应度采取罚函数处理，适应度均为正数，对其进行一定比例的缩小即可有效地排除掉不符合约束条件个体，防止出现错误的最优解，罚函数表达式为：
[0038]
j
c
＝cj
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0039]
式中j、j
c
分别表示运用罚函数处理前后的适应度，c为比例系数。
[0040]
遗传算法的具体执行步骤为：
[0041]
(1)生成新一代种群：确定初代种群规模n，并确定母体个数m、进化代数d等初始条件；
[0042]
(2)选择母体：依据计算得到的离散控制量个体及其对应的适应度函数值，从初代种群n个个体中，运用选择算子选择出m/2对母体；
[0043]
(3)交叉运算：对所选择的m/2对母体，依据交叉概率p
c
，运用十进制交叉算子执行交叉运算，生成m个中间个体；
[0044]
(4)变异运算：对m个中间个体，依据交叉概率p
m
，运用十进制变异算子执行变异运算，生成m个候选个体；
[0045]
(5)对m个候选个体进行适应度计算及选择，生成新一代母体，再次执行进化过程，直至满足终止条件；
[0046]
(6)计算输出最优个体及其对应的运动参数与弹道。
[0047]
本发明的优点及有益效果在于：优化后的弹道能够在满足过程约束、终端约束条件下，飞行能量消耗较小，在末制导阶段能有充足的能量对目标进行追击，同时满足过程约束可以防止弹体失稳，终端约束限定了弹体在终端的系统状态，以更有利于末制导进行的初始条件进入末制导阶段，保证末制导的顺利进行以及终端攻击角约束的实现。因此，本发明可以通过优化算法得到有利于末制导进行的低能耗弹道。
附图说明
[0048]
图1为制导火箭弹飞行过程示意图。
[0049]
图2为弹道优化方法流程示意图。
[0050]
图3为优化后弹道轨迹图。
[0051]
图4为攻角变化曲线图。
[0052]
图5为弹道倾角变化曲线图。
[0053]
图6为速度变化曲线图。
具体实施方式：
[0054]
步骤一、建立弹道运动学方程组
[0055]
优选地，方案飞行一般选取一个固定铅垂平面飞行，因此在纵平面建立弹道运动学方程组，经过推导化简得到如下方程组：
[0056][0057]
式(1)中式中m为弹体质量、v为弹体速度、α为攻角、θ为弹道倾角、g为重力加速度、t为飞行时刻；x、y分别为弹体的水平位置和海拔高度；c
x0
为零升阻力、c
x1
为诱导阻力；为升力系数导数、s为特征面积、q为动压q＝ρv2q＝ρv2，ρ为大气密度。
[0058]
步骤二、求解微分方程组
[0059]
在本发明中选取攻角α为控制变量，为便于优化算法，对攻角离散化处理，将弹道优化问题转化为非线性规划问题，对时间进行r等分，得到u＝[α1,α2,...α
r
,α
r 1
]作为待优化的控制变量。为保证实际控制中的连续性，在解算弹道参数前需要对控制量进行平滑处理，本发明采用五次多项式插值，其特点为被插值量的速度、加速度均为连续，保证插值后的控制过程连续，不会出现突变导致系统失稳。具体计算方法如下：
[0060][0061]
式中α、α
v
、α
a
分别为攻角、攻角变化的速率、攻角变化的加速度大小；t0为初始时刻、t为飞行时刻；b0、b1、b2、b3、b4、b5为中间函数，其具体计算方式如下：
[0062][0063]
式(3)中t0、t1、α0、α1、分别表示飞行时刻、攻角、攻角变化的速率、攻角变化的加速度大小在每一段起点和终点的状态，t为两个相邻控制点之间的时间间隔，d为攻角起点时刻和终点时刻的差值。
[0064]
对u＝[α1,α2,...α
r
,α
r 1
]中所有相邻的两个点应用上述方法计算，可得到最终的平滑连续的控制函数α(t)。
[0065]
通过α(t)以及初始状态采用四阶龙格
‑
库塔法求解微分方程组，得到弹道的位置、速度、弹道倾角等参数信息。具体求解方法如下：
[0066]
已知一个微分方程及其初值如下：
[0067]
y'＝f(t,y),y(t0)＝y0ꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0068]
式中y为一个关于t的函数，f(t,y)为包含t和y的函数，y(t0)＝y0表示在初始时刻t0的初值。
[0069]
令t
n
表示在初始时刻后经过n个时间间隔h的时刻，y
n
表示在t
n
时刻y的值。
[0070]
下一个时间间隔h后的y
n 1
值求解如下：
[0071][0072]
其中y
n 1
由现在时刻的值y
n
加上h乘以一个估算的斜率，该斜率是以下四个斜率的加权平均值：
[0073][0074]
式中
[0075]
k1是时间开始时的斜率；
[0076]
k2是时间段中点的斜率，通过欧拉法采用斜率k1来计算y在点t
n
h/2的值；
[0077]
k3是采用k2来计算y值的中点斜率；
[0078]
k4是时间段终点的斜率，采用k3计算y值
[0079]
在计算机中可以选择合适的步长h，计算初始时间后任意时刻的y值，若需要使用某一时刻的值，直接选取即可。
[0080]
步骤三、设定约束条件及目标适应度函数
[0081]
将过程约束、终端约束、优化目标及初始条件整理建立如下约束状态函数：
[0082][0083]
式(7)中y(t0)、v(t0)分别为海拔和速度的初始值，第一项为j趋于最大值，即消耗能量最小，可保证弹道在满足约束条件下能耗最小；
[0084]
式(7)中α
max
为允许控制攻角绝对值的最大值、v
min
为终端速度允许的最小值、x(t
f
)、x
min
、x
max
分别为终端水平飞行位置及其约束的最小值和最大值、y(t
f
)、y
min
、y
max
分别为终端海拔及其约束的最小值和最大值、θ(t
f
)、θ
min
、θ
max
分别为终端弹道倾角及其约束的最小值和最大值。
[0085]
式中x0、y0、θ0分别表示起控时刻初值为，x
min
、x
max、
、y
min
、y
max
、θ
min
、θ
max
分别为对系
统状态的期望约束值。
[0086]
步骤四、使用自适应遗传算法对弹道进行优化
[0087]
自适应遗传算法对弹道进行优化的流程如附图所示，具体包括如下步骤:
[0088]
算法流程由弹道参数与适应度计算和遗传算法两部分组成，其中弹道参数与适应度计算具体步骤为
[0089]
(1)采用上述五次多项式插值法对初代种群进行连续处理，即可将离散的控制变量近似为连续的控制函数；；
[0090]
(2)将控制函数代入到系统状态微分方程组即式(1)中，使用龙格
‑
库塔法计算弹体每一时刻的运动参数；
[0091]
(3)依据求解的运动参数计算目标适应度函数，对于不满足步骤三中设定的约束条件的弹道个体的适应度采取罚函数处理，适应度均为正数，对其进行一定比例的缩小即可有效地排除掉不符合约束条件个体，防止出现错误的最优解，罚函数表达式为：
[0092]
j
c
＝cj
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0093]
式中j、j
c
分别表示运用罚函数处理前后的适应度，c为比例系数。
[0094]
遗传算法的具体执行步骤为：
[0095]
(1)生成新一代种群：确定初代种群规模n，并确定母体个数m、进化代数d等初始条件；
[0096]
(2)选择母体：依据计算得到的控制量及其适应度，从初代种群n个个体中，运用选择算子选择出m/2对母体；本发明采取适应度轮盘赌法，可以依据个体被选择的概率多次选择较优的个体。其中每个个体每次被选择概率为：
[0097][0098]
式(9)中n为个体总数，i为个体的序号，j
i
为第i个个体的适应度函数值。
[0099]
轮盘赌法类似转盘，转盘面积按照适应度概率分布，选择n次旋转n次，留下选中的个体，具体选择方法实现流程如下：
[0100]
1)生成一组长度为n的0到1之间的随机数并从小到大排序；
[0101]
2)对每个个体的概率进行求和排序，相当于依次累加；
[0102]
3)依次对随机数与个体概率求和排序后结果比较，若随机数小于求和概率，即选中求和概率对应的个体，相当于每一个个体对应一个概率区间，若随机数在这个区间中即选中个体。
[0103]
(3)交叉运算：对所选择的m/2对母体，依据交叉概率p
c
，运用十进制交叉算子执行交叉运算，生成m个中间个体。为提高算法收敛速度，对交叉概率进行自适应处理，减少无效的交叉，交叉概率为：
[0104][0105]
式(10)中p
c1
为交叉概率的下界，p
c2
为交叉概率的上界；j
m
为互相交叉的两个母体
中较大的适应度值，j
avg
为所有母体适应度值的平均值，j
max
为所有母体适应度值的最大值。
[0106]
(4)变异运算：对m个中间个体，依据交叉概率p
m
，运用十进制变异算子执行变异运算，生成m个候选个体。为提高算法收敛速度，对变异概率进行自适应处理，减少无效的变异，变异概率为：
[0107][0108]
式(11)中p
m1
为交叉概率的下界，p
m2
为交叉概率的上界；j
m
为互相交叉的两个母体中较大的适应度，j
avg
为所有母体适应度的平均值，j
max
为所有母体适应度的最大值。
[0109]
(5)对m个候选个体进行适应度计算及选择，生成新一代母体，再次执行交叉、变异选择的进化过程，直至满足终止条件；
[0110]
(6)计算输出最优个体及其对应的运动参数与弹道。
[0111]
在本发明的实际应用时，可针对不同目标位置及弹体初始状态进行方案弹道规划，优化后的弹道满足所选取的约束。选取仿真条件如下：
[0112]
(1)初始状态：弹道倾角θ0＝0
°
，x0＝20000m，y0＝15000m，攻角α0＝0
°
，速度v0＝350m/s；
[0113]
(2)控制作用容许范围：|α|≤10
°
；
[0114]
(3)终端约束参数：终端速度v
f
≥450m/s，终端弹道倾角
‑
15
°
≤θ
f
≤10
°
，终端海拔3000m≤y
f
≤4000m；
[0115]
(4)初代种群规模n＝1000，母体个数m＝100，进化代数d＝50，p
c1
＝0.9，p
c2
＝0.4，p
m1
＝0.0001，p
m2
＝0.3。
[0116]
利用本发明所述的弹道优化方法对上述弹道方案进行优化，仿真结果如图3～图6所示。
[0117]
由图3～图6所示仿真结果可以看出，优化后的弹道采取先维持高空飞行，再迅速下降接近目标，最后降低弹道倾角较平缓飞行的飞行方式。
[0118]
最终终端水平位置为x
f
＝52.01km，海拔y
f
＝4.05km，弹道倾角θ
f
＝
‑
11.5
°
，终端速度v
f
＝558.2m/s，均满足终端约束要求。
[0119]
同时攻角没有频繁较大的突变，没有超过攻角控制容许范围，满足|α|≤10
°
，不会导致弹体失稳。前期弹体所处海拔较高，与低空相比，高空中空气密度较低，由于q＝ρv2，动压q较小，从弹体运动方程组式(1)中的qs(c
x0
c
x1
α2)和可以看出，需要较大的攻角实现等同效果。由于在高空中飞行动压小，消耗能量小，因此弹道前期攻角较大，降低弹道倾角变化速率，维持高空飞行并起增程作用；中期减小攻角，增大弹道倾角变化速率，以小攻角飞行消耗能量较小，海拔迅速下降接近目标；最后增大攻角控制，增大弹道倾角变化率，弹道倾角绝对值较快速减小，使其趋近于目标值。以较小的弹道倾角进入末制导阶段，在制导中弹道倾角将会有较大的变化空间，更有利于实现制导律中固定视线角打击的约束条件。
[0120]
在满足过程约束和终端约束条件下，弹道的目标适应度函数值为最大，即为能量消耗最小，经本发明方法优化后的弹道为满足系统状态约束的低能耗弹道。
[0121]
以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。