高超声速飞行器再入轨迹快速优化研究

董春云，郭 志，赵培博，蔡远利，于振华

(1.西安交通大学 电子与信息工程学院,西安 710049；2.空军工程大学 信息与导航学院,西安 710077)



高超声速飞行器再入轨迹快速优化研究

董春云1，郭 志1，赵培博1，蔡远利1，于振华2

(1.西安交通大学 电子与信息工程学院,西安 710049；2.空军工程大学 信息与导航学院,西安 710077)

基于求解最优控制问题的Chebyshev伪谱法(Chebyshev Pseudospectral Method, CPM)，研究了高超声速飞行器再入轨迹快速优化问题。针对远程多约束条件下再入轨迹优化问题的难点，提出了一种线性初值与节点更新相结合的优化策略，将攻角与倾侧角同时作为控制变量，以再入飞行时间最短为优化目标，利用CPM将轨迹优化问题转化为非线性规划问题，并使用SNOPT软件包求解，使CPM成为一种再入轨迹快速优化的通用算法。以某类高超声速再入飞行器为对象进行轨迹优化计算，并对比相同仿真条件下粒子群(PSO)算法的优化效果，仿真结果验证了该算法具有较高的求解效率和快速收敛性。

高超声速飞行器；轨迹优化；再入；Chebyshev伪谱方法

0 引言

近年来，以通用航空飞行器(Common Aero Vehicle，CAV)[1]为代表的一类高超声速再入飞行器受到各军事强国的高度重视。该类飞行器自身不携带发动机，且具有较大的升阻比，从轨道或亚轨道高度再入，依靠气动力控制实现远距离滑翔再入飞行，突破了常规弹道式再入模式，具有机动能力强、弹道灵活多变难以拦截等诸多优点，在军事上具有广泛的应用前景[1]。

轨迹优化设计是高超声速再入飞行器的关键技术之一，也是目前研究的热点问题。高超声速飞行器再入过程具有飞行速度快、空间跨度大、气动热力环境恶劣、再入轨迹对控制变量高度敏感等特点，为了满足安全飞行的需要，飞行器需满足严格的过载、动压、驻点热流密度等诸多非线性约束。同时，考虑到与末制导段交班的要求，再入终端状态也需要满足一定的约束。诸多约束条件使得再入轨迹的可行域被限制在较为狭窄的范围内，给轨迹优化设计带来了一定的困难。此外，由于升力式飞行器升阻比大，飞行时间过长极易造成再入过程的总吸热量过高，给飞行器的热防护系统设计带来困难，也不利于提高有效载荷的质量。因此，通过合适的优化策略和方法得到满足以上约束条件，且飞行时间最短的再入轨迹是非常有意义的。

上述再入轨迹优化问题是一类典型的具有路径约束和终端约束的最优控制问题，其数值求解方法主要分为间接法和直接法两类。相对于间接法，直接法在解决实际问题的适应性和收敛的鲁棒性上更具优势。依据参数化方法的不同，直接法又分为仅离散控制变量的打靶法和同时离散状态变量和控制变量的配点法。由于再入轨迹对控制变量高度敏感，仅离散控制变量的方法容易陷入局部解，甚至收敛不到可行解，而配点法以其求解精度较高、收敛性好、易实现等优点，在轨迹优化问题的求解中得到了更多的应用。近年来，配点法中的伪谱法以其高精度和高效率等优势逐渐受到重视，其收敛性得到了理论证明[3-4]。根据离散节点和插值基函数选取的不同，伪谱法可分为Gauss伪谱法、Legendre伪谱法、Radau伪谱法和Chebyshev伪谱法(Chebyshev Pseudospectral Method, CPM)等。

Fahroo和Ross于2002年提出了一种用于求解最优控制问题的Chebyeshev伪谱法[5]，但当时并未引起广泛关注。直到2008年，Trefethen证明了Clensshaw-Curtis数值求积公式与Gauss求积公式精度相当[6]；Gong等于2009年～2010年进一步完善了Chebyeshev伪谱法的理论体系，证明了其离散数值解会一致地趋近于原问题的最优解[7]，并建立了协状态估计理论[8]。此后，国内外学者开始重新关注Chebyeshev伪谱法。

