CN111897214A - 一种基于序列凸优化的高超声速飞行器轨迹规划方法 - Google Patents

Abstract

一种基于序列凸优化的高超声速飞行器轨迹规划方法，它属于高超声速飞行器轨迹规划技术领域。本发明解决了传统序列凸优化方法存在的可行性问题和收敛性问题。本发明的序列凸优化部分针对高超滑翔飞行段展开设计，提出了带罚函数的置信域加速算法。算法分为两步，第一步对非线性约束引入松弛变量，放弃置信域约束，目的是能够在更大的解空间中寻找可行解。待微分方程约束误差足够小后，转入下一步规划。第二步将目标函数重设为最小化置信域误差，主要解决子问题与原问题不等价的问题。基于这种方式能够在较差初值下，准确而迅速地完成多约束轨迹规划工作，具有极大实用性。本发明可以应用于高超声速飞行器轨迹规划。

Description

一种基于序列凸优化的高超声速飞行器轨迹规划方法

技术领域

本发明属于高超声速飞行器轨迹规划技术领域，具体涉及一种基于序列凸优化的高超声速飞行器轨迹规划方法。

背景技术

高超声速飞行器是一类飞行马赫数大于5，依靠气动升力控制在临近空间内长时间远距离飞行的飞行器。从进攻方角度而言，空基发射的高超声速飞行器航速快、航程远，在未来战场上可以实现先发制人和全球打击，具有极高的战略价值。然而，其飞行环境长时间处于临近空间的稠密大气层内，热力学环境复杂恶劣，为了保证弹载设备以及机身结构正常，飞行器必须满足热流、动压、过载等多方面的约束。此外，飞行环境等的强不确定性，更对飞行器飞行路径规划提出了更高更严格的要求。因此，多约束、强不确定性条件下，高超声速飞行器轨迹规划问题已成为现阶段高超声速飞行器研究的热点和难点。

飞行器轨迹优化设计问题本质上是一类非线性最优控制问题，受到各种强约束条件限制。求解该类问题的方法可分为两类：直接法与间接法。间接法基于庞德里亚金极大值原理，在状态变量基础上扩展协态变量转化为哈密尔顿边值问题，然后利用打靶法、有限差分法、梯度复位算法等对边值问题求解。直接法，则是将连续的弹道优化问题直接离散并参数化，转化成参数优化问题，无需求最优解的必要条件。这种方法本质上是将无限维的动态优化问题转化为有限维的非线性规划问题，然后通过迭代规划求解原问题近似解。

凸优化理论兴起于上世纪70年代，从数学上而言，对于一个凸问题局部最优解即为全局最优解，算法解集的最优性得到了保证。接着发展出了用于求解一般凸问题的自对偶内点法算法，不需要提供初值即可在指定精度内收敛。然而，凸优化方法在一般非线性问题中的应用还不够完备，其主要难点就在于问题的转换。凸优化应用时要求原问题保持凸性，而一般的非线性问题具有高度的非凸性，如何将实际问题近似为凸问题需要极大技巧。一种应用广泛的转化方式称之为序列凸优化，是基于线性化理论得来。利用泰勒展开近似后的函数一定具有凸性，因此只要利用凸优化求解算法就可以获得最优解。序列凸优化的一般性思路可以概括为，迭代中每一步运算利用上一步结果作为展开基准值，在若干次迭代后如果两次解足够接近，则认为找到了原问题解。但这种方法应用存在极大局限性，其根源在于转化并不是等价的，会产生以下问题：

1、可行性问题。子问题与原问题的近似程度与展开点选择休戚相关，在较差基准值下，当前子问题往往极大背离了原问题，并且找不到可行解；

2、收敛性问题。迭代中每一步运算都是在对不同问题求解，若不对问题进行处理，这一过程将是无穷尽且无意义的；

3、局部性问题。即使找到了问题的解，只能说明这是当前子问题的最优解，丧失了凸优化方法应有的全局最优性。

发明内容

本发明的目的是为了解决传统序列凸优化方法存在的可行性问题和收敛性问题，而提出了一种基于序列凸优化的高超声速飞行器轨迹规划方法，可应用于离线、在线高超声速飞行器轨迹规划。

本发明为解决上述技术问题所采取的技术方案是：一种基于序列凸优化的高超声速飞行器轨迹规划方法，所述方法具体包括以下步骤：

步骤一、将半速度系下的高超声速飞行器动力学模型的控制变量增广到状态变量中，根据扩维后的状态变量和设计的控制变量建立扩维高超声速飞行器动力学模型；

再根据建立的扩维高超声速飞行器动力学模型确定考虑约束的轨迹优化问题；

步骤二、利用高斯伪谱法将扩维高超声速飞行器动力学模型进行分区间段离散，建立新的轨迹优化问题；

步骤三、将步骤二建立的新的轨迹优化问题分解为若干个子问题，每个子问题的非线性约束在上一子问题解基础上采用泰勒展开近似，获得新的高超声速飞行器动力学模型；

通过考虑凸化后的禁飞区约束，速度变化约束，航迹角约束，控制变量约束以及攻角、倾侧角平滑约束，将步骤二建立的新的轨迹优化问题转化为有限维的线性规划问题；

步骤四、分步对线性规划问题的目标函数进行规划，通过对目标函数的求解，获得高超声速飞行器的轨迹规划结果。

本发明的有益效果是：本发明提出了一种基于序列凸优化的高超声速飞行器轨迹规划方法，本发明的序列凸优化部分针对高超滑翔飞行段展开设计，提出了带罚函数的置信域加速算法。算法分为两步，第一步对非线性约束引入松弛变量，放弃置信域约束，目的是能够在更大的解空间中寻找可行解。待微分方程约束误差足够小后，转入下一步规划。第二步将目标函数重设为最小化置信域误差，主要解决子问题与原问题不等价的问题。基于这种方式能够在较差初值下，准确而迅速地完成多约束轨迹规划工作，具有极大实用性。

