CN111580535B - 一种基于凸优化的再入轨迹三维剖面规划方法及*** - Google Patents

Abstract

本发明公开了一种基于凸优化的再入三维剖面规划方法及***，针对再入飞行器在滑翔段的轨迹规划问题，提出了一种基于凸优化的三维剖面规划方法。首先，建立飞行器的三自由度运动模型，并选择合适的中间变量表征三维剖面，以此作为轨迹规划的控制量。其次，通过合适的凸化与离散方法，将三维剖面规划问题建模为凸优化问题，利用凸优化求解器进行求解。最后，利用凸优化问题的最优解反算得到攻角和倾侧角指令。本发明能有效解决滑翔轨迹的三维剖面规划求解效率低的难题。

Description

一种基于凸优化的再入轨迹三维剖面规划方法及***

技术领域

本发明涉及飞行器制导控制技术领域，可用于再入飞行器滑翔段轨迹规划与制导，尤其涉及一种基于凸优化的再入轨迹三维剖面规划方法及***。

背景技术

高超声速滑翔飞行器是一种主要在临近空间飞行，且飞行速度大于5马赫的飞行器。相比于传统弹道式飞行器，此类飞行器具有机动能力强、飞行速域宽、弹道全程可控等优点，在军事领域具有广阔应用价值。该类飞行器的再入轨迹可分为初始下降段和滑翔段，滑翔段由于飞行时间长、航程远、过程约束复杂等，是轨迹规划问题研究的重点和难点。传统轨迹规划方法的不足主要包括：1)基于极大值原理的间接法需要复杂公式推导，难以处理复杂约束；2)以伪谱法为代表的直接法，计算效率低，对初值比较敏感；3)传统滑翔段轨迹规划方法多是基于二维剖面，即固定攻角剖面，将倾侧角作为轨迹规划的唯一控制量，由于参数优化范围小，所以限制了飞行器的可用机动能力。在基于三维剖面的轨迹规划中，将攻角和倾侧角同时作为控制量，增大了参数优化范围，从而充分释放了飞行器的可用机动能力。但传统基于三维剖面的规划方法为了降低计算量采用了大量模型简化方法，降低了解的精度，未能很好解决精度与计算效率之间的矛盾。为了克服该不足，研究能够充分释放飞行器机动能力，且满足在线快速计算要求的滑翔段轨迹规划方法尤为重要。

凸优化是一类具备良好理论性质和高效求解效率的数学规划问题。由于该问题能够在多项式时间内收敛到全局最优解，目前已成为在线求解轨迹规划等实时最优控制问题的常用方法之一。基于凸优化求解轨迹规划问题主要包括三个步骤：1)将原始轨迹规划问题转化为一个连续的凸优化问题，包括对目标函数、非线性动力学、控制量约束和状态量约束等原问题中非凸项的凸化处理；2)将上述连续凸优化问题进行离散化，得到一个以有限个参数为优化变量的凸优化问题；3)利用凸优化求解器进行问题求解，得到控制量。凸优化方法的难点主要在于第一步，如何对原问题进行凸化处理，目前常用的方法包括变量替换、松弛等，针对实际问题需要灵活采用多种凸化方法进行处理。特别是对于一些具有强非线性动力学或状态约束的轨迹规划问题，就目前所知，连续线性化是处理这类约束的仅有方法，这种情况下由于需要求解一个凸优化子问题序列，以这个序列的收敛解作为原问题的最优解，因此该方法又称为序列凸优化。由于再入轨迹规划问题具有强非线性动力学与非线性状态约束，因此序列凸优化是求解这类问题的主要手段。原始对偶内点法是目前求解凸优化问题最有效的算法之一，很多凸优化求解器均以该算法为基础进行二次开发，该算法对凸优化问题的求解效率明显优于其它轨迹规划方法，为轨迹规划问题的在线求解提供了技术保证。

参考文献

[1]Liu X,Shen Z,Lu P.Entry Trajectory optimization by second-ordercone programming.Journal of guidance control and dynamics,2015.DOI:10.2514/1.G001210.

