基于线性伪谱模型预测控制的自主全射向再入制导方法与流程

本发明涉及再入制导领域，更具体的说是涉及基于线性伪谱模型预测控制的自主全射向再入制导方法。

背景技术

自从“阿波罗计划”和“航天飞机计划”发起以来，近50年来，再入制导受到了广泛的关注，出现了大量的再入制导算法。虽然特殊的算法根据任务和飞行器的类型而有所不同，但是那些再入制导算法本质上分为两大类：传统的参考跟踪算法和最近开发的预测校正算法。毫无疑问，航天飞机制导是世界上最具代表性的参考跟踪算法，并被100多个发射任务成功验证。参考跟踪算法由三部分组成：第一部分是在加热率、动态压力和负载因子限制的特定安全通道内预先规划可行的纵向参考轮廓；第二部分是研究一种跟踪算法，通过调节攻角(aoa)和倾侧角来跟踪该参考轮廓；第三部分通过反复改变倾侧角符号，设计倾侧角反转逻辑来控制横向误差。随后，许多工程师和科学家致力于提高参考跟踪算法的性能；这些努力包括简化规划方法，提高规划精度，以及研究各种跟踪规律，但是参考跟踪算法需要依赖于任务的调整，鲁棒性差、计算效率低。

因此，如何提供一种不依赖于任务的调整、鲁棒性好、计算效率高的基于线性伪谱模型预测控制的自主全射向再入制导方法是本领域技术人员亟需解决的问题。

技术实现要素：

有鉴于此，本发明提供了一种基于线性伪谱模型预测控制的自主全射向再入制导方法，该方法不仅可以处理耦合的经纬方向的运动，而且通过对向心力与科氏力的技术处理，获得了考虑地球自转的高精度解；不依赖于任务的调整、鲁棒性好、计算效率高。

为了实现上述目的，本发明提供如下技术方案：

基于线性伪谱模型预测控制的自主全射向再入制导方法，包括如下步骤：

s1初始化：设置初始、终端的计算仿真参数，通过离线弹道优化获得合理的初始控制参数和分段时间初值；

s2下落段制导：以最大攻角和零倾侧角下落，直到高度的变化率为零，进入滑翔段阶段；

s3滑翔段制导倾侧反转判断：根据当前能量状态与分段时间的关系，判断是否进行倾侧反转，当当前能量大于分段时刻的能量时，不进行反转，进入s4，当当前能量小于分段时刻的能量时，且不为最后一个反转点时，重新建立非线性参数控制问题，进入s4，当为最后一个反转点时，进入s9；

s4滑翔段制导预测弹道积分：使用当前控制参数u0与倾侧反转点，积分预测终端偏差并存储所有数据点，进入s5；

s5精度判定：终端状态偏差在误差范围内，进入s6，否则，进入步骤s8；

s6过程边界控制：判断所述的u0是否满足多种过程约束界限，当所述的u0超过边界控制umax时，所述的u0等于umax，否则，所述的u0等于u0，进入s7；

s7执行控制指令：将所述的u0作用到弹道阻尼控制方程中，获得平稳滑翔的控制指令uc，并将其作用到实际的系统中去，返回s3，并且将当前的分段时间作为倾侧反转的判断基准，当前的u0和分段时间作为下一步弹道积分的控制输入；

s8滑翔段制导控制参数和分段时间更新：利用考虑地球自转的描述经纬方向耦合运动关系的降阶动力学模型进行后续制导律的推导，分别对预测的每段积分弹道进行线性化处理，结合多段伪谱法和变分法原理，获得满足修正的控制参数和分段时间解析修正关系，进而获得控制参数和分段时间的更新，并进入s6；

s9末制导段显示制导：分别根据当前状态和终端约束条件在纵、横平面采用比例导引和多项式制导，获得能够导引飞行器安全地飞向目标且满足多种终端约束的制导指令。

优选的，在上述的基于线性伪谱模型预测控制的自主全射向再入制导方法中，所述s1具体包括：

基于圆球假设条件，并考虑地球自转，获得无动力滑翔飞行器的三自由度运动方程：