通过实验证明，本发明的序列凸优化算法经过8次迭代运算解收敛，耗时123.0秒，而GPOPS方法需要经过18次迭代调整网格，需要耗时169.3秒。可见，本发明方法的收敛速度更快。

附图说明

图1为初值弹道水平面轨迹图；

图2为初值弹道空间轨迹图；

图3为初值弹道攻角规律图；

图4为初值弹道倾侧角规律图；

图5为序列凸优化迭代中水平面轨迹图；

图6为序列凸优化迭代中目标函数值曲线图；

图7为可行解水平面轨迹图；

图8为可行解空间轨迹图；

图9为可行解攻角规律图；

图10为可行解倾侧角规律图；

图11为可行解速度变化图；

图12为可行解航迹角变化图。

具体实施方式

具体实施方式一、本实施方式所述的一种基于序列凸优化的高超声速飞行器轨迹规划方法，所述方法具体包括以下步骤：

步骤一、将半速度系下的高超声速飞行器动力学模型的控制变量增广到状态变量中，根据扩维后的状态变量和设计的控制变量建立扩维高超声速飞行器动力学模型；

再根据建立的扩维高超声速飞行器动力学模型确定考虑约束的轨迹优化问题；

步骤二、利用高斯伪谱法将扩维高超声速飞行器动力学模型进行分区间段离散，建立新的轨迹优化问题；

步骤三、将步骤二建立的新的轨迹优化问题分解为若干个子问题，每个子问题的非线性约束在上一子问题解基础上采用泰勒展开近似，获得新的高超声速飞行器动力学模型；

通过考虑凸化后的禁飞区约束，速度变化约束，航迹角约束，控制变量约束以及攻角、倾侧角平滑约束，将步骤二建立的新的轨迹优化问题转化为有限维的线性规划问题；

步骤四、分步对线性规划问题的目标函数进行规划，通过对目标函数的求解，获得高超声飞行器的轨迹规划结果。

本发明主要包括三部分内容：高超声速飞行器动力学模型与规划参数设计、分段高斯伪谱参数化方法、序列凸优化算法。前两部分内容主要完成对原问题的描述与转化，为序列凸优化算法奠定基础。

高超声速飞行器动力学模型与规划参数设计部分主要介绍半速度系下飞行器动力学微分模型和规划中状态参数、控制参数选择问题。动力学微分模型描述目标受力情况及各个方向的加速度特性，建模了地心距、经度、纬度、速度、航迹角、航向角六个状态参量在攻角、倾侧角作用下变化趋势。对于规划问题而言，若直接选择攻角、倾侧角作为控制量会导致规划结果高频振荡，本发明因此扩维了状态矩阵，将攻角、倾侧角加入状态变量中，选择攻角导数、倾侧角导数作为控制变量。这种处理方式，解耦了状态变量与控制变量，增加了问题的健壮性。

分段高斯伪谱参数化方法是将原问题离散为有限维规划问题的一个过程，会在转化过程中，将控制量与状态量参数一并离散。采用高斯伪谱法优势在于，随着选择的节点增加，非线性问题会以谱精度收敛；同时，相较于一般插值方法，能避免Runge现象；并且，在求解过程中相关系数阵具有稀疏结构，能减少计算量。考虑到高斯节点具有两端密集中间稀疏的特点，会存在部分时间段函数近似效果差，同时过多数量的节点也难以计算。因此选择了分段高斯伪谱的方式，将原问题先等分若干个区间，再每段区间上取较小高斯节点。这对提高轨迹规划的精确程度起到了关键性作用。

序列凸优化部分针对高超滑翔飞行段展开了算法设计，提出了带罚函数的置信域加速算法。算法分为两步，第一步对非线性约束引入松弛变量，放弃置信域约束，目的是能够在更大的解空间中寻找可行解。待微分方程约束误差足够小后，转入下一步规划。第二步将目标函数重设为最小化置信域误差，主要解决子问题与原问题不等价的问题。基于这种方式能够在较差初值下，准确而迅速地完成多约束轨迹规划工作，具有极大实用性。

具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤一中，将半速度系下的高超声速飞行器动力学模型的控制变量增广到状态变量中，根据扩维后的状态变量和设计的控制变量建立扩维高超声速飞行器动力学模型；其具体过程为：

为了避免状态量间数量级差距过大引发矩阵奇异。采用以下形式进行无量纲化：

式中

分别代表真实的地心距、速度、角速度、与时间，具有物理量单位，而r、v、ω、t代表单位化后的地心距、速度、角速度、与时间，无量纲。

在半速度系下，建立无量纲的高超声速飞行器动力学模型为：

式中，r为单位地心距，

代表单位地心距r关于单位时间的导数；θ、φ分别为经度和纬度，

和

分别代表经度和纬度关于单位时间的导数；v为单位速度，

代表v对单位时间的导数；γ和ψ分别为航迹角和航向角，

代表γ关于单位时间的导数，

代表ψ关于单位时间的导数；σ是倾侧角，ω为地球单位自转角速度，D为归一化后阻力，L为归一化后升力；

D和L具有如下形式：

式中，R₀为地球半径，R₀＝6378km；ρ取大气密度指数模型，

ρ₀为海平面标准大气密度，ρ₀＝1.225kg/m³，h_s为高度常数，h_s＝7110m，h为当前真实的高度物理量；S为飞行器横截面积；m为飞行器质量；C_L和C_D均为气动系数；

C_L、C_D可以表示为攻角、马赫数的函数；

式(1)是关于6个状态变量与2个控制变量的高度非线性方程，当采用线性化方法进行凸化时，由于控制量与状态量间耦合关系会使得规划得到的控制变量高频振荡。因此，本发明中通过以下方式进行解耦处理；

把控制变量增广到状态变量中，则扩维后的状态变量x如下：