[2]Phillips T H.A common aero vehicle model,description,andemployment guide.http://www.dtic.Mil/matris/sbir041/srch/af03 1a.doc,2003.

发明内容

本发明所要解决的技术问题是，针对现有技术不足，提供一种基于凸优化的再入轨迹三维剖面规划方法及***，克服传统方法在求解三维剖面规划问题时计算效率低、精度低的问题，实现飞行器在大机动任务要求下的轨迹在线规划。

为解决上述技术问题，本发明所采用的技术方案是：一种基于凸优化的再入轨迹三维剖面规划方法，包括以下步骤：

1)给定飞行器的当前归一化机械能E₀、终端归一化机械能E_f、当前运动状态x₀＝[r₀λ₀φ₀θ₀σ₀]和预期终端状态x_f＝[r_fλ_fφ_fθ_fσ_f]，其中，r为地心距，λ为经度，φ为纬度，θ为速度倾角，σ为航迹偏航角，下标为0表示当前状态的取值，下标为f表示终端状态的取值；给定离散点数目N，N为正整数；

2)对滑翔段参考轨迹进行初始化，具体方法为：将地心距从当前值r₀到终端值r_f的变化规律假设为线性的，将从当前值r₀到终端值r_f的变化范围均分为N个区间，得到每个离散点上的初始地心距，采用同样方式得到经度、纬度、速度倾角、航迹偏航角、归一化机械能在每个离散点上的初始值，将这些初始值联合在一起记为初始轨迹x⁽⁰⁾，并令参考轨迹x^(k)＝x⁽⁰⁾，上标k表示当前迭代次数；

3)利用所述参考轨迹x^(k)，结合凸优化求解器，得到第k+1次迭代的参考轨迹x^(k+1)(i)；

4)以x^(k+1)(i)作为输出的最优轨迹，以第k+1次迭代的攻角指令和倾侧角指令(α^(k+1)(i),υ^(k+1)(i))作为输出的最优三维剖面。

借由上述过程，将三维剖面规划问题转化为一个凸优化问题，克服传统方法在求解三维剖面规划问题时计算效率低、精度低的问题。

上述步骤3)的具体实现过程包括：

3a)给定最小攻角α_min、最大攻角α_max、最小倾侧角υ_min和最大倾侧角υ_max；根据参考轨迹x^(k)计算升阻比下界(u₃)_min和升阻比上界(u₃)_max；基于参考轨迹x^(k)，计算动力学方程的线性化矩阵(A,B,C)、热流、过载、动压三个约束线性化后的系数

r_q、r_n；

3b)将上述基于当前参考轨迹x^(k)所得的连续变量(u₃)_min、(u₃)_max、(A,B,C)、

r_q、r_n进行等间距离散，得到每个离散点处的取值，表示为(u₃)_min(i)、(u₃)_max(i)、(A(i),B(i),C(i))、/>

r_q(i)、r_n(i)，其中，i表示离散点序号，变化范围从1到N；

3c)将离散化后的参数：(u₃)_min(i)、(u₃)_max(i)、(A(i),B(i),C(i))、

r_q(i)、r_n(i)以及预期终端状态x_f作为凸优化求解器的输入并进行求解，从而得到当前迭代的优化解(x^(k+1)(i),u^(k+1)(i))；

3d)根据(x^(k+1)(i),u^(k+1)(i))计算当前攻角和倾侧角指令，其中，

分别表示在离散点i处的三个控制量分量；/>

3e)令x^(k)＝x^(k+1)(i)，重复上述步骤3a)～3d)，直到下列收敛准则成立：

其中，ε为一个常值矢量，维数为5×1。

上述求解过程计算量相对较小，求解效率高。

本发明所述升阻比下界(u₃)_min和升阻比上界(u₃)_max分别表示的计算公式表示如下：

其中，

表示对括号中函数在攻角α的取值范围[α_min,α_max]内求极小，/>