与Gauss伪谱法相比，Chebyeshev伪谱法的离散节点Chebyshev-Gauss-Lobatto(CGL)具有显式表达式[9]，使用起来非常方便。而Gauss伪谱法只能通过大量的代数求根运算来获得节点的位置，且不能直接得到端点处的控制量，终端状态必须通过积分获得，这在一定程度上加剧了计算难度，同时也增加了计算耗时。近年来，在最短奔跑时间求解[10]、月球软着陆轨迹优化[11]、桥式吊车最优控制[12]和制导炮弹弹道优化[13]等领域，Chebyeshev伪谱法都以其高效的求解效率和良好的优化性能得到了成功应用。而对于高动态、强约束和非线性的高超声速飞行器再入过程轨迹优化问题，Chebyeshev伪谱法是否适用，文中进行了详细的讨论。

文中研究了基于Chebyshev伪谱法的高超声速飞行器再入轨迹快速优化问题。首先，给出了再入飞行器归一化动力学模型和相关约束条件，建立了以最短再入时间为优化目标的多约束轨迹优化模型；在此基础上，详细阐述了利用Chebyshev伪谱法对该轨迹优化问题进行离散化求解的具体步骤，针对设计变量初值给定的难点，给出一种兼顾求解精度与实时性的串行求解策略；最后，通过仿真对本文算法的有效性与求解效率进行了验证分析，并与基于另一种直接法-粒子群算法的轨迹优化结果进行了比较。

1 轨迹优化模型

1.1 动力学模型及归一化处理

假设地球为旋转圆球，高超声速飞行器无动力再入的非线性动力学模型可用6个状态变量来描述。由于变量间存在不可公度性，取值范围和量纲不尽相同，为了提高后续优化计算效率与求解精度，需要对模型进行归一化处理。极坐标系下无量纲运动方程[14]为

(1)

(2)

(3)

(sinγcosφ-cosγsinφcosΨ)

(4)

Ω2rcosφ(cosγcosφ+sinγsinφcosΨ)]

(5)

2ΩV(tanγcosφcosΨ-sinφ)+

(6)

式中m、S分别为再入飞行器质量和气动参考面积；ρ为大气密度；CL、CD分别为升力和阻力系数。

气动力的计算采用1976年美国标准大气模型USSA76[15]，在0～100 km高度范围内大气划分为8层，根据飞行高度计算空气密度ρ。CL、CD通过气动力数据拟合得到[16]。

1.2 控制变量

在高超声速飞行条件下，气动力系数可近似表示为攻角的函数，同时将攻角α和倾侧角σ作为优化控制变量[17]，能够充分体现出攻角在再入飞行中的控制调节作用。因此，控制变量为

u=(ασ)T

1.3 约束条件

(1)过程约束

飞行器再入飞行是一个高动态过程，要求严格满足过程约束，以保证飞行器在结构和热防护上的可靠性，主要包括热流密度、动压、过载和控制变量约束，即

(7)

(8)

(9)

(10)

式中R为飞行器头部曲率半径；ρs为海平面大气密度；C1为常数。

(2)终端约束

飞行器再入终端状态对末制导段的飞行有重要的影响，终端约束与再入飞行的任务相关，一般包括速度、飞行高度和地理位置约束等：

|V(τf)-Vf|≤εv，|r(τf)-rf|≤εr

|θ(τf)-θf|≤εθ，|φ(τf)-φf|≤εφ

(11)

其中，εV，εr，εθ，εφ分别是终端状态的误差界，考虑精确制导时，可设置εV=0，εr=0，εθ=0，εφ=0。

1.4 性能指标

优化目标按照设计要求可以有不同的形式。由于再入环境恶劣，最短的飞行时间能够减小气动加热，降低能量损耗，缩短敌方响应时间。因此，以最短飞行时间为性能指标：

J=min(τf)

(12)

1.5 优化模型

高超声速飞行器无动力再入轨迹优化问题可具体描述为：在[τ0,τf]时间内(τf未知)，确定控制变量u=(ασ)T和终端时刻τf，使得目标函数式(12)最小，并满足系统状态方程约束式(1)～式(6)、过程约束式(7)～式(10)以及终端约束式(11)。

