航空制导炸弹具有有效载荷比高、命中精度高、成本相对低廉等优点, 是作战飞机近距空中支援、执行纵深打击、压制和摧毁敌方防控系统等任务最有效的内埋式武器^[1]。由于无发动机提供动力, 其射程主要取决于投放时的高度、速度, 因此需要进行轨迹优化, 一般采用“滑翔+俯冲攻击”的方案, 但这种做法并不能得到最优的轨迹。

从原理上讲, 轨迹优化问题可归结为多约束下的非线性最优控制问题, 有间接法和直接法^[2]2大类求解方法。间接法能够求得精确解, 但存在初值难以估计, 收敛性无法保证的问题^[3]。而直接法是将最优控制问题转化为有限维的非线性规划问题(NLP), 应用NLP求解器进行求解。直接法中的配点法, 尤其是近期发展起来的伪谱法, 因其计算效率上的优势, 逐渐成为求解NLP问题的研究热点。

Elnagar等^[4]首次对伪谱法的转换原理, 收敛速度、求解效率等方面做了深入研究, 奠定了伪谱法求解最优控制问题的基础。伪谱法是将时间离散化为若干离散点, 指标函数、微分方程和约束条件由离散点上的状态变量、控制变量表示, 再利用NLP求解器求解, 相比于传统的智能优化算法在计算时间、控制律的平滑度等方面具有明显的优越性^[5]。常见的伪谱法包括:Legendre伪谱法(LPM)、Gauss伪谱法(GPM)、Radau伪谱法(RPM)。Huntington等^[6]指出上述3种方法在计算精度方面, GPM、RPM相当, 均高于LPM; Gandhi等^[7]研究表明在给定求解精度的条件下, RPM比GPM的计算代价更小; Rexius等^[8]基于RPM研究了在气动加热、禁飞区等多种约束条件下的高超声速飞行器的轨迹优化问题; 袁宴波等^[9]研究了Radau伪谱法在最优滑翔弹道应用方面的问题, 给出了基于协态映射原理的数值解验证方法。上述方法由于其节点分布形状固定, 容易导致需要较高维的插值多项式才能得到理想的近似解, 甚至在采用高维多项式后仍然得不到理想的近似解。而自适应hp伪谱法具有计算稀疏性、快速收敛性的特点, 对于需要改进求解精度的区间, 可细化网格或提高多项式的阶次, 这种方法能够大大减小计算代价, 在多约束下的轨迹快速优化问题中, 具有较大的优势。本文利用基于hp自适应Radau伪谱法(hp-RPM)的优点, 采用hp-RPM求解策略对航空制导炸弹在多约束条件下的滑翔轨迹进行快速优化。

1 轨迹优化问题 1.1 航空制导炸弹纵向运动模型

轨迹优化任务中, 一般只考虑纵向运动的情况, 为加快算法的收敛速度, 宜采用无量纲化的制导炸弹纵向运动模型。

引入参考量L_r、V_r, 则无量纲纵向运动模型为:

(1)

式中, 无量纲射程、高度、速度; 大气密度、参考面积、重力加速度g=为弹道倾角; C_D为阻力系数; C_L为升力系数; m为质量。为描述方便, 下文省去上标。

升力系数C_L、阻力系数C_D可近似为攻角的函数:

(2)

式中, Ma为飞行马赫数; C_L1、C_D0、C_D2为在不同马赫数下根据攻角α拟合得到的气动力系数。气动参数根据某航空制导炸弹的启动数据拟合获得, 拟合结果如图 1所示。

大气密度、音速采用标准大气拟合模型, 其精确度能够满足航空制导炸弹的飞行空间。

(3)

1.2 约束条件

根据制导炸弹攻击目标的特点, 对边界条件有一定的约束, 如大的弹着角、末端速度, 这样能够保证制导炸弹具有一定的侵彻能力。

1)边值约束

初值约束:

(4)

终端约束:

航空制导炸弹对地攻击对弹道倾角和碰撞速度有一定的要求, 以保持一定的动能侵彻目标, 即

(5)

2)过程约束

航空制导炸弹因无推力, 操纵舵效率有限等因素, 其需要满足一些诸如过载、攻角变化率、攻角幅值等的限制, 建立如下过程约束条件:

(6)

式中, 为最大攻角变化率; n_max为制导炸弹所能承受的最大法向过载。

1.3目标函数

对于轨迹优化, 在满足约束条件的前提下, 射程越远意味着攻击范围更广, 因此, 取目标函数为

(7)

2 hp-RPM求解最优轨迹 2.1 Radau基本原理

Radau伪谱法将控制变量和状态变量离散化为时间轴上子区间[t_k-1, t_k]的一系列Legendre-Gauss-Radau (LGR)点, 通过构造Lagrange插值多项式对其进行逼近, 并对该多项式的求导以代替动力学微分方程^[10]。

1)时域变换

Radau伪谱法的定义域为[-1,  1], 而实际问题的定义域为[t₀, t_f], 将[t₀, t_f]分成K个网格, 且∀t∈[t_k-1, t_k], k=1, …, K, t₀=t₁ < … < t_K=t_f, 令

(8)

2)状态变量与控制变量的插值近似

令(τ₁, τ₂, …, τ^(k)_Nk)表示区间[t_k-1, t_k]的LGR配点, 它与τ^(k)_Nk+1=+1构成了区间内的N_k+1个Lagrange插值点, 其中τ₁^(k)=-1。采用Lagrange插值多项式作为基函数来近似状态变量和控制变量, 即

(9)

(10)

式中, L_i(τ)和为Lagrange插值基函数。

3)动力学微分方程约束转换为代数约束

对(9)式求导, 有

(11)

则N_k个LGR点上的动力学微分方程约束变为代数方程约束, 即

(12)

式中，为第k个子区间的Radau伪谱微分矩阵[17]。

4)约束条件的离散化

子区间内N_k个配点上的路径不等式约束为

(13)

边界条件

(14)

网点约束, 有

(15)

5) LGR积分的性能指标函数

用LGR积分近似性能指标函数

(16)

式中, w_j^(k)为第k个子区间内的LGR权重, 形式如下

(17)

经过上述变换后, 最优控制问题转化为求解离散节点上的状态X_i和控制变量U_k, (i, k=1, …, K), 以及初末时刻t₀、t_f(未给定), 使得性能指标(7)最小, 并满足动力学方程约束(1), 终端条件约束(5), 以及过程约束条件(6)的最优控制问题。

2.2 hp-RPM自适应策略设计

hp自适应方法根据网格的曲率和约束方程的误差判定自动调整网格数量和多项式阶次^[11]。通过判断路径约束误差决定是否新增区间。在区间内通过判断轨迹的曲率来决定采用增加配点数的策略, 还是增加多项式阶次的策略。设允许的最大误差为ε_max, 最大允许曲率为r_max, 则误差和曲率比判据为:

1)误差判据

设[t_k-1, t_k]在第k个网格子区间, 其配点数为N_k, 令t_q^(k)为相邻配点间等间距分布的3个点, 即

(18)

式中, i=1, …, N_k为第k个子区间的第i个配点。

则该区间上的采样点在约束方程上的残差为:

(19)

(20)

式中:l=1, …, n为状态变量个数; j=1, …, s为路径约束个数; p=1, …, N_k-1为采样点; ε_lp^(k)为第l个状态在第p个采样点处的绝对值残差。ε_jp^(k)为第j个状态在第p个采样点处的绝对值残差。

对于任意, l∈(1, 2, …, n), j=(1, 2, …, s), p=(1, …, N_k-1), e_max^(k)=max{ε_lp^(k), ε_jp^(k)}=ε_lp^(k), 误差判据为ε_lp^(k)≤ε_max。

2)曲率比判据

设(x_lp)_n×3(Nk-1)表示由n个状态变量在p个采样点上的离散值所构成的状态矩阵, X_l^(k)表示矩阵中相应于ε_max^(k)的状态变量所在的行向量, 即

(21)

于是, 第l个状态在第p个采样点处的曲率为

(22)

令表示曲率的算术平均值, 即

(23)

令e为采样点上的曲率与曲率平均值的比所构成的向量。则e的表达式如下:

(24)

因此, 曲率判据为e < r_max。

2.3 算法步骤

