技术领域
[0001] 本
发明涉及高超声速飞行器的有限时间轨迹快速生成方法。
背景技术
[0002] 目前高超声速技术是一项集航天、航空技术优势的新技术领域，是21世纪航空航天技术领域的战略制高点，其发展将会对未来军事发展战略、空间技术、武器体系构建乃至整个科学技术进步产生重大影响。大
力开展高超声速技术研究对于巩固和提高我国的综合实力和国际地位具有重要的战略意义。我国已经把高超声速飞行器科技工程重大专项列入中长期发展规划。近十年在高技术方面的突飞猛进，带动了高超声速推进技术、精确导航与控制技术等的发展，特别是在高超声速巡航武器和天地往返运输系统高投入研究的带动下，涌现出体现各国研究特色的组合推进系统及高超声速飞行器新概念，为高超声速飞行器从构想到实现奠定了
基础。
[0003] 高超声速飞行器一般是指飞行
马赫数大于5，目前提出的高超声速飞行器概念，按功能分主要包括
水平起降航天运载器、空天飞机、高超声速飞机和高超声速巡航导弹等。该飞行器与常规飞行器的最大不同之处在于其飞行环境的不确定性和对飞行器自身在高速飞行条件下力学特性预测的不准确性。以X-43A为代表的吸气式高超声速飞行器所采用的超声速燃烧
冲压发动机性能同马赫数、高度等飞行状态密切相关，且发动机性能余量较小，为了获得足够的动压和发动机进气流量，这类飞行器必须始终在较稠密的大气层内飞行，而飞行器
气动性能在高超声速条件下存在许多难以预测的特性。因此大气内高超声速飞行闭环制导方式将是必然的选择，其制导系统必须具有自主、实时生成飞行轨迹和控制系统可以实现的制导指令的能力。具体而言，制导系统不但需要给控制系统提供实时的、完整的未来飞行轨迹信息，而且当控制律是基于预先规划的飞行轨迹信息时，要求必须能够根据飞行器的当前状态和飞行任务，在线实时更新轨迹设计，给出新的制导指令。由于闭环制导是基于对未来飞行轨迹的预测结果给出当前的轨迹控制指令，而实际飞行过程中，飞行器模型和大气环境条件都存在各种难以准确预测的不确定因素。因此如何实现飞行轨迹的快速优
化成为闭环制导的关键问题。
[0004] 闭环制导指令是指飞行过程中根据当前状态和任务生成的引导飞行器到达目标的飞行轨迹控制指令，它是决定高超声速飞行器能否实现其飞行方式和应用价值的关键问题。因此找到一种既能够解决这类多约束、大不确定性的非线性优化问题又能保证求解的在线性、快速性、实时性的优化方法是十分必要的。
[0005] 传统的轨迹优化方法包括间接法、直接法、动态规划法、微分方法等。间接法是基于Pontryagin极小值原理将最优控制问题转换为Hamilton边值问题(Hamilton Boundary Value Problem-HBVP)。该方法虽然有求解
精度较高、最优解满足一阶最优性条件的优点。但由于其推导最优解的过程较为复杂和繁琐并且不适宜求解有路径约束的最优控制问题，故而不能用来解决高超声速飞行器的轨迹快速优化问题。
[0006] 直接法无需求解最优性必要条件，而是将连续最优控制问题离散并参数化，直接应用数值方法对性能指标寻优。其优点包括不需要推导一阶最优性条件、收敛域相对间接法更宽广、对初值估计经度要求不高。其不足表现在不能保证获得的非线性规划解是原最优控制问题的解。然而，高超声速飞行器的轨迹优化问题本质上是有约束的非线性最优控制问题，具有高度的非线性和严格的约束条件，因此直接法仍无法胜任。
[0007] 动态规划法是解决多阶段决策过程最优化的一种数学方法，其最优性原理为：最优策略的子策略总是最优的。该方法的优点是计算原理简单，计算精度相对较高，有严格的理论支持，但对于高超声速飞行器轨迹优化这样的复杂优化问题，其计算结果需要大量的存储量，而且寻求全局最优解的过程也十分繁琐，因此无法满足轨迹优化的快速性和实时性。
[0008] 高超声速飞行器快速轨迹优化是完成飞行任务的重要保证，也是实现机动飞行的必要条件，运动模型的非线性和复杂约束条件的存在，使得高超声速飞行器轨迹优化问题变得非常复杂，优化过程中需要考虑的问题包括：再入环境及不确定性，高超声速飞行器再入经历高度、马赫数的大范围变化，将产生强热流、过载以及动压，大气环境参数变化剧烈，形成复杂的再入环境干扰和不确定性；约束条件复杂，不仅要考虑再入走廊约束，还要考虑控制约束、终端约束等；对计算时间提出更高的要求，高超声速飞行器以大于5Ma的速度飞行，因此对在线计算的时间约束很强，特别是在需要改变飞行轨迹的情况下，需要在线优化轨迹，因此需要提出高效率的
算法将优化时间控制在毫秒级；要求较高的制导精度，对于执行高精度的特殊任务(如打击动态目标点)，需要进行高精度制导。
发明内容
[0009] 本发明的目的是为了解决高超声速飞行器的轨迹快速优化问题时，间接法的推导最优解的过程较为复杂和繁琐不适宜求解有路径约束的最优控制问题、直接法不能保证获得的非线性规划解是原最优控制问题以及动态规划法无法满足轨迹优化的快速性和实时性形成复杂的再入环境干扰和环境不确定性等问题，而提出的一种高超声速飞行器的有限时间轨迹快速生成方法。
[0010] 上述的发明目的是通过以下技术方案实现的：
[0011] 步骤一、将高超声速飞行器再入末段运动分析及建模得到高超声速飞行器的运动模型；
[0012] 步骤二、将得到的高超声速飞行器再入末段轨迹的运动模型进行优化问题描述，形成非线性优化问题；
[0013] 步骤三、对高超声速飞行器再入末段非线性优化问题凸化处理，将优化问题的优化指标和优化问题的约束描述为二次型凸问题；
[0014] 步骤四：根据CVXGEN代码编写规范描述二次型凸优化问题，利用在线编译器对代码进行编译生成高速解算器；
[0015] 步骤五：将高速解算器嵌入到MATALB环境下，实现二次型凸优化问题的数值求解，并对数值求解结果进行分析；即完成了一种高超声速飞行器的有限时间轨迹快速生成方法。
[0016] 发明效果
[0017] 本发明的整个设计过程以高超声速飞行器闭环制导中的关键问题——快速轨迹优化为设计目标；建模过程中以其再入段为例进行动力学建模；优化问题描述过程中考虑了飞行器的
状态方程约束、再入走廊约束、终端约束、控制量约束并将
能量函数作为优化指标；问题转换过程中以明确的理论前提为基础实现了非线性优化问题的凸化处理并引入了滚动时域
优化算法以解决机动目标点的打击问题；最终利用凸优化求解工具CVXGEN实现高超声速飞行器再入段轨迹的快速优化。
[0018] 随着计算机性能的不断提升，一种新的优化思想进入人们的
视野。它借鉴预测控制的思想，对飞行器轨迹进行在线的有限时段优化，称为滚动时域优化(Receding Horizon-RH)。实际轨迹优化问题往往存在一定程度的不确定性，比如飞行器到达目的地之前或许目标已发生变化，且环境参数的详细信息也很难预先准确获知。而且轨迹越长，不确定性就越大。滚动优化在有限时段优化轨迹，在应对不确定性上具有优势，同时有利于实现轨迹的实时规划。
[0019] 此种高超声速飞行器的有限时间轨迹快速生成方法，实质是基于凸优化求解工具箱CVXGEN实现的多次
迭代滚动时域优化算法。该算法将多约束下的非线性轨迹优化问题转化为一系列凸优化问题，借助CVXGEN工具将凸优化问题以一种高级语言描述，并编译形成该类凸优化问题的可靠的高速解算器；与已有的直接法、间接法、动态规划法等优化方法相比具备以下几个优点：一是本发明提出一种多次迭代的滚动时域优化算法对一类具有大不确定性、多约束的非线性优化问题进行求解，通过非线性问题的凸优化转化确保了算法的快速性、实时性(单步优化计算周期能够控制在毫秒或微秒级)及精度，并且在仅考虑计算误差的前提下，精度能达到1m以内，克服了现有轨迹优化算法运算量大实时性差的问题；二是利用优化工具CVXGEN生成的解算器能实现高速求解且稳定可靠、占用内存少，且问题的转化编译过程完全由
软件自动进行，无需专业技巧，使用简单易上手，利于工程实现；三是该方法紧贴高超声速飞行器轨迹快速优化背景，针对性强、目标明确，通过实现高超声速飞行器的在线轨迹优化可以进一步解决其闭环制导问题。
附图说明
[0020] 图1是具体实施方式一中提出的一种高超声速飞行器的有限时间轨迹快速生成方法
流程图；
[0021] 图2是具体实施方式五提出的编程界面示意图；
[0022] 图3是具体实施方式五提出的代码编译界面示意图；
[0023] 图4是具体实施方式六提出的对高超声速飞行器再入末段轨迹快速优化问题的仿真结果分析流程图；
[0024] 图5是具体实施方式一提出的高超声速飞行器弹道示意图；
[0025] 图6是
实施例提出的静态单目标点优化轨迹示意图，其中，横坐标为航程，纵坐标为飞行器的飞行高度；
[0026] 图7是实施例提出的动态单目标点优化轨迹示意图，其中，横坐标为航程，纵坐标为飞行器的飞行高度。
具体实施方式
[0027] 具体实施方式一：本实施方式的一种高超声速飞行器的有限时间轨迹快速生成方法由具体是按照以下步骤制备的：
[0028] 步骤一、将高超声速飞行器再入末段运动分析及建模得到高超声速飞行器的运动模型如图5；
[0029] 步骤二、将得到的高超声速飞行器再入末段轨迹的运动模型进行优化问题描述，形成非线性优化问题；
[0030] 步骤三、对高超声速飞行器再入末段非线性优化问题凸化处理(离散线性化)，将优化问题的优化指标和优化问题的约束描述为二次型凸问题；
[0031] 步骤四：根据CVXGEN代码编写规范描述二次型凸优化问题，利用在线编译器对代码进行编译生成高速解算器；
[0032] 步骤五：将高速解算器嵌入到MATALB环境下，实现二次型凸优化问题的数值求解，并对数值求解结果进行分析，如图1所示；即完成了一种高超声速飞行器的有限时间轨迹快速生成方法。
[0033] 以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉
本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。
[0034] 本实施方式效果：
[0035] 对于高超声速飞行器而言，由于其飞行环境复杂，对象模型和参数都存在很大不确定性，飞行任务也可能根据当前局势发生变化，这要求高超声速飞行器在飞行过程中要具备一定的在线轨迹优化能力。而在线轨迹优化恰是闭环制导问题的核心，由于高超声速飞行器的高速特点和复杂动力学耦合特点，为保持系统
稳定性和制导精度，一般要求制导指令更新周期时间1s。本发明所提出的高超声速飞行器的有限时间轨迹快速生成方法实现了其飞行轨迹的快速、在线优化，并且将其优化时间控制在毫秒级，达到了技术指标的要求，同时实现了对机动目标点的精确打击。
[0036] 本实施方式的整个设计过程以高超声速飞行器闭环制导中的关键问题——快速轨迹优化为设计目标；建模过程中以其再入段为例进行动力学建模；优化问题描述过程中考虑了飞行器的状态方程约束、再入走廊约束、终端约束、控制量约束并将能量函数作为优化指标；问题转换过程中以明确的理论前提为基础实现了非线性优化问题的凸化处理并引入了滚动时域优化算法以解决机动目标点的打击问题；最终利用凸优化求解工具CVXGEN实现高超声速飞行器再入段轨迹的快速优化。
[0037] 随着计算机性能的不断提升，一种新的优化思想进入人们的视野。它借鉴预测控制的思想，对飞行器轨迹进行在线的有限时段优化，称为滚动时域优化(Receding Horizon-RH)。实际轨迹优化问题往往存在一定程度的不确定性，比如飞行器到达目的地之前或许目标已发生变化，且环境参数的详细信息也很难预先准确获知。而且轨迹越长，不确定性就越大。滚动优化在有限时段优化轨迹，因此在应对不确定性上具有优势，同时有利于实现轨迹的实时规划。
[0038] 此种高超声速飞行器的有限时间轨迹快速生成方法，实质是基于凸优化求解工具箱CVXGEN实现的多次迭代滚动时域优化算法。该算法将多约束下的非线性轨迹优化问题转化为一系列凸优化问题，借助CVXGEN工具将凸优化问题以一种高级语言描述，并编译形成该类凸优化问题的可靠的高速解算器；与已有的直接法、间接法、动态规划法等优化方法相比具备以下几个优点：一是本实施方式提出一种多次迭代的滚动时域优化算法对一类具有大不确定性、多约束的非线性优化问题进行求解，通过非线性问题的凸优化转化确保了算法的快速性、实时性(单步优化计算周期能够控制在毫秒或微秒级)及精度，并且在仅考虑计算误差的前提下，精度能达到1m以内，克服了现有轨迹优化算法运算量大实时性差的问题；二是利用优化工具CVXGEN生成的解算器能实现高速求解且稳定可靠、占用内存少，且问题的转化编译过程完全由软件自动进行，无需专业技巧，使用简单易上手，利于工程实现；三是该方法紧贴高超声速飞行器轨迹快速优化背景，针对性强、目标明确，通过实现高超声速飞行器的在线轨迹优化可以进一步解决其闭环制导问题。
[0039] 具体实施方式二：本实施方式与具体实施方式一不同的是：步骤一所述取土场的选址的过程为：步骤一中将高超声速飞行器再入末段运动建模得到高超声速飞行器的运动模型具体过程为：
[0040] (1)将高超声速飞行器再入过程发动机停止工作，不再受推力和控制力，主要受到重力和
空气动力的作用，进行远距离滑翔；此外，本发明不考虑飞行器的绕质心运动，认为绕质心运动处于瞬时平衡或配平状态，将飞行器视为质点仅研究其质心运动；简化成质心动力学方程：
[0041]
[0042] 其中，r——飞行器质心相对于地心高度，矢量表示；
[0044] R——飞行器所受空气动力，矢量表示；
[0045] mg——飞行器所受万有引力，矢量表示；
[0046] t∈[t0,tf]表示高超声速飞行器再入过程用于轨迹优化的时间；
[0047] t0表示高超声速飞行器再入过程初始时间；
[0048] tf表示高超声速飞行器再入过程终端时间；
[0049] (2)在此基础上，暂不考虑地球的自转、假设地球为一圆球，再入过程无侧滑力，即
侧滑角始终为零，得到高超声速飞行器再入末段的运动模型(即三
自由度的高超的再入运动方程是)：
[0050]
[0051] 其中，V——高超声速飞行器质心在半速度
坐标系中的速度，m/s；
[0052] ——高超声速飞行器质心在半速度坐标系中的速度相对时间的导数；
[0053] r——高超声速飞行器质心相对于地心的高度，m；
[0054] ——高超声速飞行器质心相对于地心的高度相对时间的导数；
[0055] γ——航迹角，rad；
[0056] ——航迹角相对时间的导数；
[0057] λ——经度，rad；
[0058] ——经度相对时间的导数；
[0059] φ——纬度，rad；
[0060] ——纬度相对时间的导数；
[0061] ψ——航向角，rad；
[0062] ——航向角相对时间的导数；
[0063] σ——速度斜角，rad；
[0064] D——阻力，N；
[0065] L——升力，N；
[0066] μ——地球引力常数；
[0067] m——再入飞行器质量；
[0068] 其中，状态变量x＝[r,γ,V,λ,φ,ψ]T，控制量u为
攻角和速度的倾角，即u＝T[α,σ] ；
[0069] 所谓纵向运动，是指飞行器运动参数侧滑角β、滚动角γ、速度倾斜角γv、
偏航角ψ、弹道偏角ψv、绕x,y轴的旋转
角速度ωx,ωy、z轴上的位移等恒等于零的运动；即飞行器的质心在飞行平面(或弹体坐标系纵向铅垂面)内的平移运动和绕z轴的转动运动所组成；所以，在纵向运动中，飞行器飞行速度V在弹道倾角θ、
俯仰角 及攻角α绕z轴的旋转角速度 x,y上的位移是随时间变化的，其中，飞行速度V、弹道倾角θ、俯仰角 攻角α及绕z轴的旋转角速度 和x,y上的位移通常称为纵向运动的运动学参数，简称为纵向运动参数；根据(2)式的高超声速飞行器再入末段运动模型，得到高超声速飞行器再入末段(再入飞行器攻击目标的末段)的纵向运动模型：
[0070]
[0071] 其中，h——飞行器距地球表面的距离；
[0072] ——飞行器距地球表面的距离相对时间的导数；
[0073] L——飞行器所受的升力；
[0074] D——飞行器所受的阻力；
[0075] 由于只考虑高超声速飞行器再入末段(再入飞行器攻击目标的末段)的纵向运动，不考虑其侧向运动的影响；故可取速度倾角σ＝0，控制量只有攻角即u＝α，状态量T取为x＝[h γ V] 得到升力L和阻力D的表达式如下：
[0076]
[0077]
[0078] 其中，CL0，CL1，CL2，CL3和CD0，CD1，CD2，CD3均为常数；ρ为大气层空气
密度，ρ的指-βh数形式为ρ＝ρ0e ，其中ρ0是海平面处的大气密度，β是常数；S中文含义为高超声速再入飞行器翼面参考面积，e为自然对数的底数。其它步骤及参数与具体实施方式一相同。
[0079] 具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤一所述取土场的选址的过程为：步骤二中将得到的高超声速飞行器再入末段轨迹的运动模型进行优化问题描述，形成非线性优化问题的具体过程为：
[0080] (1)、将(2)式所示高超声速飞行器再入末段运动模型表示为微分方程的形式：
[0081]
[0082] 其中，(3)式为非线性方程即状态方程约束，x＝[r,γ,V,λ,φ,ψ]T为状态变T量，u＝[α,σ] 为控制量；这是优化过程中必须要满足的状态方程约束； 表示状态变量和控制量均属于的实数域；m表示控制量的维数，n表示状态变量的维数；
[0083] 在建立高超声速飞行器再入末段的运动模型(三自由度的高超的再入运动方程)中，没有推力作用，控制量为攻角α和速度倾角σ；速度倾角σ主要控制高超的侧向机动飞行，攻角α则主要影响高超的气动系数，从而改变飞行器的纵向轨迹；满足一定的控制量约束：
[0084] umin≤u(t)≤umax,t∈[t0,tf](7)
[0085] (2)飞行轨迹需要满足的端点约束包括起始点约束和终端点约束，设初始时间t0和终端时间tf为条件，固定轨迹优化的初始时刻t0的状态为x0，固定终端时刻tf对应的状态量为xf，
[0086] 则对应的终端约束为：
[0087]
[0088] 其中，q表示终端约束的维数；ψ终端约束函数；
[0089] (3)高超声速飞行器再入末段为再入飞行器攻击目标的末段即高超声速飞行器返回过程最具特征的飞行段，是因为飞行器利用
地球大气层这种天热条件，在返回过程中减速下降，同时伴随巨大的能量消耗；因此，高超的再入过程满足严格的热流率 动压M和过载N等约束条件：
[0090]
[0091] 其中，ρ表示再入飞行器质心高度对应的大气密度；KQ表示热流密度常数；
[0092] KQ＝7.9686×10-5Js2m3.5kg0.5
[0093] 即 则(5)式可表示为(6)式，称为再入走廊约束即路径约束：此外，由于再入末段(再入飞行器攻击目标的末段)属于大气层内运动，因此无需考虑再入走廊约束；
[0094] ζ(x(t),u(t))≤0,t∈[t0,tf](6)，其中，ζ为再入走廊约束；
[0095] (4)高超声速飞行器再入末段轨迹的运动模型进行优化问题表示为如下形式：
[0096]
[0097] 其中，φ(xf,tf)和L(x,u,t)为优化指标函数，取能量函数作为优化指标；J为优化指标，(xf,tf)为终端时刻的状态量。其它步骤及参数与具体实施方式一或二相同。
[0098] 具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤三中对高超声速飞行器再入末段非线性优化问题凸化处理(离散线性化)，将优化问题的优化指标和优化问题的约束描述为二次型凸问题具体过程为：
[0099] 如前所述，本发明欲利用基于CVXGEN的凸优化算法来求解高超声速飞行器的再入末段轨迹优化问题；然而(8)式所示的优化问题并不是“凸”问题，因此需要进行优化问题的转换即非线性优化问题凸化处理：
[0100] (1)理论前提，定义非线性系统：
[0101]
[0102] 转化成状态相关系数(StateDependentCoefficient-SDC)的形式：
[0103]
[0104] 其中， 为状态量x相对时间的导数，A(x)为状态量x决定的系数矩阵，原点x＝0是一个平衡点并且假设A(x)满足局部Lipschitz条件，这是原非线性系统存在并且解唯一的最低假设；(10)式描述的非线性系统由线性时变系统来逼近；该方法使得原本应用于线性系统分析的经典方法移植到非线性系统中使用：
[0105][0]
[0106] 当j≥1时，x (t)的初值通常被取为初始状态x0；满足线性时变系统的解收敛到(10)所示的非线性系统的解，如下所示：
[0107] limj→∞{x[j](t)}→x(t)(12)
[0108] 其中，需要注意的是第一个线性估计式 x[1](0)＝x0是第一个线性时不变系统，因为其中的x(t)被初始状态x0所替换，为第一个方程提供了一个常系数阵；{x[j](t)}j≥1为迭代计算j次后得到的第j个线性时变系统，而
[0109] 为线性估计式则均为线性时变方程；
[0110] 基于这种理论，将非线性优化问题转化为一系列的线性优化问题,求得的最优解将收敛到原非线性优化问题的解；
[0111] (2)高超声速飞行器再入末段(再入飞行器攻击目标的末段)的轨迹快速优化问题仍是一非线性优化问题，欲利用凸优化算法求解则必须进行问题的转换；基于(10)、(11)、(12)式所提供的理论前提分别对优化指标和约束进行凸化处理；即优化问题的“凸”化处理；基于上述理论前提，将(8)式所描述的非线性优化问题转化为CVXGEN处理的凸优化问题：
[0112] 步骤A、将状态约束方程 t∈[t0,tf]描述为状态相关系数(SDC)的形式并给定一个初始状态：
[0113]
[0114] 其非线性在于系数矩阵A′(x),B′(x)的非线性；根据公式(11)所示的理论前提，利用线性系统来逼近该非线性系统：
[0115]
[0116] 简记为：
[0117]
[0118] 其中，A′(x),B′(x)为公式(15)的系数矩阵，从而实现了线性时变系统对原非线性系统的迭代逼近；
[0119] 步骤B、由于滚动时域优化采用离散形式，对线性系统(15)进行离散化处理，
采样周期取为T，离散化步长为P，采用差分方法得到最终的凸化形式如下：
[0120] x[j](k+1)＝A[j-1](k)x[j](k)+B[j-1](k)u[j](k)+C[j-1](k),k＝1,2,...,P(16)[0121] 其中，A[j-1](k)＝T×A′[j-1](k)+I,B[j-1](k)＝T×B′[j-1](k)C[j-1](k)表示代j次得到的第j个线性时变系统离散线性化误差；
[0122] 步骤C、对于(8)式中的再入走廊约束ζ(x(t),u(t))≤0,t∈[t0,tf]采用相同的处理方法进行线性离散化得到如下的凸化形式：
[0123] Ca[j](k)x[j](k)+Cb[j](k)u[j](k)+Cc[j](k)≤0,k＝1,2...P
[0124] Ca(k)、Cb(k)、Cc(k)为对ζ(x(t),u(t))≤0进行离散线性化得到的系数项；
[0125] 步骤D、对(8)式中的终端约束和控制量约束凸化处理，分别得到如下形式：
[0126] x[j](P+1)-xtar＝0
[0127] umin≤u[j](k)≤umax,k＝1,2...P其中，xtar为目标
位置点信息；
[0128] 步骤E、给出优化指标；以能量函数作为优化指标函数，需要说明的是，为了得到更好的优化结果，将终端硬约束转换为软约束(即公式4给出的终端约束描述为优化指标的一部分作为优化指标的一部分)；最终的优化指标形式为：
[0129]
[0130] 步骤F、高超声速飞行器再入末段轨迹优化问题由非线性优化问题转换成了凸优化工具CVXGEN可解的凸问题；可解的凸问题的最终形式为：
[0131]
[0132] (3)凸优化问题是静态优化问题，然而针对本方案中的动态目标点，为了实现高超声速飞行器轨迹的快速动态优化，引入滚动时域优化(Receding Horizon-RH)算法，最终得到基于凸优化求解工具CVXGEN的滚动时域优化算法，用来求解高超声速飞行器再入飞行器攻击目标的再入末段即高超声速飞行器再入过程的末段快速轨迹优化问题以实现对动态单目标点的打击；将当前的状态x(i)作为初始状态，在当前i时刻处预测优化时域(优化步长)P内的模型数据，并据此进行优化问题的求解；根据(19)式对优化问题的描述可得到高超声速飞行器再入末段快速轨迹优化问题具体形式如下：
[0133]
[0134] 其中，i+k|i表示在i时刻对i+k步变量的预测值，j表示迭代次数；P称为控制时域，表示预测未来的时间长度，u是控制量，Ca、Cb和Cc为再入走廊约束的三个系数阵。其它步骤及参数与具体实施方式一至三之一相同。
[0135] 具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤四中根据CVXGEN代码编写规范描述二次型凸优化问题，利用在线编译器对代码进行编译生成高速解算器过程为：
[0136] 针对(20)所描述的高超声速飞行器再入末段快速轨迹优化问题，在凸优化求解工具CVXGEN中进行相应的单步迭代优化问题代码编写；利用CVXGEN提供的在线编译功能，编译生成基于代码的高速解算器；
[0137] 由(8)式可以看出高超声速飞行器再入过程的轨迹优化是一个多约束的非线性优化问题；基于高超的飞行特性，需要我们寻求一种快速求解算法；此外，由于任务目标点的机动性，更需要该优化算法具有在线性和实时性；这便是本发明所提供的基于凸优化算法求解工具CVXGEN实现的滚动时域优化方法；利用CVXGEN工具将从代码编写到求解器过程为：
[0138] (1)代码编写，界面如图2所示，编写的代码包括维数(Dimensions)、参数(Parameters)、变量(Variables)、目标函数(Minimize)和约束(Subjectto)；
[0139] (2)代码编写后，如果编写的代码没有问题时，进行代码的编译，编译后的界面如图3所示；生成5个主要的C源文件、用于MATLAB环境中的*.m文件和用于结果验证的CVX文件；其中*.m文件可直接嵌入到MATLAB下，作为优化问题的解算器；其中，5个主要的C源文件包括求解函数及核心程序的solver.c、KKT矩阵的因式分解求解的ldl.c、矩阵向量的赋值及矢量相乘运算的matrix_support.c、简单的驱动代码的testsolver.c和测试函数的util.c；
[0140] (3)使用MATLAB
接口；使问题的仿真分析更加高效，将优化问题的解算器即自定义解算器嵌入到MATLAB环境下，十分便于优化问题的求解。其它步骤及参数与具体实施方式一至四之一相同。
[0141] 具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤五中将高速解算器嵌入到MATALB环境下，实现二次型凸优化问题的数值求解，并对数值求解结果进行分析的具体过程如下：
[0142] (1)代码编译成功后将自定义高速解算器CVXGEN嵌入到MATLAB环境下，对再入末段轨迹快速优化问题中的参数进行相应的赋值；
[0143] (2)求解高超声速飞行器再入末段轨迹快速优化问题，得到最终高超声速飞行器再入末段飞行轨迹的仿真结果；
[0144] (3)对高超声速飞行器再入末段轨迹快速优化问题的仿真结果进行分析：利用CVXGEN生成优化问题(20)的高速求解器，将求解代码嵌入到MATLAB中，按照(11)式进行多次迭代求解，得到P时域内的最优控制量u*(i|i),...,u*(i+P|i)，根据滚动优化将u*(i|i)应用于纵向非线性运动模型公式(21)，计算得到下一时刻的状态x(i+1)；若没达到
指定位置转到步骤五开始计算，直至飞行器击中目标或飞行器到达指
定位置；求解过程的原理如图4所示；其中，u*(i+P|i)表示第i次优化得到的i+P个控制变量(即在i时刻对时刻最优控制量的预测值)；指定位置为满足各项优化指标和约束的高超声速飞行器再入末段飞行轨迹。其它步骤及参数与具体实施方式一至五之一相同。
[0145] 采用以下实施例验证本发明的有益效果：
[0146] 实施例一：
[0147] 本实施例一种高超声速飞行器的有限时间轨迹快速生成方法，具体是按照以下步骤制备的：
[0148] 步骤一、将高超声速飞行器再入末段运动分析及建模得到高超声速飞行器再入末段轨迹的运动模型；
[0149] 高超声速飞行器弹道如图5所示；研究其再入过程的末端，即图中所示的再入末段的轨迹优化问题；其运动模型如下所示：
[0150]
[0151] 式中，h——飞行器距地球表面的距离；
[0152] ——飞行器距地球表面的距离相对时间的导数；
[0153] γ——航迹角，rad；
[0154] ——航迹角相对时间的导数；
[0155] V——高超声速飞行器质心在半速度坐标系中的速度，m/s；
[0156] ——高超声速飞行器质心在半速度坐标系中的速度相对时间的导数；
[0157] L——飞行器所受的升力；
[0158] D——飞行器所受的阻力；
[0159] 升力和阻力的表达式分别为：
[0160]
[0161]
[0162] 其中，CL0，CL1，CL2，CL3和CD0，CD1，CD2，CD3均为常数；ρ为大气层空气密度，ρ的指-βh数形式为ρ＝ρ0e ，其中ρ0是海平面处的大气密度，β是常数；S中文含义为高超声速再入飞行器翼面参考面积；可以看出，升力和阻力由攻角决定，因此将控制量取为攻角即uT
＝α，状态量取为飞行高度、航迹角、飞行速度即x＝[h γ V] ；其它参数的取值如下表所示：
[0163]
[0164] 步骤二、将得到的高超声速飞行器再入末段轨迹的运动模型进行优化问题描述，形成非线性优化问题；
[0165] 高超声速飞行器再入末段已处于地球大气层内，因此不用再考虑由再入大气层引起的再入走廊约束；考虑状态方程约束、终端约束、控制量约束及优化指标得到如下形式的再入末段轨迹优化问题
[0166]
[0167] 该步骤是对轨迹优化问题的数学描述，因此没有需要赋值的参数；
[0168] 步骤三：对高超声速飞行器再入末段非线性优化问题凸化处理(离散线性化)，将优化指标和约束描述为二次型凸问题；
[0169] 对公式(4)描述的多约束非线性轨迹优化问题进行凸化处理，并引入滚动时域优化策略得到CVXGEN可解的凸优化问题，具体描述如下：
[0170]
[0171] 其中，i+k|i表示在i时刻对i+k步变量的预测值，j表示迭代次数；P称为控制时域，表示预测未来的时间长度，u是控制量；参数选取结果如下：
[0172]
[0173] 步骤四：根据CVXGEN代码编写规范描述二次型凸优化问题，利用在线编译器对代码进行编译生成高速解算器；
[0174] 利用凸优化求解工具CVXGEN的规范语言完成对公式(5)描述的二次型凸优化问题的代码编写，编译生成高速解算器；由于整个代码的编写过程都是参数化的，故该步并未涉及赋值问题；
[0175] 步骤五：将高速解算器嵌入到MATALB环境下，实现二次型凸优化问题的数值求解，并对数值求解结果进行分析；
[0176] 将CVXGEN编译生成的高速解算器嵌入到MATLAB环境下，完成轨迹优化问题的最终求解，参数赋值情况如下：
[0177]
[0178] 最终得到以静态单目标点为靶点的优化轨迹，如图6所示，其优化时间为0.060272s；利用滚动时域优化方法实现对机动单目标点的打击，其优化结果如图7所示，期间共滚动优化了5次，累计优化时间为0.162109s。
