高超声速滑翔飞行器利用气动力在大气层内作远距离滑翔，具有速度快、航程远、机动性强等许多优势。但其飞行过程收到复杂的非线性气动力和动压、热流、过载、禁飞区等非线性过程约束，使其轨迹规划成为一个富有挑战性的问题，而当试图在线求解时就更加困难。
传统的直接法或者间接法通常难以收敛，或者需要巨大的计算量^[1-4]。航天飞机所采用的基于阻力加速度剖面的轨迹规划^[5]和一些衍生方法^[6-7], 曾取得巨大的成功，但是对于滑翔飞行器的机动能力往往有较大的限制。还有一些方法将轨迹规定为特定的分段解析形式^[8-9]，并在线调整参数，这将能明显提高计算效率，但可能意味着许多的可行解因为偏离这类特定形式而被排除了。总的来说，传统方法的主要不足在于：
1) 完备性不好。有许多可行解无法用这些方法找到。这相当于浪费了飞行器本来的飞行能力。
2) 计算复杂。传统数值方法收敛范围小，实际计算中经常难以收敛，或收敛到不可行解。
高超声速再入轨迹规划，由于其强非线性难以直接求解，通常要采用一系列子问题逼近原问题。这些子问题本身必须是易解的，并且尽可能接近原问题。而数学上所有易解的规划问题几乎都指向了凸优化^[10]，它是多项式复杂的，即使规模很大也能可靠求解^[11]。
近年来有****尝试将凸优化引入再入轨迹规划。Liu等^[12]基于以能量为自变量的方程在速度-攻角剖面给定的条件下，使用序列凸优化求解了再入轨迹优化问题。Wang和Grant^[13]则以时间为自变量并以倾侧角的导数为控制量，求解了终端时间固定的轨迹规划。现有的基于凸优化的滑翔再入轨迹规划方法基本上都直接使用了小扰动线性化来处理微分方程的非线性同时相应地添加了信赖域约束来保证小扰动，并需要一段特定的模拟飞行轨迹作为初值。这些处理方法简单易行，为再入轨迹规划的求解提供了新的思想，但是实际计算过程中还可能遇到以下问题：
1) 采用小扰动线性化处理动力学方程和约束，本质上类似于牛顿法，收敛范围小。
2) 子问题可行性将是一个严峻的问题。为了保证小扰动采用了信赖域约束，同时每个子问题都要求满足终端和其他路径约束。但是当初值不够精确时，在初值的信赖域内就不存在满足其他约束的解。这将使子问题不可行，迭代过程将直接失败。
为了解决上述问题，本文首先对动力学方程进行处理，以轨迹弧长为自变量，选择合适的控制变量以后，动力学方程的非线性大为减弱。以对数速度代替速度作为状态量，使最主要的路径约束——动压和热流约束完全成为线性约束。对于禁飞区约束，采用了使用类似于混合整数规划中割平面^[14-15]的思想处理，这比以往的处理方法更加适合凸优化，并尽可能避免了不必要的计算。数值仿真显示了本文算法的可靠性和快速性。
1 问题描述 1.1 再入动力学方程 考虑忽略地球自转的再入动力学方程：

(1)

(2)

(3)

(4)

(5)

(6)

式中：h为高度；r为地心距，按地球平均半径R₀归一化；θ和?分别为经纬度，可以通过选择适当的坐标系使?在飞行过程中保持较小以避免式(2)出现奇异并减弱其非线性；V为对地速度，按归一化，V_u、g₀和R₀分别为第一宇宙速度、海平面平均重力加速度和地球平均半径；γ为当地速度倾角；ψ为航向角；L和D分别为升力和阻力, 按g₀归一化；σ为侧倾角。

(7)

(8)

式中：ρ为大气密度，采用指数形式的大气模型ρ=ρ₀exp(-βh)，ρ₀为海平面上的标准大气密度，β=1.378 5×10^-4为大气模型中的常数；A_ref为飞行器的气动参考面积；C_L和C_D分别为升力系数和阻力系数；m为滑翔体的质量。
定义归一化速度在时间上的积分为归一化弧长：

(9)

式中：dτ为时间的微分。
注意到弧微分ds/dt=V, 以z=ln V代替V，并且滑翔再入时γ通常很小，高度和地心距相比是小量，可以作近似sin γ≈γ, cos γ≈1, r≈1。得到关于s的动力学方程：

(10)

(11)

