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

2020-12-01周祥张洪波何睿智汤国建包为民

航空学报 2020年11期

周祥，张洪波,*，何睿智，汤国建，包为民

1. 国防科技大学空天科学学院，长沙 410073 2. 中国航天科技集团有限公司科技委，北京 100048

可重复使用飞行器能够实现天地间的多次往返，从空间再入大气层后，依靠大升阻比气动外形实现在临近空间的远距离滑翔飞行，从而抵达预定目标[1]。再入滑翔段轨迹规划的准确性和快速性对于改善该类飞行器的制导效果具有重要意义。

为了降低轨迹规划的计算量，以航天飞机轨道器为代表的可重复使用飞行器采用固定攻角剖面，以倾侧角为唯一控制量的方式进行轨迹规划，进而实现制导任务。这种基于二维剖面的规划方法具有计算负担小、方法简单、容易实现等优点，在航天飞机等飞行器上得到了成功应用。然而，这种计算量小的规划方法难以充分发挥可重复使用飞行器的大范围机动能力。飞行器固有的机动能力主要由飞行器的总体布局和气动性能决定，上述因素确定后，如何选择合适的控制策略从而充分发挥出其固有的机动能力就显得尤为重要。1998年，文献[2]提出如果将传统二维剖面规划方法中固定攻角的条件放开，同时采用攻角和倾侧角作为控制量，将有助于充分发挥飞行器的机动能力。文献[3]将同时考虑攻角和倾侧角作为控制量的轨迹规划问题描述为三维剖面规划问题，使用阻力加速度和侧向加速度作为中间控制变量(即剖面设计量)，引入降阶动力学系统降低问题求解的计算量，这里的三维剖面从控制量的角度讲即是对攻角和倾侧角变化规律的一种表征，而传统的二维剖面规划方法只能获得对倾侧角变化规律的表征。文献[4]在加速度剖面规划中同时考虑纵向和侧向运动，首先基于二维剖面规划一条可行轨迹，其次基于三维剖面获得最优轨迹。文献[5]提出了一种分层策略的三维剖面规划方法，该方法采用侧向降阶运动模型对纵向、侧向剖面进行反复迭代，直至满足终端纵程和横程要求。文献[6]将基于分层策略的三维剖面规划方法应用到再入覆盖区求解中，获得了相比基于二维剖面规划方法更大的覆盖区，验证了机动能力得到充分发挥后的效果。整体来看，目前对二维剖面规划方法的研究较多[7-9]，而关于三维剖面规划方法的研究还比较少，还有很多技术问题有待解决，其中，计算效率是限制三维剖面规划技术的主要瓶颈之一。

近年来，内点法的快速发展使得凸优化在求解效率方面相比传统方法如伪谱法的优势愈发明显[10-11]。除此之外，凸优化问题的局部极小点就是全局最小点，这一良好的理论性质也为该方法的广泛应用奠定了基础[12]。由于实际中大量的优化问题都是非凸的，如何将轨迹规划等这类具有多种约束的原始非凸问题转化为凸问题是应用凸优化方法的关键，包括对非线性动力学、过程约束、控制量约束、优化目标的凸化。截止目前，凸优化已经在再入轨迹规划、行星着陆、无人机路径规划等方面得到了应用，一系列凸化方法被广泛使用在这些问题中，包括约束松弛和连续线性化等技术[13]。凸优化早期在航天中的应用主要是解决动力下降段的轨迹规划问题，在这类问题的凸化过程中提出了无损凸化的概念[14-15]。