其中，ωe表示地球自转角速度，r表示为飞行器质心到地球圆心的地心距离，θ和φ分别表示经度和纬度，v表示飞行器相对于地球的速度，γ表示为飞行器相对于地球速度矢量与当地水平面的夹角，称之为弹道倾角，ψ表示为飞行器相对于地球速度矢量在当地水平面的投影与正北方向的夹角，并以顺时针旋转为正，称之为航向角；m表示飞行器的质量；g＝μ/r²表示为飞行器所受的重力加速度，其中，μ＝3.9860047×10¹⁴m³/s²是地球引力常数；σ表示飞行器沿速度方向旋转的角度，称之为倾侧角；l和d分别表示飞行器所承受的升力和阻力。

优选的，在上述的基于线性伪谱模型预测控制的自主全射向再入制导方法中，所述s8中，降阶动力学方程的建立

s81在平稳滑翔段，定义弹道倾角γ及弹道倾角相对时间的导数恒为零值，代入所述无动力滑翔飞行器的三自由度运动方程，得：

引入能量e取代时间作为自变量，其中代入上式整理得到cosγ与纵向升力l的解析关系:

其中，

s82对上式两侧求导，获得能量e对时间的导数，从而在地球自转条件下将总升阻比和纵向升阻比与经纬度、航向角建立联系，得到以能量e为自变量的降阶动力学方程：

其中，表示总升阻比，f在弹道面内的分量表示纵向升阻比；

s83将控制指令参数化，将纵向升阻比曲线设计成分为三段的参数函数，其具体表达式：

其中，ere1和ere2分别表示第一次倾侧反转和第二次倾侧反转发生时的绝对能量，u为一个常数，表示两次倾侧反转之前，纵向升阻比都按常值设计；在第二次倾侧反转之后，纵向升阻比的曲线设计为一个典型的二次曲线，并且它的终端值等于总升阻比曲线的终端值，k1和k2是二次曲线的参数，k1和k2的取值通过分段函数光滑导数连续两个条件确定；

s84带有强终端约束的降阶再入动力学模型表示为

θ(ef)＝θf,φ(ef)＝φf(7)

其中，x＝[θφψ]^t表示降阶动力学系统的状态矢量，u是控制参数；通过欧拉法或龙格库塔数值积分方法预测全局的状态量，获得预测的终端状态与实际所需的终端状态之间的终端误差表示为δxf＝x(ef)-xf，将带有强终端约束的降阶再入动力学模型表达式在预测弹道周围进行高阶泰勒展开，获得一组以状态偏差为自变量的误差传播动力学方程：

其中，x＝xref-δx,u＝uref-δu；a1，a2，和a3是函数f在不同区间对状态量的偏导；b1，b2，和是函数f在不同区间对控制量的偏导；δx是状态量的偏差；δu是控制量的偏差；

s85使用gauss伪谱法求解控制量的偏差，修正控制序列；

先推导一般的线性伪谱校正公式，将线性动力学方程在lg节点上离散为：

其中，d是微分逼近矩阵，下标k、l分别代表第k个与第l个插值点；以所述公式(9)的第一段代数约束关系，第一段的终端状态偏差δxrel通过gauss型积分公式表示为初始状态偏差和在lg节点上的状态偏差的显示函数关系如公式(10)：

其中，ω是gauss型积分的权函数，所述δxrel为关于初始状态偏差δx0和控制偏差δu的显式解析表达式，其表达形式如公式(11)所示，

其中，是一个3×3的系数矩阵，而是一个3×1的系数矩阵，二者分别表示终端状态偏差对初始状态偏差以及控制修正量的偏导数；同理得到和和根据链式法则建立起关于初始状态偏差、控制偏差以及终端状态偏差之间的函数关系如公式(12)：

终端状态写为

这里，g是一个将三个控制参数，u,ere1,ere2映射至终端状态的矢量函数；在控制变量微小改变的条件下，即假设δu,δere1,δere2都是小量，终端状态表示为

忽略高阶小量，终端状态的变分表示为对不同控制参数偏微分的线性组合

其中，

制导过程中一般不考虑航向角约束；在第一次反转之前，控制参数u和第一反转点被用于消除终端经纬度偏差；与此同时，第二个反转点在计算时被认为是一个定值；在第一次反转之后，通过改变控制参数u和第二反转点消除终端经纬度偏差；相应的控制变量的增量通过下式获得

利用求解的控制修正量更新控制参数

控制参数u与倾侧角指令σcom之间的变换关系

其中，(l/d)imu是惯导组件测得的升阻比；