(12)

(13)

(14)

(15)

式中：k_f为常数，定义为

(16)

近似带来的误差，将在多次在线规划中不断的得到修正。
1.2 基本约束 再入过程要满足热流、动压和过载等路径约束，可分别表示为

(17)

(18)

(19)

式中：k_Q为取决于头部形状的系数；、q_max和n_max分别为热流、动压和过载的最大值。
对式(17)、式(18)两边取对数，就得到完全线性的形式为

(20)

(21)

定义如下控制量：

(22)

(23)

(24)

记控制向量u=[u₁, u₂, u₃]^T，则过载约束重新表述为

(25)

由于内在的气动关系，升力系数C_L和阻力系数C_D是耦合的，所以控制量u₁、u₂、u₃是耦合的，但是可以暂时将它们看作独立的。虽然这样得到的解是松弛解，不能满足准确的气动关系，但之后可以通过在迭代过程中逐步改进使其逼近真实的气动关系。
具体说来，由于升阻比λ变化通常不大。在线求解可以采用以下约束代替具体的气动模型，其带来的误差可以通过飞行过程中多次重规划来修正，还可以在优化指标中加入对应的项来弥补。

(26)

(27)

(28)

2 凸化处理 2.1 动力学方程进一步处理 动力学方程式(11)、式(12)含有三角函数，仍然具有较强的非线性。为此定义

(29)

(30)

从而

(31)

(32)

其中：

(33)

(34)

以上定义中隐含一个附加约束：

(35)

可以松弛为

(36)

在适当的条件下最优解将能够在式(36)的边界上取到，这就使得松弛约束成为原约束，否则松弛问题的解将不再是原问题的解。具体地说，注意到球面两点间大圆弧最短，由于有

(37)

式中：R_min为连接起点和终点的大圆弧的长度。式(37)是终点约束所隐含的，它倾向于使式(36)靠近边界。如果总弧长s_f估计比较准，在终端位置约束的作用下，可以保证式(36)几乎总在边界上取到，否则适当减小s_f即可。事实上，由于R_min是s_f的低估，并且通常误差不大，实际计算表明以R_min作为s_f的初值就能基本保证松弛问题的解也同时是原问题的解。
至此，新的动力学方程可以记为

(38)

(39)

(40)

(41)

式中：新的状态和新的控制分别为

(42)

可以看出，B(x)为常矩阵，通过选择合适的坐标网格可以使A(x)接近常数。C(x)相对其他两项为小量，其变化对可行性和收敛性几乎没有影响。所以第k+1迭代求解时可以使用^[16-17]

(43)

2.2 过程约束进一步处理 动压和热流约束已经是完全线性的，因而是凸的。过载约束不是凸的，但是可以松弛成：

(44)

式中：M_z为Lipschitz常数，满足

(45)

M_z估值过大可以造成约束作用不明显，降低收敛速度。注意到成立以下不等式(a≥b):

(46)

可以取

(47)

其中：z_min(s)≤z(s)，是一个下界估计，可以离线给出。注意到不等式(45)右边含有绝对值符号，这将会导致混合整数规划而不是凸规划，严重降低计算效率。为此可以选择z(s)的一个上界曲线z_max(s)≥z(s), 作为初值，则绝对值符号可取消。

(48)

再次迭代时z(s)已经得到矫正，比较接近真值，所以可以继续采用

(49)

这样相比将比基于一阶泰勒展开的形式有两个主要优势：一是式(48)相对原不等式是松弛，保证了初次迭代的成功性；二是M_z是导数值的适当高估，提高了收敛性和收敛范围。
约束式(26)可以进一步改写为

(50)

(51)

式中：C_Dminexp(-βh)/ρ₀为指数锥，因而是凸约束。而C_Dmaxexp(-βh)/ρ₀可以采用类似式(48)的处理技巧，即

(52)

(53)

(54)

式中：h_max和h_min分别为h(s)的一对上下界，它们可以离线地给出，并且估计得越精确求解速度越快。
不等式(27)可以改写为

(55)

(56)

暂时略去式(55), 将式(56)处理成

(57)

而式(28)现在成为

(58)