文献[16]利用序列凸优化的方法求解再入轨迹规划问题，相比伪谱法显示出了很好的求解效率，文中所提的序列凸化方法是目前处理强非线性动力学约束与非线性过程约束的主要方法。文献[17]中讨论了利用序列凸优化求解再入轨迹规划问题的收敛性问题，序列凸优化的收敛性与解的等价性紧密相关。文献[18]实现了基于凸优化在线轨迹规划的闭环跟踪制导，将轨迹规划子问题和轨迹跟踪制导指令求解子问题均转化为凸问题，从而可实现在线快速求解，取得了很好的制导效果。文献[19-21]对于改进序列凸优化方法的性能提供了多种技巧。文献[22]利用凸优化求解再入轨迹规划问题时将攻角与倾侧角同时作为控制量，通过定义中间变量实现了对控制量约束的凸化处理。但并未给出规划后的指令计算方法，限制了该方法在实际问题中的应用。

总体来看，目前基于三维剖面的再入轨迹规划方法研究还比较少，主要存在的问题包括：① 动力学模型简化较多，导致问题求解精度较低；② 为了减低计算量，将加速度剖面进行简单参数化，人为限制了剖面求解的搜索空间，不利于充分发挥飞行器的能力。总的来看，已有三维剖面规划方法并未很好解决轨迹规划的求解精度与计算效率之间的矛盾。

基于上述分析，本文采用完整动力学模型，利用凸优化方法求解再入三维剖面规划问题。主要贡献是：① 将原始三维剖面规划问题描述为最优控制问题，采用合适的凸化技巧将考虑多种约束条件的原问题转化为凸优化问题，包括动力学约束、路径约束和控制量约束的凸化，其次对连续凸问题进行离散化；② 引入指令反解步骤，通过优化后的中间控制量计算攻角和倾侧角指令，并建立了相应的序列迭代算法，通过迭代求解凸优化子问题，从而获得原问题的可行解。通过数值仿真证明了本文方法的可行性和收敛性。

本文的框架结构：引言部分介绍研究背景和研究现状；第1部分介绍问题描述；第2部分介绍如何对原问题进行凸化与离散处理，将其转化为一个凸的数学规划问题；第3部分介绍如何利用凸优化的解(即中间变量和状态量)反解攻角和倾侧角指令，并给出了序列凸化算法的步骤；第4部分进行了数值仿真；第5部分给出了结论。

1 问题描述

1.1 三维剖面规划问题

首先建立轨迹规划中使用的动力学模型。在建模过程中主要引入以下假设： ① 考虑地球自转；② 仅考虑球形引力；③ 认为地球外形是一个标准圆球。可重复使用飞行器在再入过程中进行无动力飞行，由于大气阻力的作用，飞行器能量是单调递减的，为了降低优化变量的数目，本文采用能量E作为自变量，定义为

E=V2/2-μ/r

(1)

式中：μ为地球引力常数；r为地心距；V为速度。

进而得到以能量为自变量的再入质心动力学模型为

(2)

式中：λ为经度；φ为纬度；θ为当地速度倾角；σ为航迹偏航角；L、D分别为升力和阻力加速度；υ为倾侧角；由地球自转引起的哥氏加速度项Cθ和Cσ以及牵连加速度项C′θ和C′σ可分别表示为

(3)

上述所有变量均可采用如下归一化方法进行无量纲化处理。

(4)

(5)

式中:R0为地球平均半径;g0为R0处的引力;ωe为地球旋转角速度。升力和阻力加速度的定义为

(6)

式中：Mr为飞行器质量；Sr为参考面积；CL、CD分别为升力系数和阻力系数，通常为攻角α和马赫数Ma的函数；ρ为大气密度，有

ρ=ρ0e-β h

(7)

其中：β=1/hs，hs=7 110 m；ρ0=1.225 kg/m3，h为高度。如无特殊说明，下文所有变量均为无量纲化处理后的值。

以E为自变量的动力学方程包括5个状态参数和2个控制指令。前者指地心距r、经度λ、地心纬度φ、当地速度倾角θ和航迹偏航角σ，而速度与自变量有关，可以通过式(8)确定

(8)

将独立描述飞行器运动状态的参数构成状态矢量x，可表示为

x=[r,λ,φ,θ,σ]T

(9)

2个控制指令是攻角和倾侧角，利用它们实现对气动力大小和方向的控制，所以再入轨迹规划中的控制矢量U′可表示为

