CN103728976B - 一种基于广义标控脱靶量概念的多过程约束和多终端约束末制导律 - Google Patents

Abstract

本发明公开了一种基于广义标控脱靶量概念的多过程约束和多终端约束末制导律，包括以下几个步骤：步骤1：加载存储的标准控制及飞行器状态初值；步骤2：弹道积分预测；步骤3：过程约束修正；步骤4：标控弹道一阶变分模型求解；步骤5：更新标准控制求解；本发明相较于传统的末制导律，该制导律不仅可以满足脱靶量、落角、落点攻角等传统的终端约束的要求，还能够使强非线性状态变量V收敛到所需要的值，同时还能闭环的处理过程约束。

Description

一种基于广义标控脱靶量概念的多过程约束和多终端约束末制导律

技术领域

本发明涉及一种基于广义标控脱靶量概念的多过程约束和多终端约束末制导律，属于航天技术、武器技术领域。

背景技术

高超声速飞行器是指飞行马赫数大于或者等于5的飞行器，它具有飞行速度快、突防能力强、全球到达、毁伤威力大等独特的有时，已经成为当今世界各国武器研制的热点和焦点。而下压段多约束末制导律则是决定其作战效能的一个重要因素。

在高超声速飞行器的飞行末段，为了增大突防能力，往往要求以较大的弹道倾角下压，这就意味着末段的需用过载较大；而飞行器的速度大、高度高，则意味着它的机动能力差；即容易出现控制饱和的情况。这就要求在末制导中不仅要考虑脱靶量、碰撞角和落速等终端约束，还需考虑最大攻角约束和最大过载约束等过程约束。

现在末制导常用的显示导律(如比例导引，广义显示制导等)，往往基于忽略速度变化的线性化几何弹道模型，通过控制方法使得弹道最终收敛到所需要的目标，但是这些制导律没法满足落速精确控制的要求。虽然在再入制导律中存在对末速的考虑，但末制导段的近垂直打击要求使得再入制导律不适用。并且，这些制导律在考虑过程约束时往往没有预先的闭环考虑，在控制能力不足的情况下很容易出现不能命中目标。因此，需要一种能够闭环考虑落速要求和过程约束的制导律。

发明内容

本发明的目的是为了解决上述问题，提出一种同时满足脱靶量、落角、落速和落点攻角等终端约束以及最大攻角约束和最大过载约束的过程约束的高超声速飞行器末制导律。

广义标控脱靶量是指从t₀时刻开始飞行器按照标准控制u^*(t)飞行，在t_f时刻飞行器的飞行状态与给定末状态的偏差。广义标控脱靶量不仅包含传统意义上的脱靶量，也可以包含终端碰撞角偏差，甚至还可以包含终端碰撞速度偏差。在获得了广义标控脱靶量后，进一步采用梯度法为制导修正算法，对非线性问题取广义梯度，将其化为线性时变优化问题，利用多次线性时变优化问题来逼近非线性优化问题。最终通过求解基于当前弹道的线性化模型来获取对当前弹道的改进，从而不断的使标控弹道不断逼近最优解。

对于在制导中需要考虑的最大攻角、最大法向过载等过程约束，本发明提出了一种基于控制走廊概念的处理方法。在进行在线制导的时候，由于每一个周期积累的误差都为一小量，因此每一个新的迭代周期的规划弹道相对于上一周期的规划弹道的改变也会是一个小量，从而使得局部线性化模型能够很好的逼近真实的非线性模型。

本发明的优点在于：

(1)相较于传统的末制导律，该制导律不仅可以满足脱靶量、落角、落点攻角等传统的终端约束的要求，还能够使强非线性状态变量V收敛到所需要的值，同时还能闭环的处理过程约束；

(2)该制导律基于真实的非线性运动模型，所获得的制导律更加接近真实的最优解。

(3)该制导律还可以根据实际的制导需求减少或者增加相应的标控脱靶量个数，从而应对不同的需求，具有广泛的适用性。

附图说明

图1是本发明制导律制导指令生成的流程图；

图2是本发明制导律的使用流程图；

图3是弹道迭代收敛图；

图4是无拉偏制导仿真结果；

图5是大气风模型；

图6是风拉偏制导仿真结果；

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种基于广义标控脱靶量概念的多过程约束和多终端约束末制导律，在运用基于广义标控脱靶量概念的多过程约束和多终端约束末制导律，首先需要给出标准控制的初值估计，在这里，选用广义显示制导律生成控制初值猜测。接下来，需要计算标控弹道和标控脱靶量。所谓标控弹道，是指从当前状态出发依据给定的标准控制飞行所得到的弹道，标控弹道终端状态与需要终端状态的差为标控脱靶量。最后采用梯度法求解基于标准弹道线性化，得到的最优控制修正，从而实现对标准控制的更新，整个制导过程如图1所示，包括以下几个步骤：

步骤1：加载存储的标准控制及飞行器状态初值。

所谓标准控制，即为飞行器事先设定的控制指令规律，设为u_b(x)，飞行器状态初值则是当前时刻飞行器的飞行状态，包括初始速度V₀，初始弹道倾角γ₀，初始高度h₀以及初始纵程x₀。需要注意的是，在制导初始时刻，标准控制可能存在较大的偏差，但随着制导迭代修正的进行，标准控制会迅速收敛到所需要的最优控制上。

步骤2：弹道积分预测

首先要建立高超声速飞行器末制导段动力学模型。由于高超声速飞行器具有飞行距离远、突防能力强的特点，因此一般用于攻击地面高价值目标。同时，由于其高速性导致末端机动能力有限，因此，可以假定在中制导阶段已经使飞行器对准目标，在本发明中，仅考虑纵向平面内的制导问题，飞行器的运动方程为，

$\stackrel{\·}{V}=-g\sin \mathrm{\γ}-\frac{D}{m}$

$\stackrel{\·}{\mathrm{\γ}}=-\frac{gcos\mathrm{\γ}}{V}+\frac{L}{mV}---\left(1\right)$

$\stackrel{\·}{h}=V\sin \mathrm{\γ}$