hp-RPM求解航空制导炸弹轨迹优化设计问题的步骤如下所示, 其迭代流程如图 2所示。

1)初始化, 给定初始状态x₀, 误差门限值ε_max, 选取K个区间, 每个区间设定N个点;

2)求解NLP, 若ε_max^(k)≤ε_max, 迭代停止; 否则进入步骤3);

3)若ε_max^(k)≤ε_jp^(k), 则转步骤4), 否则转步骤5);

4)在局部残差最大值所对应的采样点上将区间进一步细化, 新增区间配点数为N;

5)若e < r_max, 则增加该区间配点数, 否则以局部最大值所对应的采样点作为新的网点, 新增区间配点数为N;

6)所有子区间更新后, 返回步骤2)。

3 仿真算例与结果分析 3.1 仿真条件

本文以某航空制导炸弹为例, 仿真条件为:制导炸弹质量m=50 kg、参考面积S=0.020 1 m², 迭代误差门限值ε_max=1×10^-6, 当ε_max^(k)≤ε_max时, 停止迭代。

初始条件:

终端条件:

过程约束:

3.2 仿真结果

根据上述条件, 选取V_r=240 m/s, L_r=7 400 m代入(1)式的无量纲化模型, 采用第2节的hp-RPM优化求解策略对制导炸弹轨迹进行优化, 优化后的状态变量和控制变量曲线如图 3所示。

其中, CM代表常规方法, Hamilton函数曲线如图 4所示。为验证该增程效果, 对其他不同初始条件下的轨迹快速优化问题进行了仿真, 表 1为在不同初始条件下经hp-RPM优化后和与常规方案相比的仿真结果。

如图 3a)所示, 经hp-RPM优化后的曲线显得更加平缓光滑, 特别是弹道的后半段, 这样的安排, 使得制导炸弹飞的更远。图 3b)显示, 两者在速度方面, 常规方案的波动幅度要大, 图 3c)展示, 两者的弹道倾角基本变化趋势基本一致, 但优化后的方案速度变化更加平滑, 图 3d)~图 3f)所示两者在弹道初始段和弹道末段的控制量相差较大, 在30~100 s的中段两者基本一致, 两者都满足了攻角、过载、攻角速率、终端速度、终端弹道倾角的约束, 而优化后, 制导炸弹的飞行时间由原来的180 s增加至195 s。另外, 若终端时刻自由, 由极小值原理可知, 哈密尔顿函数应该保持为零, 图 4可知, 哈密顿函数在10^-5量级, 在一定程度上说明了解的最优性。

该优化计算在CPU为3.3GHz/i5-4590, 内存8G, Windows 7系统, Matlab环境下编程实现。算例中, 优化的时间开销为1.14 s, 而在Matlab环境下求解一个二维优化问题, 一般需要30分钟。因此, hp-RPM在计算速度上具有很大的优势。

从表 1可知, 经hp-RPM优化后制导炸弹的射程得到了较大的改善, 均有不同程度的提高。其中在序号1和2的条件下, 增程小于10%, 这主要在于这2个初始条件的初始高度较高, 制导炸弹滑翔段时间较长, 占比整个弹道的比例较大, 常规方案中俯冲攻击段使得飞行时间减少。在7 000 m高度, 240 m/s速度投放, 能够增程10.94%, 而在180 m/s时能够增程8.42%, 这是因为220 m/s的末端速度约束, 迫使俯冲攻击段需要调整较大的角度, 导致飞行时间减少。在5 000 m高度投放, 最大可以增加17.37%的射程, 增程比较明显。

综上, hp-RPM能在多约束条件下快速生成符合要求的轨迹, 计算量小, 计算精度高, 适合于工程应用。

4 结论

本文采用了一种基于hp自适应的Radau伪谱法对航空制导炸弹进行快速轨迹优化, 首先建立了优化模型, 将优化问题转换为NLP问题后, 利用hp-RPM对其进行求解。仿真结果表明, 在多约束条件下, 该方法能够快速的生成符合要求的轨迹, Hamilton函数曲线表明优化结果满足一阶最优性条件, 计算速度快, 射程增加幅度平均达到11.46%。对多约束条件下的轨迹优化问题, 具有一定的参考价值。