U′=[α,υ]T

(10)

在轨迹规划中需要对控制指令的幅值进行约束，可以表示为

αmin≤α≤αmax

(11)

υmin≤|υ|≤υmax

(12)

式中：|·|表示绝对值，υmax>υmin≥0。

过程约束主要有：热流密度、动压、法向过载和气动力，前3个约束用于避免飞行器结构受到破坏，第4个约束用于确保飞行器有足够气动力用于轨迹控制，分别表示为

(13)

(14)

(15)

(16)

轨迹规划的边界条件包括初始状态约束和终端状态约束，表示为

(17)

式中：x0和xf由具体任务给定;E0、Ef分别为初始状态和终端状态对应的能量。

优化目标可以表示为一般形式

(18)

式中：Φ(·)为终端目标函数；L(·)为过程目标函数。

至此，将再入轨迹三维剖面规划问题记为P0，表示如下

(19)

在问题P0中，动力学方程和过程约束均是非线性的，因此这是一个连续的非线性最优控制问题。同时，由于控制指令攻角并不显含在方程式(2)中，这也增加了该问题的求解难度。

1.2 新控制量的选择

再入轨迹三维剖面规划问题的控制量是攻角和倾侧角，但是，由于攻角在再入动力学方程式(2) 中是隐含的，攻角通过影响气动系数从而改变气动力，倾侧角则通过与三角函数进行复合而显含于式(2)中，因此，原动力学方程式(2)关于攻角和倾侧角是强非线性的。由于凸优化问题中等式约束必须是线性的，若继续使用攻角和倾侧角作为轨迹规划的控制量，则显著增加了对原方程式(2) 进行线性化的难度，并且线性化后的动力学方程也容易引起再入轨迹的“抖振”现象[16]。

为了避免上述问题的出现，可以定义新控制量，从而获得新的线性化方程。新控制量的选择应满足以下2个条件：① 新控制量可以比较方便地将动力学方程中与之相关的一项表示为关于它的线性表达式；② 新控制量与原始控制量之间应该存在唯一映射关系。因此，本文定义如下新的控制变量:

(20)

(21)

这里引入u3的意义在于建立u1和u2之间的关系，方便后续进行控制量约束的等价转换。

定义了新的控制量后，式(2)所示的动力学方程可以重新写为

(22)

令U=[u1,u2,u3]Τ，将式(22)简写为

(23)

式中：

(24)

(25)

(26)

在动力学方程式(23)中，新控制量矢量U关于矩阵G(x)是线性的；g(x)为关于状态量的非线性项，但其中包含阻力加速度项D,这与原始控制量攻角有关，关于该问题对求解精度的影响将在数值仿真中进行分析；h(x)为与地球自转相关的小量。

2 凸化与离散

2.1 凸化

1)控制量约束的处理

定义了新控制量(u1,u2,u3)后，需将如式(11) 和式(12)所示的原始控制量约束转化为对新控制量的约束。由于控制量u3仅与攻角和马赫数有关，由参考轨迹对应的高度和速度可得到参考马赫数Mar。基于攻角取值范围Cα，采用一维优化方法线搜索确定控制量u3的取值范围，表示为

(27)

(28)

式中：Cα=[>αmin,αmax]，则攻角幅值约束可以转化为控制量u3的幅值约束，表示为

(u3)min≤u3≤(u3)max

(29)

倾侧角幅值约束可以转化为对控制量u1的约束，表示为

u3cos(υmax)≤u1≤u3cos(υmin)

(30)

至此，实现了对控制量幅值约束的转化处理，且转化后的约束(29)和(30)均是凸的。

除了上述对控制量幅值约束进行转化外，对式(21)所示的非凸附加约束也要进行凸化处理。这里采用松弛技术对其进行凸化[16]，得到如下凸约束为

(31)

上述松弛过程的合理性证明将在后面进行介绍。

2) 动力学方程的线性化