$\stackrel{\·}{x}=Vcos\mathrm{\γ}$

式中，V为飞行器相对地球速度，γ为弹道倾角，h为高度，x为纵程，和分别表示飞行器相对地球速度、弹道倾角、高度和纵程对时间的导数。L和D分别为升力和阻力，其表达式分别为L＝0.5ρV²S_refC_L和D＝0.5ρV²S_refC_D。其中ρ为大气密度，一般与高度相关；S_ref为飞行器的气动参考面积；C_L和C_D分别为升力系数和阻力系数，一般与攻角α相关。g为当地重力加速度，与高度相关，一般变化较小；m为飞行器质量，由于在末制导段发动机关机，故质量为常量。观察运动方程(1)，可以发现其中不显含时间相关项，并且在整个末制导段纵程x单调递增，因此可以将所有状态变量转为对x进行微分，这样运动方程可降维变为：

$\frac{dV}{dx}=\frac{-g}{V}tan\mathrm{\γ}-\frac{D}{mVcos\mathrm{\γ}}$

$\frac{d\mathrm{\γ}}{dx}=\frac{u}{cos\mathrm{\γ}}---\left(2\right)$

$\frac{dh}{dx}=tan\mathrm{\γ}$

其中，u＝-gcosγ/V²+L/(mV²)，设为控制变量。

进一步，利用上一步中的标准控制u_b(x)，根据式(2)以及式(6)(步骤2和步骤3同时进行)，可以获得过程约束影响下的控制修正量权重因子c(x)；经过过程约束修正的标准控制标准速度曲线V_b(x)、标准弹道倾角曲线γ_b(x)，标准高度曲线h_b(x)，以及标准的阻力曲线D_b(x)。(注：公式(2)为飞行器运动的动力学方程，需要根据给定的控制u(x)才能得到对应的弹道V(x),γ(x),h(x)，这里给定的控制为u_b(x)，故得到的弹道为上述的V_b(x)，γ_b(x)，h_b(x)，之所以加下标为与控制对应)

步骤3：过程约束修正

在进行步骤2的同时，同时开展步骤3。

在高超声速飞行器的末制导段，考虑的过程约束主要是最大攻角约束和最大过载约束。

(1)攻角约束：

假设攻角约束为，

α_min≤α≤α_max(3)

其中，α_min为最小攻角；α_max为最大攻角。将其转化为关于升力系数的约束，

C_Lmin(Ma)≤C_L(Ma)≤C_Lmax(Ma)(4)

其中，C_Lmin(Ma)为在当前马赫数下最小攻角对应的升力系数；C_Lmax(Ma)为当前马赫数下最大攻角对应的升力系数。若标准弹道的升力系数满足方程的约束关系，则不做任何处理。若不满足，假定(其中，为根据标准控制积分弹道所对应的升力系数；为根据标准控制积分弹道所对应的最小升力系数；为在标准控制u_b(x)下不满足该约束的纵程x序列，在标控弹道计算的时候生成)，则令同时，令控制指令更新时的指令修正量Δu(x^*)＝0(如式(22)所示)，这样修正弹道在处的升力系数也满足C′_L(Ma，x^*)＝C′_Lmin(Ma，x^*)。同样，对于(其中，为根据标准控制积分弹道所对应的最小升力系数；为在标准控制u_b(x)下不满足该约束的纵程x序列，在标控弹道计算的时候生成)，则令同时，令控制指令更新时的指令修正量(如式(22)所示)。

(2)过载约束处理

假定过载约束为：

n_y≤n_ymax(5)

其中，n_y为法向过载，n_ymax为最大法向过载。若标准弹道的过载满足方程满足约束关系，则不做任何处理。若不满足，假定(其中，为根据标准控制积分弹道所对应的法向过载；为在标准控制u_b(x)下不满足该约束的纵程x序列，在标控弹道计算的时候生成)则令同时，令控制指令更新时的指令修正量(如式(22)所示)。

为了防止在迭代过程中偶然出现的升力系数超出范围导致整个迭代过程升力系数被锁定的情况出现，可令Δu(x^*)＝Δu₀(x^*)/C_limit，其中，为不满足约束的纵程序列集合；Δu₀(x^*)为不经过过程约束修正的控制量修正；Δu(x^*)为经过过程约束修正的控制修正量；C_limit为一足够大的常数，目的是使得在约束碰边界之后，控制修正量为零，即Δu(x^*)≈0同样实现了对攻角的约束。

总的来说，过程约束对标准控制及控制修正的影响如下式所示。

$u_b^{old}\left(x\right)=\left\{\begin{array}{cc}\mathrm{\α}\left(u_b^{old}\right)={\mathrm{\α}}_{\min }& \mathrm{\α}={\mathrm{\α}}_{\min }\\ \mathrm{\α}\left(u_b^{old}\right)={\mathrm{\α}}_{\max }& \mathrm{\α}={\mathrm{\α}}_{\max }\\ n_y\left(u_b^{old}\right)=n_{y\max }& n_y=n_{y\max }\\ u_b\left(x\right)& {\mathrm{\α}}_{\min }\≤\mathrm{\α}\≤{\mathrm{\α}}_{\max }\∪n_y\≤n_{y\max }\end{array}\right.---\left(6\right)$

其中，c(x)为过程约束影响下的控制修正量权重因子；经过过程约束修正的标准控制；u_b(x)为弹上存储的标准控制量；表示选择使得α＝α_min时的控制量作为当前控制量；表示选择使得α＝α_max时的控制量作为当前控制量；表示选择使得n_y＝n_ymax时的控制量作为当前控制量。。

步骤4：标控弹道一阶变分模型求解

在标准弹道附近求解一阶变分可得：

$\left[\begin{array}{c}\frac{d\mathrm{\Δ}V}{dx}\\ \frac{d\mathrm{\Δ}\mathrm{\γ}}{dx}\\ \frac{d\mathrm{\Δ}h}{dx}\end{array}\right]=\left[\begin{array}{ccc}{f}_1& f_2& f_3\\ 0& f_4& 0\\ 0& f_5& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}V\\ \mathrm{\Δ}\mathrm{\γ}\\ \mathrm{\Δ}h\end{array}\right]+\left[\begin{array}{c}{g}_1\mathrm{\Δ}u\\ g_2\mathrm{\Δ}u\\ 0\end{array}\right]---\left(7\right)$

