高超声速飞行器飞行速度高，机动范围大，具备入轨和再入大气的能力，可实现远程快速物资运送或精确打击，具有较高的经济军事价值，近年来已成为各国的研究热点^[1]，其中具有代表性的为美国通用空天飞行器(Common Aero Vehicle, CAV)。再入飞行段作为此类飞行器重要的飞行阶段，由于飞行器的力热载荷环境严酷，且对于控制量高度敏感，因此再入轨迹优化设计也面临着一定挑战。
此类问题可以描述为具有多种复杂约束的非线性连续时间最优控制问题。在近20年里，直接配点法已经成为数值求解非线性最优控制问题的有力工具之一。直接配点法的基本思想是在所研究时域内合理的选取一系列不同的点，对状态变量和控制变量进行离散并建立一定的约束关系，将连续时间最优控制问题转化为有限维的非线性规划(Nonlinear Programming, NLP)问题，从而通过著名的NLP求解器，例如SNOPT^[2]、IPOPT^[3]、KNITRO^[4]等进行有效求解。相比于间接法，直接配点法无需推导Pontryagin极小值原理的最优控制一阶必要条件^[5]，降低了初值敏感性，减小了求解的难度，在飞行器轨迹优化等领域得到了广泛应用。
起初，直接配点法中发展了h方法，例如Runge-Kutta方法等，其通过将时间区间划分为多个子区间，在每个子区间上利用固定阶数的多项式来近似状态变量和控制，通过增加子区间数量来实现h方法的收敛^[6]。
近年来，大量的研究都聚焦于采用一系列Gauss正交配点的方法^[7-10]，称为伪谱法或正交配点法。由于高精度Gauss积分规则，LG(Legendre- Gauss)、LGR(Legendre-Gauss-Radau)以及LGL(Legendre-Gauss-Lobatto)这3种配点类型被广泛采用，它们分别是Legendre多项式的根、Legendre多项式线性组合的根以及Legendre多项式导数的根。伪谱法属于p方法，通常基于Gauss正交规则，在整个时间区间选取一系列配点，通常采用Lagrange全局插值多项式近似状态变量和控制变量，通过配置微分代数方程来近似微分约束，将连续时间最优控制问题转换为NLP问题。通过增加多项式的阶数来实现p方法的收敛^[11]。
自适应伪谱法^[12-14]结合了h方法和p方法的特点，通过估计当前解的相对误差，按照一定规则调整区间数量和多项式阶数，相比单独使用h方法或p方法，提高了求解的精度及效率，并能够有效应对状态控制等不连续的问题。
本文针对高超声速飞行器再入轨迹优化问题，建立了三自由度再入运动方程，以CAV-H飞行器为对象，考虑主要约束条件，基于Legendre多项式近似理论，建立了LGR伪谱法误差估计关系式，采用自适应网格重构技术，实现了3种典型性能指标再入轨迹优化问题的高效求解。本文算法相对传统自适应伪谱法具有配点和区间分配合理，收敛速度快及对人工参数不敏感等优点。
1 再入轨迹优化问题 1.1 三自由度再入运动方程 本文考虑地球为旋转圆球，飞行器无侧滑情况下，建立无动力高超声速飞行器再入运动方程：

(1)

式中：h、V、θ、?、γ、ψ、m、ω_e、σ和r分别为高度、速度、经度、纬度、航迹角、航向角、飞行器质量、地球自转角速度、倾侧角和地心距；引力常数μ=3.986×10¹⁴ m³/s²；地球半径R_e=6 378 137 m；L和D分别为升力和阻力，其表达式为

(2)

其中：ρ为大气密度；海平面大气密度ρ₀=1.2256kg/m³；归一化高度H₀=7 110 m；升力系数C_L和阻力系数C_D由飞行器外形和飞行状态决定；S为参考面积。本文以CAV-H为例，相关总体和气动参数来自文献[15]，将气动系数拟合为马赫数Ma和迎角α的函数。

(3)

1.2 约束条件 1) 端点约束
为满足再入段末端能量管理要求，对末端高度h_f、速度V_f及航迹角γ_f提出限制：