将Chebyshev伪谱法用于该轨迹优化问题的求解，力求使其成为一种再入轨迹快速优化的通用算法。

2 Chebyshev伪谱法与轨迹优化

2.1 方法描述

对于上述连续时间最优控制问题，Chebyshev伪谱法的基本求解思路为：将连续时间状态变量和控制变量在一系列CGL点上离散，并以这些离散点为节点构造Largrange插值多项式来逼近真实状态和控制；通过对全局插值多项式求导来近似状态变量对时间的导数，将微分方程约束转换为代数约束；性能指标中的积分项由Clenshaw-Curtis数值积分计算。通过上述方法，可将最优控制问题转化为具有一系列代数约束的NLP问题[5,8]。

CPM的具体计算步骤简述如下：

(1)时间区间转换

CPM的CGL离散点定义在[-1,1]区间上，需要将原最优控制问题的时间域τ从[τ0,τf]线性转换到t∈[-1,1]区间，转换式为

τ=[(τf-τ0)t+(τf+τ0)]/2

(2)离散节点计算

CPM离散节点选取为N阶Chebyshev多项式TN(t)=cos(Ncos-1t)的极值点，即CGL点，它们不均匀的分布在[-1,1]区间上，显式计算式为

显然，离散节点tk满足t0=-1，tN=1。

(3)状态变量和控制变量插值近似

取上述N+1个离散点处的状态变量和控制变量，分别构造Lagrange插值多项式作为连续状态和控制的近似。真实状态变量x(t)与控制变量u(t)的近似表达式为

(13)

(14)

其中，Lagrange插值基函数：

由Lagrange插值的性质可知，离散节点处的状态近似值与实际状态相等，控制近似值与实际控制相等。

(4)动态约束处理

对式(13)求导，得到状态向量在tk点处的导数近似表达式为

(15)

式中Dkj为(N+1)×(N+1)微分矩阵D的第k行第j列元素；D的计算公式见文献[5]。

(16)

式中f为式(1)～式(6)所示的状态方程。

对于式(7)～式(10)所示的过程约束，需在离散节点处严格满足。

(5)性能指标积分近似

对于优化性能指标中存在积分项的情况，利用Clenshaw-Curtis数值积分对其进行近似。对于 [-1,1]区间上的任一连续函数p(t)，其积分可用N+1个CGL离散点处的函数值累加和近似，即

(17)

式中，ωk(k=0,1,…,N)为Clenshaw-Curtis加权，其计算式参见文献[5]。

通过以上对最优控制问题在CGL节点的近似，将原连续时间最优控制问题转化为以下离散化的非线性规划问题，即确定CGL节点处的离散状态变量X、离散控制变量U和终端时间τf，使得性能指标J最小，并满足离散化后的状态方程约束、过程约束与终端约束等。

2.2 轨迹优化求解策略

利用SNOPT软件包[18]求解经Chebyshev伪谱法转化获得的非线性规划问题，该软件包基于序列二次规划算法，在处理大规模NLP的求解上具有良好的性能，目前已得到广泛应用。

根据数值验证结果，CGL节点的数量越大，求解精度越高，但计算量与计算耗时也会增加。另一方面，高超声速飞行器再入轨迹优化模型复杂，考虑的约束条件较多，根据本文模型，当CGL离散节点数为N+1时，设计变量数目为6(N+1)+2(N+1)+1，约束总数为6(N+1)+4(N+1)+4，N的数值较大时，设计变量初值的设置会非常繁杂，不恰当的初值会增加计算复杂度与程序耗时，甚至使问题无法收敛到可行解。因此，如何确定CGL节点的数目与设计变量的初值，是优化过程中首先需要解决的问题。为此，提出了一种线性初值与节点更新相结合的优化策略。

如图1所示，首先采用线性化方法获得设计变量初值，求解较少CGL节点下的初步优化数值解；然后，增加CGL节点，利用初步优化数值解插值获得需要的设计变量初值，从而提高求解效率。

所谓线性化方法即根据设计变量初始位置x0与终端约束xf构造简单线性函数，具体可表述为