表示对括号中函数在攻角α的取值范围[α_min,α_max]内求极大；G表示升阻比计算函数；Ma_r是由参考轨迹x^(k)根据大气模型得到的参考马赫数。

本发明矩阵A、B和矩阵C的计算公式表示如下：

C＝f(x^(k),u^(k),E)-Ax^(k)+h(x^(k))

其中，E表示当前离散点处的归一化机械能，所述当前离散点处的归一化机械能通过对归一化机械能变化域[E₀,E_f]进行等间距离散获得；

/>

上标k表示当前迭代次数，r^(k)表示当前参考轨迹x^(k)对应的地心距，φ^(k)表示当前参考轨迹x^(k)对应的纬度，V^(k)表示当前参考轨迹x^(k)对应的速度，θ^(k)表示当前参考轨迹x^(k)对应的速度倾角，σ^(k)表示当前参考轨迹x^(k)对应的航迹偏航角，D^(k)表示当前参考轨迹x^(k)对应的阻力加速度，

表示当前参考轨迹x^(k)对应的阻力加速度关于地心距的偏导数，且

R₀＝6371004m，u^(k)表示当前参考轨迹x^(k)对应的中间控制量，f(x^(k),u^(k),E)和h(x^(k))表示如下：

ω_e为地球自转角速度，ω_e＝7.292115×10^-5rad/s。

本发明中，热流、过载、动压三个约束线性化后的系数

r_q、r_n计算公式如下：

其中，

q_max和n_max分别为最大热流、最大动压、最大过载；/>

q_k和n_k分别为当前参考轨迹x^(k)确定的热流、动压和过载；r_k为当前参考轨迹x^(k)确定的地心距；β为大气密度指数模型中的常数。

所述攻角指令α^(k+1)(i)由

进行升阻比反插值得到，倾侧角指令由下式得到：

可见，攻角指令和倾侧角指令的计算简单，容易实现。

本发明中，为了进一步简化计算，常值矢量ε中的元素对应小于[5000.20.20.20.5]^T中的元素，其中[5000.20.20.20.5]^T第一项单位为米，后面四项单位为度。

相应地，本发明还提供了一种基于凸优化的再入轨迹三维剖面规划***，包括：

初始化单元，对滑翔段参考轨迹进行初始化，具体执行如下操作：将地心距从当前值r₀到终端值r_f的变化规律假设为线性的，将从当前值r₀到终端值r_f的变化范围均分为N个区间，得到每个离散点上的初始地心距，采用同样方式得到经度、纬度、速度倾角、航迹偏航角、归一化机械能在每个离散点上的初始值，将这些初始值联合在一起记为初始轨迹x⁽⁰⁾，并令参考轨迹x^(k)＝x⁽⁰⁾，上标k表示当前迭代次数；其中，所述归一化机械能包括飞行器的当前归一化机械能E₀、终端归一化机械能E_f；飞行器的当前运动状态x₀＝[r₀λ₀φ₀θ₀σ₀]和预期终端状态x_f＝[r_fλ_fφ_fθ_fσ_f]，r为地心距，λ为经度，φ为纬度，θ为速度倾角，σ为航迹偏航角，下标为0表示当前状态的取值，下标为f表示终端状态的取值；给定离散点数目N，N为正整数；

计算单元，用于根据所述参考轨迹x^(k)，结合凸优化求解器，得到第k+1次迭代的参考轨迹x^(k+1)(i)；

输出单元，用于输出最优轨迹x^(k+1)(i)和最优三维剖面；所述最优三维剖面为第k+1次迭代的攻角指令和倾侧角指令(α^(k+1)(i),υ^(k+1)(i))。

计算单元包括：

第一计算模块，用于根据参考轨迹x^(k)计算升阻比下界(u₃)_min和升阻比上界(u₃)_max；基于参考轨迹x^(k)，计算动力学方程的线性化矩阵(A,B,C)、热流、过载、动压三个约束线性化后的系数

r_q、r_n；