x＝[r；θ；φ；v；γ；ψ；α；σ] (4)

式中，α为攻角；

设计控制变量u为：

式中，

为攻角α的一阶导数，

为倾侧角σ的一阶导数；

则扩维高超声速飞行器动力学模型如下：

式中，

为状态变量x的一阶导数，f_p(x)为简略项(不考虑地球自转)的动力学列向量，f_p(x)∈R⁸，B为控制变量系数矩阵，f_ω(x)为角速度相关列向量，f_ω(x)∈R⁸，R代表实数域；

在式(6)中将缓变的地球自转扰动项在动力学模型中分离，式(6)中各部分具体形式如下：

其中，f_p(x)的表达式为：

f_ω(x)的表达式为：

B的表达式为：

具体实施方式三：本实施方式与具体实施方式二不同的是，所述步骤一中，根据建立的扩维高超声速飞行器动力学模型确定考虑约束的轨迹优化问题，其具体过程为：

禁飞区是飞行器为了避免威胁或出于地缘限制而必须避免飞入的区域。禁飞区通常被描述为指定位置和半径，无限高度的圆柱。

若高超声速飞行器到第i′个禁飞区中心(θ_i′,φ_i′)的距离大于禁飞区半径R_i′，则认为满足禁飞区约束：

L((θ,φ),(θ_i′,φ_i′))-R_i′≥0 (10)

其中，(θ,φ)为高超声速飞行器的位置，(θ_i′,φ_i′)为第i′个禁飞区中心的位置，L((θ,φ),(θ_i′,φ_i′))为高超声速飞行器到第i′个禁飞区中心的距离；

记高超声速飞行器运动过程中的初值约束和末值约束为：

x(t₀)＝x₀,x(t_f)＝x_f (11)

式中，t₀代表初始时刻，x₀为初始时刻t₀对应的状态变量，t_f代表终止时刻，x_f为终止时刻t_f对应的状态变量；

考虑控制系统性能，控制量攻角、倾侧角会存在幅值约束、变化率约束，同时飞行过程中对于飞行器状态也存在约束要求。

状态变量和控制变量存在的约束为：

x∈[x_min,x_max]，|u≤u_max (12)

其中，x_max和x_min分别为状态变量约束的上界和下界，u_max为控制变量的上界；

综上，轨迹优化问题可表述为：

问题1:

其中，J是目标函数。J根据飞行任务等要求有多种不同形式。

公式(13)的轨迹优化问题是一个连续非线性优化问题，具有非凸性，求解该优化问题的关键是将上述非凸问题转化为连续凸优化问题。

对于无限维的规划问题，需先采取一种参数化方法离散原问题。本发明采用的是高斯伪谱法。本发明基本思路可以概括为：将未知的状态时间历程与控制时间历程在一系列高斯点上离散，之后用这些离散的状态与控制分别构造拉格朗日插值多项式去逼近真实的状态与控制的时间历程，并借助简单的数学知识将运动学微分方程约束转化为一系列代数约束，经过以上步骤，最优控制问题最终转化为受一系列代数约束的参数优化问题。

高斯节点是在[-1,1]区间上非均匀分布的一些点，具有两边密集中间稀疏的性质。选择节点数量越多，非线性函数会以谱精度收敛。

具体实施方式四：本实施方式与具体实施方式三不同的是，所述步骤二的具体过程为：

利用高斯伪谱法，将扩维的高超声速飞行器动力学模型进行分区间段离散，将公式(13)的原问题时域均分为N个区间，在每段区间上均选择N₁个高斯点来近似原问题，并在N个区间内引入约束；

每一区间时间范围Δt＝(t_f-t₀)·N^-1，各区间端点可以表示为{t₀,t₁,...,t_N}。对于规划问题中时域范围[t_i-1,t_i]段，首先映射到时间区间[-1,1]上：

转换后，高斯节点τ取代实际时间t成为独立变量，τ＝-1时表示为τ₀对应t_i-1，τ＝1时表示为

对应t_i。依据选择的节点数量N₁，每段区间上高斯节点值τ_j,(j＝1,...,N₁)是确定的。

对于区间段[t_i-1,t_i]，将区间段映射到时间区间[-1,1]上，区间段[t_i-1,t_i]内的扩维高超声速飞行器动力学模型离散为：

式中，N₁为区间段[t_i-1,t_i]内的高斯节点数量，

为微分近似阵，τ_k为区间段[t_i-1,t_i]内的第k个高斯节点值，x(τ_k)为第k个高斯节点处的状态变量值，当k取0值时，x(τ₀)代表区间起点的状态变量值。f(x(τ_j),u(τ_j))为在高斯节点τ_j处状态变量的一阶导数；

对于微分近似阵的推导过程简述如下：

基于高斯节点τ_k，可构造拉格朗日插值多项式来近似状态量的时间历程。在高斯伪谱法下，插值需要同时使用N₁个高斯点τ₁,…,τ_N1和初始端点τ₀＝-1上的离散状态。任意时间点τ处的近似状态值x(τ)形式如下：

L_k(τ)为拉格朗日插值基函数，k＝0,1,2,…,N₁，形式如下：

将上式对转换后的时间τ求导，

微分近似阵D在每一个高斯节点处(j＝1,...,N₁)的取值可以表示为：

引入新的函数f_aug(x(τ_j),u(τ_j),t)＝(t_i-t_i-1)·f(x(τ_j),u(τ_j))，这步意义在于把飞行时间视为了一个新的输入参数。因此，对于全程飞行时间未知的自由飞行问题，采用该算法也能求解，区别仅在于多了一维新控制量——时间。本发明仅研究约束更强的固定时间飞行问题。高斯伪谱法下系统微分方程只在高斯点处进行离散近似。

对于区间段端点x₁与x_-1，利用高斯求积公式约束：

