试验飞行器Max-Max类再入轨迹优化的序列差凸规划方法

邓泽晓，李骜，刘鲁华，欧朝; DENG Zexiao; LI Ao; LIU Luhua; OU Chao

doi: 10.11887/j.issn.1001-2486.25010026

邓泽晓¹ ，李骜¹ ，刘鲁华¹ ，欧朝^2,3

1. 中山大学航空航天学院, 广东深圳 518107

2. 华中科技大学航空航天学院, 湖北武汉 430074

3. 中国空气动力研究与发展中心空天技术研究所，四川绵阳 621000

基金项目: 跨域飞行交叉技术实验室开放课题资助项目（2024-KF02201）

详细信息

作者简介

邓泽晓（1996—），男，广西北流人，博士研究生，E-mail: dengzx286@gmail.com

通讯作者

刘鲁华（1977—），男，陕西西安人，教授，博士，博士生导师，E-mail: flightdc@mail.sysu.edu.cn；

欧朝（1988—），男，湖南衡阳人，高级工程师，博士，E-mail: ouchao_jg@126.com

中图分类号: V249

文献标识码: A

文章编号: 1001-2486(2026)01-183-13

Max-Max type entry trajectory optimization for testing vehicle by successive difference-of-convex programming

DENG Zexiao¹ ， LI Ao¹ ， LIU Luhua¹ ， OU Chao^2,3

1. School of Aeronautics and Astronautics, Sun Yat-sen University, Shenzhen 518107 , China

2. School of Aerospace Engineering, Huazhong University of Science and Technology, Wuhan 430074 , China

3. Aerospace Technology Institute, China Aerodynamics Research and Development Center, Mianyang 621000 , China

摘要

靶场试验前，利用轨迹优化方法对飞行器的极端边界进行模拟是一项有效仿真手段。采用差凸（difference-of-convex，DC）规划方法，对试验再入飞行器在峰值热流密度方面的极端性能进行了研究。通过DC分解方法对热流、动压和过载约束进行处理，将其分解思路推广到Max-Max类代价函数，如峰值热流、峰值动压和峰值过载等。采用大-M法将原问题转为混合整数非线性规划子问题，将凹凸分解和罚函数结合，解决迭代过程中代价函数的振荡不收敛问题，并提出基于DC松弛模型的改进序列DC规划算法。数值实验表明基于DC松弛模型比传统直接线性化的近似精度要高，且提出的算法具有较高的数值稳定性、算法鲁棒性和代价函数最优性。

关键词

再入轨迹优化 / 试验弹道 / Max-Max类轨迹优化 / DC规划

Abstract

Simulating the extreme boundaries of the entry vehicle using trajectory optimization methods is an effective simulation method before the range test. The DC (difference-of-convex) programming method was employed to study the extreme performance of the testing entry vehicle in terms of peak heat flux. The DC decomposition method was utilized to handle constraints such as heat flux, dynamic pressure, and normal load, and this method was extended to Max-Max type cost functions, such as peak heat flux, peak dynamic pressure, and peak normal load. The Big-M method was adopted to transform the primal problem into a mixed-integer nonlinear programming sub-problem, combining concave-convex decomposition with penalty function technique to address the oscillation and non-convergence issues for the cost function during the iteration process. An improved successive DC programming algorithm based on the DC relaxation model was proposed. Numerical experiments show that the DC relaxation model-based approach has higher approximation accuracy than traditional direct linearization methods, and the proposed algorithm demonstrates high numerical stability, robustness, and optimality of the cost function.

Keywords

entry trajectory optimization / testing trajectory / Max-Max type trajectory optimization / DC programming

1 再入轨迹优化模型 2 模型离散与DC分解 2.1 动力学方程线性化和离散化 2.2 约束的DC分解 2.3 Max-Max类代价函数的DC分解 3 序列差凸规划方法 4 数值算例 4.1 最小终端误差再入轨迹 4.2 最大峰值热流再入轨迹 5 结论

在靶场试验中，为测试飞行器在材料、气动和结构等方面的特性，需要在有限的条件下，如有限射程和有限空域，通过进行力、热和大范围机动等破坏性试验，探索获得飞行器的极端性能边界。利用轨迹优化设计方法对极端边界进行模拟逼近是其中一项有效仿真手段。分析飞行器的最大/最小热流密度、动压和过载轨迹能为新防热材料、气动/结构、大范围机动能力提供设计和测试参考依据。

凸优化是其中一种常见的飞行器轨迹优化方法，包含丰富的研究成果、严密的理论证明以及多项式时间收敛速度，被认为具备在线轨迹生成应用的潜力^[1-4]。Liu等^[4-7]首先将凸优化方法应用于高超声速飞行器再入轨迹优化问题，通过定义新“虚拟”控制量的方式将微分方程转换为仿射控制形式，并给出了附加控制量等式约束的无损凸松弛方法与等价证明。Wang^[8-9]在研究再入飞行器机动性能和结构承载能力指标时发现，再入最大峰值过载轨迹优化不同于传统的轨迹优化问题，也不同于微分对策论的Min-Max/Max-Min代价函数类型的最优控制问题，是一类特殊的Max-Max代价函数类型的最优控制问题，无法直接采用基于梯度的优化方法，进而提出了求解离散事件Max-Max优化问题的序列混合整数凸规划方法，优化得到了最大峰值过载轨迹。祝静娴等^[10]研究传统序列凸优化方法代价函数后期振荡不收敛问题，提出了一种基于罚函数的改进保守凸近似的序列凸优化方法。陈扬等^[11]针对复杂多禁绕飞区绕飞问题，使用整数变量方法重新对其进行建模，提出了一种初值依赖性小、适应能力更广的序列混合整数二阶锥规划方法。更多关于凸优化方法和高超声速飞行器轨迹优化方面应用可以参考文献或综述^[12-15]。

差凸（difference-of-convex，DC）规划理论是优化领域的一个重要分支，它通过变量替代、等价变换和DC分解等数学变换，将代价函数及其约束条件转化为两个凸函数之差的形式，即DC结构，利用这些函数结构为解决航空航天领域中非线性、非凸优化轨迹优化问题提供了新的视角和方法。Misra等^[16]在文献^[17]的基础上，提出了一种基于凹凸程序、平方和规划的序列凸优化的方法，并在包含羽流撞击、阻挡区和避障约束和Clohessy-Hill-Witshire动力学的轨迹优化模型上，验证了该方法可用于相对航天器燃料最优制导问题。Misra的方法本质上就是序列凹凸程序方法在航空航天方面的一种应用。Lu^[18]提出使用凸-凹分解来处理一类非线性等式约束，在松弛问题中该类非线性等式约束由三个或凸或凹的不等式约束表示，并通过在燃料优化的有限推力航天器环绕问题中的应用证明了该方法的有效性。Doff-Sotta等^[19]将凸优化、DC分解和管状非线性模型预测控制等方法结合提出了序列鲁棒优化方法，并将其应用于倾斜旋翼垂直起降飞行器过渡轨迹的优化，获得了鲁棒的最优轨迹。Zhao和Shang^[20]结合饱和函数和罚函数技术来研究航天器编队重构问题，扩展了DC优化问题的可行域，避免了传统序列凸优化方法常见的人为不可行性问题，提高了算法的鲁棒性和收敛性。