其中，ΔV、Δγ和Δh分别为在标控弹道附近的速度、弹道倾角和高度的微小变化量；Δu为控制量的小增量。

$\begin{array}{cc}{f}_1=\frac{g}{V_b^2}{\mathrm{tan\γ}}_b-\frac{D_b}{{\mathrm{mV}}_b^2{\mathrm{cos\γ}}_b}& f_4=-u_b^{old}\frac{{\mathrm{sin\γ}}_b}{{\cos }^2{\mathrm{\γ}}_b}\end{array}$

$\begin{array}{cc}{f}_2=-\frac{g}{V_b{\cos }^2{\mathrm{\γ}}_b}-\frac{D_b{\mathrm{sin\γ}}_b}{{\mathrm{mV}}_b{\cos }^2{\mathrm{\γ}}_b}& f_5=\frac{1}{{\cos }^2{\mathrm{\γ}}_b}\end{array}$

$\begin{array}{cc}{f}_3=\frac{D_b}{{\mathrm{mV}}_b{\mathrm{cos\γ}}_b}{\mathrm{\β}}_h& g_1=-\frac{S_{ref}{V}_b}{2m{\mathrm{cos\γ}}_b}\frac{\∂C_D}{\∂C_L}\end{array}$

$g_2=\frac{1}{{\mathrm{cos\γ}}_b}$

其中，β_h为指数大气模型常数，一般为1/7200；为阻力系数对升力系数的偏导数，满足为阻力系数对攻角的变化率，为升力系数对攻角的变化率。

定义纵向弹道优化的目标函数为：

$\min J=K_V{(V_f-V_{object})}^2+{\∫}_{x_0}^{x_f}R\left(x\right)u^{2}dx---\left(8\right)$

其中，K_V为末速度约束松弛因子，用以调整末速度约束的强度；R(x)为控制权重因子；V_f为末速度大小(非标控弹道末速)；V_object目标末速度大小。对目标函数进行线性化处理，可得，

$\min J=K_V{({\mathrm{\ΔV}}_f+V_b(x_f)-V_{object})}^2+{\∫}_{x_0}^{x_f}2{\mathrm{Ru}}_b^{old}\mathrm{\Δ}u+{\mathrm{R\Δu}}^{2}dx---\left(9\right)$

其中x_f为目标所在处的纵程，ΔV_f为摄动弹道的终端速度增量。为了使得标控弹道收敛到目标值，定义如下边界条件，

$\left.\begin{array}{c}\mathrm{\Δ}h\left(x_0\right)=0,\mathrm{\Δ}V\left(x_0\right)=0,\mathrm{\Δ}\mathrm{\γ}\left(x_0\right)=0\\ {\mathrm{\Δh}}_f=h_{object}-h_b\left(x_f\right),{\mathrm{\Δ\γ}}_f={\mathrm{\γ}}_{object}-{\mathrm{\γ}}_b\left(x_f\right)\end{array}\right\}---\left(10\right)$

其中，Δh(x₀)、ΔV(x₀)和Δγ(x₀)分别为高度、速度和弹道倾角在弹道起始处的摄动值；Δh_f和Δγ_f分别为高度、弹道倾角在弹道终点处的摄动值(即为标控脱靶量)；h_object和γ_object分别为目标高度和目标弹道倾角。从而得到了纵平面内的标控弹道修正模型，包括末状态偏差(即标控脱靶量)、线性化模型系数矩阵(即由f₁、f₂、f₃、f₄、f₅、f₆、g₁和g₂组成的矩阵，见式(7))。

步骤5：更新标准控制求解

标控弹道修正模型的哈米尔顿函数为：

H＝λ₁(f₁ΔV+f₂Δγ+f₃Δh+g₁Δu)+λ₂(f₄Δγ+g₂Δu)+λ₃f₅Δγ+2Ru_bΔu+RΔu²(11)

其中，λ₁、λ₂和λ₃分别为ΔV、Δγ和Δh对应的协态变量。为将H对状态变量ΔV、Δγ和Δh求偏导数，可得协态方程：

${\stackrel{\·}{\mathrm{\λ}}}_1=-\frac{\∂H}{\∂\mathrm{\Δ}V}=-{\mathrm{\λ}}_1{f}_1---\left(12\right)$

${\stackrel{\·}{\mathrm{\λ}}}_2=-\frac{\∂H}{\∂\mathrm{\Δ}q}=-({\mathrm{\λ}}_1{f}_2+{\mathrm{\λ}}_2{f}_4+{\mathrm{\λ}}_3{f}_5)---\left(13\right)$

${\stackrel{\·}{\mathrm{\λ}}}_3=-\frac{\∂H}{\∂\mathrm{\Δ}h}=-{\mathrm{\λ}}_1{f}_3---\left(14\right)$

迈耶函数为：

Φ＝K_V(ΔV(x_f)+V_b(x_f)-V_object)²

+ν₁(Δγ(x_f)-Δγ_f)+ν₂(Δh(x_f)-Δh_f)

其中，x_f为弹道终点；ν₁和ν₂分别为拉格朗日因子(对应每一个终端约束)；Δγ(x_f)、Δh(x_f)为和ΔV(x_f)由标控弹道一阶变分模型求解获得的终端弹道倾角、高度和速度摄动量；Δγ_f和Δh_f为标控脱靶量(见式(10))。

由式(12)可得：

λ₁＝λ₁₀f_λ11(15)

其中，λ₁₀为λ₁的协态初值；f_λ11为积分表达式，如下，

$f_{\mathrm{\λ}11}=\exp (-{\∫}_0^x{f}_{1}dx)$

由式(14)和式(15)可得：

${\mathrm{\λ}}_3={\mathrm{\λ}}_{30}-{\mathrm{\λ}}_{10}{\∫}_0^x{f}_{\mathrm{\λ}11}{f}_3{\mathrm{dx}}_1---\left(16\right)$