s86路径约束控制策略，通过将驻点热流、动压或过载等路径约束将转化为倾侧角的幅值约束，通过限制实际倾侧角的幅值以达到对路径约束的严格控制；倾侧角的最大值如下

其中，rmin为当前速度给定情况下，通过驻点热流约束函数所确定的最小高度；其中，l为加速度计测量获得的当前升力大小，为通过惯组单元估计的当前升力大小，dr/dv和分别表示动力学方程与路径约束中高度对速度的偏导；

s87弹道阻尼抑制策略，弹道倾角γm构造负反馈抑制弹道震动，令弹道倾角的二阶导数为零。

优选的，在上述的基于线性伪谱模型预测控制的自主全射向再入制导方法中，所述s82中得到能量e对时间的导数为：

代入所述无动力滑翔飞行器的三自由度运动方程，并忽略关于速度的运动方程将所述无动力滑翔飞行器的三自由度运动方程的状态空间由六维压缩至五维：

在平衡滑翔条件下，省略上式中地心距r和弹道倾角γ的运动方程，并利用所述cosγ与纵向升力l的解析关系，获得所述降阶动力学方程。

优选的，在上述的基于线性伪谱模型预测控制的自主全射向再入制导方法中，所述s82中预先确定攻角与所述能量e的变化规律，并根据所述总升阻比f随攻角的变化规律，离线计算所述总升阻比f随能量e变化的函数。

优选的，在上述的基于线性伪谱模型预测控制的自主全射向再入制导方法中，所述s9中具体的步骤如下：设计纵、横向制导平面的指令分配算法，通过调整攻角和倾侧角实现两个制导平面的制导指令；

其中，横向比例制导平面采用比例导引：

上式中，nlateral是横向制导指令，npc是比例导引常数，是视线角速度，v是投影在横向平面的速度；

纵向多项式制导平面采用多项式制导方法，其制导指令表示如下：

其中，nvertical是纵向制导指令，sf为taem交班点与目标之间的距离，并且制导增益mg和ng分别为满足mg＞ng＞0的常值

经由上述的技术方案可知，与现有技术相比，本发明公开提供了基于线性伪谱模型预测控制的自主全射向再入制导方法，其中，降阶动力学方程的建立通过选取能量e为自变量并引入平稳滑翔条件，将六组运动方程压缩为三组方程，并充分考虑地球自转的对降阶动力学方程的影响，不仅可以处理耦合的经纬方向的运动，而且通过对向心力与科氏力的技术处理，获得了考虑地球自转的高精度解。另外，分别设计了作用于球面地球假设下的比例导引和多项式导引，在最后一次反转后，分别在纵横向方向保证各自平面内的多种终端约束，此外，该方法的适用性，有效性和制导性能将通过高敏的数值算例进行验证，该方法的鲁棒性也进行验证，仿真结果表明该算法不仅具有很快的计算效率，而且具有很高的计算精度，即使在恶劣的再入环境下也能提供稳定和鲁棒的再入制导性能。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。

图1附图为本发明的流程图；

图2附图为本发明的实施例中不同算例的地面轨迹；

图3附图为本发明的实施例中不同算例的高度变化曲线；

图4附图为本发明的实施例中不同算例的弹道倾角曲线；

图5附图为本发明的实施例中不同算例的攻角曲线；

图6附图为本发明的实施例中不同算例的倾斜角曲线；

图7附图为本发明的实施例中不同算例的cpu计算耗时曲线；

图8附图为本发明的实施例中不同算例的热流密度曲线；

图9附图为本发明的实施例中不同算例的动压曲线；

图10附图为本发明的实施例中不同算例的过载曲线。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例公开了一种基于线性伪谱模型预测控制的自主全射向再入制导方法，其中，降阶动力学方程的建立通过选取能量e为自变量并引入平稳滑翔条件，将六组运动方程压缩为三组方程，并充分考虑地球自转的对降阶动力学方程的影响，不仅可以处理耦合的经纬方向的运动，而且通过对向心力与科氏力的技术处理，获得了考虑地球自转的高精度解。

实施例1

基于线性伪谱模型预测控制的自主全射向再入制导方法，具体步骤如下：

s1初始化：设置初始再入参数、终端约束条件，并通过离线弹道优化设置合适的初始控制参数。

即初始化建模，具体如下：

s11再入动力学方程

基于圆球假设条件，并考虑地球自转，无动力滑翔飞行器的三自由度运动方程可以描述为如下形式