式中，x_-1为t_i-1时刻对应的状态变量值，x₁为t_i时刻对应的状态变量值，w_j为在高斯节点τ_j处的高斯权重；

式中，τ_k′为区间段[t_i-1,t_i]内的第k′个高斯节点值，k′＝1,2，...，N₁且k′≠j；

采用这种方式，在每一区间段处将原问题中连续非线性动力学约束转化为了有限维的两个代数等式约束，虽然仍是非线性的，但已经能够利用常见的非线性规划算法求解，下一步将讨论如何利用序列凸优化方法线性化该问题。

则新的轨迹优化问题的形式如下：

问题2:

采用凸优化方法求解问题，要求目标函数为凸函数并且约束方程具备凸性，问题2中离散后的动力学约束、禁飞区约束显然仍是非线性并且非凸的，因此序列凸优化的第一步便是采用连续线性化方法转化这些约束。

连续线性化过程可以概括为：将原问题分解为若干个子问题，每个子问题的非线性约束在上一子问题解基础上采用泰勒展开近似。这样处理后子问题会具有二阶锥规划的形式，满足凸约束条件。在迭代求解这些子问题过程中，如果解是收敛的，那么最后一个子问题的解能满足原问题，认为找到了一组可行解。

在高超飞行器飞行任务中，目标函数一般设计为简单线性函数形式，如最小飞行时间、最大再入速度等，因此本发明中不对目标函数线性化处理。

具体实施方式五：本实施方式与具体实施方式四不同的是，所述步骤三的具体过程为：

对于公式(6)的扩维高超声速飞行器动力学模型，以x^*表示上一步子问题解，则一阶泰勒展开后新的动力学模型形式如下：

式中，A(x^*)为f_p(x^*)的偏导系数阵；f_p(x^*)表示上一步子问题计算得到的动力学列向量，f_ω(x^*)表示上一步子问题计算得到的角速度相关列向量；

动力学方程中未展开地球自转项而是直接利用上一步解，使f_ω(x)≈f_ω(x^*)，原因在于该项在动力学方程中量级较小，近似处理不会产生较大误差。

式中，a₂₁、a₃₁、a₄₁、a₅₁、a₆₁、a₂₃、a₆₃、a₂₄、a₃₄、a₄₄、a₅₄、a₆₄、a₂₅、a₃₅、a₄₅、a₅₅、a₆₅、a₂₆、a₃₆、a₆₆、a₄₇、a₅₇、a₆₇、a₅₈和a₆₈均为中间变量；

禁飞区约束要求解集位于圆的外侧，是该问题中另一个非凸约束，采用相似的处理方式。用f₁(θ,φ)＝L((θ,φ),(θ_i′,φ_i′))-R_i′来简记该约束形式，则

凸化后的禁飞区约束为：

f₁(θ,φ)≈f₁(θ^*,φ^*)+f₁′(θ^*,φ^*)·[θ-θ^*；φ-φ^*]

f₁(θ^*,φ^*)+f₁′(θ^*,φ^*)·[θ-θ^*；φ-φ^*]≥0 (21)

式中，L((θ^*,φ^*),(θ_i′,φ_i′))代表上一步子问题计算得到的高超声速飞行器到第i′个禁飞区中心的距离，中间变量f₁(θ^*,φ^*)＝L((θ^*,φ^*),(θ_i′,φ_i′))-R_i′，f₁′(θ^*,φ^*)为f₁(θ^*,φ^*)的偏导系数，θ^*与φ^*是泰勒展开基点对应的经度与纬度，通过上一步迭代求解得到，f₁′(θ^*,φ^*)的形式如下：

f₁′(θ^*,φ^*)＝[2(θ^*-θ_i′)2(φ^*-φ_i′)] (22)

由于采用了泰勒展开的线性化形式，只有在基准值附近才能保证近似精度，因此需要在原问题中附加一个置信域约束。

|x-x^*|≤δ (23)

δ表示各归一化后状态量的极限偏差取值范围：

R₀、v₀分别代表实际的地心距、速度初值。

若能引入额外的满足真实条件的约束，能使得原问题更加完备，在一定程度上加速求解过程。

对于无动力滑翔飞行器而言，飞行过程伴随着能量耗散，因而速度曲线呈现出递减的趋势。

速度变化约束为：

Δv≤0 (24)

式中，Δv代表后一时刻飞行器速度与前一时刻飞行器速度的差值；

其次，本发明针对的是高超飞行器平衡滑翔飞行问题，因此飞行中航迹角存在约束要求，可以采用拟平衡滑翔约束，限制航迹角在接近0的小量δ_γ范围内。

航迹角约束为：

δ_γ≤γ≤0 (25)

δ_γ为任意小量；

接下来，考虑对规划结果平滑程度的要求。一般而言，对于同一个规划问题会存在若干个满足要求的可行解，规划算法只能获得其中一种解的形态，因此想要获得控制量曲线更平滑的解，最直接的方式就是引入对平滑程度的控制。从数学角度上，控制量曲线抖动表征

这两个量对应的变化率变化剧烈。但该问题中的规划变量没有攻角、倾侧角二阶导相关项，无法直接约束。于是思路变为，尝试对现有的控制量施加额外约束条件，以使控制量变化率不过分剧烈。该规划问题由于采用了分段高斯伪谱法，每一段区间上状态控制量与其他区间相对独立，因此可以做如下处理：

对于每个区间段，除端点外的控制变量u_j在区间内同号，j＝1,...,N₁；

u_j′·u_j＞0,j′≠j,j′＝1,...,N₁,j＝1,...,N₁ (26)

式中，u_j′代表在高斯节点τ_j′处的控制变量，j′＝1,...,N₁且j′≠j，u_j代表在高斯节点τ_j处的控制变量；

该约束显然也是非凸的，并且上一步凸化方式泰勒展开，此时约束形式无法使用。但由于规划过程是迭代进行的，随着目标函数逐渐减小，迭代前后两次状态量也逐渐趋近，因此可以利用迭代中上一步状态量符号约束当前待规划状态。