是一个二阶锥约束，因而是凸的。
2.3 禁飞区约束的处理 禁飞区作为一个有限空间，无论其具体形式为何，它所代表的约束总是非凸的，这将给子问题的求解带来困难。如果禁飞区很密，以至于可行解存在的范围很狭窄而未知，此时很难有找到可行解的有效方法。本文针对的是禁飞区相对系稀疏的情况，这时可以采用“割平面法”处理。
步骤1??首先忽略禁飞区约束规划轨迹，如果得到的轨迹恰好不经过禁飞区，则求解完成，否则转入下一步。
步骤2??选择轨迹上离禁飞区最近的一点p, 设它对应的轨迹弧长为s_c, 在禁飞区的边界上取一个平面(即割平面)，添加约束使弧长为s_c处对应的点和禁飞区不在平面的同一侧。然后转入步骤1。典型的平面圆形禁飞区对应的割平面如图 1所示。图中:θ_p和?_p分别为p点的经纬度，k_θ?为斜率。

图 1 禁飞区和割平面 Fig. 1 No-fly zone and cutting plane

图选项

步骤3??迭代进行直到所有的禁飞区约束得到满足。
2.4 优化指标的选择 对于在线轨迹规划来说，快速求出可行解比花费长时间去寻求某种最优解更加重要。因而可以选择特定的有利于求解的优化指标。注意到整个问题的非线性主要来源于大气密度随着高度的剧烈变化，为此可以选择u₂在轨迹上的平方积分为指标，以抑制高度的振荡。

(59)

为了克服气动模型不准确带来的问题，可以在优化指标中加入气动模型的拟合项。考虑C_D、C_L都是迎角α和马赫数Ma的函数:

(60)

(61)

式(60)、式(61)消去迎角，设得到的方程为

(62)

用控制量u表示可得

(63)

考虑采用F_aero的平方积分评估气动模型的误差。由于F_aero是非线性的，通常也不具备凸性，需要对其进行局部线性化，以便于采用凸优化求解。

(64)

从而加入气动模型拟合项并适当加权之后的优化指标可以写成

(65)

式中：ε为权系数，是事先给定的正数。因为希望通过改变控制量而不是高度来匹配气动模型，所以只加入了关于控制的偏导数。另外Ma对方程的影响较小，并且Ma(s)曲线可以粗略地事先估计出来，故可以在迭代中保持为固定。
在气动拟合项和粗略的气动约束共同作用下，求出的实际解将非常接近给定的气动模型。虽然存在一定误差，但是考虑到气动模型本身就有不小的误差，因而这样做不会影响最终解的实用性。另外如果F_aero可以通过在线气动参数辨识来实时估计，则本文的方法就可以几乎不经修改地与之结合以提高制导精度。
2.5 轨迹弧长初值和迭代矫正 由于高度变化相对于射程是小量，弧长主要取决于射程，后者可以通过球面三角理论估计：

(66)

式中：s_f⁰为轨迹弧长的初始估计；L_Rf为归一化待飞航程；?_f和θ_f分别为终点的纬度和经度；?₀和θ₀分别为起点的纬度和经度。
设真实弧长为s_f。为了得到其矫正估计，设第k次迭代使用的弧长为s_f^k, 考虑弧长重参数化

(67)

则终端约束无论部分还是全部给定，都可以写成

(68)

式中：x_f和x₀分别为状态变量的末值和初值。
从而

(69)

所以第k次迭代时可以引入弧长松弛量ζ，并使用如下终点约束和优化指标：

(70)

(71)

那么下一次迭代就可以取矫正弧长：

(72)

3 轨迹规划算法 3.1 问题离散 为了求得无限维规划子问题的数值解，有必要对其进行离散。注意到通常的数值积分格式都可以用来构造微分方程的离散形式，但是为了简单起见，这里仍然采用梯形积分。与通常不同的是，控制量由样条曲线上控制点的表达式表示, 而状态量则像通常一样离散成等距的值型点。计算表明这样做能十分有效的抑制数值求解中的锯齿化现象，同时减少需要求解的变量个数。具体来说，取m + 3个样条控制点代表控制曲线, 取N=m×n个等距的状态量节点, 比率N控制了状态节点的相对稠密程度。在第k次迭代中，对于每连续4个样条控制控制点构造如下约束：

(73)

式中：

(74)

(75)

(76)

(77)

同时相应的路径约束式(20)、(21)、(36)、(48)、(50)、(52)、(57)、(58)和优化指标式(71)也做相应离散。而边界条件则成为

(78)