离散模块，用于将基于当前参考轨迹x^(k)所得的连续变量(u₃)_min、(u₃)_max、(A,B,C)、

、r_q、r_n进行等间距离散，得到每个离散点处的取值，表示为(u₃)_min(i)、(u₃)_max(i)、(A(i),B(i),C(i))、/>

r_q(i)、r_n(i)，其中，i表示离散点序号，变化范围从1到N；

第二计算模块，用于将离散化后的参数：(u₃)_min(i)、(u₃)_max(i)、(A(i),B(i),C(i))、

r_q(i)、r_n(i)以及预期终端状态x_f作为凸优化求解器的输入并进行求解，从而得到当前迭代的优化解(x^(k+1)(i),u^(k+1)(i))；

第三计算模块，用于根据(x^(k+1)(i),u^(k+1)(i))计算当前攻角和倾侧角指令，其中，

分别表示在离散点i处的三个控制量分量；/>

迭代模块，用于判断

是否成立，若成立，则输出x^(k+1)(i)；否则，重复第一计算模块、离散模块、第二计算模块和第三计算模块的操作。

本发明相对现有技术的有益效果为：本发明将三维剖面规划问题转化为一个凸优化问题，利用这类问题的良好理论性质，实现对原问题高效精确的求解，既充分释放了飞行器的机动能力，又克服了传统方法求解三维剖面效率较低的难题，实现了飞行器在大机动任务要求下的轨迹在线规划。

附图说明

图1为本发明地面轨迹对比结果；

图2为本发明空间三维轨迹曲线；

图3为本发明速度倾角和航迹偏航角随时间的变化曲线；

图4为本发明攻角和倾侧角随时间的变化曲线；

图5为本发明热流、过载和动压随时间的变化曲线；

图6为本发明高度、经纬度、速度倾角和航迹偏航角随迭代次数的收敛曲线；

图7为本发明方法流程图；

图8为本发明***结构框图；

图9为本发明计算单元结构框图。

具体实施方式

下面将结合附图和具体实例对本发明做进一步详细说明。

该实例中以美国CAV-H为仿真对象[2]，质量m＝907kg，参考面积Sm＝0.4839m²，仿真初始条件设置为：初始高度h₀＝80km，速度V₀＝7000m/s，速度倾角θ₀＝-2°，初始经纬度均为0°，初始航迹偏航角为90°。预期终端条件设置为：高度h_f＝25km，速度V_f＝2000m/s，速度倾角θ_f＝-4°，终端经度为100°，终端纬度为50°，终端航迹偏航角为30°。攻角取值范围为10到20度，倾侧角取值范围为-85到85度。最大热流

q_max＝90kPa，n_max＝4.5。计算机硬件条件为Inter Core i5-3470 [email protected]，内存RAM为8G。

如图7所示，具体包括如下步骤：

1)步骤1：由公式E＝V²/2-μ/r计算飞行器当前归一化机械能

E₀＝-0.5955和终端归一化机械能E_f＝-0.9641，其中，V表示飞行器速度，μ表示引力常数，μ＝3.986005×10¹⁴m³/s²。根据所给参数确定当前运动状态

x₀＝[1.012600-0.03491.5708]和预期终端状态

x_f＝[1.00391.74530.8727-0.06980.5236]；确定离散点数目N＝300；

2)步骤2：对滑翔段参考轨迹进行初始化，具体方法为：将地心距

r₀＝1.0126到r_f＝1.0039的变化范围均匀离散为300个区间，则初始化的地心距可以表示为

采用同样的方法得到初始化经度

初始化纬度/>

初始化速度倾角/>

初始化航迹偏航角

和离散后的归一化机械能{E₀,E₁,…,E_N}，其中E_N＝E_f；其中，上标0表示初始化参考轨迹的状态值，下标表示对应的离散点顺序。将这些初始化状态量联合在一起记为初始轨迹x⁽⁰⁾，并令

x^(k)＝x⁽⁰⁾，上标k表示当前迭代次数；

步骤3：迭代执行以下步骤：