动力学方程作为等式约束，在凸优化问题中必须是线性的。这里介绍如何对式(23)所示的动力学方程进行线性化处理。基于参考轨迹(x(k),α(k),υ(k))，其中上标k表示参考轨迹序号，由一阶泰勒展开可得

(32)

将式(32)化简为

(33)

其中：

(34)

B=G|x=x(k)

(35)

C=g(x(k))-Ax(k)+h(x(k))

(36)

根据推导，线性化矩阵A和B的表达式为

(37)

B(xk)=

(38)

其中：

(39)

(40)

(41)

(42)

为了降低线性化误差，需添加信赖域约束，表示为

|x-x(k)|≤δ

(43)

式中：δ表示信赖域半径，与x维数相同。

3) 过程约束的线性化

在式(13)～式(16)所示的4项过程约束中，对应函数均具有强非线性，且既不是凸函数也不是凹函数，因此仍需要对这些函数进行线性化。本文采用文献[16]提出的处理方法，利用一阶泰勒展开，将4项过程约束转化为关于地心距的线性约束，统一记为

(44)

4) 目标函数的选择

在凸优化问题中，目标函数必须是凸的或线性的。为了简化处理，这里以终端位置误差最小为主要指标。由于终端位置由经度和纬度共同描述，这也是一个终端性能指标，可以表示为

(45)

(46)

式中：C1和C2均为权重系数，用于平衡不同项之间的量级差异。

由于引入的附加项是线性的，因此式(46)所示的目标函数仍然是凸的。

基于上述凸化处理，可以得到具有凸目标函数、凸或线性约束的连续优化问题，记为问题P1，表示为

(47)

由于问题P1中优化目标是终端位置误差最小，因此相应的终端约束应去除对经纬度的限制，将终端约束调整为

xf=[rf, ～, ～,θf,σf]Τ

(48)

式中：～表示该项无约束。从工程实现的角度考虑，在终端约束中对速度倾角和航迹偏航角施加了严格等式约束。

为了证明问题P1中松弛过程的合理性，这里引入如下2个假设：

假设1式(43)所示的信赖域约束|x-x(k)|≤δ在定义域的大部分区间内不会取到边界，即除了少数有限个连接点外，在[>E0,Ef]上|x-x(k)|<δ是恒成立的。

假设2式(44)所示的过程约束rmin≤r≤rmax在定义域的大部分区间内不会取到边界，即除了少数有限个连接点外，在[>E0,Ef]上rmin

对于假设1，只要在求解过程中设置一个充分大的δ，即可保证在[>E0,Ef]上|x-x(k)|<δ是恒成立的。对于假设2，由于求解过程采用了梯度下降算法，当迭代解接近过程约束边界时，控制量将迫使当前解远离边界，所以能够保证在大部分连续定义域内rmin

2.2 离散

问题P1是一个连续凸优化问题，其目标函数是凸的，约束是凸或线性的，但状态量和控制量是连续的，这为直接求解该问题带来了困难。因此，将其转化为对有限参数进行寻优的数学规划问题，应用数值方法进行求解是一种可行的方法。这里介绍对问题P1进行离散化，从而便于求解。

将自变量变化域[>E0,Ef]等间距离散为N个间隔，则自变量离散为E0,E1,E2,…,EN-1,EN，且EN=Ef；状态量离散为x0,x1,x2,…,xN-1,xN，且xN=xf；控制量离散为U0,U1,U2,…,UN-1,UN。

采用梯形法则对问题P1中的线性动力学方程进行近似离散，则离散后的动力学方程为

Bj+1Uj+1+Cj+1)

(49)

式中：ΔE=(Ef-E0)/N，下标j表示离散点编号，j=0,1,…,N-1。

整理式(49)，可得到更简洁的动力学方程形式，表示为

LjXj-Lj+1Xj+1+MjUj+Mj+1Uj+1=Gj

(50)

完成动力学方程离散化后，再对问题P1中的其他约束和目标函数积分项进行离散化，进而得到一个凸数学规划问题，记为问题P2，表示为