节点更新过程中，通过循环依次增加节点的数目，计算优化结果的数值积分终端约束偏差，当最大相对误差小δ时，认为精度符合要求，停止计算。图1中，变量nx、ny可根据仿真效果灵活调整。

2.3 CPM数值验证

在应用CPM求解再入轨迹优化问题之前，首先选择标准最速降线问题对本文算法的有效性进行验证，并分析其求解精度。

图1 轨迹优化求解策略Fig.1 Combination strategy for trajectory optimization

最速降线问题即寻找一条最优曲线，使得质点在重力作用下，沿该曲线从给定点到不在它垂直下方的另一点的时间最短。其数学模型可以概括为

式中，控制量θ是曲线斜率随时间变化的函数。

该问题的解析形式为

x(τ)=(gτf/π){τ-(τf/π)sin[π(1-τ/τf)]}

将CPM得到的控制变量插值后代入原运动方程，采用四阶龙格库塔法积分求解运动轨迹，验证方法可行性。CGL节点数目K=N+1，设N=10时对比CPM数值解、积分求解结果及原问题解析解曲线图，如图2与图3所示。结果表明，三者的位置坐标x、y轨迹均能较好吻合，CPM和SNOPT软件包的计算时间短，求解精度高。

图2 最速降线控制变量对比图Fig.2 Time history of control for the brachistochrone problem

为了分析CGL节点数目对CPM求解精度与收敛速度的影响，分别取N=5,10,20，对比CPM数值解与原问题解析解在离散节点处的求解误差以及程序仿真运行时间，结果如表1所示。

从表1看出，3种情况下CPM获得的数值解精度均较高，且随着离散点数的增加，求解精度越来越高，一致趋向于原问题最优解，与文献[7]结论一致。需注意的是，仿真程序耗时随着离散点数增加而增加。

图3 最速降线状态变量对比图Fig.3 Time history of states for the brachistochrone problem

表1 CPM求解精度与仿真时间比较Table 1 Comparison of solution accuracy and computation cost of CPM

3 仿真校验与结果分析

以远程高超声速滑翔式再入飞行为例，研究Chebyshev伪谱法在再入轨迹快速优化问题中的应用。相同仿真条件下，将本文算法的优化结果与粒子群(Particle Swarm Optimization, PSO)算法的优化结果进行对比分析。

PSO算法的详细介绍与具体步骤可参考文献[19-20]，在此不展开讨论。

3.1 参数设定

飞行器结构参数和气动参数参考1998年美国洛克希勒-马丁公司设计的高超声速飞行器CAV-H[16]，该气动参考面积为0.483 9 m2，质量为907.2 kg，长度为2.717 8～3.627 6 m。

3.2 结果及分析

仿真平台是2.33G主频CPU、2G内存的PC机，软件环境是Matlab 2009。初始节点个数取为K=5，采用线性化方法设置初值。节点更新过程中，将获得的控制变量带入动力学方程，利用四阶龙格库塔法进行数值积分，积分步长取为0.05，积分轨迹终端约束最大相对误差阈值取为0.05。

初始优化结果获得的积分轨迹中，终端最大相对误差0.496。更新CGL节点数目，重新优化并积分，当节点个数K=15时，终端最大相对误差为0.117。继续更新节点数目，当K=28时，获得终端高度20.54 km，速度1 018.21 m/s，经度236.36°，纬度36.95°，最大相对误差为0.027，小于最大误差阈值，停止更新。

至此，CPM优化获得的再入飞行时间为1 486.46 s，初始优化耗时约1～2 s，末次优化耗时约为75～80 s。相同仿真条件下，PSO优化程序耗时约16 min，优化得到的再入飞行时间为1 485.11 s。Chebyshev伪谱法较一般直接法(如本文PSO算法)计算效率更高，另一方面，作者在实际应用过程中发现，良好的初值条件能够有效改善优化效率，减少程序耗时。

飞行器的主要轨迹参数，即高度、速度及航迹的优化结果如图4～图6所示。图中，实线为Chebyeshev伪谱法优化结果，双划线为其数值积分结果，虚线表示粒子群算法。可以发现，CPM及PSO 2种方法获得的飞行器高度、速度及航迹曲线变化趋势较为一致，均能较好地收敛到指定终端状态，且优化获得的再入飞行时间非常接近。从程序耗时角度考虑，本文算法运行速度更快，能够满足再入轨迹快速优化的要求。