3a)根据参考轨迹x^(k)计算升阻比下界(u₃)_min和升阻比上界(u₃)_max；

3b)基于参考轨迹x^(k)，计算动力学方程的线性化矩阵(A,B,C)，矩阵A、B和C的计算公式表示如下：

C＝f(x^(k),u^(k),E)-Ax^(k)+h(x^(k))

其中，E表示当前离散点对应的归一化机械能，当前离散点处的归一化机械能通过对归一化机械能变化域[E₀,E_f]进行等间距离散获得，对归一化机械能变化域[E₀,E_f]进行等间距离散得到的能量依次取值为{E₀,E₁,…,E_N}，

/>

上标k表示当前迭代次数，r^(k)表示当前参考轨迹x^(k)对应的地心距，φ^(k)表示当前参考轨迹x^(k)对应的纬度，V^(k)表示当前参考轨迹x^(k)对应的速度，θ^(k)表示当前参考轨迹x^(k)对应的速度倾角，σ^(k)表示当前参考轨迹x^(k)对应的航迹偏航角，D^(k)表示当前参考轨迹x^(k)对应的阻力加速度，

表示当前参考轨迹x^(k)对应的阻力加速度关于地心距的偏导数，且

R₀＝6371004m，u^(k)表示当前参考轨迹x^(k)对应的中间控制量，f(x^(k),u^(k),E)和h(x^(k))的表达式如下所示

ω_e为地球自转角速度，ω_e＝7.292115×10^-5rad/s。

3c)基于参考轨迹x^(k)，计算热流、过载、动压三个约束线性化后的系数

r_q、r_n；

3d)将上述基于当前参考轨迹x^(k)所得的连续变量(u₃)_min、(u₃)_max、(A,B,C)、

r_q、r_n进行均匀离散，分别得到(u₃)_min、(u₃)_max、(A,B,C)、/>

r_q、r_n在每个离散点处的取值，表示为(u₃)_min(i)、(u₃)_max(i)、(A(i),B(i),C(i))、/>

r_q(i)、r_n(i)，其中，i表示离散点序号，变化范围从1到N。

3e)将离散化后的参数：(u₃)_min、(u₃)_max(i)、(A(i),B(i),C(i))、

r_q(i)、r_n(i)以及预期终端状态x_f作为凸优化求解器的输入并进行求解，得到当前迭代的优化解(x^(k+1)(i),u^(k+1)(i))，其中，凸优化求解器可以通过开源方式获得或购买商用软件；/>

3f)根据(x^(k+1)(i),u^(k+1)(i))计算当前攻角和倾侧角指令，其中，

攻角指令α^(k+1)(i)由u₃ ^(k+1)(i)进行升阻比反插值得到，倾侧角指令/>

3g)令x^(k)＝x^(k+1)(i)，重复上述步骤3a)～3f)，直到下式收敛准则成立

其中，ε为一个常值矢量，本实施例中取值为ε＝[10001115]^T。

步骤4：以x^(k+1)(i)作为输出的最优轨迹，以(α^(k+1)(i),υ^(k+1)(i))作为输出的最优三维剖面。

根据上述步骤，本发明所提方法规划一条满足各项约束的可行轨迹计算耗时为12.6秒，而参考方法(即伪谱法)在同样参数设置和硬件条件下的计算耗时为235.00秒，由此证明了本发明中所提方法在同等计算条件下具有求解效率高的有益效果。

如图8，本发明的规划***包括：

初始化单元，具体执行如下操作：将地心距从当前值r₀到终端值r_f的变化规律假设为线性的，将从当前值r₀到终端值r_f的变化范围均分为N个区间，得到每个离散点上的初始地心距，采用同样方式得到经度、纬度、速度倾角、航迹偏航角、归一化机械能在每个离散点上的初始值，将这些初始值联合在一起记为初始轨迹x⁽⁰⁾，并令参考轨迹x^(k)＝x⁽⁰⁾，上标k表示当前迭代次数；其中，所述归一化机械能包括飞行器的当前归一化机械能E₀、终端归一化机械能E_f；飞行器的当前运动状态x₀＝[r₀λ₀φ₀θ₀σ₀]和预期终端状态x_f＝[r_fλ_fφ_fθ_fσ_f]，r为地心距，λ为经度，φ为纬度，θ为速度倾角，σ为航迹偏航角，下标为0表示当前状态的取值，下标为f表示终端状态的取值；给定离散点数目N，N为正整数；