(51)

问题P2中的优化变量包括每个离散点处的状态量xj和控制量Uj。由于N个离散间隔对应N+1个离散点，因此优化变量数目为8(N+1)。

3 指令反解与迭代算法

3.1 指令反解

由于问题P2中采用状态量x和新控制量U作为优化变量，经过凸优化求解后不能直接获得用于制导的攻角和倾侧角指令，因此需要利用凸优化解进行指令反解，获得每个离散点处相应的攻角和倾侧角指令。

假设优化后第j个离散点处的状态量为xj，新控制量为Uj=[u1,j,u2,j,u3,j]Τ，则对应攻角指令αj可利用牛顿法或反插值方法求解下列非线性方程确定。

F(αj,Maj)=u3,j

(52)

式中：Maj为由xj所确定的当前马赫数。

倾侧角指令可利用式(53)确定

(53)

在问题P2的求解过程中，控制量附加约束很难保证在每个离散点处都满足等式成立。这里定义附加约束违反程度函数，用于定量描述该约束的违反情况，表示为

(54)

当违反程度VU超出给定容许值εU时，则由式(52)和式(53)所得的反解指令(αj,υj)是不可信的，反之，则是可信的。本文将在仿真中分析这种情况对解精度的影响。

违反程度容许值的取值可以根据如下经验公式确定:

(55)

式中：Z为经验常数，一般可以取100～500之间。如果Z取值过大，则容许值εU较小，使得较多反解指令是不可信的；反之，若容许值较大，则大部分反解指令是可信的，则失去了判断反解指令可信度的意义。常数Ζ的选择应根据具体问题进行权衡确定。

3.2 序列凸化算法

由于原最优控制问题P0具有强非线性，如果仅采用一次线性化来处理非线性动力学与过程约束，即仅求解一次问题P2，则所得凸优化解对于原问题P0基本是不可行的。因此需采用连续线性化方法，对动力学与过程约束进行连续线性化，从而产生一个迭代求解序列。

在这个序列中，每次迭代都会求解一个凸优化子问题(即问题P2)，利用该子问题的解再进行下一次的线性化，直至这个解序列最终满足收敛条件。为了保证该求解序列的收敛性[15]，在信赖域约束中引入松弛系数κ，并将其作为一个附加项添加至目标函数中，从而得到新问题P3，表示如下

(56)

式中：C3为松弛系数项的权重系数。

动力学约束是轨迹规划问题中必须满足的一个约束，问题P3所得的解仅满足线性动力学约束，相比问题P0中的非线性约束必然存在一定误差，定量描述这个误差可用于表示凸优化解相对原问题P0的可行性。因此，下面定义非线性约束违反程度，用于定量描述凸优化解对原非线性动力学约束的违反情况。

求解问题P3得到的凸优化解可表示为(x(k),U(k))，通过指令反解得到的控制指令为(α(k),υ(k))。以该控制指令为输入，对式(2)所示的原始动力学方程进行积分，得到相对于凸优化解的积分验证解，表示为

(57)

由此可定义非线性约束违反程度，也可称为凸优化解的近似误差，用于表征优化后的状态量x(k)与控制指令(α(k),υ(k))对原始非线性动力学约束的违反程度，计算公式为

(58)

该迭代算法以相邻两次凸优化解对应的状态量最大偏差作为收敛准则，表示为

(59)

式中：收敛误差限定义为ε=[εr,ελ,εφ,εθ,εσ]Τ，上述不等式左边的表达式可定义为每个状态量的相邻2次凸优化解最大偏差。

基于上述分析，这里给出基于凸优化的序列凸化算法求解原最优控制问题P0的步骤：

步骤1给出一条初始参考轨迹x(0)，初始参考指令(α(0),υ(0))，令k=0。

步骤2基于参考轨迹x(k)和参考指令(α(k),υ(k))计算控制量约束边界、过程约束边界以及线性化矩阵。