图4 飞行高度对比Fig.4 Comparison of altitude profile

图5 飞行速度对比Fig.5 Comparison of velocity profile

图6 地面航迹对比Fig.6 Comparison of ground track

从图4～图6的积分曲线可看出，CPM优化结果与数值积分结果基本一致，具有较高的精度，说明本文给出的节点更新策略有效保证了算法优化结果的可行性与有效性。其中，飞行高度数值积分结果在约500～1 000 s的中间阶段略有偏差，而在初始与末尾阶段均保持了较高的积分精度。这是因为CGL节点分布不等距，两端较密集中间较稀疏，且实际飞行中弹道中段高度变化剧烈。此外，由高度轨迹可看出，为了增大射程，从而减小飞行时间，飞行器进行了若干次跳跃滑翔，开始跳跃幅度较大，但在滑翔末段，轨迹趋于平缓，为与末制导的交班提供了良好的条件。

图7给出了控制变量优化结果的对比曲线，包括攻角与倾侧角曲线。从图7可看出，Chebyshev伪谱法与粒子群算法获得的控制变量变化规律有一定的差别。这是由于高超声速飞行器再入环境恶劣、再入过程复杂，不同优化方法得到的飞行轨迹并不一定是严格意义上的最优，而有可能是“次优控制”[21]，这也导致了2种方法在飞行轨迹上存在差异。

Chebyshev伪谱法获得的飞行器再入热流密度、动压和过载严格满足过程约束且数值较小，航迹角与航向角曲线变化均较为平缓，满足设计指标要求，结果如图8～图12所示。从而可看出，再入初期，飞行器速度较高，热流密度约束起主要作用；在随后的滑翔阶段，随着再入高度的降低以及大气密度的增加，动压和过载逐渐起主要作用。

图7 控制变量对比Fig.7 Comparison of control profile

图8 CPM热流密度变化图Fig.8 Heating rate profile

图9 CPM动压变化图Fig.9 Dynamic pressure profile

图10 CPM过载变化图Fig.10 Overload profile

图11 CPM航迹角变化图Fig.11 Flight path angle profile

图12 CPM航向角变化图Fig.12 Heading angle profile

4 结论

将Chebyshev伪谱法用于高超声速无动力再入飞行器的轨迹优化领域，提出了一种线性初值与节点更新相结合的优化策略。仿真结果表明，与一般直接法比较，该方法能够利用较少的离散节点，在短时间内生成满足约束条件的再入飞行轨迹，求解效率高，收敛速度快，且步骤简单，使用方便，是一种有效的高超声速飞行器再入轨迹快速优化算法。

[1] George Richie.The common aero vehicle:space delivery system of the future[C] //AIAA Space Technology Conference and Exposition.Albuquerque,New Mexico,1999.

[2] 谢愈,刘鲁华,汤国建,等.多约束条件下高超声速滑翔飞行器轨迹优化[J].宇航学报,2011,32(12):2499- 2504.

[3] Rao A V.A survey of numerical methods for optimal control[J].Advances in the Astronautical Sciences,2009,135(1):497-528.

[4] Begum Senses,Rao A V.A preliminary analysis of small spacecraft finite-thrust aeroassisted orbital transfer[C] //AIAA/AAS Astrodynamics Specialist Conference.Minneapolis,Minnesota,2012.

[5] Fahroo Fariba,Ross I Michael.Direct trajectory optimization by a Chebyshev pseudospectral method[J].Journal of Guidance,Control,and Dynamics,2002,25(1):160-166.

[6] Trefethen Lloyd N.Is Gauss quadrature better than Clenshaw-Curtis[J].SIAM review,2008,50(1):67-87.

[7] Qi Gong,Ross I Michael,Fahroo Fariba.A Chebyshev pseudospectral method for nonlinear constrained optimal control problems[C]//The 48th IEEE Conference on Decision and Control.Shanghai,China:IEEE,2009:16-18.

[8] Gong Qi,Michael Ross I,Fahroo Fariba.Costate computation by a Chebyshev pseudospectral method[J].Journal of Guidance,Control,and Dynamics,2010,33(2):623-628.