其中，λ₃₀为λ₃的协态初值。

将式(15)和式(16)带入式(13)可得：

λ₂＝λ₁₀f_λ21+λ₂₀f_λ22+λ₃₀f_λ23(17)

其中，λ₂₀为λ₂的协态初值；f_λ21、f_λ22和f_λ23均为积分表达式，

$f_{\mathrm{\λ}21}=\exp (-{\∫}_0^x{f}_{4}dx){\∫}_0^x(-(f_2{f}_{\mathrm{\λ}11}-f_5{\∫}_0^x{f}_{\mathrm{\λ}11}{f}_{3}dx\left)\exp \right({\∫}_0^x{f}_{4}dsx\left)\right)dx$

$f_{\mathrm{\λ}22}=\exp (-{\∫}_0^x{f}_{4}dx)$

$f_{\mathrm{\λ}23}=\exp (-{\∫}_0^x{f}_{4}dx){\∫}_0^x(-f_5\exp ({\∫}_0^x{f}_{4}dx\left)\right)dx$

最优轨线满足：

$\frac{\∂H}{\∂\mathrm{\Δ}u}=2\mathrm{\Δ}u+{\mathrm{\λ}}_2{g}_2+{\mathrm{\λ}}_1{g}_1+2u_b^{old}=0---\left(18\right)$

将式(17)代入式(18)，并考虑过程约束修正，可得：

$\begin{array}{c}\mathrm{\Δ}u=-\frac{({\mathrm{\λ}}_{10}{f}_{\mathrm{\λ}21}+{\mathrm{\λ}}_{20}{f}_{\mathrm{\λ}22}+{\mathrm{\λ}}_{30}{f}_{\mathrm{\λ}23})g_2+{\mathrm{\λ}}_{10}{f}_{\mathrm{\λ}11}{g}_1+2u_b^{old}}{2c}\\ =-{\mathrm{\λ}}_{10}\frac{f_{\mathrm{\λ}11}{g}_1+f_{\mathrm{\λ}21}{g}_2}{2c}-{\mathrm{\λ}}_{20}\frac{f_{\mathrm{\λ}22}{g}_2}{2c}-{\mathrm{\λ}}_{30}\frac{f_{\mathrm{\λ}23}{g}_2}{2c}-\frac{u_b^{old}}{c}\end{array}---\left(19\right)$

其中，c即为c(x)，由式(6)确定。将式(19)带入式(7)，积分后带入边界条件，可得关于λ₁₀、λ₂₀和λ₃₀的三元一次方程组：

$\left\{\begin{array}{c}{V}_b\left(x_f\right)-V_{object}=f_{V0}\left(x_f\right)+{\mathrm{\λ}}_{10}\left(f_{V1}\right(x_f)-f_{\mathrm{\λ}11}(x_f)/2K_V)\\ +{\mathrm{\λ}}_{20}{f}_{V2}\left(x_f\right)+{\mathrm{\λ}}_{30}{f}_{V3}\left(x_f\right)\\ {\mathrm{\Δ\γ}}_f=f_{\mathrm{\γ}0}\left(x_f\right)+{\mathrm{\λ}}_{10}{f}_{\mathrm{\γ}1}\left(x_f\right)+{\mathrm{\λ}}_{20}{f}_{\mathrm{\γ}2}\left(x_f\right)+{\mathrm{\λ}}_{30}{f}_{\mathrm{\γ}3}\left(x_f\right)\\ {\mathrm{\Δh}}_f=f_{h0}\left(x_f\right)+{\mathrm{\λ}}_{10}{f}_{h1}\left(x_f\right)+{\mathrm{\λ}}_{20}{f}_{h2}\left(x_f\right)+{\mathrm{\λ}}_{30}{f}_{h3}\left(x_f\right)\end{array}\right.---\left(20\right)$

其中，f_V0(x_f)、f_γ0(x_f)与f_h0(x_f)为由式中的带入式(7)得到的结果；f_V1(x_f)、f_γ1(x_f)与f_h1(x_f)为由式中的-(f_λ11g₁+f_λ21g₂)/(2c)带入式(7)得到的结果；f_V2(x_f)、f_γ2(x_f)与f_h2(x_f)为由式中的-f_λ22g₂/(2c)带入式(7)得到的结果；f_V3(x_f)、f_γ3(x_f)与f_h3(x_f)为由式中的-f_λ23g₂/(2c)带入式(7)得到的结果。求解方程即可获得λ₁₀、λ₂₀和λ₃₀，如下所示，

$\left[\begin{array}{c}{\mathrm{\λ}}_{10}\\ {\mathrm{\λ}}_{20}\\ {\mathrm{\λ}}_{30}\end{array}\right]=F^{-1}\left[\begin{array}{c}{V}_b\left(x_f\right)-V_{object}\\ {\mathrm{\Δ\γ}}_f\\ {\mathrm{\Δh}}_f\end{array}\right]---\left(21\right)$

其中，F为影响函数矩阵，表达式如下。

$F=\left[\begin{array}{ccc}{f}_{V1}\left(x_f\right)-f_{\mathrm{\λ}11}\left(x_f\right)/2K_V& f_{V2}\left(x_f\right)& f_{V3}\left(x_f\right)\\ f_{\mathrm{\γ}1}\left(x_f\right)& f_{\mathrm{\γ}2}\left(x_f\right)& f_{\mathrm{\γ}3}\left(x_f\right)\\ f_{h1}\left(x_f\right)& f_{h2}\left(x_f\right)& f_{h3}\left(x_f\right)\end{array}\right]$

从而可得修正后的控制规律为，

$u_b^{new}=\mathrm{\Δ}u+u_b^{old}=-\[\begin{array}{ccc}\frac{f_{\mathrm{\λ}11}{g}_1+f_{\mathrm{\λ}21}{g}_2}{2c}& \frac{f_{\mathrm{\λ}22}{g}_2}{2c}& \frac{f_{\mathrm{\λ}23}{g}_2}{2c}\end{array}\]F^{-1}\left[\begin{array}{c}{V}_b\left(x_f\right)-V_{object}\\ {\mathrm{\Δ\γ}}_f\\ {\mathrm{\Δh}}_f\end{array}\right]---\left(22\right)$

上式中，Δu为由式(19)给出的最优控制修正量；为经过过程约束修正的初始标准控制量，在步骤2中给出；为最终得到的修正控制。式(22)即为基于广义标控脱靶量概念的多过程约束和多终端约束末制导律。在下一个制导周期中，新的重复步骤1至5。在制导过程中，该制导律将根据设定的更新频率以及飞行器当前的状态，不断生成新的参考弹道，保证飞行器一直沿着参考弹道飞行直至命中目标(如图2所示)。

相较于传统的制导律，本发明制导律不仅仅考虑了脱靶量、碰撞角等传统约束的影响，还闭环考虑了终端速度脱靶量以及控制能力约束的影响。特别适用于高超声速巡航导弹末制导或者再入飞行器下压段的制导问题。同时，该制导律还可以根据实际的制导需求减少或者增加相应的标控脱靶量个数，从而应对不同的需求。

为了进一步说明该制导方法的原理及优势，下面结合仿真算例来具体说明。

整个制导过程如图1所示，飞行器在飞行的过程中采用跟踪标准控制指令的方式，而积累的偏差则由按照指定周期更新标准控制来修正。标准控制的更新步骤如图1所示，具体实施如下。

步骤1：加载存储的标准控制及飞行器状态初值。在制导初始时刻，标准控制可能存在较大的偏差，但随着制导迭代修正的进行，标准控制会迅速收敛到所需要的最优控制上。

步骤2：弹道积分预测。这部分是利用标准控制以及飞行器的动力学模型，进行积分弹道预测，得到所需要的标控弹道。

步骤3：过程约束修正。这一部分主要讨论控制过程约束的处理办法，提出了两个处理过程约束的简单函数：过程约束修正的标准控制函数以及控制修正的影响函数。

步骤4：标控弹道一阶变分模型求解。这一步是利用上一步得到的标控弹道，计算在标控弹道附近的一阶变分模型，包括线化模型系数以及标控脱靶量。

$\left[\begin{array}{c}\frac{d\mathrm{\Δ}V}{dx}\\ \frac{d\mathrm{\Δ}\mathrm{\γ}}{dx}\\ \frac{d\mathrm{\Δ}h}{dx}\end{array}\right]=\left[\begin{array}{ccc}{f}_1& f_2& f_3\\ 0& f_4& 0\\ 0& f_5& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}V\\ \mathrm{\Δ}\mathrm{\γ}\\ \mathrm{\Δ}h\end{array}\right]+\left[\begin{array}{c}{g}_1\mathrm{\Δ}u\\ g_2\mathrm{\Δ}u\\ 0\end{array}\right]$