计算单元，用于根据所述参考轨迹x^(k)，结合凸优化求解器，得到第k+1次迭代的参考轨迹x^(k+1)(i)；

输出单元，用于输出最优轨迹x^(k+1)(i)和最优三维剖面；所述最优三维剖面为第k+1次迭代的攻角指令和倾侧角指令(α^(k+1)(i),υ^(k+1)(i))。

如图9，本发明的计算单元包括三个计算模块和一个离散模块，分别执行如下操作：

第一计算模块，用于根据参考轨迹x^(k)计算升阻比下界(u₃)_min和升阻比上界(u₃)_max；基于参考轨迹x^(k)，计算动力学方程的线性化矩阵(A,B,C)、热流、过载、动压三个约束线性化后的系数

r_q、r_n；

离散模块，用于将基于当前参考轨迹x^(k)所得的连续变量(u₃)_min、(u₃)_max、(A,B,C)、

r_q、r_n进行等间距离散，分别得到(u₃)_min、(u₃)_max、(A,B,C)、/>

r_q、r_n在每个离散点处的取值，表示为(u₃)_min(i)、(u₃)_max(i)、(A(i),B(i),C(i))、/>

r_q(i)、r_n(i)，其中，i表示离散点序号，变化范围从1到N；

第二计算模块，用于将离散化后的参数：(u₃)_min(i)、(u₃)_max(i)、(A(i),B(i),C(i))、

r_q(i)、r_n(i)以及预期终端状态x_f作为凸优化求解器的输入并进行求解，从而得到当前迭代的优化解(x^(k+1)(i),u^(k+1)(i))；

第三计算模块，用于根据(x^(k+1)(i),u^(k+1)(i))计算当前攻角和倾侧角指令，其中，

分别表示在离散点i处的三个控制量分量；

迭代模块，用于判断

是否成立，若成立，则输出x^(k+1)(i)；否则，重复第一计算模块、离散模块、第二计算模块和第三计算模块的操作。

为了突出本发明中所提方法可以充分释放飞行器机动能力的有益效果，将所提方法得到的轨迹与基于二维剖面方法得到的轨迹进行对比。图1给出了地面轨迹的对比结果，可以看出，在给定再入点和预期目标点后，基于三维剖面的方法可以得到一条到达该目标的可行轨迹，但基于二维剖面的方法无法得到一条到达该目标的可行轨迹，因此基于三维剖面的规划方法可以使得飞行器的可达目标“更远”。由于任务参数设置是相同的，只有控制量选择的差异，说明本发明所提的三维剖面规划方法可以充分释放飞行器的机动能力，证明了本方法的在释放飞行器机动能力方面的有益效果。图2至图4给出了由本发明所提方法得到轨迹对应的空间飞行轨迹及其它参数与控制量的变化规律，均满足给定约束。图5说明了所得轨迹满足热流、动压和过载三项过程约束，最大值均未超过约束要求。图6说明了本发明所提方法是有效收敛的。

Claims (8)

1.一种基于凸优化的再入轨迹三维剖面规划方法，其特征在于，包括以下步骤：

1)给定飞行器的当前归一化机械能E₀、终端归一化机械能E_f、当前运动状态x₀＝[r₀λ₀φ₀θ₀σ₀]和预期终端状态x_f＝[r_fλ_fφ_fθ_fσ_f]，其中，r为地心距，λ为经度，φ为纬度，θ为速度倾角，σ为航迹偏航角，下标为0表示当前状态的取值，下标为f表示终端状态的取值；给定离散点数目N，N为正整数；