其中，r表示为飞行器质心到地球圆心的地心距离，θ和φ分别表示经度和纬度，v表示飞行器相对于地球的速度，γ表示为飞行器相对于地球速度矢量与当地水平面的夹角，称之为弹道倾角，ψ表示为飞行器相对于地球速度矢量在当地水平面的投影与正北方向的夹角，并以顺时针旋转为正，称之为航向角；m表示飞行器的质量；g＝μ/r²表示为飞行器所受的重力加速度，其中，μ＝3.9860047×10¹⁴m³/s²是地球引力常数；ωe表示地球自转角速度，其值等于7.2921×10-5rad/s；σ表示飞行器沿速度方向旋转的角度，称之为倾侧角；l和d分别表示飞行器所承受的升力和阻力，其表达式表示为：

其中，ρ表示当地大气密度，sref表示飞行器的特征面积，cl和cd分别表示飞行器的升力系数和阻力系数，它们只与马赫数(ma)和攻角(aoa)有关。

能量参数e被定义为

显然，随着飞行器在大气中滑行，能量e严格单调递减。因此，它常作为制导律设计中的自变量，如此不仅可以降低动力学模型的维度，而且更方便设置包括终端高度与终端速度在内的截止条件。

s12再入弹道约束

再入弹道约束包括路径约束与终端约束。

对于再入飞行器，典型的路径约束包括驻点热流密度、动压以及法向过载，它们通常可以表示为如下的形式：

q＝0.5ρv²≤qmax

其中，qmax和nmax分别为飞行器所能承受的最大热流、动压和过载，kq是用于计算飞行器表面驻点处受热率的常用参数。

对于水平着陆的再入飞行器，再入飞行段结束时，需要将飞行器交班至终端区域能量管理(taem)制导系统中，因此，这需要飞行器在taem交班时刻满足正确的终端状态约束，除此之外，还需保证飞行器在taem交班时刻的速度矢量方向要与目标视线方向的夹角尽量小，具体的终端约束条件如下：

h(ef)＝hf，v(ef)＝vf，s(ef)＝sf

γ(ef)＝0,σ(ef)＝0,

δψ＝ψ(ef)-ψlos(ef)＝0

其中，ef表示在taem交班时刻的截至能量：s表示飞行器的剩余距离；δψ表示航向误差角，它是航向角与目标视线角之间的夹角。

s13飞行器模型参数

通用航空飞行器(cav)是迄今为止最具代表性的高升阻比再入飞行器。从phillips的研究报告中获知，高升力cav的最大升阻比接近3.5。因此，采用这类飞行器用于评估本发明所提出的制导律。将cav-h的升力系数和阻力系数采用高阶的多项式函数进行拟合，获得以马赫数和攻角为输入的升力系数函数和阻力系数函数。除此之外，因为最大升阻比所在的攻角为10度附近，将攻角的变化范围拓展到至5～20度。cav-h在飞行过程中所能承受的最大热流、动压和过载分别为600w/cm2、100kpa和2g。在再入飞行过程中，将飞行器的攻角变化设计为以能量为自变量的分段函数。

其中，emid定义为飞行器最后一次倾侧反转时刻的能量值，aα1和aα2分别为多项式系数项，用于平滑过渡攻角曲线。

s2下落段制导：以最大攻角和零倾侧角下落，直到下落段结束，开始滑翔段。

初始下降段是指飞行器从初始高度下落到高度变化率为零的一段高度范围，该段结束时，飞行器所受到的气动升力已经充分建立起来并能够支撑飞行器进行滑翔飞行。在这一阶段，由于飞行器速度很高并且大气密度急剧增加，攻角和倾侧角的调整将对飞行器受到的局部最大驻点热流和之后的滑翔弹道产生巨大的影响。一般情况下，局部最大热流将发生在该段结束时刻，因此，为了保证飞行器的局部驻点热流最小且飞行器具有最大的能量调整裕度，飞行器在该段将按最大攻角和零度倾侧角下落。

αcmd＝αmax

σcmd＝0deg

如果飞行器在稠密大气层中的弹道倾角增加至0度，则结束下落段制导，开始滑翔段制导过程。