子问题中，每一区间段上攻角、倾侧角平滑约束表述为：

式中，

代表上一步子问题计算中区间段上首个高斯节点处的攻角导数值，

代表上一步子问题计算中区间段上首个高斯节点处的倾侧角导数值；

代表区间段上高斯节点τ_j处的攻角导数值，

代表区间段上高斯节点τ_j处的倾侧角导数值；

在全区间段[t₀,t_f]上，将全区间段分为N个区间，各区间端点分别表示为{t₀,t₁,...,t_N}，其中，t_N＝t_f，共需要考虑N次该约束条件，攻角、倾侧角平滑约束的完整形式为：

式中，

代表在区间段[t_i-1,t_i]上第j个高斯节点的攻角导数值，

代表在区间段[t_i-1,t_i]上第j个高斯节点的倾侧角导数值；

代表在上一步子问题中计算得到的在区间段[t_i-1,t_i]上首个高斯节点的攻角导数值；

代表在上一步子问题中计算得到的在区间段[t_i-1,t_i]上首个高斯节点的倾侧角导数值；

至此，将问题2已经转化为了有限维的线性规划问题；

问题3:

式中，Δt代表区间段[t_i-1,t_i]的时间范围，w代表高斯权重集合，

具体实施方式六：本实施方式与具体实施方式五不同的是，所述步骤四的具体过程为：

第一步规划目标函数：

其中，ξ₁为动力学约束的松弛因子(即式28第一行的动力学约束)，ξ₂为区间端点约束的松弛因子(即式28第二行的区间端点约束)，ξ₃为禁飞区约束下的松弛因子，上角标T代表转置，||·||₁为向量的1范数；c₁、c₂、p₁为各项的系数，在仿真时分别取值为100、200、100；

在第n-1步子问题的解基础上，求解第n步线性规划问题：

式中，

代表第n步子问题的动力学约束的松弛因子，

代表第n步子问题的区间端点约束的松弛因子，

代表第n步子问题的禁飞区约束下的松弛因子，公式(29)中，上角标n代表在第n步子问题下的值，上角标n-1代表在第n-1步子问题下的值；Δvⁿ代表第n步子问题与第n-1步子问题的单位速度差，γⁿ代表第n步子问题与第n-1步子问题的航迹角差；

代表第i个区间段内动力学约束的松弛因子，

代表第i个区间段内区间端点约束的松弛因子，i＝1,2,…,N，令中间变量

中间变量

中间变量

不断的进行迭代求解，直至

时停止迭代；

转入第二步规划目标函数：min 0.5qΔx_p ^TΔx_p，其中，Δx_p为置信域误差，q为系数，在仿真时取1即可；Δx_p具体形式如下：设置为前后两次计算部分状态参数差；

Δx_p＝[αⁿ-α^n-1 σⁿ-σ^n-1]^T (30)

式中，αⁿ代表第n步子问题下的攻角值，α^n-1代表第n-1步子问题下的攻角值，σⁿ代表第n步子问题下的倾侧角值，σ^n-1代表第n-1步子问题下的倾侧角值；

在第n-1步子问题解的基础上，求解第n步线性规划问题：

式中，

代表第n步子问题下的置信域误差，

代表第n步子问题下第i个区间段内的置信域误差，i＝1,2,…,N，当|Δx|≤ε时达到迭代收敛条件，将当前步子问题的计算结果作为高超声速飞行器的轨迹规划结果，其中：|Δx|代表前后两步子问题的状态变量的差值。

本实施方式是序列凸优化算法设计

若直接采用凸优化算法迭代求解问题3会产生以下问题：

1)可行性问题。凸优化算法的优势之一在于采用了自对偶内点法求解，对于一个凸问题不需要给初值即可得到最优解。然而在采用了连续线性化方法凸化问题之后，这一优势不再存在。泰勒展开使得每步运算的系数阵均来自于上次运算，初值直接决定了子问题解的性质。此外，子问题与原问题并非等价，子问题无解更是常见情况。现有的很多凸优化文章都是在已经可行的轨迹上，开展轨迹优化工作，无法推广到工程应用中。

2)收敛性问题。基于泰勒展开得到的子问题本质上是一系列不等价的问题，直接求解这些子问题，得到的解一般意义上是不会收敛的。前文已经提及只有在解收敛的情况下，最后一个子问题解才是原问题的可行解。因此，解的收敛性需要设计算法来保证。

3)局部性问题。凸优化算法的另外一个优势在于局部极值点即为全局极值点。在序列凸优化算法下，这一优势也失去了，其根源即在于设计的子问题与原问题不可能等价。在这种算法下规划出的结果只具有局部极值性。

为此，本发明设计了带罚函数的置信域加速算法，来解决可行性问题与收敛性问题；算法分为两大部分，第一部分对转化后的非线性约束引入松弛变量，放弃置信域约束。该部分目标函数选择为松弛变量内积，随迭代进行松弛变量逐渐减小，代表着微分方程约束误差也在逐渐减小。其值小于一定阈值后，转入下一步规划。第二部分将目标函数重设为最小化置信域误差，主要解决子问题与原问题不等价的问题。置信域误差越小，子问题越接近原问题。

第一步规划目标函数：

式中ξ₁、ξ₂为微分方程约束下的松弛变量，ξ₃为禁飞区约束下的松弛变量。理论上，随迭代进行以上所有松弛变量应该趋于0。

第二步规划目标函数：

min 0.5qΔx_p ^TΔx_p

式中Δx_p对应泰勒展开下两次状态变量部分项之差，该项整体称为置信域误差。对于问题3，置信域加速算法具体形式如下：

第一阶段在第k-1次值基础上求解如下第k次线性规划问题：

Δv^k≤0,δ_γ≤γ^k≤0