(4)

2) 过程约束
高超声速条件下载荷环境严酷，考虑飞行器防热、结构强度、控制能力等因素，对飞行过程驻点热流密度、法向过载n_y、动压q提出限制：

(5)

式中：K_Q为与飞行器外形和材料相关的常数，本文取K_Q=2.5820×10^-7; g₀为海平面重力加速度。
3) 控制约束
为减轻姿态控制系统压力，要求迎角和倾侧角变化较为平缓，同时对其大小提出限制：

(6)

为了方便描述控制约束，本文控制变量选择为，并引入迎角α和倾侧角σ为状态变量，则状态变量为[h??θ?????V??γ??ψα??σ]^T。
1.3 目标函数 考虑3种典型再入段轨迹优化目标，即最大纵程，最大横程及最小航迹角变化。通过引入权重系数表征设计偏好，则一般性目标函数可以表示为

(7)

式中：w₁、w₂、w₃为权重系数；t₀为初始时刻；t_f为末端时刻。
2 Bolza型连续时间最优控制问题 考虑如下Bolza型连续时间最优控制问题：

(8)

式中：x(t)∈R^N_x为状态变量，N_x为状态变量维数；u(t)∈R^N_u为控制变量，N_u为控制变量维数。Mayer型性能指标Φ，Lagrange型性能指标L，动力学微分约束f，路径约束C及端点约束φ定义为

(9)

其中：N_C和N_φ分别为路径约束和端点约束维数。
3 多区间LGR伪谱法 多区间伪谱法的思想是将整个时间区间划分为多个小的子区间，在每个子区间内采用较少数量的配点来离散最优控制问题。全局伪谱法则可认为是只有一个区间的特例。不失一般性，将t∈[t₀, t_f]分为s(s≥1)个子区间S_k=[t_k－1, t_k]，其中:k=1, 2, …, s; t₀ < t₁ < … < t_s=t_f。引入τ∈[－1, 1]，则可将t映射至[－1, 1]，即

(10)

LGR伪谱法的配点采用τ∈[－1, 1]上的N个LGR点或反转的LGR点，本文选用前者。LGR点为(τ+1)(P_N(τ)+P_N－1(τ))的根，其中:P_N(τ)为N阶Legendre多项式。因此对t进行区间转化，可使每个子区间满足Legendre正交多项式的定义区间。
LGR配点和Gauss-Radau积分权重的计算方法主要为特征值方法和牛顿迭代法^[16]。本文采用特征值方法进行求解，N－1阶Jacobi矩阵为

(11)

式中：

(12)

A_N－1的第i个特征值为λ_i，对应单位特征向量的首个元素为v_1i。根据Radau点与Gauss点的对应关系^[16]，可得到LGR配点τ_i及积分权值ω_i为

(13)

式中：i=1, 2, …, N_k-1。
在每个时间子区间内采用Lagrange基函数近似状态变量和控制变量。状态变量的插值节点为N_k个LGR点和τ_Nk=1点，即插值节点数比配点数多1个；控制变量的插值节点为N_k个LGR点，末端控制点通过外插得到。因此第k个子区间的状态变量和控制变量可近似表示为

(14)

式中：

(15)

状态变量在配点处对τ的微分可表示为

(16)

式中：D^N_k×(N_k+1)为Radau微分矩阵，表示各Lagrange基函数在各LGR配点处的微分值，其表达式为^[16]

(17)

(18)

至此最优控制问题(8)可以离散为微分形式的NLP问题：

(19)

4 自适应网格重构技术 自适应伪谱法的基本思想：基于当前解的最大相对误差来评估解的质量，当不满足求解精度要求时，按照一定规则调整子区间内多项式阶数或调整子区间数量，以达到改善解精度的目的。
4.1 误差评估 文献[13]给出了一种近似相对误差的算法，其基本思想是通过对比2种状态近似解，从而构建相对误差的近似表达式。
在当前解的每个子区间S_k=[t_k－1, t_k]内增加1个配点，即配点数为N_k+1，定义新的配点序列为高阶网格，其中，，定义原始配点序列(τ₀^(k), τ₁^(k)…, τ_Nk^(k))为低阶网格。通过Lagrange插值获得低阶网格的状态变量和控制变量在高阶网格处的值X^(k)()和U^(k)()，i=1, 2，…, N_k+1。对于LGR伪谱法，其积分形式与微分形式等价^[9]，因此可基于积分矩阵构建高阶网格的状态近似解：