s3滑翔段制导倾侧反转判断：根据当前能量状态与分段时间的关系，判断是否进行倾侧反转，当当前能量大于分段时刻的能量时，不进行反转，进入s4，当当前能量小于分段时刻的能量时，且不为最后一个反转点时，重新建立非线性参数控制问题，进入s4，当为最后一个反转点时，进入s9；

s4滑翔段制导预测弹道积分：使用当前控制参数u0与倾侧反转点，积分预测终端偏差并存储所有数据点，进入s5；

s5精度判定：终端状态偏差在误差范围内，进入s6，否则，进入步骤s8；

s6过程边界控制：判断所述的u0是否满足多种过程约束界限，当所述的u0超过边界控制umax时，所述的u0等于umax，否则，所述的u0等于u0，进入s7；

s7执行控制指令：将所述的u0作用到弹道阻尼控制方程中，获得平稳滑翔的控制指令uc，并将其作用到实际的系统中去，返回s3，并且将当前的分段时间作为倾侧反转的判断基准，当前的u0和分段时间作为下一步弹道积分的控制输入；

s8滑翔段制导控制参数和分段时间更新：利用考虑地球自转的描述经纬方向耦合运动关系的降阶动力学模型进行后续制导律的推导，分别对预测的每段积分弹道进行线性化处理，结合多段伪谱法和变分法原理，获得满足修正的控制参数和分段时间解析修正关系，进而获得控制参数和分段时间的更新，并进入s6；

s81非线性降阶模型

考虑地球自转的描述经纬方向耦合运动关系的降阶动力学模型，该降阶动力学模型将被用于后续制导律的推导。模型的规划变量的是可以直接被观测的侧向与纵向升阻比。首先假设飞行器的滑翔段的轨迹位于与地心重合的球面上，也就是说，在平稳滑翔段，飞行器的高度不变，弹道倾角恒为零值，所以有

而后，向上式添加零阶项2ωesinψcosφcosγ-2ωesinψcosφcosγ＝0，并整理为如下形式

仔细分析上式的第三项和第四项的变化范围，定义变量

ε＝2ωesinψcosφ(1-cosγ)

因为滑翔段的弹道倾角在0.5°的范围内，最大地心距为6438245m，而且最小速度为2500m/s，所以变量ε和ξ的范围可以保守地估计为

由于变量ε的最大值远小于6e-9，而且比变量ξ小四个量级，忽略变量ε，并得到如下等式

进一步地，通过引入能量e取代时间作为自变量，利用严格单调递减的能量e将公式(1)另外五个变量统一表述为该变量的函数，状态空间进而从六维被压缩至五维。

e被定义为

其中，v²可以视为e与r的函数，即v²＝2(e+μ/r)。

显然，将e带入无动力滑翔飞行器的三自由度运动方程，整理后可以得到cosγ的表达式

其中，

需要强调的是，获得cosγ与纵向升力的解析关系对于构建经纬方向运动与纵向升阻比之间的数学关系至关重要。

进一步地，s82通过引入能量e取代时间作为自变量。在考虑地球自转的条件下，对公式两侧求导，获得能量e对时间的导数

由于公式中的第二项包含ωe的二次方，是可以忽略的小量，所以公式可以近似写为

将公式代入公式中，可以得到以能量为自变量的考虑地球自转的飞行器运动方程组，其中速度可以由能量和地心距确定。并利用拟平衡滑翔条件，省略地心距r和弹道倾角γ的表达式，一组将纵向和侧向升阻比与经纬度、航向角建立联系的降阶动力学方程可以表述为如下形式

其中，f表示总升阻比，u表示纵向升阻比。由于升阻比随攻角的变化规律可以通过飞行器的气动资料获得，因此，一旦攻角随能量的变化规律事前确定，则与之相关的升阻比曲线也可以离线计算。也就是说，f是随能量e变化的一组可以离线确定的函数。

需要强调的是，航向角表达式第一项的分母中的cosγ直接近似为1，如此可以使降阶动力学系统仅作为纵向和侧向升阻比的函数。此外，该项的符号与倾侧角符号相反。

补充两组后文线化过程中会用到的关于函数θ的偏微分公式

s83控制指令的参数化