|x^k-x^k-1|≤δ_x

|u^k|≤u_max

该子问题约束中，前两项为高斯伪谱法下转化的微分约束；第三项为禁飞区约束；第四项代表初值约束与末态约束；第五项为滑翔过程状态约束；第六项为置信域约束；第七项为控制量约束，该处可以为攻角、倾侧角的幅值约束或者攻角变化率、倾侧角变化率约束。

该步目的在于，在一个较差的初值下能寻找到合适可行域以初步满足约束条件。这种问题与约束设计方式，旨在消除前文提及的因连续线性化而引起的可行性难题。首先，对原问题中所有非线性约束引入了松弛变量，以减小因不合适的线性化产生的问题不可行程度。同时，新的子问题目标函数设置为了使所有松弛因子最小化，这样在若干步迭代后，松弛因子足够小，可认为当前子问题的解初步满足了原问题的所有约束条件。

令

判别条件：当

时，认为当前解对应的子问题与原问题在约束意义上足够接近。在已经找到较好满足原问题约束条件的解的情况下，可以在更严格约束下规划问题。接下来转入第二步规划算法：

f₁(θ^k-1,φ^k-1)+f₁′·(θ^k-θ^k-1,φ^k-φ^k-1)^T≤0

Δv^k≤0,δ_γ≤γ^k≤0

|u^k|≤u_max

第二步规划问题中目标函数设置为前后两次计算部分状态参数差。式中，

dx_p＝[α^k-α^k-1 σ^k-σ^k-1]^T

该步解决前文提到的收敛性问题。当dx_p足够小时，线性化误差可以忽略，认为当前子问题是原问题的一个可行解。目标函数中状态参数差只选择了攻角增量与倾侧角增量，一方面是因为，作为实际控制量，攻角序列与倾侧角序列能完备的表示整个问题其他的过程参数；另一方面在于，减少对其他状态量的限制，能最大程度上自由地寻找可行解；同时，去除冗余项还能加速计算过程。在约束方面与上一步相比，舍弃了所有松弛因子与置信域约束。前者是因为，第一步完成后的规划结果已经能初步满足约束条件，以此为基准生成的子问题基本可行；后者是因为，现在置信域约束已经以目标函数的形式重新包含在了问题中。此外，前文提及的平滑约束仅在第二步规划时引入，这是考虑到第一步算法只提供合适初值即可。

与线性规划问题LP1相比，无论从目标函数复杂程度还是规划变量个数上而言，线性规划问题LP2都得到了极大简化，因此计算效率更快。仿真中，LP1每次迭代求解时间一般在17秒左右，而LP2的求解时间一般在10秒左右。

整个问题迭代收敛条件设置为：

|Δx|≤ε

迭代求解问题中，若某一次计算结果使上式成立，则认为当前规划结果为满足原问题的可行解。算例中收敛邻域取如下值：

仿真部分

本发明涉及的仿真结果是在以下仿真环境中得到：MATLAB2014a版本，IntelCorei7-7700HQ，2.80GHz CPU。

该部分设计了一个仿真算例来验证本发明序列凸优化算法的可行性。某型号高超飞行器启滑段状态参数固定，根据飞行任务要求在1500秒后，需要抵达西南方向1500公里外的目标点，飞行中途存在一个禁飞区。相关参数如表1和表2所示：

表1仿真参数

表2控制约束

算法第一步置信域约束取为：

为了展现本发明设计算法的性能，将与非线性规划领域已广泛应用的求解工具包GPOPS-2作对比。在数值仿真阶段，GPOPS软件与序列凸优化方法输入的原问题相同，并且求解动力学模型均采用式(6)所述的扩维模型。在参数设置阶段，由于该软件在参数化过程中也采用了伪谱法，并同样对整个区间进行了等分处理。因此，配点个数与区间长度选择应与本发明保持一致。GPOPS软件的一大优点在于采用了自适应网格细化方法，针对给定的容许误差e_r，若当前配点参数下对应的解满足不了要求，会自动调整配点以降低误差。因而，容许误差e_r取值是影响求解精度与速度的关键。在该软件中，e_r被定义为“等式约束近似时最大相对误差”。在仿真时将先采用序列凸优化方法，确定出迭代收敛时容许误差值e_r，再设置GPOPS中参数，最后对比二者结果，如表3所示，注：配点个数差异是由于伪谱方法不同所造成的。

表3伪谱法参数

首先，采用本发明设计的置信域加速算法求解原问题P1。在对迭代初值设置时，该算法已能很好消除初值敏感性，不需要提供较好初值即可规划出可行解。如图3、图4所示，初值弹道按常值攻角、常值倾侧的控制规律提供。其水平面弹道形态如图1所示，可见未满足该问题禁飞区约束、末态约束；三维弹道如图2所示，呈现出跳跃弹道的趋势，而本问题中需要以平衡滑翔抵达目标点。

在该初值下，序列凸优化算法经过8次迭代运算解收敛，耗时123.0秒。迭代中每次运算，轨迹变化如图5所示。运算中目标函数变化如图6所示，随迭代进行趋于0。在最后一次计算中，LP2问题中等式约束近似时产生的最大相对误差为1.527×10^-8，因而在GPOPS中设置容许误差极值e_r＝1×10^-7。

两种方法下的可行解见图7至图12，均满足所有约束要求。由于未额外约束可行解过程参数，因此从结果上二者有差异，凸优化方法找到的可行解是从禁飞区上侧绕飞后抵达目标，而GPOPS方法的解是从禁飞区下侧绕飞后抵达目标，如图7、8所示。因此两种方法下，控制量规律也会产生相应变化。具体而言，凸优化方法下前500秒飞行，倾侧角保持为负值以能够从上侧绕过禁飞区；后500秒飞行倾侧角保持为正值以纠正航向能抵达目标。而GPOPS解与之相反，前500秒保持倾侧角为正值，后500秒保持为负值。速度变化方面如图11所示，虽然下降趋势不同，但末速度都位于10Ma至11Ma之间。而对于航迹角，由于采用了拟平衡滑翔约束，因此飞行期间一直保持在0附近，两种方式下变化趋势也不尽相同。总体而言，序列凸优化方法与GPOPS软件均找到了合适可行解，但就计算速度而言，GPOPS方法经过了18次迭代调整网格，耗时169.3秒，慢于本发明设计的序列凸优化算法，本发明方法与GPOPS方法的仿真结果对比如表4所示。