(20)

式中：为Radau积分矩阵，其表达式为

(21)

其中：为微分矩阵去除第1列元素后的方阵。
此时可以构建第i维状态变量绝对误差和相对误差的近似表达式为

(22)

则S_k内最大相对误差为

(23)

4.2 Legendre多项式近似 根据收敛性理论^[17]，对于光滑问题，全局伪谱法以指数形式收敛。文献[17]指出，对于分段光滑函数，采用Legendre多项式的近似误差与多项式系数具有相同的衰减率，且衰减率与多项式的阶数相关。因此可以基于Legendre多项式系数构建多项式阶数与近似误差的关系。
对于分段光滑有界函数y(τ)可以展开为Legendre级数：

(24)

式中：p_i为Legendre系数；P_i(τ)为Legendre多项式。取N阶近似为

(25)

则截断误差为

(26)

其中：为N阶近似的Legendre系数。
文献[17]表明p_i和通常非常相近，因此截断误差的第2项相比于第1项可近似忽略。
根据Legendre多项式的正交性：

(27)

式中：δ_mn为Kronecker函数。则式(26)化简为

(28)

根据文献[17]分段光滑有界函数y(τ)的Legendre系数近似为

(29)

式中：c为一常数；为衰减系数。
式(29)的形式不利于构建近似误差关系式，注意到存在一个略小于的υ使得

(30)

因此可根据离散Legendre近似形式：

(31)

采用最小二乘法求解c和υ。
同时考虑到级数关系式：

(32)

因此可构建近似误差关系式：

(33)

4.3 光滑性判定 假设子区间S_k内e_max^(k)不满足给定容许误差ε。采用4.2节算法计算υ，当υ较大时, 表明函数光滑性较好，采用单段多项式近似程度良好；当υ较小时, 表明函数光滑性较差，有必要采用分段多项式来改善近似程度。因此给定一正参数υ_ε，当υ>υ_ε时，需要增加多项式阶数；当υ≤υ_ε时，需要对当前子区间进行分段。
4.4 多项式阶数更新 当υ>υ_ε时，根据式(34)可计算当前网格的求解误差：

(34)

假设通过增加多项式阶数至N_k^(new)使得最大相对误差满足给定容许误差：

(35)

根据式(26)和(27)可求解N_k^(new)：

(36)

为了严格增加配点数量，采用向上取整算子:

(37)

4.5 区间更新 当υ≤υ_ε时，考虑到υ接近于0时，由式(33)确定区间分段数会产生很大的区间数量，从而导致NLP问题的规模过大，使得求解效率降低。因此采用υ_ε代替υ，根据式(36)估计所需的多项式阶数：

(38)

令每个新子区间与当前子区间的配点数相同，则区间分段数H_k为

(39)

4.6 网格缩减 对于子区间S_k内e_max^(k)满足给定容许误差ε，文献[14]提出了一种减少子区间配点数及合并相邻等配点数子区间的算法。
减少区间配点数的基本思想是将Lagrange插值多项式展开为τ的幂级数形式，并计算各阶次系数。为了减少该多项式阶数同时保证相对误差仍满足要求，考虑到τ∈[－1, 1]，令τ=τ_max=1，因此仅需从最高阶次系数向最低阶次系数依次累加，直至累加值超过容许误差ε，则可减少的配点数不超过累加次数。
合并相邻子区间的基本思想与上述算法类似，只是将相邻子区间的值在当前区间对应多项式上进行外插，然后类似地判断容许误差是否满足。
网格缩减算法的具体推导过程见文献[14]，本文不再赘述。
4.7 算法流程 本文自适应伪谱法的计算流程如下：
步骤1??给定求解容许误差ε和初始网格。
步骤2??求解NLP问题(19)。
步骤3??采用式(23)计算每个子区间的最大相对误差e_max^(k)，若max(e_max^(k))≤ε，则求解完成；否则转至步骤4。
步骤4??对每一个子区间，若e_max^(k)>ε，采用4.4节和4.5节算法更新网格；否则采用4.6节算法缩减网格。
步骤5??重复步骤2。
综上，本文算法的计算流程图见图 1。