步骤5：更新标准控制。这一部分将采用线性最优控制的理论对上述一阶变分模型进行优化求解，获得最优的控制修正，从而完成标准控制的修正。

$u_{new}=\mathrm{\Δ}u+u_b=-\[\begin{array}{ccc}\frac{f_{\mathrm{\λ}11}{g}_1+f_{\mathrm{\λ}21}{g}_2}{2}& \frac{f_{\mathrm{\λ}22}{g}_2}{2}& \frac{f_{\mathrm{\λ}23}{g}_2}{2}\end{array}\]F^{-1}\left[\begin{array}{c}\mathrm{\Δ}V\left(x_f\right)-f_{V0}\left(x_f\right)\\ \mathrm{\Δ}\mathrm{\γ}\left(x_f\right)-f_{\mathrm{\γ}0}\left(x_f\right)\\ \mathrm{\Δ}h\left(x_f\right)-f_{h0}\left(x_f\right)\end{array}\right]$

相较于传统的末制导方法，该制导方法不仅可以满足脱靶量、落角、落点攻角等传统的终端约束的要求，还能够使强非线性状态变量V收敛到所需要的值，同时还能闭环的处理过程约束。

实施例：

为了验证上述制导律的效果，这里选用国外某飞行器模型进行数值仿真效验。飞行器的仿真参数如下：初始速度V₀＝1900m/s、初始高度h₀＝18.5km、初始弹道倾角γ₀＝0°、终端高度h_f＝0m、终端弹道倾角γ_f＝-80°、终端速度为V_f＝1000m/s、纵向飞行距离50km、攻角范围-8°≤α≤6°、最大法向过载n≤13。

(1)修正算法收敛性验证

为了使基于广义梯度法的弹道规划算法能够用于在线弹道优化，首先需要检测算法的收敛性。图3给出了弹道随着迭代的进行逐渐收敛的过程示意图。由图可见初始弹道偏差很大，经过三次迭代修正，即获得了最优解。由此可知，基于广义梯度法的标控脱靶量修正算法能够很快的收敛到最优解。

(2)无拉偏情况下的制导仿真

无拉偏情况下仿真得到的结果脱靶量为0.045m，落角误差为0.0059°，末速度误差为2.66m/s。得到的结果如图4所示。由图4(d)可知，控制规律非常光滑，并且弹道也非常光滑，这有助于在弹上的实现；由图4(e)可知最大法向过载小于限制；由图4(f)可以看出，迭代耗时分布非常有规律，散布非常小，这说明该弹道规划算法稳定性好。每一次更新耗时均小于0.1s，并且随着迭代的进行呈逐渐下降趋势，这是由于该算法的每一次迭代相当于进行一次微分方程求解，越接近弹道末段，所需的积分步数越少，因而计算速度越快。

(3)风拉偏的影响

为了进一步分析制导律的性能，现考虑风拉偏的影响。由标准大气模型可得不同高度情况下最大风速的分布如图5所示。

风拉偏的结果如表格1和图6所示，可以看出，在最大风拉偏的影响下，制导律仍然能够使得飞行器准确命中目标，同时正碰撞角的精度。但是风拉偏的情况下会对落速造成一定的影响，会使得落速误差有所增大，这是由于速度的调节属于长周期，仅通过弹道调节无法完全消除速度脱靶量。