离散后的子问题是完全凸的，即使规模很大也能高效的求解。
3.2 算法迭代流程 先分别给出轨迹高度和对数速度的一对上下界, h_min、h_max、z_min和z_max。
步骤1??k=0, 取h⁰=h_max, z⁰=z_max, 按式(66)计算弧长初值s_f⁰, 求解离散子问题。
步骤2??以上次计算的结果为初值，并按照式(72)更新弧长，迭代计算直到收敛。判定收敛的条件为：。其中, δ、δ_aero为事先给定的小整数。
步骤3??检验计算结果。主要检验松弛问题是否收敛于原问题，即|1-(ψ_x²+ψ_y²)| ≤δ_relax。其中, δ_relax为事先给定的小整数。
4 数值仿真 为了检验本文算法的有效性，以X-33再入任务为例进行数值测试。其质量m=104 305 kg, 气动参考面积A_ref=391.22 m², 仿真采用的任务参数见表 1，气动模型可以写成如式(79)、式(80)形式^[13]。计算得到的轨迹如图 2~图 9所示。其中横坐标s是归一化轨迹弧长。

(79)

(80)

表 1 主要任务参数 Table 1 Main task parameters

参数	数值
初始高度h0/km	70
初始经度θ0/(°)	0
初始纬度?0/(°)	0
初始速度V0/(m·s-1)	7 000
初始速度倾角γ0/(°)	-1.5
初始航向角ψ0/(°)	0
初始倾侧角σ0/(°)	0
终点高度hf/km	25
终点经度θf/(°)	12
终点纬度?f/(°)	70
终点速度Vf/(m·s-1)	700
最大热流/(kW·m-2)	1 500
最大动压qmax/kPa	18
最大过载nmax/g0	2.5

表选项

图 2 高度曲线 Fig. 2 Altitude curves

图选项

图 3 速度曲线 Fig. 3 Speed curves

图选项

图 4 水平轨迹 Fig. 4 Horizontal trajectory

图选项

图 5 动压曲线 Fig. 5 Dynamic pressure curves

图选项

图 6 热流曲线 Fig. 6 Heat flux curves

图选项

图 7 过载曲线 Fig. 7 Overload curves

图选项

图 8 迎角曲线 Fig. 8 Angle of attack curves

图选项

图 9 倾侧角曲线 Fig. 9 Angle of bank curves

图选项

仿真时采用MATLAB编程环境，以yalmip^[18]+ Mosek9^[19]作为凸子问题求解器，直接以水平线为初值，大约迭代3~5次就基本收敛，表现出很好的潜力。尽管每次迭代仍然需要约0.3 s(i5+8GB RAM), 但是由于凸优化是多项式复杂的，通过有针对性地对该问题编程调优，有望达到远高于此的效率。
从图 10、图 11可以看出，最终的解和气动模型符合的较好, 并且松弛问题收敛于原问题。

图 10 气动模型匹配误差 Fig. 10 Aerodynamic model matching error

图选项

图 11 松弛误差 Fig. 11 Relaxation error

图选项

上述算例采用了常值函数作为初值。由于算法本身由良好的收敛性，实际计算中更精确的初值对收敛速度的影响并不显著。相比之下，下界估计h_min、z_min, 对收敛速度的影响更为明显，因为更好的下界估计会得到更精确的M_z、M_h，从而使不等式(48)和(52)更加准确。算例中取h_min=h_f、z_min=z_f，如果取h_min=(1-s/s_f)(h₀-dh)+(s/s_f)h_f，将只需要迭代2~4次即可收敛。(h₀-dh)是一个略低于初始高度的值。
5 结论 1) 基于轨迹弧长和对数速度的动力学方程相对以往的再入动力学模型，其非线性减弱且相应过程约束大为简化，具有明显优势。
2) 以气动系数与大气密度的乘积作为控制量，绕开了复杂的特定气动模型，再通过在优化指标中加入气动匹配项，这样做在计算上更为简单，更具有普适性。
3) 通过引入弧长自适应参数，避免了弧长估计不准造成求解失败，保证迭代过程能够继续下去，同时又能自动更新得到更准确的弧长。
4) 通过割平面法处理禁飞区约束，建立了方便使用的统一框架，同时兼顾了计算效率。

参考文献