表4仿真结果对比

本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (6)

1.一种基于序列凸优化的高超声速飞行器轨迹规划方法，其特征在于，所述方法具体包括以下步骤：

步骤一、将半速度系下的高超声速飞行器动力学模型的控制变量增广到状态变量中，根据扩维后的状态变量和设计的控制变量建立扩维高超声速飞行器动力学模型；

再根据建立的扩维高超声速飞行器动力学模型确定考虑约束的轨迹优化问题；

步骤二、利用高斯伪谱法将扩维高超声速飞行器动力学模型进行分区间段离散，建立新的轨迹优化问题；

步骤三、将步骤二建立的新的轨迹优化问题分解为若干个子问题，每个子问题的非线性约束在上一子问题解基础上采用泰勒展开近似，获得新的高超声速飞行器动力学模型；

通过考虑凸化后的禁飞区约束，速度变化约束，航迹角约束，控制变量约束以及攻角、倾侧角平滑约束，将步骤二建立的新的轨迹优化问题转化为有限维的线性规划问题；

步骤四、分步对线性规划问题的目标函数进行规划，通过对目标函数的求解，获得高超声速飞行器的轨迹规划结果。

2.根据权利要求1所述的一种基于序列凸优化的高超声速飞行器轨迹规划方法，其特征在于，所述步骤一中，将半速度系下的高超声速飞行器动力学模型的控制变量增广到状态变量中，根据扩维后的状态变量和设计的控制变量建立扩维高超声速飞行器动力学模型；其具体过程为：

在半速度系下，建立无量纲的高超声速飞行器动力学模型为：

式中，r为单位地心距，

代表单位地心距r关于单位时间的导数；θ、φ分别为经度和纬度，

和

分别代表经度和纬度关于单位时间的导数；v为单位速度，

代表v对单位时间的导数；γ和ψ分别为航迹角和航向角，

代表γ关于单位时间的导数，

代表ψ关于单位时间的导数；σ是倾侧角，ω为地球单位自转角速度，D为归一化后阻力，L为归一化后升力；

D和L具有如下形式：

式中，R₀为地球半径，R₀＝6378km；ρ取大气密度指数模型，

ρ₀为海平面标准大气密度，ρ₀＝1.225kg/m³，h_s为高度常数，h_s＝7110m，h为当前真实的高度物理量；S为飞行器横截面积；m为飞行器质量；C_L和C_D均为气动系数；

把控制变量增广到状态变量中，则扩维后的状态变量x如下：

x＝[r；θ；φ；v；γ；ψ；α；σ] (4)

式中，α为攻角；

设计控制变量u为：

式中，

为攻角α的一阶导数，

为倾侧角σ的一阶导数；

则扩维高超声速飞行器动力学模型如下：

式中，

为状态变量x的一阶导数，f_p(x)为简略项的动力学列向量，f_p(x)∈R⁸，B为控制变量系数矩阵，f_ω(x)为角速度相关列向量，f_ω(x)∈R⁸，R代表实数域；

其中，f_p(x)的表达式为：

f_ω(x)的表达式为：

B的表达式为：

3.根据权利要求2所述的一种基于序列凸优化的高超声速飞行器轨迹规划方法，其特征在于，所述步骤一中，根据建立的扩维高超声速飞行器动力学模型确定考虑约束的轨迹优化问题，其具体过程为：

若高超声速飞行器到第i′个禁飞区中心(θ_i′,φ_i′)的距离大于禁飞区半径R_i′，则认为满足禁飞区约束：

L((θ,φ),(θ_i′,φ_i′))-R_i′≥0 (10)

其中，(θ,φ)为高超声速飞行器的位置，(θ_i′,φ_i′)为第i′个禁飞区中心的位置，L((θ,φ),(θ_i′,φ_i′))为高超声速飞行器到第i′个禁飞区中心的距离；

记高超声速飞行器运动过程中的初值约束和末值约束为：

x(t₀)＝x₀,x(t_f)＝x_f (11)

式中，t₀代表初始时刻，x₀为初始时刻t₀对应的状态变量，t_f代表终止时刻，x_f为终止时刻t_f对应的状态变量；

状态变量和控制变量存在的约束为：

x∈[x_min,x_max]，|u|≤u_max (12)

其中，x_max和x_min分别为状态变量约束的上界和下界，u_max为控制变量的上界；

轨迹优化问题表述为：

问题1:

其中，J是目标函数。

4.根据权利要求3所述的一种基于序列凸优化的高超声速飞行器轨迹规划方法，其特征在于，所述步骤二的具体过程为：

利用高斯伪谱法，将扩维的高超声速飞行器动力学模型进行分区间段离散，将公式(13)的原问题时域均分为N个区间，在每段区间上均选择N₁个高斯点来近似原问题，并在N个区间内引入约束；

对于区间段[t_i-1,t_i]，将区间段映射到时间区间[-1,1]上，区间段[t_i-1,t_i]内的扩维高超声速飞行器动力学模型离散为：

式中，N₁为区间段[t_i-1,t_i]内的高斯节点数量，

为微分近似阵，τ_k为区间段[t_i-1,t_i]内的第k个高斯节点值，x(τ_k)为第k个高斯节点处的状态变量值，f(x(τ_j),u(τ_j))为在高斯节点τ_j处状态变量的一阶导数；