图 1 自适应伪谱法流程图 Fig. 1 Flowchart of adaptive pseudospectral method

图选项

5 算例分析 本文以CAV-H飞行器为例，飞行器质量为m=907.2 kg，参考面积为S=0.484 m²，飞行过程中最大热流密度=6 000 kW/m²，最大法向过载n_y?max=4，最大动压q_max=50 kPa，迎角0°≤α≤25°，倾侧角σ≤85°，迎角和倾侧角变化率分别为≤5(°)/s，≤5(°)/s。
再入轨迹优化算例设置见表 1，边界条件设置见表 2。计算环境为Windows XP 2.83 GHz，3.25 GB内存，软件环境为MATLAB R2014a，NLP求解器采用SNOPT v7.0，梯度计算采用有限差分算法。
表 1 算例设置 Table 1 Setup of calculation examples

算例	状态	权重系数
1	最大纵程	w1=1，w2=0，w3=0
2	最大横程	w1=0，w2=1，w3=0
3	最小航迹角变化	w1=0，w2=0，w3=1

表选项






表 2 边界条件设置 Table 2 Setup of boundary conditions

参数	初始	末端
高度h/km	60	20
经度θ/(°)	0	自由
纬度?/(°)	0	自由
速度V/(km·s-1)	5	1
航迹角γ/(°)	0	5
航向角ψ/(°)	90	自由

表选项






为了验证本文自适应伪谱法求解的有效性与快速性，采用本文算法(记为hpL)求解最优控制变量，并采用变步长Runge-Kutta-Fehlberg法(记为RKF45)进行积分验证，积分相对误差为10^-13。对比算法分别采用文献[12]算法(记为hp)以及文献[13]算法(记为ph)。本文算法轨迹优化结果见图 2~图 6，不同算法对比结果见表 3。

图 2 3D轨迹 Fig. 2 Three-dimensional trajectory

图选项

图 3 航迹角变化 Fig. 3 Path angle versus time

图选项

图 4 再入走廊边界 Fig. 4 Boundary of reentry corridor

图选项

图 5 迎角与倾侧角变化 Fig. 5 Angle of attack and heeling angle versus time

图选项

图 6 迎角与倾侧角变化率变化 Fig. 6 Change rate of angle of attack and heeling angle versus time

图选项

表 3 不同算法轨迹优化结果 Table 3 Trajectory optimization results by different methods

算法	算例	误差容限ε	衰减容限υε	区间配点限制	配点总数	区间总数	网格迭代次数	平均求解时间(5次)/s	性能指标
hpL	1	10-4	0.25	[3, 10]	130	33	3	3.520 1	-0.112 9
hpL	1	10-4	0.5	[3, 10]	122	37	4	3.878 5	-0.112 8
hpL	1	10-5	0.25	[3, 10]	209	30	5	7.467 0	-0.112 8
hp	1	10-4		[3, 10]	132	39	5	4.462 8	-0.112 8
hp	1	10-5		[3, 10]	250	58	13	56.272 2	-0.112 8
ph	1	10-4		[3, 10]	105	20	4	42.716 9	-0.112 6
ph	1	10-4		[3, 6]	114	28	4	4.461 6	-0.112 6
ph	1	10-5		[3, 6]	273	67	7	10.378 1	-0.112 9
hpL	2	10-4	0.25	[3, 10]	79	20	4	3.270 8	-0.014 8
hpL	2	10-4	0.5	[3, 10]	96	26	4	3.784 4	-0.014 8
hpL	2	10-5	0.25	[3, 10]	126	23	3	2.997 0	-0.014 8
hp	2	10-4		[3, 10]	82	23	4	3.059 3	-0.014 8
hp	2	10-5		[3, 10]	155	40	7	10.297 3	-0.014 8
ph	2	10-4		[3, 10]	79	20	4	6.972 8	-0.014 8
ph	2	10-4		[3, 6]	81	22	4	5.901 6	-0.014 8
ph	2	10-5		[3, 6]	152	34	5	5.333 2	-0.014 8
hpL	3	10-4	0.25	[3, 10]	65	20	4	11.745 5	0.467 8×10-4
hpL	3	10-4	0.5	[3, 10]	65	20	4	12.306 4	0.4678×10-4
hpL	3	10-5	0.25	[3, 10]	91	21	6	15.380 1	0.467 2×10-4
hp	3	10-4		[3, 10]	70	22	4	34.250 6	0.467 2×10-4
hp	3	10-5		[3, 10]	86	26	8	19.334 9	0.467 2×10-4
ph	3	10-4		[3, 10]	66	20	3	83.696 5	0.467 8×10-4
ph	3	10-4		[3, 6]	69	22	4	21.790 3	0.467 8×10-4
ph	3	10-5		[3, 6]	89	26	5	17.796 7	0.467 2×10-4