[1]	ISTRATIE V.Optimal profound entry into atmosphere with minimum heat and constraints: AIAA-2001-4069[R].Reston: AIAA, 2001.
[2]	BOLLINO K, OPPENHEIMER M, DOMAN D.Optimal guidance command generation and tracking for reusable launch vehicle reentry: AIAA-2006-6691[R].Reston: AIAA, 2006.
[3]	REA J.Launch vehicle trajectory optimization using a Legendre pseudospectral method: AIAA-2003-5640[R].Reston: AIAA, 2003.
[4]	黄长强, 国海峰, 丁达理. 高超声速滑翔飞行器轨迹优化与制导综述[J]. 宇航学报, 2014, 35(4): 369-379. HUANG C Q, GUO H F, DING D L. A survey of trajectory optimization and guidance for hypersonic gliding vehicle[J]. Journal of Astronautics, 2014, 35(4): 369-379. DOI:10.3873/j.issn.1000-1328.2014.04.001 (in Chinese)
[5]	HARPOLD J C, GAVERT D E. Space shuttle entry guidance performance results[J]. Journal of Guidance, Control, and Dynamics, 1983, 6(6): 442-447. DOI:10.2514/3.8523
[6]	ROENNEKE A J, MARKL A. Re-entry control to a drag-vs-energy profile[J]. Journal of Guidance, Control, and Dynamics, 1994, 17(5): 916-920. DOI:10.2514/3.21290
[7]	LEAVITT J A, MEASE K D. Feasible trajectory generation for atmospheric entry guidance[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(2): 473-481.
[8]	李昭莹, 黄兴李, 李惠峰. 基于闭环解析解的可重复使用运载器轨迹在线生成方法[J]. 宇航学报, 2013, 34(6): 755-762. LI Z Y, HUANG X L, LI H F. On-line trajectory generation method based on closed-form analytical solution for reusable launch vehicle[J]. Journal of Astronautics, 2013, 34(6): 755-762. DOI:10.3873/j.issn.1000-1328.2013.06.003 (in Chinese)
[9]	胡锦川, 张晶, 陈万春. 高超声速飞行器平稳滑翔弹道解析解及其应用[J]. 北京航空航天大学学报, 2016, 42(5): 961-968. HU J C, ZHANG J, CHEN W C. Analytical solutions of steady glide trajectory for hypersonic vehicle and planning application[J]. Journal of Beijing University of Aeronautics and Astronautics, 2016, 42(5): 961-968. DOI:10.13700/j.bh.1001-5965.2015.0330 (in Chinese)
[10]	BOYD S, VANDENBERGHE L. Convex optimization[M]. Cambridge: Cambridge University Press, 2004: 157-160.
[11]	NESTEROV Y, NEMIROVSKII A. Interior-point polynomial algorithms in convex programming[M]. Philadelphia: Society for Industrial and Applied Mathematics Press, 1994: 102-137.
[12]	LIU X, SHEN Z, LU P. Entry trajectory optimization by second-order cone programming[J]. Journal of Guidance, Control, and Dynamics, 2016, 39(2): 227-241.
[13]	WANG Z, GRANT M J.Constrained trajectory optimization for planetary entry via sequential convex programming: AIAA-2016-3241[R].Reston: AIAA, 2016.
[14]	BALAS E, CERIA S, CORNUéJOLS G. A lift-and-project cutting plane algorithm for mixed 0-1 programs[J]. Mathematical Programming, 1993, 58(1-3): 295-324. DOI:10.1007/BF01581273
[15]	MARCHAND H, MARTIN A, WEISMANTEL R, et al. Cutting planes in integer and mixed integer programming[J]. Discrete Applied Mathematics, 2002, 123(1-3): 397-446. DOI:10.1016/S0166-218X(01)00348-1
[16]	TOMáS-RODRíGUEZ M, BANKS S P. Linear, time-varying approximations to nonlinear dynamical systems:With applications in control and optimization[M]. Berlin: Springer, 2010: 112-131.
[17]	BANKS S P, DINESH K. Approximate optimal control and stability of nonlinear finite-and infinite-dimensional systems[J]. Annals of Operations Research, 2000, 98: 19-44. DOI:10.1023/A:1019279617898
[18]	LOFBERG J.YALMIP: A toolbox for modeling and optimization in MATLAB[C]//2004 IEEE International Conference on Robotics and Automation.Piscataway, NJ: IEEE Press, 2004: 284-289.
[19]	ANDERSEN E D, ANDERSEN K D.The MOSEK interior point optimizer for linear programming: An implementation of the homogeneous algorithm[M]//HANS F, KEES R, TAMáS T.High performance optimization.Berlin: Springer, 2000: 197-232.