控制指令将根据平稳滑翔弹道的特点直接进行规划处理，因为攻角是通过离线规划获得，控制指令的规划也就是指将纵向升阻比曲线离散为几个参数的函数。实际上，因为纵向升阻比和倾侧角之间具有简单的一一映射函数关系，离散纵向升阻比曲线也就是离散倾侧角曲线。在滑翔过程中，将安排两次倾侧反转以保证飞行器具有足够的能力消除积累起来的横向误差，第一次倾侧反转安排发生在飞行中段，第二次倾侧反转安排发生在离taem段不远处。因此，纵向升阻比曲线就可以设计成如下一个分为三段的参数函数，其具体表达形式如下所示：

其中，ere1和ere2分别表示第一次倾侧反转和第二次倾侧反转发生时的绝对能量，u为一个常数，表示两次倾侧反转之前，纵向升阻比都按常值设计；在第二次倾侧反转之后，纵向升阻比的曲线设计为一个典型的二次曲线，并且它的终端值等于总升阻比曲线的终端值，k1和k2是二次曲线的参数，k1和k2的取值通过分段函数光滑导数连续两个条件确定；其中解析表达式为：

因而，当u,ere1,ere2给定，纵向升阻比曲线也就确定，当初始的状态已知时，就能够通过对降阶的非线性运动方程进行积分而确定全部再入滑翔弹道。

s84带有强终端约束的降阶再入动力学模型可以一般地表示为

θ(ef)＝θf,φ(ef)＝φf

x＝[θφψ]^t表示降阶动力学系统的状态矢量，u是控制参数；通过欧拉法或龙格库塔数值积分方法预测全局的状态量，获得预测的终端状态与实际所需的终端状态之间的终端误差表示为δxf＝x(ef)-xf，将带有强终端约束的降阶再入动力学模型表达式在预测弹道周围进行高阶泰勒展开，获得一组以状态偏差为自变量的误差传播动力学方程：

其中，x＝xref-δx,u＝uref-δu；a1，a2，和a3是函数f在不同区间对状态量的偏导；b1，b2，和是函数f在不同区间对控制量的偏导；δx是状态量的偏差；δu是控制量的偏差；其中，

通过推到获得矩阵ai和bi可以表示为

其中，各分量的形式如下

s85线性伪谱模型预测控制

应用gauss伪谱法的将每一段的自变量区间映射至[-1,1]，因为用于lagrange插值多项式的正交配点位于[-1,1]区间，区间变换公式如下

将τ作为自变量，于是线性动力学系统可以表示为

之后，定义ln(τ)为n阶的拉格朗日插值多项式，τi表示n阶legendre多项式的根，也就是之前提到的lg节点。值得注意的是，这里没有解析求解n阶legendre多项式根的方法，必需通过数值算法才能获得lg节点。于是，就可以将线性动力学方程中的状态量表示为一组由lg节点为支撑点所形成的拉格朗日插值多项式基来线性组合。

根据拉格朗日插值多项式的基本理论，ln(τ)满足以下的基本性质

并且

δxⁿ(τl)＝δx(τl)

之后，lg节点上的状态导数值就能通过微分逼近矩阵转化为一组代数方程，并且建立起与各个lg节点上状态量之间的显示关系。其中，微分逼近矩阵d是一个n×n+1的矩阵，通过对拉格朗日插值多项式的各个元素分别求导获得，d矩阵各元素的具体表达形式以及状态量与其导数的代数关系如下所示：

其中，

假定每一段的状态量在lg节点上表示为如下形式

将方程带入到方程中，那么，微分动力学约束不仅能够转化为一组代数约束，并且能够表示为lg节点上各个状态的函数。

其中，k＝1,2,...，n。以方程中的第一段代数约束关系式为列，重新组合微分逼近矩阵d，一个以状态矢量形式表示且更为简单的关系式如下所示：

其中，d1和d2:n的具体表达形式如下：

其中，s表示为状态变量的个数，矩阵a1和b1分别表示如下：

再次对方程进行重新组合，在lg节点上的各个状态量可以解析地表示为：

因为gauss正交配点不包括终点，终端值可以通过gauss积分获得，以公式中的第一个公式为例

gauss积分具有最高的数值精度，可以精确拟合不高于2n-1阶的多项式。所以有

显然与gauss配点τi相匹配的gauss积分权重ωi是lagrange插值多项式的积分，可以表示为

其中，是n阶legendre多项式的导数。因此，第一段的终端状态偏差可以表示为

整理上式，终端状态量的变分可以表示为矢量形式

其中，1是单位列向量，w1由高斯积分权重组成