表选项






由图 2~图 4可以发现，对于3种典型算例，采用hpL算法的求解结果与采用变步长RKF45算法积分的结果均保持一致。优化获得的飞行器再入轨迹末端均严格满足末端高度约束，速度约束以及航迹角约束。由图 4可以发现，h-V平面内, 3种算例的最优轨迹均位于再入走廊下边界的上方，即全程满足飞行过程中的热流率约束、法向过载约束以及动压约束要求。可以注意到，对于本文CAV-H飞行器，法向过载在再入段过程中均不是主导约束，再入初始阶段，飞行器飞行高度高，大气非常稀薄，即使飞行速度较高，相比于动压，热流率为主导约束，由图 4可以发现，此时飞行器采用较大迎角以获得较大升力，从而有效减缓飞行高度的快速下降，以满足热流率要求。当飞行器能量下降至一定时，动压为主导约束。对于最小航迹角变化率问题，整个飞行轨迹平滑无跳跃，初始以较大迎角飞行，飞行速度快速下降以匹配平稳飞行的能量要求。由图 5和图 6可以发现，整个飞行过程中，3种算例的迎角和倾侧角均变化缓慢，其变化率均不超过1(°)/s。从而验证了本文算法的有效性。
由表 3可以发现，本文hpL算法具有较好的稳健性，人工参数衰减容限υ_ε一般取值为(0, 1]之间，本文设置为0.25和0.5时，3种算例的求解时间和网格迭代次数差别不大，原因是采用Legendre衰减系数算法进行光滑性判定时，衰减系数υ的区分度较大，使得人工参数在一定区间内取值时，对增加配点数和分割区间数的选择影响较小，因此轨迹优化效率对于人工参数的选择不敏感。对于hpL算法，增加的配点数取决于衰减系数υ，即与状态平滑性相关，而ph算法可认为是υ≡lg N的算法，当N < 10时，ph算法估计的多项式阶数偏大，其求解效率受到区间最大配点数的影响较大，对于最大射程问题，最大配点数为6的求解时间仅为最大配点数为10的0.1倍，此时ph算法更倾向于分段，从而避免稠密的约束Jacobian矩阵。类似地，对于hp算法可认为是υ≡1的算法。本文phL算法对于3种算例均展现了较好的求解快速性，并且相比于hp算法和ph算法，hpL算法的网格迭代次数更少，收敛速度更快，同时采用网格缩减技术，使得配点和区间总数均相对更少。验证了基于Legendre衰减系数的误差估计的有效性以及本文算法的快速性。
6 结论 本文针对高超声速飞行器再入段轨迹优化问题，以CAV-H飞行器对象，基于多区间Radau伪谱法，对3种典型性能指标的连续时间最优控制问题进行离散，基于Legendre多项式近似理论构建相对误差估计关系式，采用自适应网格重构技术实现了再入轨迹优化的高效求解。仿真结果表明：
1) 本文算法结果与变步长Runge-Kutta-Fehlberg法积分结果一致。
2) 本文算法较2种传统自适应伪谱法，配点与区间分配更合理，计算效率更高，且对于人工参数不敏感，具有较好的工程应用价值。