步骤3求解凸优化子问题P3，获得凸优化解(x(k+1),U(k+1))。

步骤4利用凸优化解进行指令反解，得到新的指令序列(α(k+1),υ(k+1))，以此作为积分输入，计算非线性动力学违反程度函数Rnonlinear。

步骤5判断是否满足如式(59)所示的收敛条件，若满足，则算法结束；否则，令k=k+1，转入步骤2。

4 数值仿真

本文采用CAV-H模型进行仿真分析[23]，初始化轨迹生成方式可参考文献[16]，初始攻角指令取为20°常值，初始倾侧角指令取为0°常值。边界条件设置如表1所示，主要约束设置如表2所示。

表1 边界条件取值Table 1 Values of boundary constraints

表2 其他约束条件取值Table 2 Values of other constraints

仿真中采用MATLAB平台下CVX作为建模语言[24]，计算精度设为1×10-6，即当相邻2次迭代的目标函数相对变化量小于预设计算精度时，求解器迭代终止。默认SDPT3求解器求解凸优化子问题P3，运行计算机搭载Intel Core i5-3470 3.20 GHz处理器。

由于整个滑翔飞行过程中能量E一直为负值，为了显示方便，将其转换至[0,1]范围内，公式为

(60)

此时由再入点开始，能量E′从0到1单调递增。

4.1 可行性分析

本节主要分析所提三维剖面规划方法优化解的可行性，并验证该方法相比二维剖面规划方法在发挥飞行器固有机动能力方面的优势。这里的可行主要包括2个方面，一是优化解的状态轨迹与控制指令(包括攻角和倾侧角)之间是否满足非线性动力学约束，二是由控制指令积分获得的过程量是否满足非线性约束要求。

图1～图5给出了本文仿真条件下利用所提方法得到的再入轨迹。其中，蓝色粗实线代表的“优化解”为迭代过程收敛时最后一次求解问题P3获得的解，以红色点划线代表的“积分解”由式(57) 定义可得。这里以积分解的结果作为精确值，验证“优化解”的可行性。

图1中对轨迹规划的起始点进行了标注，以该点对应的状态作为轨迹规划的起始状态，在此之前的轨迹状态按照常值控制量进行积分获得。由图1～图3可以看出，优化解对应的状态量与积分解的结果基本一致，说明优化解中的状态量与控制指令之间基本满足原始非线性动力学约束，没有出现较大偏差。由图3可知，优化解对应的轨迹满足期望目标点位置约束。

图4给出了热流密度、动压和过载3项过程约束的变化规律。可以看出，优化解与积分解对应的轨迹均满足给定的过程约束，说明相应的状态量基本一致。同时，由凸优化模型得到的状态量仍然满足原始过程约束，即相对于原问题是可行的，说明应用线性化方法处理非线性过程约束在该问题中是可行的。

图1 高度和速度变化规律Fig.1 History of altitude and velocity

图2 速度倾角和航迹偏航角变化规律Fig.2 History of flight path angle and heading angle

图3 地面轨迹对比示意图Fig.3 Comparison of footprints

图5给出了利用优化解进行指令反解得到的攻角和倾侧角曲线，从控制量的角度讲它们即是对三维剖面的一种表征。在规划起点之前，攻角设置为20°常值，倾侧角设置为0°常值。图1～图4 中积分解均是以图5所示的控制指令为输入进行积分获得的。

下面定量分析优化解的可行性，主要指优化解相对于原最优控制问题的非线性约束违反程度。图6给出了优化解对应的控制量附加约束违反程度V在所有离散点上的分布情况。可以看出，在第153～157个离散点之间，违反程度异常增大，说明这几个离散点处的反解指令是不可信的，而其他离散点处的反解指令是可信的。

图7给出了优化解在每个离散点处的近似误差分布情况。可以明显发现，在反解指令异常情况出现前后，单个离散点的近似误差出现明显变化，由此也导致后续离散点处的近似误差都显著增大。这个现象与图2中，在归一化能量取值0.8以后，优化解与积分解开始出现微小差异是一致的。这些误差的产生都是由于优化解在极个别离散点处不严格满足如式(21)所示的控制量附加约束引起的，而这种个别离散点处的现象一般是由于数值求解不稳定造成的。