2)对滑翔段参考轨迹进行初始化，具体方法为：将地心距从当前值r₀到终端值r_f的变化规律假设为线性的，将从当前值r₀到终端值r_f的变化范围均分为N个区间，得到每个离散点上的初始地心距，采用同样方式得到经度、纬度、速度倾角、航迹偏航角、归一化机械能在每个离散点上的初始值，将这些初始值联合在一起记为初始轨迹x⁽⁰⁾，并令参考轨迹x^(k)＝x⁽⁰⁾，上标k表示当前迭代次数；

3)利用所述参考轨迹x^(k)，结合凸优化求解器，得到第k+1次迭代的参考轨迹x^(k+1)(i)；

4)以x^(k+1)(i)作为输出的最优轨迹，以第k+1次迭代的攻角指令和倾侧角指令(α^(k+1)(i),υ^(k+1)(i))作为输出的最优三维剖面；

步骤3)的具体实现过程包括：

3a)给定最小攻角α_min、最大攻角α_max、最小倾侧角υ_min和最大倾侧角υ_max；根据参考轨迹x^(k)计算升阻比下界(u₃)_min和升阻比上界(u₃)_max；基于参考轨迹x^(k)，计算动力学方程的线性化矩阵(A,B,C)、热流、过载、动压三个约束线性化后的系数

r_q、r_n；

3b)将上述基于当前参考轨迹x^(k)所得的连续变量(u₃)_min、(u₃)_max、(A,B,C)、

r_q、r_n进行等间距离散，得到(u₃)_min、(u₃)_max、(A,B,C)、/>

r_q、r_n在每个离散点处的取值，分别表示为(u₃)_min(i)、(u₃)_max(i)、(A(i),B(i),C(i))、/>

r_q(i)、r_n(i)，其中，i表示离散点序号，变化范围从1到N；

3c)将离散化后的参数：(u₃)_min(i)、(u₃)_max(i)、(A(i),B(i),C(i))、

r_q(i)、r_n(i)以及预期终端状态x_f作为凸优化求解器的输入并进行求解，从而得到当前迭代的优化解(x^{(k +1)}(i),u^(k+1)(i))；

3d)根据(x^(k+1)(i),u^(k+1)(i))计算当前攻角和倾侧角指令，其中，

分别表示在离散点i处的三个控制量分量；

3e)令x^(k)＝x^(k+1)(i)，重复上述步骤3a)～3d)，直到下列收敛准则成立：

其中，ε为一个常值矢量，维数为5×1。

2.根据权利要求1所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，所述升阻比下界(u₃)_min和升阻比上界(u₃)_max分别表示的计算公式表示如下：

其中，

表示对括号中函数在攻角α的取值范围[α_min,α_max]内求极小，/>

表示对括号中函数在攻角α的取值范围[α_min,α_max]内求极大；G表示升阻比计算函数；Ma_r是由参考轨迹x^(k)根据大气模型得到的参考马赫数。

3.根据权利要求1所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，矩阵A、B、C的计算公式表示如下：

C＝f(x^(k),u^(k),E)-Ax^(k)+h(x^(k))；

其中，E表示当前离散点处的归一化机械能，所述当前离散点处的归一化机械能通过对归一化机械能变化域[E₀,E_f]进行等间距离散获得；

上标k表示当前迭代次数，r^(k)表示当前参考轨迹x^(k)对应的地心距，φ^(k)表示当前参考轨迹x^(k)对应的纬度，V^(k)表示当前参考轨迹x^(k)对应的速度，θ^(k)表示当前参考轨迹x^(k)对应的速度倾角，σ^(k)表示当前参考轨迹x^(k)对应的航迹偏航角，D^(k)表示当前参考轨迹x^(k)对应的阻力加速度，

表示当前参考轨迹x^(k)对应的阻力加速度关于地心距的偏导数，且

R₀＝6371004m，β为大气密度指数模型中的常数；u^(k)表示当前参考轨迹x^(k)对应的中间控制量，f(x^(k),u^(k),E)和h(x^(k))，表示如下：