参考文献

[1]	雍恩米, 陈磊, 唐国金. 高超声速无动力远程滑翔飞行器多约束条件下的轨迹快速生成[J]. 宇航学报, 2008, 29(2): 397-406. YONG E M, CHEN L, TANG G J. A survey of numerical methods for trajectory optimization of spacecraft[J]. Journal of Astronautics, 2008, 29(2): 397-406. DOI:10.3873/j.issn.1000-1328.2008.02.002 (in Chinese)
[2]	GILL P E, MURRAY W, SAUNDERS M A. SNOPT:An SQP algorithm for large-scale constrained optimization[J]. SIAM review, 2005, 47(1): 99-131. DOI:10.1137/S0036144504446096
[3]	W?CHTER A, BIEGLER L T. On the implementation of an interior-point filter line-search algorithmfor large-scale nonlinear programming[J]. Mathematical Programming, 2006, 106(1): 25-57. DOI:10.1007/s10107-004-0559-y
[4]	BYRD R H, NOCEDAL J, WALTZ R A, et al. Knitro:An integrated package for nonlinear optimization[J]. Large-Scale Nonlinear Optimization, 2006, 83: 35-59.
[5]	PONTRYAGIN L S. Mathematical theory of optimal processes[M]. New York: Routledge, 2018.
[6]	BETTS J T. Practical methods for optimal control and estimation using nonlinear programming[M]. Philadelphia: SIAM, 2010.
[7]	BENSON D A, HUNTINGTON G T, THORVALDSEN T P, et al. Direct trajectory optimization and costate estimation via an orthogonal collocation method[J]. Journal of Guidance, Control, and Dynamics, 2006, 29(6): 1435-1440. DOI:10.2514/1.20478
[8]	DARBY C L, GARG D, RAO A V. Costate estimation using multiple-interval pseudospectral methods[J]. Journal of Spacecraft Rockets, 2011, 48(5): 856-866. DOI:10.2514/1.A32040
[9]	GARG D, PATTERSON M, HAGER W W, et al. A unified framework for the numerical solution of optimal control problems using pseudospectral methods[J]. Automatica, 2010, 46(11): 1843-1851. DOI:10.1016/j.automatica.2010.06.048
[10]	PATTERSON M A, RAO A V. Exploiting sparsity in direct collocation pseudospectral methods for solving optimal control problems[J]. Journal of Spacecraft Rockets, 2012, 49(2): 364-377.
[11]	CANUTO C, HUSSAINI M Y, QUARTERONI A, et al. Spectral methods in fluid dynamics[M]. New York: Springer Science & Business Media, 2012.
[12]	DARBY C L, HAGER W W, RAO A V. An hp-adaptive pseudospectral method for solving optimal control problems[J]. Optimal Control Applications and Methods, 2011, 32(4): 476-502. DOI:10.1002/oca.957
[13]	PATTERSON M A, HAGER W W, RAO A V. A ph mesh refinement method for optimal control[J]. Optimal Control Applications and Methods, 2015, 36(4): 398-421. DOI:10.1002/oca.2114
[14]	LIU F J, HAGER W W, RAO A V. Adaptive mesh refinement for optimal control using nonsmoothness detection and mesh size reduction[J]. Journal of the Franklin Institude, 2015, 352(10): 4081-4106. DOI:10.1016/j.jfranklin.2015.05.028
[15]	PHILLIPS T H.A common aero vehicle (CAV) model, description, and employment guide[R].[S.L.]: Schafer Corporation for AFRL and AFSPC, 2003.
[16]	SHEN J, TANG T, WANG L L. Spectral methods:algorithms, analysis and applications[M]. New York: Springer Science & Business Media, 2011: 55-85.
[17]	WANG H, XIANG S. On the convergence rates of legendre approximation[J]. Mathematicas of Computation, 2012, 81(278): 861-877.