表格1风拉偏对制导精度的影响

弹道序列	最大风速拉偏	脱靶量(m)	碰撞角误差(deg)	落速误差(m/s)
					1	顺风	0.011	-0.0014	5.42
2	无风	0.001	-0.0005	-2.81
					3	逆风	0.02	0.0079	-8.43

Claims (1)

1.一种基于广义标控脱靶量概念的多过程约束和多终端约束末制导律，包括以下几个步骤：

步骤1：加载存储的标准控制及飞行器状态初值；

标准控制为飞行器事先设定的控制指令规律，设为u_b(x)，飞行器状态初值则是当前时刻飞行器的飞行状态，包括初始速度V₀，初始弹道倾角γ₀，初始高度h₀以及初始纵程x₀；

步骤2：弹道积分预测；

建立高超声速飞行器末制导段动力学模型，飞行器的运动方程为，

$\stackrel{\·}{V}=-gsin\mathrm{\γ}-\frac{D}{m}$

$\stackrel{\·}{\mathrm{\γ}}=-\frac{gcos\mathrm{\γ}}{V}+\frac{L}{mV}---\left(1\right)$

$\stackrel{\·}{h}=Vsin\mathrm{\γ}$

$\stackrel{\·}{x}=Vcos\mathrm{\γ}$

式中，V为飞行器相对地球速度，γ为弹道倾角，h为高度，x为纵程，和分别表示飞行器相对地球速度、弹道倾角、高度和纵程对时间的导数；L和D分别为升力和阻力，其表达式分别为L＝0.5ρV²S_refC_L和D＝0.5ρV²S_refC_D；其中ρ为大气密度，S_ref为飞行器的气动参考面积；C_L和C_D分别为升力系数和阻力系数，g为当地重力加速度，m为飞行器质量，将所有状态变量转为对x进行微分，运动方程降维变为：

$\frac{dV}{dx}=\frac{-g}{V}tan\mathrm{\γ}-\frac{D}{mVcos\mathrm{\γ}}$

$\frac{d\mathrm{\γ}}{dx}=\frac{u}{cos\mathrm{\γ}}---\left(2\right)$

$\frac{dh}{dx}=tan\mathrm{\γ}$

其中，u＝-gcosγ/V²+L/(mV²)，设为控制变量；

步骤3：过程约束修正

在步骤2与步骤3同时进行；

过程约束对标准控制及控制修正的公式如下；

$u_b^{old}\left(x\right)=\left\{\begin{array}{cc}\mathrm{\α}\left(u_b^{old}\right)={\mathrm{\α}}_{\min }& \mathrm{\α}={\mathrm{\α}}_{\min }\\ \mathrm{\α}\left(u_b^{old}\right)={\mathrm{\α}}_{\max }& \mathrm{\α}={\mathrm{\α}}_{max}\\ n_y\left(u_b^{old}\right)=n_{y\max }& n_y=n_{y\max }\\ u_b\left(x\right)& {\mathrm{\α}}_{\min }\≤\mathrm{\α}\≤{\mathrm{\α}}_{\max }\∪n_y\≤n_{ymax}\end{array}\right.---\left(3\right)$

其中，α表示攻角，α_min为最小攻角；α_max为最大攻角，n_y为法向过载，n_ymax为最大法向过载，C_limit为过程约束修正因子，为一足够大的常数；c(x)为过程约束影响下的控制修正量权重因子；经过过程约束修正的标准控制；u_b(x)为弹上存储的标准控制量；表示选择使得α＝α_min时的控制量作为当前控制量；表示选择使得α＝α_max时的控制量作为当前控制量；表示选择使得n_y＝n_ymax时的控制量作为当前控制量；

利用标准控制u_b(x)，根据式(2)以及式(3)，获得过程约束影响下的控制修正量权重因子c(x)；经过过程约束修正的标准控制标准速度曲线V_b(x)、标准弹道倾角曲线γ_b(x)，标准高度曲线h_b(x)，以及标准的阻力曲线D_b(x)；

步骤4：标控弹道一阶变分模型求解；

在标准弹道附近求解一阶变分可得：

$\left[\begin{array}{c}\frac{d\mathrm{\Δ}V}{dx}\\ \frac{d\mathrm{\Δ}\mathrm{\γ}}{dx}\\ \frac{d\mathrm{\Δ}h}{dx}\end{array}\right]=\left[\begin{array}{ccc}{f}_1& f_2& f_3\\ 0& f_4& 0\\ 0& f_5& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}V\\ \mathrm{\Δ}\mathrm{\γ}\\ \mathrm{\Δ}h\end{array}\right]+\left[\begin{array}{c}{g}_1\\ g_2\\ 0\end{array}\right]\mathrm{\Δ}u---\left(4\right)$

其中，ΔV、Δγ和Δh分别为在标控弹道附近的速度、弹道倾角和高度的微小变化量；Δu为控制量的小增量；

$\begin{array}{cc}{f}_1=\frac{g}{V_b^2}{\mathrm{tan\γ}}_b-\frac{D_b}{{\mathrm{mV}}_b^2{\mathrm{cos\γ}}_b}& f_4=-u_b^{old}\frac{{\mathrm{sin\γ}}_b}{{\cos }^2{\mathrm{\γ}}_b}\end{array}$

$\begin{array}{cc}{f}_2=-\frac{g}{V_b{\cos }^2{\mathrm{\γ}}_b}-\frac{D_b{\mathrm{sin\γ}}_b}{{\mathrm{mV}}_b{\cos }^2{\mathrm{\γ}}_b}& f_5=\frac{1}{{\cos }^2{\mathrm{\γ}}_b}\end{array}$

$\begin{array}{cc}{f}_3=\frac{D_b}{{\mathrm{mV}}_b{\mathrm{cos\γ}}_b}{\mathrm{\β}}_h& g_1=-\frac{S_{ref}{V}_b}{2m{\mathrm{cos\γ}}_b}\frac{\∂C_D}{\∂C_L}\end{array}$