与传统的再入轨迹优化问题不同，试验再入飞行器的轨迹设计有着极高的力、热等极端环境模拟需求，以确保试验考核的充分性。本文采用DC规划方法，对试验再入飞行器在力热方面的固有边界进行深入研究，以提高飞行器探索复杂环境的能力。为提高对热流密度、动压和过载等过程约束边界逼近精度，采用凹凸程序、变量替代和罚函数方法对热流、动压、过载、阻力等复杂约束进行DC松弛逼近。引入二值变量，运用DC松弛、大-M法等方法将最大峰值热流问题转化为混合整数非线性规划子问题形式，并提出基于DC松弛改进的序列DC规划方法，用于Max-Max类再入轨迹优化问题的求解。最后以数值算例验证提出方法的精度、稳定性以及最优性。

1 再入轨迹优化模型

由于大气阻力的作用，飞行器能量是单调递减的，为降低优化变量的数目，本文采用E作为自变量，忽略地球自转。

E = \frac{V^{2}}{2} - \frac{μ}{r}

(1)

\frac{d E}{d τ} = V \frac{d V}{d τ} + \frac{μ d r}{r^{2} d τ} \approx - D V

(2)

其中：μ为地球引力常数；r为地心距；V为速度；τ为无量纲时间。以能量为自变量的再入质心动力学模型^[21]为