通过将上式代入终端状态量的变分的矢量形式中，可以将第一段的终端状态偏差表示为初始状态偏差δx0和控制参数的修正量δu的函数

其中，是一个3×3的系数矩阵，而是一个3×1的系数矩阵，二者分别表示终端状态偏差对初始状态偏差以及控制修正量的偏导数。

公式l^x，l^u可以具体展开为

类似地，可以给出第二段的偏导矩阵

同理，给出第三段的偏导矩阵

然后，将上一段的终端状态偏差作为下一段的初始偏差，并以递归的方式，可以将全段的终端状态偏差表示为

显然，上式建立了终端状态偏差与初始状态偏差以及控制修正量之间的数学关系，但是，有必要强调的是，公式中只考虑了一个控制参数。但是这个表达式中包含了初始状态偏差对终端误差的影响，后文将从初始状态偏差的角度推导倾侧角反转时刻对终端状态的影响。

终端状态可以一般写为

这里g是一个将三个控制参数，u,ere1,ere2映射至终端状态的矢量函数；在控制变量微小改变的条件下，即假设δu,δere1,δere2都是小量，终端状态表示为

忽略高阶小量，终端状态的变分可以表示为对不同控制参数偏微分的线性组合

其中，

以下逐一推导各项偏导数的解析表达式。

(a)

表示终端状态变量对控制变量u的偏导数，而这部分关系已经在公式中表示，所以有

(b)

表示终端状态偏差相对第一个反转点的偏导数。必须指出的是，从传动的动力学方程中，难以建立二者关联。通过分析反转点的改变对初始反转点处状态量的影响，进而再分析该状态量的改变是如何体现在终端状态偏差中。我们获得了的表达式。

首先，我们考虑第一个反转点的变分对初始反转点处状态量的影响，由变分法则

δx(ere1)＝[f1(u,ere1,ere2)-f2(u,ere1,ere2)]δere1

其中，f1(u,ere1,ere2)δere1是状态量在反转点增加δere1后的变分，-f2(u,ere1,ere2)δere1是状态量在反转点缩减δere1后的变分。二者的合表示受反转点变分的影响，初始反转点处状态量的改变量。一旦确定了初始改变量，我们可以通过假定其他控制参数不变的线性动力学方程计算终端状态改变量，即

与公式结构相似的，不难得到反转点的改变对终端状态的影响

显然，可以被定义为

(c)

表示终端状态偏差相对第二个反转点的偏导数。推导过程同上类似，不同点在于只需要考虑最后一段，所以的表达式为

当公式中的各项参数得以确定后，相应的用以消除终端误差的控制参数的修正量可以表示为

其中，dxf通过积分前文所述的降阶动力学模型获得。

制导过程中一般不考虑航向角约束；在第一次反转之前，控制参数u和第一反转点被用于消除终端经纬度偏差；与此同时，第二个反转点在计算时被认为是一个定值；在第一次反转之后，通过改变控制参数u和第二反转点消除终端经纬度偏差；相应的控制变量的增量通过下式获得

利用求解的控制修正量更新控制参数

有必要说明的是，对于再入飞行器，倾侧角是一个基本的控制变量，因此，给出控制参数u与倾侧角指令σcom之间的变换关系

其中，(l/d)imu是惯导组件测得的升阻比。

s86路径约束控制策略

通过将驻点热流、动压或过载等路径约束将转化为倾侧角的幅值约束，通过限制实际倾侧角的幅值以达到对路径约束的严格控制。倾侧角的最大值如下

其中，rmin为当前速度给定情况下，通过驻点热流约束函数所确定的最小高度；其中，l为加速度计测量获得的当前升力大小，为通过惯组单元估计的当前升力大小，dr/dv和分别表示动力学方程与路径约束中高度对速度的偏导；

s87弹道阻尼抑制策略弹道倾角γm构造负反馈抑制弹道震动，令弹道倾角的二阶导数为零，并进一步求解了γm的表达式

其中，

其中，上标为“*”的项表示由攻角指令所确定的升力与阻力或升力系数与阻力系数，之后，弹道跳跃的抑制就通过在纵向升力方向加入该特殊弹道倾角的负反馈，并保持横向升力不变而实现，实际作用的倾侧角和升力系数就通过式确定。

其中，下标“1”表示由制导律给出的倾侧角和攻角指令，下标为“2”的项表示实际需要作用的倾侧角和以实际作用的攻角所确定的升力系数，kγ表示为确定弹道抑制效果的负反馈增益。显然，倾侧角指令可以解析地表示为