图4 三项过程约束变化规律Fig.4 History of three path constraints

图5 控制指令曲线Fig.5 History of control command

表3给出了优化解与积分解的终端状态对比结果。可以看出终端状态偏差是比较小的，即解精度在一定范围内是可以保证的。这也说明根据本文定义的新控制量得到的线性动力学相对于原始非线性动力学的约束误差是有限的。所以，在一定精度条件下，由本文方法所得的凸优化解可以作为原问题的一个可行解。

图6 控制量附加约束违反程度分布Fig.6 Violation distribution of additional constraint

图7 逐点近似误差分布Fig.7 Distribution of approximation error on each discretization point

表3 终端状态对比Table 3 Comparison of terminal states

为了比较三维剖面规划方法相较于二维剖面规划方法的优势，以文献[16]中的二维剖面规划方法为参考，仿真条件与本文相同，2种方法所得的地面轨迹对比结果如图8所示，终端位置偏差对比结果如表4所示。可以看出，在相同条件下，三维剖面规划方法的优化解对应的终端位置可以到达给定目标点，距离偏差为0，但二维剖面规划方法的优化解对应的终端位置与给定目标点相距1 678.57 km。由于仿真设置的优化目标是终端位置偏差最小，此时二维剖面规划方法的优化结果代表了其所能发挥的最大机动能力，对比而言，三维剖面规划方法所能发挥的机动能力更大。由此可以得出，三维剖面规划方法在再入滑翔轨迹规划时可以充分发挥飞行器固有的机动能力。

4.2 收敛性分析

序列凸化算法的收敛性对于该方法的可行性至关重要，本节对4.1节算例求解过程的收敛性进行分析。

根据式(59)所示的收敛准则，当序列凸化算法迭代至第15次时，各个状态量的相邻两次凸优化解最大偏差分别为

|Δrj|max=78.286 6 m，|Δλj|max=0.044 5°，|Δφj|max=0.020 0°，|Δθj|max=0.045 4°，|Δσj|max=0.198 4°，满足仿真中所给的收敛误差限，迭代终止。

为了体现凸优化的求解效率，在相同计算条件和计算精度下，利用GPOPS 5.0伪谱法工具包求解该问题，则本文所提方法与伪谱法的比较结果如表5所示，其中，凸优化方法的求解时间由每次子问题求解输出的CPU时间累加获得，伪谱法的求解时间由每次迭代输出的实际求解时间累加获得。

由表5可以看出，凸优化和伪谱法均可获得满足各项约束条件的再入轨迹，实际计算精度均小于预设计算精度，2种方法获得的性能指标差异小于10%。相比伪谱法，凸优化方法迭代次数更少，在计算效率方面具有明显优势。

图9给出了序列凸化迭代过程中得到的历次轨迹与初始化轨迹的对比结果，其中，蓝色实线表示由线性假设得到的初始化轨迹，红色点划线表示迭代收敛后输出的优化轨迹，其他细实线表示中间迭代解对应的轨迹。可以看出，随着迭代次数的增加，相邻2次轨迹的偏差逐渐减小，优化轨迹逐渐收敛。同时，收敛后的优化轨迹与初始化轨迹相差较大，说明本文所提方法的收敛性对于初始化轨迹的选择具有一定的适应能力。

图10给出了5个状态量的收敛情况，所有子图的横轴取值从1开始。图10(a)中高度收敛图的纵坐标是对数坐标。本文中状态量的收敛过程主要由相邻2次凸优化解的最大偏差来衡量，由式(59)计算可得。总的来看，5个状态量的收敛趋势是比较一致的，最大高度偏差在最后2次迭代时趋于收敛，而最大速度倾角偏差在第10次迭代后趋于收敛，剩余3个最大状态偏差则在第6次迭代后趋于收敛。这说明不同状态量的收敛速率是不一致的。