\{\begin{array}{l} \frac{d r}{d E} = - \frac{\sin θ}{D} \\ \frac{d λ}{d E} = - \frac{\cos θ \sin σ}{D r \cos ϕ} \\ \frac{d ϕ}{d E} = - \frac{\cos θ \cos σ}{D r} \\ \frac{d θ}{d E} = - \frac{L \cos ν}{D V^{2}} + (\frac{1}{r^{2}} - \frac{V^{2}}{r}) \frac{\cos θ}{D V^{2}} + C_{θ} + C_{θ}^{'} \\ \frac{d σ}{d E} = - \frac{L \sin ν}{D V^{2} \cos θ} - \frac{\tan ϕ \cos θ \sin σ}{D r} + C_{σ} + C_{σ}^{'} \end{array}

(3)

其中：λ为经度；φ为地心纬度；θ为速度倾角；σ为速度方位角；ν为倾侧角。L和D分别为升力和阻力加速度，

\{\begin{matrix} L = R_{0} ρ V^{2} S_{r e f} C_{L} / (2 M_{r e f}) \\ D = R_{0} ρ V^{2} S_{r e f} C_{D} / (2 M_{r e f}) \end{matrix}

式中，M_ref为质量，S_ref为参考面积，C_L和C_D分别为升力和阻力系数。大气模型采用指数大气，ρ=ρ₀exp（-βh），h=rR₀-R₀为飞行高度，R₀为地球平均半径，β=1/h_Scale，h_Scale为标称高度，ρ₀为参考密度。C_σ、C_θ、

C_{σ}^{'}

和

C_{θ}^{'}

分别为哥氏加速度项及牵连加速度项。

\begin{matrix} C_{θ} = - \frac{2 ω_{e}}{D V} s i n σ c o s ϕ \\ C_{θ}^{'} = - \frac{ω_{e}^{2} r}{D V^{2}} c o s ϕ (c o s ϕ c o s θ + s i n ϕ c o s σ s i n θ) \\ C_{σ} = - \frac{2 ω_{e}}{D V} (s i n ϕ - c o s σ t a n θ c o s ϕ) \\ C_{σ}^{'} = - \frac{ω_{e}^{2} r c o s ϕ s i n ϕ s i n σ}{D V^{2} c o s θ} \end{matrix}

其中，ω_e为地球自转角速度。

直接采用攻角α和倾侧角ν作为控制变量会增加后续动力学方程线性化的复杂性，因此定义新的控制变量如下：

\begin{matrix} u = {[\begin{matrix} u_{1} & u_{2} & u_{3} \end{matrix}]}^{T} \\ u_{1} = χ c o s ν, u_{2} = χ s i n ν, u_{3} = χ \end{matrix}

式中，χ=L/D是升阻比，为攻角和马赫数的拟合函数，χ=f（α，Ma）。

攻角和倾侧角的幅值约束转换为

χ_{m i n} ⩽ u_{3} ⩽ χ_{m a x}

(4)

u_{3} s i n ν_{m i n} ⩽ u_{2} ⩽ u_{3} s i n ν_{m a x}

(5)

其中：χ_min和χ_max分别表示最小和最大升阻比；ν_min和ν_max分别表示最小和最大倾侧角。

在三个虚拟控制变量u₁、u₂和u₃之间隐含着非线性等式关系，为附加控制量等式约束：

u_{1}^{2} + u_{2}^{2} = u_{3}^{2}

(6)

过程约束包括热流密度、动压、法向过载和最小阻力加速度：

\dot{Q} = k_{Q} \sqrt{ρ} V^{3.15} V_{Scale}^{3.15} ⩽ {\dot{Q}}_{m a x}

(7)

q = \frac{1}{2} ρ V^{2} V_{Scale}^{2} ⩽ q_{m a x}

(8)

n = \sqrt{L^{2} + D^{2}} ⩽ n_{m a x}

(9)

D ⩾ D_{m i n}

(10)

其中：

\dot{Q}

、q和n分别是热流密度、动压和法向过载；k_Q为飞行器材料参数；V_Scale为速度无量纲因子。值得说明的是，最小阻力加速度是一个“软约束”，在飞行器初始下降时可以不满足条件。

以E为自变量的初始点和终端点边界条件为

\{\begin{matrix} x (E_{0}) = x_{0} \\ x (E_{f}) = x_{f} \end{matrix}

(11)

再入轨迹优化问题的代价函数是为满足特定的任务目标而人为设定的，可以分为三个不同的类别：积分型、终端型和复合型。代价函数的一般形式如下所示：

J (x, u, t) = Φ (x_{f}, t_{f}) + \int_{t_{0}}^{t_{f}} L (x, u, t) d t

(12)

以能量为自变量时式中时间t可替换为能量E。最小化终端误差代价函数为

J : = |λ_{f} - λ_{f}^{*}| + |ϕ_{f} - ϕ_{f}^{*}|

（13

为了评估升力体飞行器的力和热性能，并为热和结构设计部门提供设计指导，代价函数可以定制为Max-Max类型，如代价函数定义为最大头部驻点处的峰值热流密度、最大峰值动压、最大峰值过载：

J : = {m a x m a x \dot{Q}, m a x m a x q, m a x m a x n}

(14)

由于无法提前确定最大峰值点发生的能量状态，Max-Max类型代价函数与传统的积分型或终端型代价函数不同，不能直接使用变分法或庞特里亚金最大原理求解此类问题。

2 模型离散与DC分解

附加控制量等式约束和热流、动压、过载、最小阻力等约束都是非凸非线性的，给再入轨迹优化问题的求解带来了极大的挑战。利用约束函数的DC性质，通过变量替代、等价变换和DC分解等数学变换，对非凸非线性方程约束进行处理，得到易于序列求解的形式。

2.1 动力学方程线性化和离散化

动力学方程（3）关于新的控制变量重新表述为如下仿射形式：

x^{'} = f (x, u) + h (x)

(15)

式中：x为状态向量，x=[r λ φ θ σ]^T；

f (x, u) = [\begin{matrix} - \frac{\sin θ}{D} \\ - \frac{c o s θ s i n σ}{D r c o s ϕ} \\ - \frac{c o s θ c o s σ}{D r} \\ - \frac{u_{1}}{D V^{2}} + (\frac{μ}{r^{2}} - \frac{V^{2}}{r}) \frac{\cos θ}{D V^{2}} \\ - \frac{u_{2}}{D V^{2} c o s θ} - \frac{t a n ϕ c o s θ s i n σ}{D r} \end{matrix}]

h (x) = {[\begin{matrix} 0 0 0 C_{θ} + C_{θ}^{'} C_{σ} + C_{σ}^{'} \end{matrix}]}^{T}

在参考轨迹x^（^k^）处作泰勒展开：

x^{'} = A x + B u + C

(16)

式中，

\begin{matrix} A = {\frac{\partial f (x, u)}{\partial x}|}_{x = x^{(k)}} \\ B = {\frac{\partial f (x, u)}{\partial u}|}_{u = u^{(k)}} \\ C = f (x^{(k)}, u^{(k)}) - {\frac{\partial f (x, u)}{\partial x}|}_{x = x^{(k)}} x^{(k)} - \\ {\frac{\partial f (x, u)}{\partial u}|}_{u = u^{(k)}} u^{(k)} + h (x^{(k)}) \end{matrix}

将能量、状态变量、控制变量分别等距离散为E₁，E₂，···，E_N_-1，E_N；x1，x2，···，x_N_-1，x_N；u₁，u₂，···，u_N_-1，u_N。对式（16）应用梯形规则：

\begin{matrix} x_{j + 1} - x_{j} = \frac{Δ E}{2} (A_{j} x_{j} + B_{j} u_{j} + C_{j} + \\ A_{j + 1} x_{j + 1} + B_{j + 1} u_{j + 1} + C_{j + 1}) \end{matrix}

(17)

式中：ΔE=（E_f-E₀）/N；下标j表示离散点的索引号，j=1，2，···，N-1。整理上述方程得

L_{j} x_{j} - L_{j + 1} x_{j + 1} + M_{j} u_{j} + M_{j + 1} u_{j + 1} + F_{j} = 0

(18)

式中，

L_{j} = I + \frac{Δ E}{2} A_{j} ， L_{j + 1} = I - \frac{Δ E}{2} A_{j + 1} ， M_{j} = \frac{Δ E}{2} B_{j} ， M_{j + 1} = \frac{Δ E}{2} B_{j + 1} ， F_{j} = - \frac{Δ E}{2} (C_{j} + C_{j + 1})

。

写成线性形式表达如下：

P z + f = 0

(19)

式中，z为决策变量，

z = {[x_{1}^{T} ， \dots ， x_{N}^{T} ， u_{1}^{T} ， \dots ， u_{N}^{T}]}^{T}

。

为了降低线性化带来的误差和防止“人为无界性”^[6，22]，引入信赖域约束：

‖x - x^{(k)}‖ ⩽ \tilde{δ}

(20)

式中，为信赖域半径。

热流密度、动压、法向过载和最小阻力加速度约束是非线性的，可通过反解方法^[6]或线性化方法^[8-9]得到线性约束形式：

f_{D} ⩽ r ⩽ m a x \{f_{\dot{Q}}, f_{q}, f_{n}\}

(21)

式中，

f_{\dot{Q}} ， f_{q} ， f_{n}

和

f_{D}

由热流密度、动压、法向载荷和最小阻力加速度表达式反解得到。

以最大峰值热流为代价函数时，放宽热流密度约束，即

f_{D} ⩽ r ⩽ m a x \{f_{q}, f_{n}\}

(22)

以最大峰值动压和过载为代价函数时，分别放宽动压和过载约束。

2.2 约束的DC分解

附加控制变量约束式（6）等价的约束形式为

u_{1}^{2} + u_{2}^{2} - u_{3}^{2} ⩽ 0

(23)

u_{3}^{2} - (u_{1}^{2} + u_{2}^{2}) ⩽ 0

(24)

引入一个松弛变量ξ，可以得到约束条件

u_{3}^{2} - ξ ⩽ 0

(25)

ξ - (u_{1}^{2} + u_{2}^{2}) ⩽ 0

(26)

其中，

ξ \in [ξ_{m i n} ， ξ_{m a x}] ， ξ_{m i n} = m i n u_{3}^{2} ， ξ_{m a x} = m a x u_{3}^{2}

。式（23）和式（25）为凸不等式约束，而式（26）为凹不等式约束。为保证凹方程的非负特性，引入非负松弛变量η，可以得到

η + ξ - (u_{1}^{2} + u_{2}^{2}) ⩽ 0

(27)

式中，

η \in [0 ， η_{m a x}] ， η_{m a x} ⩾ m a x (u_{1}^{2} + u_{2}^{2}) - ξ_{m i n}

。进一步，对其进行泰勒展开

η + ξ ⩽ {(u_{1}^{(k)})}^{2} + {(u_{2}^{(k)})}^{2} + 2 [\begin{matrix} u_{1}^{(k)} u_{2}^{(k)} \end{matrix}] [\begin{matrix} u_{1} - u_{1}^{(k)} \\ u_{2} - u_{2}^{(k)} \end{matrix}]

(28)

在代价函数中添加松弛变量对应的罚项，

J^{'} (z) = J (x, u, E) + C_{η} η

(29)

式中，C_η为一正常数。上述分解过程称为凹凸程序分解^[23]或DC分解。

热流、动压、过载和阻力等约束为非凸非线性过程约束，直接使用序列线性化容易引起代价函数的振荡不收敛问题，为此本文提出使用DC分解的方式对热流、动压、过载和阻力等约束作DC松弛逼近。

将热流约束表达式（7）作以下变换：

e x p [- \frac{β R_{0} r}{2} + \frac{β R_{0}}{2} + l n (\frac{k_{Q} V_{Scale}^{3.15} \sqrt{ρ_{0}}}{{\dot{Q}}_{max}})] V^{3.15} ⩽ 1

(30)

引入替代变量

κ_{\dot{Q}}

满足：

e x p (κ_{\dot{Q}}) ⩾ V = \sqrt{2 (E + \frac{μ}{r})}

(31)

式中，

2 (E + \frac{μ}{r}) \in （ 0，1 ） ， κ_{\dot{Q}}

表示飞行速度对数的上界，

κ_{\dot{Q}} \in [l n V_{m a x} ， l n V_{m i n}] ， κ_{\dot{Q}} < 0

。热流约束不等式（30）等效于如下形式：

\begin{matrix} e x p [- \frac{β R_{0} r}{2} + \frac{β R_{0}}{2} + l n (\frac{k_{Q} V_{Scale}^{3.15} \sqrt{ρ_{0}}}{{\dot{Q}}_{m a x}})] e x p (3.15 κ_{\dot{Q}}) ⩽ 1 \\ \Leftrightarrow (1,1, - \frac{β R_{0} r}{2} + 3.15 κ_{\dot{Q}} + \frac{β R_{0}}{2} + l n (\frac{k_{Q} V_{Scale}^{3.15} \sqrt{ρ_{0}}}{{\dot{Q}}_{m a x}})) \in K_{e x p} \end{matrix}

(32)

式中，

K_{e x p}

表示指数锥约束^[24]。

对式（31）作进一步处理

2 E + \frac{2 μ}{r} - e x p (2 κ_{\dot{Q}}) ⩽ 0

(33)

引入替代变量

ι_{\dot{Q}}

，不等式（33）等效于

2 E + ι_{\dot{Q}} - e x p (2 κ_{\dot{Q}}) ⩽ 0

(34)

\frac{2 μ}{r} ⩽ ι_{\dot{Q}} \Leftrightarrow (r, ι_{\dot{Q}}, 2 \sqrt{μ}) \in Q_{r}^{3}

(35)

其中：

ι_{\dot{Q}}

表示2倍重力势能的上边界，

ι_{\dot{Q}} \in [\frac{2 μ}{r_{m a x}} ， \frac{2 μ}{r_{m i n}}] ， ι_{\dot{Q}} ⩾ 0; Q_{r}^{3}

为3维旋转二阶锥。对式（34）取反可得

\begin{matrix} 2 E + ι_{\dot{Q}} - e x p (2 κ_{\dot{Q}}) ⩾ 0 \\ \Leftrightarrow (2 E + ι_{\dot{Q}}, 1,2 κ_{\dot{Q}}) \in K_{e x p} \end{matrix}

(36)

显然，式（36）是一个凸锥，故式（34）为一凹不等式约束。对于再入轨迹规划问题，迭代求解困难主要来源于约束的非凸部分。通过线性化方法可得约束式（34）的凸逼近形式，但直接对非凸部分进行泰勒展开往往是迭代不可行的，因此给出一种可行的松弛近似形式：

2 E + ι_{\dot{Q}}^{(k)} - [e x p (2 κ_{\dot{Q}}^{(k)}) + 2 e x p (2 κ_{\dot{Q}}^{(k)}) (κ_{\dot{Q}} - κ_{\dot{Q}}^{(k)})] ⩽ δ

(37)

式中，δ为松弛变量，表示松弛约束式（37）与原始约束式（34）之间的线性化误差，δ≥0。并在代价函数中添加松弛变量对应的罚项

J^{''} (z) = J^{'} (z) + C_{δ} δ

(38)

式中，C_δ为一正常数。

对动压、过载和最小阻力加速度处理的方式与热流相类似，下面直接给出约束的DC松弛形式。

引入松弛变量κ_q和替代变量ι_q，动压约束式（8）的DC松弛形式为

\{\begin{array}{l} (1, 1, - β R_{0} r + 2 κ_{q} + β R_{0} + \ln (\frac{ρ_{0} V_{Scale}^{2}}{2 q_{max}})) \in K_{\exp} \\ 2 E + ι_{q}^{(k)} - [\exp (2 κ_{q}^{(k)}) + 2 \exp (2 κ_{q}^{(k)}) (κ_{q} - κ_{q}^{(k)})] ⩽ δ \\ (r, ι_{q}, 2 \sqrt{μ}) \in Q_{r}^{3} \end{array}

(39)

式中，κ_q表示飞行速度对数的上界，κ_q∈[lnV_max，lnV_min]，κ_q＜0；ι_q表示2倍重力势能的上边界，

ι_{q} \in [\frac{2 μ}{r_{m a x}} ， \frac{2 μ}{r_{m i n}}] ， ι_{q} ⩾ 0

。

由于大气和飞行器气动外形的复杂性，气动力与飞行状态紧耦合，过载约束对飞行器路径的限制是极其复杂的。假设气动力系数和升阻比随时间变化较小，可用上一步迭代结果近似代替。引入松弛变量κ_n和替代变量ι_n，过载约束式（9）的DC松弛形式为

\{\begin{array}{l} (1, 1, - β R_{0} r + 2 κ_{n} + β R_{0} + \ln (\frac{ρ_{0} R_{0} S_{ref} C_{D} \sqrt{1 + χ^{2}}}{2 M_{ref} n_{max}})) \in K_{\exp} \\ 2 E + ι_{n}^{(k)} - [\exp (2 κ_{n}^{(k)}) + 2 \exp (2 κ_{n}^{(k)}) (κ_{n} - κ_{n}^{(k)})] ⩽ δ \\ (r, ι_{n}, 2 \sqrt{μ}) \in Q_{r}^{3} \end{array}

(40)

式中，κ_n表示飞行速度对数的上界，κ_n∈[lnV_max，lnV_min]，κ_n＜0；ι_n表示2倍重力势能的上边界，

ι_{n} \in [\frac{2 μ}{r_{m a x}} ， \frac{2 μ}{r_{m i n}}] ， ι_{n} ⩾ 0

。

引入松弛变量κ_D和替代变量ι_D，最小阻力加速度约束式（10）的DC松弛形式为

\{\begin{array}{l} (1, 1, β R_{0} r - 2 κ_{D} - β R_{0} - \ln (\frac{ρ_{0} R_{0} S_{ref} C_{D}}{2 M_{ref} D_{min}})) \in K_{\exp} \\ (2 E + ι_{D}, 1, 2 κ_{D}) \in K_{\exp} \\ (r^{(k)} ι_{D}^{(k)} - 2 μ) + [ι_{D}^{(k)} r^{(k)}] {[r - r^{(k)} ι_{D} - ι_{D}^{(k)}]}^{T} ⩽ δ \end{array}

(41)

式中，κ_D表示飞行速度对数的上界，κ_D∈[lnV_max，lnV_min]，κ_D＜0；ι_D表示2倍重力势能的上边界，

ι_{D} \in [\frac{2 μ}{r_{m a x}} ， \frac{2 μ}{r_{m i n}}] ， ι_{D} ⩾ 0

。

2.3 Max-Max类代价函数的DC分解

以最大峰值热流作为代价函数进行分析，首先需要求解得到一条以最大峰值热流为代价函数的热流剖面，再从中计算出最大峰值热流的数值大小和对应的飞行器能量状态，即可求出最大峰值热流轨迹。而实际情况是事先无法确定飞行器到达峰值热流点的能量状态，无法直接给出以最大峰值热流作为代价函数的显式表示形式，更无法确定含峰值热流点的热流剖面。

换种角度思考模型的近似形式，采用离散形式代替。事先确定若干有限个离散状态点，然后求解该点出现最大热流的热流轨迹剖面，再比较若干有限个点对应的最大峰值热流数值，选择最大热流数值最大的那一条，即标记出当前离散问题最大峰值热流点和最大峰值热流点对应的最大峰值热流剖面。只要离散点足够密集，则离散问题就足够逼近原问题。

近似代价函数如下：

\begin{matrix} J = \underset{i = 1, \dots, N}{m a x} \underset{x \in ε}{m a x} \{{\dot{Q}}_{i}\} \\ = \underset{x \in ε}{m a x} \underset{i = 1, \dots, N}{m a x} \{{\dot{Q}}_{i}\} \end{matrix}

(42)

式中：

ε

为所有约束集合的交集，称为可行解集；N为离散节点个数。

引入无量纲松弛变量

Z_{\dot{Q}}

，

Z_{\dot{Q}} = \underset{i = 1, \dots, N}{m a x} \{\frac{{\dot{Q}}_{i}}{{\dot{Q}}_{Scale}}\}

(43)

式中，

Z_{\dot{Q}}

表示最大峰值热流点的峰值热流数值，

{\dot{Q}}_{Scale}

为热流密度无量纲因子。

引入二值变量y_i，当y_i=1时，表示在当前指标i的状态点上，热流取得最大值。应用大-M法^[25]，式（43）可等效于

\frac{{\dot{Q}}_{i}}{{\dot{Q}}_{Scale}} ⩽ Z_{\dot{Q}} ⩽ \frac{{\dot{Q}}_{i}}{{\dot{Q}}_{Scale}} + M (1 - y_{i})

(44)

\sum_{i = 1}^{N} y_{i} = 1

(45)

y_{i} \in {0,1}

(46)

其中：M是一个充分大的正数；i=1，···，N。

式（44）是一个非凸非线性不等式约束，需要作进一步处理。引入替代变量

κ_{\dot{Q}}

，式（44）的左边等效于

\begin{matrix} e x p [- \frac{β R_{0} r}{2} + \frac{β R_{0}}{2} + l n (\frac{k_{Q} V_{Scale}^{3.15} \sqrt{ρ_{0}}}{{\dot{Q}}_{Scale}})] e x p (3.15 κ_{\dot{Q}}) ⩽ Z_{\dot{Q}} \\ \Leftrightarrow (Z_{\dot{Q}}, 1, - \frac{β R_{0} r}{2} + 3.15 κ_{\dot{Q}} + \frac{β R_{0}}{2} + l n (\frac{k_{Q} V_{Scale}^{3.15} \sqrt{ρ_{0}}}{{\dot{Q}}_{Scale}})) \in K_{e x p} \end{matrix}

(47)

由于不等式（44）右边包含一个充分大的正实数M，对它进行DC分解操作，把一个充分大的正数值放入指数锥内，容易导致“病态”的数值问题，因此对其直接采用线性化处理：

Z_{\dot{Q}} - M (1 - y_{i}) - [\frac{{\dot{Q}}_{i}^{(k)}}{{\dot{Q}}_{Scale}} + \frac{{\dot{Q}}_{i}^{(k)}}{{\dot{Q}}_{Scale}} (\frac{0.5 ρ_{r}}{ρ} + \frac{3.15 V_{r}}{V}) (r_{i} - r_{i}^{(k)})] ⩽ δ

(48)

式中，ρ_r、V_r分别为大气密度、速度对地心距的导数。

由此得到基于DC松弛的最小终端误差再入轨迹优化问题模型：

\begin{matrix} m i n J = C_{σ} Δ E [E_{2, (N, :)} F_{σ}^{'} 0_{1 \times 3 N}] z / 2 + \\ C_{T E} (|λ_{f} - λ_{f}^{*}| + |ϕ_{f} - ϕ_{f}^{*}|) + C_{η} η + C_{δ} δ \end{matrix}

(49)

s.t. 式（4）~（5）、式（11）、式（19）~（20）、式（23）、式（25）、式（28）、式（32）、式（35）、式（37）、式（39）~（41）

式中，C_TE为一正常数。

2 个混合整数非线性形式的最大峰值热流密度再入轨迹优化问题模型：

1）基于线性化泰勒展开的最大峰值热流密度再入轨迹优化问题模型（MMP-Ⅰ）：

\begin{matrix} m i n J = C_{σ} Δ E [E_{2, (N, :)} F_{σ}^{'} 0_{1 \times 3 N}] z / 2 + \\ C_{Z_{\dot{Q}}} Z_{\dot{Q}} + C_{η} η + C_{δ} δ \end{matrix}

(50)

s.t. 式（4）~（5）、式（11）、式（19）~（20）、式（22）~（23）、式（25）、式（28）、式（39）、式（48）

2）基于DC松弛的最大峰值热流密度再入轨迹优化问题模型（MMP-Ⅱ）：

\begin{matrix} m i n J = C_{σ} Δ E [E_{2, (N, :)} F_{σ}^{'} 0_{1 \times 3 N}] z / 2 + \\ C_{Z_{\dot{Q}}} Z_{\dot{Q}} + C_{η} η + C_{δ} δ \end{matrix}

(51)

s.t. 式（4）~（5）、式（11）、式（19）~（20）、式（23）、式（25）、式（28）、式（39）~（41）、式（47）~（48）

其中，

\begin{matrix} C_{Z_{\dot{q}}} = {\dot{Q}}_{Scale} \\ E_{1} = {[\begin{matrix} 0 & 0 & 0 & \dots & 0 & 0 \\ 1 & 1 & 0 & \dots & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & 2 & 2 & \dots & 1 & 0 \\ 1 & 2 & 2 & \dots & 2 & 1 \end{matrix}]}_{N \times N} \\ E_{2} = E_{1} \otimes 1_{1 \times 5} \\ f_{σ} = {[\begin{matrix} σ_{1} & σ_{2} & \dots & σ_{N} \end{matrix}]}^{T} \\ F_{σ}^{'} = d i a g \{{(f_{σ}^{'})}_{1} {(f_{σ}^{'})}_{2} \dots {(f_{σ}^{'})}_{N}\} \end{matrix}

同理，推广到最大峰值动压和最大峰值过载等Max-Max类轨迹优化设计问题，这里不再赘述。

3 序列差凸规划方法

序列差凸规划方法（successive difference-of-convex programming method，SDCP）的算法流程如图1所示。因核心步骤为凹凸分解，故此又称序列凹凸程序方法（successive concave convex procedure method，SCCP）。

图1序列DC规划方法流程图

Fig.1Flowchart of SDCP method

输入：计算模型和求解控制参数。

输出：状态变量x，控制变量u。

步骤1：初始化再入轨迹x₀，u₀，η₀，k=1。

步骤2：如果k≤k_maxIter，转到步骤3，否则，算法因到达最大迭代步数而停止。

步骤3：求解动力学模型。

步骤4：使用文献^[21]方法获得式（4）和式（5）的边界。

步骤5：添加再入轨迹优化模型的多种约束，包括离散动力学约束、过程约束、控制量约束、热流约束和边界约束等。

步骤6：使用凸优化算法求解式（49）或使用混合整数非线性规划算法求解模型MMP-Ⅰ和MMP-Ⅱ得到（x^（^k^+1），u^（^k^+1））。

步骤7：如果

‖x^{(k+ 1)} - x^{(k)}‖ ⩽ ϵ

，算法停止，输出（x^（^k^+1），u^（^k^+1）），否则转入步骤8。

步骤8：更新迭代，x^（^k^）←x^（^k^）+Δ（x^（^k^+1）-x^（^k^）），u^（^k^）←u^（^k^+1），k←k+1，转到步骤2。

计算模型和求解控制参数包括：

1）地球模型参数R₀、g₀、ω_e和μ等；

2）指数大气模型参数ρ₀、h_Scale和β等；

3）飞行器特征参数M_ref、S_ref、k_Q、C_L和C_D等；

4）约束参数

{\dot{Q}}_{m a x} 、 q_{m a x} 、 n_{m a x}

和D_min等。

4 数值算例

使用再入飞行器CAV-H来展示所提出方法的性能，结构和气动参数如文献^[26]描述。采用中国2000大地坐标系作为地球参考坐标系。编程遵循CVX2.2的定制化规则编写，并调用锥规划求解器MOSEK解混合整数非线性规划问题。

参数设置如下：节点个数为

N = 100

，最大迭代步数

k_{m a x I t e r} = 50

，迭代步长

Δ = 0.5

。停止阈值设置为

ϵ = {[\frac{200}{R_{0}} ， \frac{0.05 π}{180} ， \frac{0.05 π}{180} ， \frac{0.05 π}{180} ， \frac{0.05 π}{180}]}^{T} ，

再入任务参数如表1所示，包括初始和终端高度、纬度、经度、速度、速度倾角和速度方位角。决策变量边界幅值如表2所示，包括地心距、经纬度等状态变量、虚拟控制变量、松弛变量和替代变量等。过程约束设置如下：

{\dot{Q}}_{m a x}

=3 900 kW/m²，q_max=100 kPa，n_max=3.5g₀，D_min=0.05g₀。控制变量的振幅边界条件设置为[α_min，α_max]=[10°，20°]，[ν_min，ν_max]=[-85°，85°]，[

{\dot{α}}_{m i n} ， {\dot{α}}_{m a x}

]=[-0.01（°）/s，0.01（°）/s]和[

{\dot{ν}}_{min} ， {\dot{ν}}_{max}

]=[-0.5（°）/s，0.5（°）/s]。其他的终端条件设置为

r_{N}^{*}

=25 km、

θ_{N}^{*}

=-4°和

σ_{N}^{*}

∈[80°，100°]。

算法初始化：在再入过程的初始阶段，再入飞行器上的热流影响显著，调整攻角对飞行器的气动控制效果被认为可以忽略不计。使用固定的大攻角α=20°和固定的倾侧角ν=0°作为飞行控制变量，直到速度倾角首次达到零为初始下降阶段（initial descending phase，IDP）。初始轨迹通过从起滑点到终端点的线性插值得到^[6]：

u_{3} = χ_{m a x} \cdot 1_{N \times 1} ， u_{1} = u_{3} c o s （ 0 ） ， u_{2} = u_{3} s i n （ 0 ） ， κ_{\dot{Q}} = κ_{\dot{Q} ， m i n} \cdot 1_{N \times 1} ， κ_{q} = κ_{q ， m i n} \cdot 1_{N \times 1} ， κ_{n} = κ_{n ， m i n} \cdot 1_{N \times 1} ， κ_{D} = κ_{D ， m i n} \cdot 1_{N \times 1} ， ι_{\dot{Q}} = ι_{\dot{Q} ， m i n} \cdot 1_{N \times 1} ， ι_{q} = ι_{q ， m i n} \cdot 1_{N \times 1} ， ι_{n} = ι_{n ， m i n} \cdot 1_{N \times 1} ， ι_{D} = ι_{D ， m i n} \cdot 1_{N \times 1} ， η = η_{m i n} ， δ = δ_{m i n} ， Z = 1

。

表1初始和终端条件

Tab.1Initial and terminal conditions

表2决策变量的边界约束

Tab.2Box constraints of the decision variables

采用文献^[6]的序列凸规划方法（successive convex programming method，SCP）作为参考对比方法。通过优化得到攻角和倾侧角控制指令，输入动力学式（3）进行仿真，作为比较验证组，记作SIMU。建立迭代步数（iter.）、CPU时间、经纬度相对误差ε_{rel_gap1}、高相对误差ε_{rel_gap2}、球面相对误差ε_gap3和终端点相对距离误差ε_f作为算法的计算效率和计算精度评价指标，其中，

\begin{matrix} ε_{rel_gap1} = \frac{d i s t [T_{SIMU} (λ_{f}, ϕ_{f}) - T_{OPT} (λ_{f}, ϕ_{f})]}{s_{SIMU}} \times 100 % \\ ε_{rel_gap 2} = \frac{‖h_{f, SIMU} - h_{f, OPT}‖}{s_{SIMU}} \times 100 % \end{matrix}

\begin{matrix} {\{ε_{{gap}_{3}}\}}_{k m} = d i s t [T_{OPT} (λ_{f}, ϕ_{f}) - T_{DES} (λ_{f}, ϕ_{f})] \\ ε_{f} = \frac{‖λ_{N}^{*} - λ_{f}‖ + ‖ϕ_{N}^{*} - ϕ_{f}‖}{‖λ_{f}‖ + ‖ϕ_{f}‖} \times 100 % \end{matrix}

式中，dist（·）表示计算球面距离的函数，T_DES（λ_f，φ_f）表示设计终端点，T_SIMU（λ_f，φ_f）表示SIMU仿真的终端点，T_OPT（λ_f，φ_f）表示通过不同优化方法得到的终端点，s_SIMU表示仿真航程。ε_{rel_gap1}和ε_{rel_gap2}反映了为物理飞行设计的控制变量与实际使用的控制变量之间的差异，值越小，表明控制变量的设计越精确，控制量可执行度越高。ε_gap3表示优化终端点与设计终端点之间的球面距离误差。

4.1 最小终端误差再入轨迹

算例一：设计最小终端误差轨迹算例验证过程约束的DC分解的有效性。设置C_TE=1×10³、C_η=1×10^-5、C_δ=1×10²和C_σ=1×10²。

数值仿真结果如表3所示。图2~11给出了由SCP和SDCP设计的最小终端误差再入轨迹剖面。图2和图3为设计的高度-速度走廊，并展示了迭代收敛的过程。图4为经度-纬度-高度曲线，图5~7为热流密度、动压和过载曲线。图8和图9为由SCP和SDCP设计的攻角和倾侧角控制量剖面。图10和图11为速度倾角和速度方位角曲线。

表3算例一的结果比较

Tab.3Comparison results in example1

图2高度-速度剖面（SCP）

Fig.2Altitude-vs-velocity histories by SCP

图3高度-速度剖面（SDCP）

Fig.3Altitude-vs-velocity histories by SDCP

图4经度-纬度-高度曲线（终端误差指标）

Fig.4Longitude-latitude-altitude footprints (terminal error index)

图5热流密度曲线（终端误差指标）

Fig.5Heat flux density histories (terminal error index)

由图2、图3及表3可知，基于DC松弛模型的SDCP迭代步数要比基于线性化模型的SCP要少，在ε_{rel_gap1}和ε_{rel_gap2}性能指标上表现要好，表明SDCP设计控制变量要比SCP更容易插值还原，可执行度更高。但在终端误差表现方面，即ε_gap3，SCP终端误差在分米级别以下，而SDCP为分米级别，SCP要比SDCP更优。由图4~11可以看出，SCP与SDCP设计的再入轨迹极其接近，而攻角剖面略微有点差异。由图9可以看出，飞行器在能量为-0.82和-0.95附近发生了倾侧反转。从计算效率方面看，由表3可知，SCP比SDCP表现要好，这是由于DC松弛逼近过程约束的方法增加了决策变量的规模和复杂凸约束（如指数锥、旋转二阶锥）的个数。

图6动压曲线（终端误差指标）

Fig.6Dynamic pressure histories (terminal error index)

图7过载曲线（终端误差指标）

Fig.7Normal load histories (terminal error index)

图8攻角曲线（终端误差指标）

Fig.8Angle of attack profiles (terminal error index)

图9倾侧角曲线（终端误差指标）

Fig.9Bank angle profiles (terminal error index)

图10速度倾角曲线（终端误差指标）

Fig.10Flight-path angle histories (terminal error index)

图11速度方位角曲线（终端误差指标）

Fig.11Heading angle histories (terminal error index)

4.2 最大峰值热流再入轨迹

算例二：设置

C_{Z \dot{Q}}

=1×10⁶、C_TE=1×10⁵、C_σ=1×10²、C_δ=1×10²、M=1×10¹⁰和C_η=1×10^-5。采用文献^[8-9]方法作为对比，并标记为MMP-Ⅰ。改变不同节点个数，以显示算法在数值稳定性、收敛特性和最优性方面的表现。

统计结果如表4所示，迭代过程如图12所示。图13为节点N=200时，运用模型MMP-Ⅱ计算得到的高度-速度再入走廊，用以展示迭代收敛的过程。图14~21为节点N=100时，运用模型MMP-Ⅰ和模型MMP-Ⅱ计算得到的再入轨迹剖面。图14为经度-纬度-高度曲线，图15~17为热流密度、动压和过载曲线，图18和图19为设计的攻角和倾侧角剖面，图20和图21描述了速度倾角和速度方位角的变化情况。

表4算例二的结果比较

Tab.4Comparison results in example2

~: 迭代停止，到达最大迭代步数或最大运行时间。

图12代价函数迭代过程

Fig.12Cost function iteration process

图13高度-速度剖面（N=200）

Fig.13Altitude-vs-velocity histories (N=200)

图14经度-纬度-高度曲线（峰值热流指标）

Fig.14Longitude-latitude-altitude footprints (peak heat flux index)

图15热流密度曲线（峰值热流指标）

Fig.15Heat flux density histories (peak heat flux index)

图16动压曲线（峰值热流指标）

Fig.16Dynamic pressure histories (peak heat flux index)

图17过载曲线（峰值热流指标）

Fig.17Normal load histories (peak heat flux index)

图18攻角曲线（峰值热流指标）

Fig.18Angle of attack profiles (peak heat flux index)

在计算精度方面，比较表4中的数据，求解基于DC松弛的模型MMP-Ⅱ在最大峰值热流数值接近或高于求解基于线性化的模型MMP-Ⅰ，在终端经纬度误差接近或小于线性化的模型MMP-Ⅰ；查看图12迭代过程，随着节点数量的增加，模型MMP-Ⅱ收敛到指定停止阈值条件，没有遇到像模型MMP-Ⅰ那样的代价函数数值跳跃现象。在计算效率方面，比较表4中的数据，并查看图12迭代过程：当节点数N为25、50时，求解基于线性化的模型MMP-Ⅰ和基于DC松弛的模型MMP-Ⅱ都能收敛到指定的停止阈值，但求解模型MMP-Ⅱ的解所需的迭代步数明显少于模型MMP-Ⅰ；而当节点数N为100、200时，求解模型MMP-Ⅱ满足停止阈值条件而停止，但求解模型MMP-Ⅰ到达最大迭代步数被停止。

图19倾侧角曲线（峰值热流指标）

Fig.19Bank angle profiles (peak heat flux index)

图20速度倾角曲线（峰值热流指标）

Fig.20Flight-path angle histories (peak heat flux index)

图21速度方位角曲线（峰值热流指标）

Fig.21Heading angle histories (peak heat flux index)

数值算例表明，基于DC松弛的模型MMP-Ⅱ要比基于线性化的模型MMP-Ⅰ对Max-Max类再入轨迹优化问题的逼近精度更高，在数值稳定性、收敛特性和代价函数的最优性方面表现更优，适合应用于高精度再入热轨迹规划任务，可为飞行器热试验提供参考。由于模型MMP-Ⅱ增加了更多的线性和非线性约束，求解模型MMP-Ⅱ的速度明显慢于求解模型MMP-Ⅰ。

5 结论

为了评估再入飞行器的力热性能，给材料、气动热和结构设计部门提供测试参考，本文采用序列DC规划方法研究了一类特殊的Max-Max类再入轨迹优化问题。主要贡献在于：

1）结合DC分解和罚函数法，推导获得了比线性化模型逼近精度更高的Max-Max类代价函数DC松弛模型，解决了线性化模型迭代过程中发生代价函数振荡或不收敛问题；

2）运用大-M法将Max-Max类轨迹优化问题转化为可解的混合整数规划子问题，提出了基于DC松弛模型的改进序列DC规划方法。

以峰值热流再入轨迹优化问题为研究对象，数值实验表明，基于DC松弛模型的峰值热流轨迹优化方法比基于线性化模型方法在最优性上要高0.578 1%，终端经纬度精度高2个量级，所提方法具有较高的稳定性、收敛性和最优性。

图1序列DC规划方法流程图

Fig.1Flowchart of SDCP method

图2高度-速度剖面（SCP）

Fig.2Altitude-vs-velocity histories by SCP

图3高度-速度剖面（SDCP）

Fig.3Altitude-vs-velocity histories by SDCP

图4经度-纬度-高度曲线（终端误差指标）

Fig.4Longitude-latitude-altitude footprints (terminal error index)

图5热流密度曲线（终端误差指标）

Fig.5Heat flux density histories (terminal error index)

图6动压曲线（终端误差指标）

Fig.6Dynamic pressure histories (terminal error index)

图7过载曲线（终端误差指标）

Fig.7Normal load histories (terminal error index)

图8攻角曲线（终端误差指标）

Fig.8Angle of attack profiles (terminal error index)

图9倾侧角曲线（终端误差指标）

Fig.9Bank angle profiles (terminal error index)

图10速度倾角曲线（终端误差指标）

Fig.10Flight-path angle histories (terminal error index)

图11速度方位角曲线（终端误差指标）

Fig.11Heading angle histories (terminal error index)

图12代价函数迭代过程

Fig.12Cost function iteration process

图13高度-速度剖面（N=200）

Fig.13Altitude-vs-velocity histories (N=200)

图14经度-纬度-高度曲线（峰值热流指标）

Fig.14Longitude-latitude-altitude footprints (peak heat flux index)

图15热流密度曲线（峰值热流指标）

Fig.15Heat flux density histories (peak heat flux index)

图16动压曲线（峰值热流指标）

Fig.16Dynamic pressure histories (peak heat flux index)

图17过载曲线（峰值热流指标）

Fig.17Normal load histories (peak heat flux index)

图18攻角曲线（峰值热流指标）

Fig.18Angle of attack profiles (peak heat flux index)

图19倾侧角曲线（峰值热流指标）

Fig.19Bank angle profiles (peak heat flux index)

图20速度倾角曲线（峰值热流指标）

Fig.20Flight-path angle histories (peak heat flux index)

图21速度方位角曲线（峰值热流指标）

Fig.21Heading angle histories (peak heat flux index)

表1初始和终端条件

Tab.1Initial and terminal conditions

表2决策变量的边界约束

Tab.2Box constraints of the decision variables

表3算例一的结果比较

Tab.3Comparison results in example1

表4算例二的结果比较

Tab.4Comparison results in example2

引用提醒

图(21) / 表(4)

引用本文

邓泽晓, 李骜, 刘鲁华, 等. 试验飞行器Max-Max类再入轨迹优化的序列差凸规划方法[J]. 国防科技大学学报, 2026, 48(1)： 183-195.

复制

DENG Z X, LI A, LIU L H, et al. Max-Max type entry trajectory optimization for testing vehicle by successive difference-of-convex programming[J]. Journal of National University of Defense Technology, 2026, 48(1)： 183-195.

Copy

计量

图1序列DC规划方法流程图

Fig.1Flowchart of SDCP method

图2高度-速度剖面（SCP）

Fig.2Altitude-vs-velocity histories by SCP

图3高度-速度剖面（SDCP）

Fig.3Altitude-vs-velocity histories by SDCP

图4经度-纬度-高度曲线（终端误差指标）

Fig.4Longitude-latitude-altitude footprints (terminal error index)

图5热流密度曲线（终端误差指标）

Fig.5Heat flux density histories (terminal error index)

图6动压曲线（终端误差指标）

Fig.6Dynamic pressure histories (terminal error index)

图7过载曲线（终端误差指标）

Fig.7Normal load histories (terminal error index)

图8攻角曲线（终端误差指标）

Fig.8Angle of attack profiles (terminal error index)

图9倾侧角曲线（终端误差指标）

Fig.9Bank angle profiles (terminal error index)

图10速度倾角曲线（终端误差指标）

Fig.10Flight-path angle histories (terminal error index)

图11速度方位角曲线（终端误差指标）

Fig.11Heading angle histories (terminal error index)

图12代价函数迭代过程

Fig.12Cost function iteration process

图13高度-速度剖面（N=200）

Fig.13Altitude-vs-velocity histories (N=200)

图14经度-纬度-高度曲线（峰值热流指标）

Fig.14Longitude-latitude-altitude footprints (peak heat flux index)

图15热流密度曲线（峰值热流指标）

Fig.15Heat flux density histories (peak heat flux index)

图16动压曲线（峰值热流指标）

Fig.16Dynamic pressure histories (peak heat flux index)

图17过载曲线（峰值热流指标）

Fig.17Normal load histories (peak heat flux index)

图18攻角曲线（峰值热流指标）

Fig.18Angle of attack profiles (peak heat flux index)

图19倾侧角曲线（峰值热流指标）

Fig.19Bank angle profiles (peak heat flux index)

图20速度倾角曲线（峰值热流指标）

Fig.20Flight-path angle histories (peak heat flux index)

图21速度方位角曲线（峰值热流指标）

Fig.21Heading angle histories (peak heat flux index)

表1初始和终端条件

Tab.1Initial and terminal conditions

表2决策变量的边界约束

Tab.2Box constraints of the decision variables

表3算例一的结果比较

Tab.3Comparison results in example1

表4算例二的结果比较

Tab.4Comparison results in example2

WANG Z B, GRANT M J. Autonomous entry guidance for hypersonic vehicles by convex optimization[J]. Journal of Spacecraft and Rockets,2018,55(4):993-1006.

WANG J B, CUI N G, WEI C Z. Optimal rocket landing guidance using convex optimization and model predictive control[J]. Journal of Guidance Control and Dynamics,2019,42(5):1078-1092.

LIU X F. Fuel-optimal rocket landing with aerodynamic controls[J]. Journal of Guidance Control and Dynamics,2019,42(1):65-77.

LIU X F, SHEN Z J. Rapid smooth entry trajectory planning for high lift/drag hypersonic glide vehicles[J]. Journal of Optimization Theory and Applications,2016,168(3):917-943.

LIU X F, SHEN Z J, LU P. Exact convex relaxation for optimal flight of aerodynamically controlled missiles[J]. IEEE Transactions on Aerospace and Electronic Systems,2016,52(4):1881-1892.

LIU X F, SHEN Z J, LU P. Entry trajectory optimization by second-order cone programming[J]. Journal of Guidance Control and Dynamics,2016,39(2):227-241.

LIU X F, SHEN Z J, LU P. Solving the maximum-crossrange problem via successive second-order cone programming with a line search[J]. Aerospace Science and Technology,2015,47:10-20.

WANG Z B. Maximum-normal-load entry trajectory optimization for hypersonic glide vehicles[C]//Proceedings of AIAA SciTech 2019 Forum,2019:1-16.

WANG Z B. Optimal trajectories and normal load analysis of hypersonic glide vehicles via convex optimization[J]. Aerospace Science and Technology,2019,87:357-368.

祝静娴, 王剑颖, 罗云皓, 等. 基于改进保守凸可分近似的高超声速飞行器轨迹规划方法[J]. 宇航学报,2024,45(6):935-947. ZHU J X, WANG J Y, LUO Y H,et al. A hypersonic vehicle trajectory planning method based on an improved conservative convex separable approximation[J]. Journal of Astronautics,2024,45(6):935-947.(in Chinese)

陈扬, 梁卓, 吕瑞, 等. 一种基于MISOCP的复杂绕飞区轨迹优化方法[J]. 宇航学报,2024,45(5):720-730.CHEN Y, LIANG Z, LYU R,et al. An MISOCP based trajectory optimization method for complex no-fly zone[J]. Journal of Astronautics,2024,45(5):720-730.(in Chinese)

LIU X F, LU P, PAN B F. Survey of convex optimization for aerospace applications[J]. Astrodynamics,2017,1(1):23-40.

MALYUTA D, YU Y, ELANGO P,et al. Advances in trajectory optimization for space vehicle control[J]. Annual Reviews in Control,2021,52:282-315.

WANG Z B. A survey on convex optimization for guidance and control of vehicular systems[J]. Annual Reviews in Control,2024,57:100957.

陈顺毅, 徐小平, 刘双喜, 等. 面向时间协同的高超声速滑翔飞行器集群再入轨迹规划[J/OL]. 国防科技大学学报.(2025-01-03)[2025-03-01].https://link.cnki.net/urlid/43.1067.t.20250103.0931.002. CHEN S Y, XU X P, LIU S X,et al. Time-coordinated re-entry trajectory planning for hypersonic glide vehicle cluster[J/OL]. Journal of National University of Defense Technology.(2025-01-03)[2025-03-01].https://link.cnki.net/urlid/43.1067.t.20250103.0931.002.(in Chinese)

MISRA G, BAI X L. Iteratively feasible optimal spacecraft guidance with non-convex path constraints using convex optimization[C]//Proceedings of AIAA SciTech 2020 Forum,2020:1-18.

AHMADI A A, HALL G. DC decomposition of nonconvex polynomials with algebraic techniques[J]. Mathematical Programming,2018,169(1):69-94.

LU P. Convex-concave decomposition of nonlinear equality constraints in optimal control[J]. Journal of Guidance Control and Dynamics,2021,44(1):4-14.

DOFF-SOTTA M, CANNON M, BACIC M. Robust trajectory optimisation for transitions of tiltwing VTOL aircraft[C]//Proceedings of 2023 IEEE Conference on Control Technology and Applications(CCTA),2023:425-431.

ZHAO Z C, SHANG H B. Trajectory planning for spacecraft formation reconfiguration using saturation function and difference-of-convex decomposition[J]. IEEE Transactions on Aerospace and Electronic Systems,2024,60(1):857-866.

DENG Z X, LIU L H, WANG Y J. Entry trajectory optimization of lifting-body vehicle by successive difference-of-convex programming[J]. Advances in Space Research,2024,74(11):5837-5859.

LUO Y H, WANG J Y, JIANG J,et al. Reentry trajectory planning for hypersonic vehicles via an improved sequential convex programming method[J]. Aerospace Science and Technology,2024,149:109130.

YUILLE A L, RANGARAJAN A. The concave-convex procedure[J]. Neural Computation,2003,15(4):915-936.

MOSEK ApS. MOSEK modeling cookbook[Z/OL].(2025-11-05)[2025-11-10].https://docs.mosek.com/MOSEKModelingCookbook.pdf.

BALAS E. Disjunctive programming[M]. Cham: Springer,2018.

PHILLIPS T. A common aero vehicle(CAV)model,description,and employment guide[EB/OL].(2003-01-27)[2025-03-10].https://www.doc88.com/p-501932383682.html.

Home

About Journal

Guide for Authors

Editorial Board

Publication Statement

Open Access Statement

Contact

Journal Subscription

Rss

AI assistant

Chinese

1 再入轨迹优化模型

2 模型离散与DC分解

3 序列差凸规划方法

4 数值算例

5 结论