$g_2=\frac{1}{{\mathrm{cos\γ}}_b}$

其中，β_h为指数大气模型常数，为阻力系数对升力系数的偏导数，满足为阻力系数对攻角的变化率，为升力系数对攻角的变化率；

定义纵向弹道优化的目标函数为：

$\min J=K_V{(V_f-V_{object})}^2+{\∫}_{x_0}^{x_f}R\left(x\right)u^{2}dx---\left(5\right)$

其中，K_V为末速度约束松弛因子；R(x)为控制权重因子；V_f为末速度大小；V_object目标末速度大小；对目标函数进行线性化处理，可得，

$\min J=K_V{({\mathrm{\ΔV}}_f+V_b(x_f)-V_{object})}^2+{\∫}_{x_0}^{x_f}2{\mathrm{Ru}}_b^{old}\mathrm{\Δ}u+{\mathrm{R\Δu}}^{2}dx---\left(6\right)$

其中x_f为目标所在处的纵程，ΔV_f为摄动弹道的终端速度增量；为了使得标控弹道收敛到目标值，定义如下边界条件：

$\left.\begin{array}{c}\mathrm{\Δ}h\left(x_0\right)=0,\mathrm{\Δ}V\left(x_0\right)=0,\mathrm{\Δ}\mathrm{\γ}\left(x_0\right)=0\\ {\mathrm{\Δh}}_f=h_{object}-h_b\left(x_f\right),{\mathrm{\Δ\γ}}_f={\mathrm{\γ}}_{object}-{\mathrm{\γ}}_b\left(x_f\right)\end{array}\right\}---\left(7\right)$

其中，Δh(x₀)、ΔV(x₀)和Δγ(x₀)分别为高度、速度和弹道倾角在弹道起始处的摄动值；Δh_f和Δγ_f分别为高度、弹道倾角在弹道终点处的摄动值；h_object和γ_object分别为目标高度和目标弹道倾角；从而得到了纵平面内的标控弹道修正模型，包括末状态偏差Δh_f、Δγ_f，线性化模型系数矩阵，如下所示：

$\left[\begin{array}{ccc}\frac{g}{V_b^2}{\mathrm{tan\γ}}_b-\frac{D_b}{{\mathrm{mV}}_b^2\cos \mathrm{\γ}}& \frac{g}{V_b{\cos }^2{\mathrm{\γ}}_b}-\frac{D_b{\mathrm{tan\γ}}_b}{{\mathrm{mV}}_b{\cos }^2{\mathrm{\γ}}_b}& \frac{D_b}{{\mathrm{mV}}_b{\mathrm{cos\γ}}_b}{\mathrm{\β}}_h\\ 0& -u_b^{old}\frac{{\mathrm{sin\γ}}_b}{{\cos }^2{\mathrm{\γ}}_b}& 0\\ 0& \frac{1}{{\cos }^2{\mathrm{\γ}}_b}& 0\end{array}\right]\left[\begin{array}{c}-\frac{S_{ref}{V}_b}{2{\mathrm{mcos\γ}}_b}\frac{\∂C_D}{\∂C_L}\\ \frac{1}{{\mathrm{cos\γ}}_b}\\ 0\end{array}\right]$

步骤5：更新标准控制求解；

标控弹道修正模型的哈米尔顿函数为：

H＝λ₁(f₁ΔV+f₂Δγ+f₃Δh+g₁Δu)+λ₂(f₄Δγ+g₂Δu)+λ₃f₅Δγ+2Ru_bΔu+RΔu²(8)

其中，λ₁、λ₂和λ₃分别为ΔV、Δγ和Δh对应的协态变量；为将H对状态变量ΔV、Δγ和Δh求偏导数，得到协态方程：

${\stackrel{\·}{\mathrm{\λ}}}_1=-\frac{\∂H}{\∂\mathrm{\Δ}V}=-{\mathrm{\λ}}_1{f}_1---\left(9\right)$

${\stackrel{\·}{\mathrm{\λ}}}_2=-\frac{\∂H}{\∂\mathrm{\Δ}q}=-({\mathrm{\λ}}_1{f}_2+{\mathrm{\λ}}_2{f}_4+{\mathrm{\λ}}_3{f}_5)---\left(10\right)$

${\stackrel{\·}{\mathrm{\λ}}}_3=-\frac{\∂H}{\∂\mathrm{\Δ}h}=-{\mathrm{\λ}}_1{f}_3---\left(11\right)$

迈耶函数为：

Φ＝K_V(ΔV(x_f)+V_b(x_f)-V_object)²

+ν₁(Δγ(x_f)-Δγ_f)+ν₂(Δh(x_f)-Δh_f)

其中，x_f为弹道终点；ν₁和ν₂分别为拉格朗日因子；Δγ(x_f)、Δh(x_f)为和ΔV(x_f)由标控弹道一阶变分模型求解获得的终端弹道倾角、高度和速度摄动量；Δγ_f和Δh_f为标控脱靶量；

由式(9)可得：

λ₁＝λ₁₀f_λ11(12)

其中，λ₁₀为λ₁的协态初值；f_λ11为积分表达式，如下，

$f_{\mathrm{\λ}11}=\exp (-{\∫}_0^x{f}_{1}dx)$

由式(11)和式(12)可得：

${\mathrm{\λ}}_3={\mathrm{\λ}}_{30}-{\mathrm{\λ}}_{10}{\∫}_0^x{f}_{\mathrm{\λ}11}{f}_{3}dx---\left(13\right)$

其中，λ₃₀为λ₃的协态初值；

将式(12)和式(13)带入式(10)可得：

λ₂＝λ₁₀f_λ21+λ₂₀f_λ22+λ₃₀f_λ23(14)

其中，λ₂₀为λ₂的协态初值；f_λ21、f_λ22和f_λ23均为积分表达式，

$f_{\mathrm{\λ}21}=\exp (-{\∫}_0^x{f}_{4}dx){\∫}_0^x(-(f_2{f}_{\mathrm{\λ}11}-f_5{\∫}_0^x{f}_{\mathrm{\λ}11}{f}_{3}dx\left)\exp \right({\∫}_0^x{f}_{4}dx\left)\right)dx$