虽然，滑翔段的攻角指令受非线性气动系数的影响不能给出解析表达式，但是，指令升力系数可以表示为

而后，newton-raphson迭代被用于求解指令攻角

当使用上一步获得攻角作为初始值代入计算，只需消耗很少的时间(几毫秒)，仅仅几步计算就能获得满足精度要求的解，完全适用于在线应用。

s9末段调整段制导：采用比例导引与多项式制导获得满足终端约束的控制指令，直至导引飞行器进入能量管理段。至此，再入段制导结束。

横向比例制导律

比例导引作为最著名的且最具鲁棒性的制导律，将在这一章节内被推广到末端调整段的横向导引中，从而导引飞行器在横向平面飞向目标，并且保证终端的横向加速度收敛到零，正如比例导引的显示表达式所示，比例导引的制导加速度指令是通过视线角变化率、投影在横向平面的速度以及比例导引常数所确定。

nlateral是横向制导指令，npc是比例导引常数，是视线角速度，v是投影在横向平面的速度；

纵向多项式制导律

在这一节中，由kim提出的多项式制导策略将用于导引飞行器在纵向平面飞向taem交班点，并且满足高度、弹道倾角等纵向状态约束。其制导指令表示如下：

其中，其中，nvertical是纵向制导指令，sf为taem交班点与目标之间的距离，并且制导增益mg和ng分别为满足mg＞ng＞0的常值。这样，能够导引飞行器在纵向平面飞向taem交班点，并且满足终端高度和弹道倾角约束。

实施例2

在本实施例中，十组不同目标点的算例将用于评估本发明所提方法的性能。这十组算例的初始条件除航向角外保持一致，而且终端与能量管理段的交班条件也相同，分别如表1和表2所示。十组算例的终端位置约束如表3所示。需要说明的是，当剩余射程满足截止条件时，仿真结束。所有仿真结果如图2至图10所示。

表1初始再入条件

表2终端能量管理段的交班条件

表3不同算例的终端位置约束

图2展示了十组算例在地面的轨迹。明显地，所有曲线从相同的发射点出发并截止于各自的目标点。所以，我们有理由认为本发明所提出的算法能够适应全球分布的目标点，其中经度范围达到20°至340°，纬度范围达到-60°至45°。由于考虑了地球自转，所以虽然目标点对称分布，但是向东飞行与向西飞行所形成的地面轨迹还是略有不同。图3中描绘了高度变化曲线，显然所有曲线都由下落段、滑翔段、末段调整段三部分组成，并最终收敛到能量管理段的交班条件。此外，由于地球自转的影响，不同算例的下降段截至高度不同。图4展示了弹道倾角变化曲线。可以看出，弹道倾角在飞行全程中都是小量，特别是滑翔段，最大值不超过0.2°。这幅图也反映出由于考虑了自动驾驶仪的一阶延迟效应，倾侧角的反转过程会导致弹道产生尖点。

图5展示了攻角的变化曲线。显然，下落段保持最大攻角飞行；滑翔段则参考事先设定的攻角曲线，当然也会略有调整进而抑制弹道振动，这种微调在倾侧角反转时会更明显；末段调整段，攻角直接由制导指令生成。倾侧角曲线如图6所示，显然，滑翔段倾侧角模值基本呈直线变化，并在能量管理段的交班点收敛至零。尽管事先没有根据特定任务修改制导律，但我们可以看到不同飞行过程的倾侧反转时刻是不同的，这些计算都可以由本发明所提算法在线完成。

图7展示了每个制导周期中用于产生制导指令的cpu计算耗时。其中，最大耗时不超过0.02s，而且普遍低于0.01s。因此，本文所提算法具有令人满意的计算效率，完全适用于在线应用。图8至图10分别展示了飞行全程中的热流、动压与过载，显然所有路径约束都得到严格满足。

表4总结了十组算例的终端状态偏差。所有的高度偏差都不超过0.1m，所有的速度偏差都不超过1m/s，而且弹道倾角和航向角的偏差也都不超过0.01°。相比之下，倾侧角偏差略大，但完全在可接受范围内。综上，本发明所提算法适用于各种各样的再入条件，并且都有非常好的表现。

表4十组算例的终端状态偏差

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。