引入新的函数f_aug(x(τ_j),u(τ_j),t)＝(t_i-t_i-1)·f(x(τ_j),u(τ_j))，对于区间段端点x₁与x_-1，利用高斯求积公式约束：

式中，x_-1为t_i-1时刻对应的状态变量值，x₁为t_i时刻对应的状态变量值，w_j为在高斯节点τ_j处的高斯权重；

式中，τ_k′为区间段[t_i-1,t_i]内的第k′个高斯节点值，k′＝1,2，…，N₁且k′≠j；

则新的轨迹优化问题的形式如下：

问题2:

5.根据权利要求4所述的一种基于序列凸优化的高超声速飞行器轨迹规划方法，其特征在于，所述步骤三的具体过程为：

对于公式(6)的扩维高超声速飞行器动力学模型，以x^*表示上一步子问题解，则一阶泰勒展开后新的动力学模型形式如下：

式中，A(x^*)为f_p(x^*)的偏导系数阵；f_p(x^*)表示上一步子问题计算得到的动力学列向量，f_ω(x^*)表示上一步子问题计算得到的角速度相关列向量；

式中，a₂₁、a₃₁、a₄₁、a₅₁、a₆₁、a₂₃、a₆₃、a₂₄、a₃₄、a₄₄、a₅₄、a₆₄、a₂₅、a₃₅、a₄₅、a₅₅、a₆₅、a₂₆、a₃₆、a₆₆、a₄₇、a₅₇、a₆₇、a₅₈和a₆₈均为中间变量；

凸化后的禁飞区约束为：

f₁(θ^*,φ^*)+f₁′(θ^*,φ^*)·[θ-θ^*；φ-φ^*]≥0 (21)

式中，L((θ^*,φ^*),(θ_i′,φ_i′))代表上一步子问题计算得到的高超声速飞行器到第i′个禁飞区中心的距离，中间变量f₁(θ^*,φ^*)＝L((θ^*,φ^*),(θ_i′,φ_i′))-R_i′，f₁′(θ^*,φ^*)为f₁(θ^*,φ^*)的偏导系数，θ^*与φ^*是泰勒展开基点对应的经度与纬度，f₁′(θ^*,φ^*)的形式如下：

f₁′(θ^*,φ^*)＝[2(θ^*-θ_i′)2(φ^*-φ_i′)] (22)

|x-x^*|≤δ (23)

δ表示各归一化后状态量的极限偏差取值范围：

速度变化约束为：

Δv≤0 (24)

式中，Δv代表后一时刻飞行器速度与前一时刻飞行器速度的差值；

航迹角约束为：

δ_γ≤γ≤0 (25)

δ_γ为任意小量；

对于每个区间段，除端点外的控制变量u_j在区间内同号，j＝1,...,N₁；

u_j′·u_j＞0,j′≠j,j′＝1,...,N₁,j＝1,...,N₁ (26)

式中，u_j′代表在高斯节点τ_j′处的控制变量，j′＝1,...,N₁且j′≠j，u_j代表在高斯节点τ_j处的控制变量；

每一区间段上攻角、倾侧角平滑约束表述为：

式中，

代表上一步子问题计算中区间段上首个高斯节点处的攻角导数值，

代表上一步子问题计算中区间段上首个高斯节点处的倾侧角导数值；

代表区间段上高斯节点τ_j处的攻角导数值，

代表区间段上高斯节点τ_j处的倾侧角导数值；

在全区间段[t₀,t_f]上，将全区间段分为N个区间，各区间端点分别表示为{t₀,t₁,...,t_N}，共需要考虑N次该约束条件，攻角、倾侧角平滑约束的完整形式为：

式中，

代表在区间段[t_i-1,t_i]上第j个高斯节点的攻角导数值，

代表在区间段[t_i-1,t_i]上第j个高斯节点的倾侧角导数值；

代表在上一步子问题中计算得到的在区间段[t_i-1,t_i]上首个高斯节点的攻角导数值；

代表在上一步子问题中计算得到的在区间段[t_i-1,t_i]上首个高斯节点的倾侧角导数值；

至此，将问题2已经转化为了有限维的线性规划问题；

问题3:

式中，Δt代表区间段[t_i-1,t_i]的时间范围，w代表高斯权重集合，

6.根据权利要求5所述的一种基于序列凸优化的高超声速飞行器轨迹规划方法，其特征在于，所述步骤四的具体过程为：

第一步规划目标函数：

其中，ξ₁为动力学约束的松弛因子，ξ₂为区间端点约束的松弛因子，ξ₃为禁飞区约束下的松弛因子，上角标T代表转置，||·||₁为向量的1范数；c₁、c₂、p₁为各项的系数；

在第n-1步子问题的解基础上，求解第n步线性规划问题：

代表第i个区间段内动力学约束的松弛因子，

代表第i个区间段内区间端点约束的松弛因子，i＝1,2,…,N，令中间变量

中间变量

中间变量

不断的进行迭代求解，直至

时停止迭代；

转入第二步规划目标函数：min 0.5qΔx_p ^TΔx_p，其中，Δx_p为置信域误差，q为系数；Δx_p具体形式如下：

Δx_p＝[αⁿ-α^n-1σⁿ-σ^n-1]^T (30)

式中，αⁿ代表第n步子问题下的攻角值，α^n-1代表第n-1步子问题下的攻角值，σⁿ代表第n步子问题下的倾侧角值，σ^n-1代表第n-1步子问题下的倾侧角值；

在第n-1步子问题解的基础上，求解第n步线性规划问题：

式中，

代表第n步子问题下的置信域误差，

代表第n步子问题下第i个区间段内的置信域误差，i＝1,2,…,N，当|Δx|≤ε时达到迭代收敛条件，将当前步子问题的计算结果作为高超声速飞行器的轨迹规划结果，其中：|Δx|代表前后两步子问题的状态变量的差值。

Images