图11给出了信赖域松弛系数的收敛情况。可以看出，除了前2次迭代时得到的松弛系数比较大以外，后续迭代过程中松弛系数始终比较小，且最后达到收敛条件时松弛系数取值为κ*=0.002 7，说明随着信赖域松弛系数的减小，序列迭代算法逐渐趋近于收敛。

图10 5个状态量收敛情况Fig.10 Convergence process of five state variables

图12给出了最优性能指标的收敛情况，该图中纵坐标为对数坐标。可以看出，优化目标的收敛是比较快的，基本在第6次迭代以后就趋于稳定。当序列凸化算法收敛时，目标函数中的终端位置偏差和信赖域松弛变量都将是一个小量，此时影响目标函数取值的主要是航迹偏航角的积分项，而该项的变化仅与侧向运动相关。所以，经度、纬度和航迹偏航角这3个与侧向运动相关的状态量，其收敛过程与优化目标基本一致。由于目标函数中缺少纵向运动参数，所以与纵向运动相关的高度和速度倾角这2个状态量在迭代过程中出现了较多震荡，一定程度上降低了算法的整体收敛速率。

以上分析从数值上验证了本文所提方法的收敛性，且收敛后的优化解仅可能是原问题的一个局部最优解，关于这种方法最优性的理论证明目前还并不充分[17]。

图12 最优性能指标收敛情况Fig.12 Convergence process of optimal objective function

由于控制量附加约束的违反程度对于优化解精度有一定影响，因此在目标函数中引入关于航迹偏航角积分的附加项，如式(46)所示，下面分析该附加项的引入对解精度的影响。图13给出了迭代过程中，每次迭代的解中，严格满足附加约束的离散点占总离散点数目的比例变化情况。可以看出，引入附加项后，除了第8次和第11次的比例很低外，其他每次迭代解对应的比例都比较高，迭代后期这个比例值基本趋于稳定，而在不引入附加项的迭代过程中，满足附加约束的离散点占比在后期一直较低，由此说明该附加项的引入可以有效确保附加控制量约束松弛过程的有效性。

图14给出了迭代过程中近似误差的收敛情况。可以看出，引入附加项后，由式(58)所定义的近似误差在迭代过程中的收敛性是比较好的，在迭代前期很快能够减小到一个很低的水平，说明随着迭代算法逐步趋于收敛，凸优化解对原始非线性约束的违反程度逐渐降低，即解的精度逐步提高。但在不引入附加项的迭代过程中，由于大部分离散点处的优化值均不满足附加约束，所以导致解精度较低。由此说明附加项的引入对于提高优化解精度是有显著意义的。

图13 附加项对离散点占比的影响Fig.13 Effect of additional term on proportions of discretization points

图14 附加项对凸优化解的近似误差影响Fig.14 Effect of additional term on approximation errors of solution by convex optimization

序列凸化算法的启动必须依赖于一条初始轨迹。本文采用线性假设进行初始轨迹的生成，即认为运动状态的滑翔段均是线性变化；作为对比的是以表5中提到的伪谱法所得优化轨迹作为初始轨迹。下面分析不同初始轨迹对本文所提算法收敛性的影响。图15～图17分别给出了信赖域松弛系数、最优性能指标和凸优化解近似误差的收敛对比情况。

由图15～图17可以看出，采用伪谱法得到>的优化轨迹作为初始轨迹后，序列凸化算法在迭代12次即可达到收敛条件，而同样条件下线性假设得到的初始轨迹需要15次迭代。由于收敛条件的参数不变，所以2组优化解的计算精度和近似误差也基本一致。在相同计算精度条件下，由于迭代次数的减少，以伪谱法优化轨迹作为初始轨迹的迭代过程耗时也更短，由于每次迭代的CPU时间基本一致，所以节约时间比例与迭代次数的减少比例也是基本一致的。