$f_{\mathrm{\λ}22}=\exp (-{\∫}_0^x{f}_{4}dx)$

$f_{\mathrm{\λ}23}=\exp (-{\∫}_0^x{f}_{4}dx){\∫}_0^x(-f_5\exp ({\∫}_0^x{f}_{4}dx\left)\right)dx$

最优轨线满足：

$\frac{\∂H}{\∂\mathrm{\Δ}u}=2\mathrm{\Δ}u+{\mathrm{\λ}}_2{g}_2+{\mathrm{\λ}}_1{g}_1+2u_b^{old}=0---\left(15\right)$

将式(14)代入式(15)，并考虑过程约束修正，可得：

$\begin{array}{c}\mathrm{\Δ}u=-\frac{\left({\mathrm{\λ}}_{10}{f}_{\mathrm{\λ}21}+{\mathrm{\λ}}_{20}{f}_{\mathrm{\λ}22}+{\mathrm{\λ}}_{30}{f}_{\mathrm{\λ}23}\right)g_2+{\mathrm{\λ}}_{10}{f}_{\mathrm{\λ}11}{g}_1+2u_b^{old}}{2}\\ =-{\mathrm{\λ}}_{10}\frac{f_{\mathrm{\λ}11}{g}_1+f_{\mathrm{\λ}21}{g}_2}{2c}-{\mathrm{\λ}}_{20}\frac{f_{\mathrm{\λ}22}{g}_2}{2c}-{\mathrm{\λ}}_{30}\frac{f_{\mathrm{\λ}23}{g}_2}{2c}-\frac{u_b^{old}}{c}\end{array}---\left(16\right)$

其中，c即为c(x)，由式(3)确定；将式(16)带入式(4)，积分后带入边界条件，可得关于λ₁₀、λ₂₀和λ₃₀的三元一次方程组：

$\left\{\begin{array}{c}{V}_b\left(x_f\right)-V_{object}=f_{V0}\left(x_f\right)+{\mathrm{\λ}}_{10}\left(f_{V1}\right(x_f)-f_{\mathrm{\λ}11}(x_f)/2K_V)\\ +{\mathrm{\λ}}_{20}{f}_{V2}\left(x_f\right)+{\mathrm{\λ}}_{30}{f}_{V3}\left(x_f\right)\\ {\mathrm{\Δ\γ}}_f=f_{\mathrm{\γ}0}\left(x_f\right)+{\mathrm{\λ}}_{10}{f}_{\mathrm{\γ}1}\left(x_f\right)+{\mathrm{\λ}}_{20}{f}_{\mathrm{\γ}2}\left(x_f\right)+{\mathrm{\λ}}_{30}{f}_{\mathrm{\γ}3}\left(x_f\right)\\ {\mathrm{\Δh}}_f=f_{h0}\left(x_f\right)+{\mathrm{\λ}}_{10}{f}_{h1}\left(x_f\right)+{\mathrm{\λ}}_{20}{f}_{h2}\left(x_f\right)+{\mathrm{\λ}}_{30}{f}_{h3}\left(x_f\right)\end{array}\right.---\left(17\right)$

其中，f_V0(x_f)、f_γ0(x_f)与f_h0(x_f)为由式中的带入式(4)得到的结果；f_V1(x_f)、f_γ1(x_f)与f_h1(x_f)为由式中的-(f_λ11g₁+f_λ21g₂)/(2c)带入式(4)得到的结果；f_V2(x_f)、f_γ2(x_f)与f_h2(x_f)为由式中的-f_λ22g₂/(2c)带入式(4)得到的结果；f_V3(x_f)、f_γ3(x_f)与f_h3(x_f)为由式中的-f_λ23g₂/(2c)带入式(4)得到的结果；求解方程即可获得λ₁₀、λ₂₀和λ₃₀，如下所示，

$\left[\begin{array}{c}{\mathrm{\λ}}_{10}\\ {\mathrm{\λ}}_{20}\\ {\mathrm{\λ}}_{30}\end{array}\right]=F^{-1}\left[\begin{array}{c}{V}_b\left(x_f\right)-V_{object}\\ {\mathrm{\Δ\γ}}_f\\ {\mathrm{\Δh}}_f\end{array}\right]---\left(18\right)$

其中，F为影响函数矩阵，表达式如下；

$F=\left[\begin{array}{ccc}{f}_{V1}\left(x_f\right)-f_{\mathrm{\λ}11}\left(x_f\right)/2K_V& f_{V2}\left(x_f\right)& f_{V3}\left(x_f\right)\\ f_{\mathrm{\γ}1}\left(x_f\right)& f_{\mathrm{\γ}2}\left(x_f\right)& f_{\mathrm{\γ}3}\left(x_f\right)\\ f_{h1}\left(x_f\right)& f_{h2}\left(x_f\right)& f_{h3}\left(x_f\right)\end{array}\right]$

从而可得修正后的控制规律为，

$u_b^{new}=\mathrm{\Δ}u+u_b^{old}=-\left[\begin{array}{ccc}\frac{f_{\mathrm{\λ}11}{g}_1+f_{\mathrm{\λ}21}{g}_2}{2c}& \frac{f_{\mathrm{\λ}22}{g}_2}{2c}& \frac{f_{\mathrm{\λ}23}{g}_2}{2c}\end{array}\right]F^{-1}\left[\begin{array}{c}{V}_b\left(x_f\right)-V_{object}\\ {\mathrm{\Δ\γ}}_f\\ {\mathrm{\Δh}}_f\end{array}\right]---\left(19\right)$

上式中，Δu为由式(16)给出的最优控制修正量；为经过过程约束修正的初始标准控制量，在步骤2中给出；为最终得到的修正控制；式(19)即为基于广义标控脱靶量概念的多过程约束和多终端约束末制导律；在下一个制导周期中，新的重复步骤1至5；在制导过程中，该制导律将根据设定的更新频率以及飞行器当前的状态，不断生成新的参考弹道，保证飞行器一直沿着参考弹道飞行直至命中目标。