ω_e为地球自转角速度，ω_e＝7.292115×10^-5rad/s。

4.根据权利要求1所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，热流、过载、动压三个约束线性化后的系数

r_q、r_n计算公式如下：

其中，

q_max和n_max分别为最大热流、最大动压、最大过载；/>

q_k和n_k分别为当前参考轨迹x^(k)确定的热流、动压和过载；r_k为当前参考轨迹x^(k)确定的地心距；β为大气密度指数模型中的常数。

5.根据权利要求4所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，β＝1/7110。

6.根据权利要求1～5之一所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，所述攻角指令α^(k+1)(i)由

进行升阻比反插值得到，倾侧角指令由下式得到：

7.根据权利要求1所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，常值矢量ε中的元素对应小于[5000.20.20.20.5]^T中的元素，其中[5000.20.20.20.5]^T第一项单位为米，后面四项单位为度。

8.一种基于凸优化的再入轨迹三维剖面规划***，其特征在于，包括：

初始化单元，对滑翔段参考轨迹进行初始化，具体执行如下操作：将地心距从当前值r₀到终端值r_f的变化规律假设为线性的，将从当前值r₀到终端值r_f的变化范围均分为N个区间，得到每个离散点上的初始地心距，采用同样方式得到经度、纬度、速度倾角、航迹偏航角、归一化机械能在每个离散点上的初始值，将这些初始值联合在一起记为初始轨迹x⁽⁰⁾，并令参考轨迹x^(k)＝x⁽⁰⁾，上标k表示当前迭代次数；其中，所述归一化机械能包括飞行器的当前归一化机械能E₀、终端归一化机械能E_f；飞行器的当前运动状态x₀＝[r₀λ₀φ₀θ₀σ₀]和预期终端状态x_f＝[r_fλ_fφ_fθ_fσ_f]，r为地心距，λ为经度，φ为纬度，θ为速度倾角，σ为航迹偏航角，下标为0表示当前状态的取值，下标为f表示终端状态的取值；给定离散点数目N，N为正整数；

计算单元，用于根据所述参考轨迹x^(k)，结合凸优化求解器，得到第k+1次迭代的参考轨迹x^(k+1)(i)；

输出单元，用于输出最优轨迹x^(k+1)(i)和最优三维剖面；所述最优三维剖面为第k+1次迭代的攻角指令和倾侧角指令(α^(k+1)(i),υ^(k+1)(i))；

所述计算单元包括：

第一计算模块，用于根据参考轨迹x^(k)计算升阻比下界(u₃)_min和升阻比上界(u₃)_max；基于参考轨迹x^(k)，计算动力学方程的线性化矩阵(A,B,C)、热流、过载、动压三个约束线性化后的系数

r_q、r_n；

离散模块，用于将基于当前参考轨迹x^(k)所得的连续变量(u₃)_min、(u₃)_max、(A,B,C)、

r_q、r_n进行等间距离散，分别得到(u₃)_min、(u₃)_max、(A,B,C)、/>

r_q、r_n在每个离散点处的取值，表示为(u₃)_min(i)、(u₃)_max(i)、(A(i),B(i),C(i))、/>

r_q(i)、r_n(i)，其中，i表示离散点序号，变化范围从1到N；

第二计算模块，用于将离散化后的参数：(u₃)_min、(u₃)_max(i)、(A(i),B(i),C(i))、r_Q(i)、r_q(i)、r_n(i)以及预期终端状态x_f作为凸优化求解器的输入并进行求解，从而得到当前迭代的优化解(x^(k+1)(i),u^(k+1)(i))；

第三计算模块，用于根据(x^(k+1)(i),u^(k+1)(i))计算当前攻角和倾侧角指令，其中，

分别表示在离散点i处的三个控制量分量；

迭代模块，用于判断

是否成立，若成立，则输出x^(k+1)(i)；否则，重复第一计算模块、离散模块、第二计算模块和第三计算模块的操作。/>

Images