技术领域
[0001] 本
发明属于自适应制导方法，具体涉及一种多级助推
滑翔飞行器主动段在线自适应制导算法。
背景技术
[0002] 多级助推滑翔飞行器在大气层内的主动段飞行是一个重要阶段，实际飞行过程中通过箭载计算机中的制导律，将飞行器从当前的飞行状态导向终端目标状态。制导律的主要功能是实时根据飞行器当前的飞行状态输出程序
角指令或其他指令信息。目前，多级助推滑翔飞行器的主动段制导任务需求被进一步提高，需要满足终端约束(主动段到滑翔段的交班条件，如速度约束、高度约束)和过程约束(如过载约束、
攻角约束)。传统的多级助推滑翔飞行器在主动段飞行时通常采用摄动制导方法作为制导律，即为了满足各种约束基于标称条件下的额定参数进行离线优化设计标准轨迹(将级间
滑行时间、程序角指令固定，不能在线调整)，而将实际飞行中各种偏差(
发动机推
力、比冲、
起飞质量以及
风干扰等)均认为是小量，基于线性摄动制导理论，在线进行偏差量的导引修正。该种摄动制导模式过于依赖于标准轨迹，并且发射前需要进行大量的发射诸元计算(导引参数，轨迹参数)，导致飞行任务不能快速装订，另外一旦实际飞行过程中的偏差(固体
火箭发动机性能受
温度影响较大)超出了摄动制导理论所允许的范围，该方法会产生较大的制导
精度误差。
[0003] 此外，现有的自适应制导方法如“
迭代制导”和“闭路制导”方法虽然能够在线计算满足多终端约束的制导指令，但均适用于
真空环境下的制导飞行，无法保障飞行器在大气层内飞行时的飞行过程约束。最后，对于多级助推-滑行-助推模式的飞行器，级间滑翔时间的自适应，一方面能够有效的抑制飞行包络的散布，另一方面能够降低制导指令的修正角幅值，能够为满足飞行过程约束及多终端交班约束提供有利条件。上述方法并未针对滑翔时间的自适应问题进行研究，通常采用离线设计好的飞行时序进行级间滑翔飞行，因此严重影响了系统的快速响应和实战化能力。
发明内容
[0004] 本发明的目的是解决上述问题，提供一种满足新特性需求下多级助推滑翔飞行器的主动段制导任务需求的多级助推滑翔飞行器主动段在线自适应制导算法。
[0005] 为实现上述目的，本发明提供如下技术方案：
[0006] 一种多级助推滑翔飞行器主动段在线自适应制导算法，首先，由于大气层内主动段飞行环境复杂，其动力学模型涉及推力曲线模型、J2引力模型以及
气动力模型等严重非线性问题，无法直接获取解析解，需要通过数值积分的方式求解微分方程组，因此采用全新的滑行时间在线数值预测迭代算法，来求解滑行点火时间。然后，在动力飞行段，以预测的基准飞行轨迹为
基础上采用模型预测静态规划算法，根据最小修正角的二次型指标计算出同时满足终端多约束条件和飞行过程约束的制导指令。最后，根据导航器件测量得到的实时状态信息，采用由模型预测静态规划算法构成的闭环迭代制导格式，在每个制导周期输出闭环指令。
[0007] 以典型的多级“助推-滑行-助推”飞行器为研究对象，均采用耗尽关机的固体发动机，主动段飞行弹道采用低弹道模式，包含一级动力飞行段，无动力滑行段、二级动力飞行段、三级动力飞行段。
[0008] 进一步的，通过对动力学模型进行无量纲化，来确定飞行轨迹数值
预测模型，其中
位置用地球平均半径R0无量纲化，速度用 无量纲化，
加速度用g0无量纲化，时间用无量纲化，质量用初始时刻的质量m0进行无量纲化；所述无量纲化的动力学方程模型表示为：
[0009]
[0010] 其中，无量纲化的阻力加速度、升力加速度以及无量纲化的推力加速度的表达式为：
[0011]
[0012] 离心加速度为：ae＝ωe×(ωe×r)，哥氏加速度为： 这两项加速度在发射系的具体表达式为：
[0013]
[0014] 对于“助推—滑行—助推”的飞行模式中的推力项，还应该区分无动力段及动力段过程，推力项的表达式为：
[0015]
[0016] 进一步的，建立点火时间与终端高度约束的一维非线性微分方程；而此时的
俯仰角指令采用动力段全程零攻角即 的模式和
偏航角指令全程为零ψ(t)＝0，尽可能的为后续满足多约束条件的指令迭代过程预留充足的空间；因此，点火时间的求解方程和终端约束方程为：
[0017]
[0018] 进一步的，将无量纲化的飞行器运动方程式，考虑为离散形式的一般非线性系统，其状态和输出动力学如下：
[0019]
[0020] 其中，X∈Rn,U∈Rm,Y∈Rp,k＝1,2,......，N，提高一组合适的控制变量Uk,k＝1,2......，N-1，使得输出的最终值YN达到预期值 YN关于用 表示，用泰勒展开得：
[0021]
[0022] 根据小偏差摄动原理，可以将残差表示为：
[0023]
[0024] 然后，根据式可以进一步将k+1步的误差表示成：
[0025]
[0026] 其中，dXk,dUk分别是状态向量和控制向量在第k步的误差项；同样地，时间步骤N-1的状态误差dXN-1可以根据迭代步骤N-2的状态误差和控制误差进行递推；依次将dXN-2展开为dXN-3和dUN-3的递推形式，直到k＝1，则模型预测静态规划中误差的递推表达式为：
[0027]
[0028] 其中，状态转移矩阵A及灵敏度矩阵B的表达式为：
[0029]
[0030] 进一步的，为了应用MPSP方法，来实现对制导指令的过程约束并保证制导指令的平滑性，引入虚拟控制量
姿态角速率：
[0031]
[0032] 称 为虚拟控制量；定义状态向量 控制向量 引入虚拟控制量后动力学方程变为：
[0033]
[0034] 对式应用Euler积分法则进行离散化，可以得到：
[0035]
[0036] 为了得到状态转移矩阵A及灵敏度矩阵B，分别对式求偏导数，得到 以及 的符号表达式；把它们代入式即可递归地求得敏感矩阵Bj,j＝1,2,…,N-1；
[0037] 其中，链式递推偏导数的表达式依次为：
[0038]
[0039]
[0040] 为了通过最小化的控制量来满足终端多约束条件，采用二次型性能指标，表达式为：
[0041]
[0042] 其中， 是上一步的控制量的解决方案，而dUk是控制量中相应的修正量；根据最优化理论可得最小的指令修正角序列为：
[0043]
[0044] 其中， 以及
[0045]
[0046] 权值函数Rj的取值定义为：
[0047]
[0048] 进一步的，在MPSP方法的推导过程中用Euler积分法则对连续系统进行了离散化，在仿真过程中应用四阶Runge-Kutta方法进行积分。
[0049] 进一步的，所述制导的流程为：
[0050] S1：滑行段开始；
[0051] S2：在线数值预测程序，生成初始飞行轨迹并输出点火条件，进入主动段；
[0052]
[0053] dYN＝YN-Yd
[0054] U＝Up-dU
[0055] S3：用四阶龙格库塔方法对动力学方程积分，dYN＝YN-Yd；
[0056] S4：判断是否满足||dY||＜ε1，不满足时回到U＝UP-dU；满足时为RK4；
[0057] S5：判断是否满足|Z-Zf|＜ε2，不满足时回到步骤S2，满足时结束。
[0058] 与
现有技术相比，本发明的有益效果在于：
[0059] 1)射前无需进行大量的计算和诸元装订，主要由制导计算机在线完成，能够适应任务灵活、对发射响应
块的任务；
[0060] 2)飞行包络抑制能力强，飞行过程中采用滑行时间动态调整的制导方式，能够有效的降低模型不确定性、参数散差以及大气环境干扰的影响，使飞行散布管道比较集中；
[0061] 3)滑行时间数值迭代算法收敛性强可靠性高，基于重力转弯零攻角飞行假设，将滑行点火时间与终端约束建立理论关系，并通过行之有效的一维迭代求解方式保证在线应用的可靠性；
[0062] 4)具有很强的飞行过程约束和终端多约束能力，采用高精度和快收敛迭代模型预测静态规划算法，在零攻角飞行假设的基准轨迹基础上，以最优二次型指标静态搜寻满足终端多约束的制导修正角，并通过权值函数来满足飞行过程约束；
[0063] 5)闭环迭代制导格式保证高精度交班点分离条件，将模型预测静态规划算法改进为闭环递推迭代格式，根据飞行器实时导航输入信息，闭环生成实时制导指令，来抵抗飞行过程中的偏差干扰和不确定性，到达高精度的班点分离条件。
附图说明
[0064] 为了更清楚地说明本发明
实施例技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是为了更清楚地说明本发明实施例或现有技术中的技术方案，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0066] 图2为本发明的闭环多约束迭代制导流程图；
[0067] 图3为滑行时间自适应计算和固定常值下飞行高度仿真对比曲线图；
[0068] 图4为重力转弯方法和本发明方法下飞行当地倾角仿真对比曲线图；
[0069] 图5为重力转弯方法和本发明方法下飞行攻角仿真对比曲线图。
具体实施方式
[0070] 为了使本领域技术人员更好地理解本发明的技术方案能予以实施，下面结合具体实施例对本发明作进一步说明，但所举实施例只作为对本发明的说明，不作为对本发明的限定。
[0071] 1.如图1-5所示的一种多级助推滑翔飞行器主动段在线自适应制导算法。
[0072] 首先，由于大气层内主动段飞行环境复杂，其动力学模型涉及推力曲线模型、J2引力模型以及气动力模型等严重非线性问题，无法直接获取解析解，需要通过数值积分的方式求解微分方程组，因此采用全新的滑行时间在线数值预测迭代算法，来求解滑行点火时间。然后，在动力飞行段，以预测的基准飞行轨迹为基础上采用模型预测静态规划算法，根据最小修正角的二次型指标计算出同时满足终端多约束条件和飞行过程约束的制导指令。最后，根据导航器件测量得到的实时状态信息，采用由模型预测静态规划算法构成的闭环迭代制导格式，在每个制导周期输出闭环指令。
[0073] 以典型的多级“助推-滑行-助推”飞行器为研究对象，均采用耗尽关机的固体发动机，主动段飞行弹道采用低弹道模式，包含一级动力飞行段，无动力滑行段、二级动力飞行段、三级动力飞行段。
[0074] 2.1无量纲化的数值预测模型
[0075] 为了保证数值积分计算的有效性及
稳定性，对动力学模型进行无量纲化，来确定精度较高且行之有效的飞行轨迹数值预测模型，其中位置用地球平均半径R0无量纲化，速度用 无量纲化，加速度用g0无量纲化，时间用 无量纲化，质量用初始时刻的质量m0进行无量纲化。则无量纲化的动力学方程模型可以表示为：
[0076]
[0077] 其中，无量纲化的阻力加速度、升力加速度以及无量纲化的推力加速度的表达式为：
[0078]
[0079] 离心加速度为：ae＝ωe×(ωe×r)，哥氏加速度为： 这两项加速度在发射系的具体表达式为：
[0080]
[0081] 此外，对于“助推—滑行—助推”的飞行模式中的推力项，还应该区分无动力段及动力段过程，推力项的表达式为：
[0082]
[0083] 2.2滑行点火时间迭代求解方程
[0084] 根据飞行轨迹数值预测方程式，为实现终端交班点参数约束，控
制模型中需要确定的变量为点火时间tig和程序指令角。对于推力大小不可调节的固体火箭，
燃料消耗秒流量不可调节且燃料必须耗尽关机，导致动力段飞行时间不可控制，并且需要在耗尽关机点处满足终端状态约束，通常为箭下点高度、当地弹道倾角以及弹道偏角等状态量。
[0085] 为了保证自适应制导在线应用的计算效率以及收敛性，将滑行点火时间的求解与多约束制导指令的求解解耦计算。由于点火时间的变化对终端高度影响更为敏感，故建立点火时间与终端高度约束的一维非线性微分方程。而此时的俯仰角指令采用动力段全程零攻角即 的模式和偏航角指令全程为零ψ(t)＝0，尽可能的为后续满足多约束条件的指令迭代过程预留充足的空间。因此，点火时间的求解方程和终端约束方程为：
[0086]
[0087] 对于一维求根问题，当方程的解在迭代区间内时，黄金分割法在迭代5～10步以内可获得高精度解，此时的点火时间迭代区间为[0，Tmax]。
[0088] 2.3多约束条件下模型预测静态规划方法
[0089] 对于连续的微分函数，采用数值积分后通常可以描述为离散形式。将无量纲化的飞行器运动方程式，考虑为离散形式的一般非线性系统，其状态和输出动力学如下：
[0090]
[0091] 其中，X∈Rn,U∈Rm,Y∈Rp,k＝1,2,......，N，提高一组合适的控制变量Uk,k＝1,2......，N-1，使得输出的最终值YN达到预期值 YN关于用 表示，用泰勒展开得：
[0092]
[0093] 根据小偏差摄动原理，可以将残差表示为：
[0094]
[0095] 然后，根据式可以进一步将k+1步的误差表示成：
[0096]
[0097] 其中dXk,dUk分别是状态向量和控制向量在第k步的误差项。同样地，时间步骤N-1的状态误差dXN-1可以根据迭代步骤N-2的状态误差和控制误差进行递推。依次将dXN-2展开为dXN-3和dUN-3的递推形式，直到k＝1，则模型预测静态规划中误差的递推表达式为：
[0098]
[0099] 其中，状态转移矩阵A及灵敏度矩阵B的表达式为：
[0100]
[0101] 由于离散模型中每个
节点的控制变量与其他节点的状态变量是线性无关的，因此可以在任何时间点做出独立的控制。另外，当离散节点N很大时，状态转移矩阵A及灵敏度矩阵B将是一个计算密集型变量，可以通过递归计算的形式来节省计算时间。
[0102] 2.4具有过程约束和终端约束的MPSP制导模型
[0103] 为了应用MPSP方法，来实现对制导指令的过程约束并保证制导指令的平滑性，引入虚拟控制量姿态角速率：
[0104]
[0105] 称 为虚拟控制量。定义状态向量 控制向量 引入虚拟控制量后动力学方程变为：
[0106]
[0107] 对式应用Euler积分法则进行离散化，可以得到：
[0108]
[0109] 为了得到状态转移矩阵A及灵敏度矩阵B，分别对式求偏导数，得到 以及 的符号表达式。把它们代入式即可递归地求得敏感矩阵Bj,j＝1,2,…,N-1。 的符号表达式见表1。
[0111]
[0112] 其中，链式递推偏导数的表达式依次为：
[0113]
[0114]
[0115] 为了通过最小化的控制量来满足终端多约束条件，采用二次型性能指标，表达式为：
[0116]
[0117] 其中， 是上一步的控制量的解决方案，而dUk是控制量中相应的修正量。根据最优化理论可得最小的指令修正角序列为：
[0118]
[0119] 其中， 以及
[0120]
[0121] 对于正对角矩阵Rj＞0，最小化目标函数式使得虚拟控制量 相对于前一次迭代得到的虚拟控制量剖面的变化量最小，进而保证满足多约束条件的制导指令保持在由数值预测模型生成的小攻角轨迹附近。另一方面，当飞行过程中迭代的攻角超过过程约束αmax时，飞行轨迹应先保证飞行过程约束条件，然后控制向量再尽可能贴近数值预测模型生成的小攻角轨迹的控制向量，响应的数学处理手段为，增大对应过程约束变量的权值。为了满足这一要求，权值函数Rj的取值定义为：
[0122]
[0123] 3.闭环多约束迭代制导流程
[0124] 数值预测模型生成的小攻角轨迹确定后，MPSP方法不断迭代更新制导指令剖面直到末端约束的误差小于
指定的
门限值，在实现终端多约束条件的同时构成闭环迭代形式。当一个制导周期结束后，以得到的制导指令剖面作为下一制导周期的初始参考轨迹继续迭代更新制导指令剖面，如此循环，直到固体火箭达到交班条件。虽然在MPSP方法的推导过程中用Euler积分法则对连续系统进行了离散化，但在仿真过程中可以应用四阶Runge-Kutta方法进行积分以提高精度。整个制导的流程如图5所示。
[0125] 本发明中未做详细描述的内容均为现有技术。
[0126] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