[9] Vlassenbroeck J,Van Dooren R.A Chebyshev technique for solving nonlinear optimal control problems[J].IEEE Transactions on Automatic Control,1988,33:333-340.

[10] Maroński R,Rogowski K.Minimum-time running:a numerical approach[J].Acta of Bioengineering and Biomechanics/Wroclaw University of Technology,2011,13(2):83-86.

[11] Jianhui Z,Peng Q.Lunar soft landing trajectory optimization by a Chebyshev pseudospectral method[C]//Computer Science and Automation Engineering (CSAE),2011 IEEE International Conference on IEEE,2011,2:425-430.

[12] 刘熔洁,李世华.桥式吊车系统的伪谱最优控制设计[J].控制理论与应用,2013,30(8):981-989.

[13] 陈琦,王中原,常思江,等.不确定飞行环境下的滑翔制导炮弹方案弹道优化[J].航空学报,2014,35(9):2593-2604.

[14] Vinh Nguyen X,Adolf Busemann,Robert D.Hypersonic and planetary entry flight mechanics[J].NASA STI/Recon Technical Report A,1980,81:16245.

[15] 张毅,肖龙旭,等.弹道导弹弹道学[M].长沙:国防科技大学出版社,2005:82-85.

[16] Jorris Timothy R.Common aero vehicle autonomous reentry trajectory optimization satisfying waypoint and no-fly zone constraints[D].Wright-Patterson:Air Force Institute of Technology,2007:88-94.

[17] 赵欣,闫循良,张金生,等.助推-滑翔导弹再入弹道快速优化[J].固体火箭技术,2012,35(4):427-433.

[18] Gill Philip E,Walter Murray,Michael A.SNOPT:An SQP algorithm for large-scale constrained optimization[J].SIAM Journal on Optimization,2002,12(4):979-1006.

[19] Rahimi A,Dev Kumar K,Alighanbari H.Particle swarm optimization applied to spacecraft reentry trajectory[J].Journal of Guidance,Control,and Dynamics,2012,36(1):307-310.

[20] Gao,Xiao-zhi,Wu Ying,Wang Xiao-lei,et al.Trajectory optimization in reentry phase for hypersonic gliding vehicles using swarm intelligence algorithms[M].Practical Applications of Intelligent Systems.Springer Berlin Heidelberg,2012:361-371.

[21] 符俊,蔡洪,李安梁.基于 Gauss 伪谱法的航天器气动力辅助平面变轨问题研究[J].国防科技大学学报,2011,33(6):95-99.

(编辑：薛永利)

Rapid trajectory optimization for hypersonic reentry vehicle

DONG Chun-yun1, GUO Zhi1, ZHAO Pei-bo1, CAI Yuan-li1, YU Zhen-hua2

(1.School of Electronic and Information Engineering,Xi'an Jiaotong University,Xi'an 710049,China;2.School of Information and Navigation, Air Force Engineering University,Xi'an 710077,China)

Trajectory optimization of a hypersonic reentry vehicle was investigated via the optimal control method-Chebyshev Pseudospectral Method (CPM). Upon the difficulties of long-range reentry process with many constraints, a combination strategy of the linear initial guess and the nodes update for optimization was adopted. Both the angle of attack and the bank angle were chosen as the control variables to minimize the flight time. The trajectory optimization problem was then translated to a nonlinear programming problem (NLP) via CPM and optimized by the SNOPT software package to obtain a general algorithm for rapid trajectory optimization. In the trajectory optimization of a hypersonic reentry vehicle，simulation results show high efficiency and rapid convergence of the presented strategy compared to the Particle Swarm Optimization (PSO) algorithm.

hypersonic vehicle; trajectory optimization; reentry; Chebyshev pseudospectral method

2014-08-26；

：2014-11-04。

国家自然科学基金(61202128)；宇航动力学国家重点实验室开放基金(2011ADL-JD0202)。

董春云(1989—)，女，博士生，研究方向为飞行器轨迹优化与方法评估。E-mail:dongdong_2007y@163.com

V412.4

A

1006-2793(2015)06-0757-07

10.7673/j.issn.1006-2793.2015.06.002