空天飞机再入轨迹的变精度序列优化方法

朱俊杰，余雄庆，罗东明，王宇

南京航空航天大学 飞行器先进设计技术国防重点实验室，南京 210016

空天飞机再入轨迹的变精度序列优化方法

朱俊杰，余雄庆*，罗东明，王宇

南京航空航天大学 飞行器先进设计技术国防重点实验室，南京 210016

为提高空天飞机再入轨迹优化方法的稳健性和精度，提出一种变精度序列优化方法。优化计算过程分为3个步骤：首先进行预优化，即在粗离散情况下，应用混合优化方法进行轨迹优化，获得各时间点设计变量的最优解；然后对预优化的结果进行一元全区间插值处理，获得细离散情况下设计变量的初值；最后进行精细优化，获得精度更高的计算结果。算例结果表明，变精度序列优化方法能稳健地求解空天飞机再入轨迹优化问题,并能以较少的计算量获得满足精度要求的优化结果。

空天飞机；再入轨迹；优化；混合优化；插值

空天飞机轨迹优化是在飞行动力学和物理学约束的条件下，寻求最优解的过程，它贯穿于整个空天飞机设计过程中，影响着总体、气动布局、制导控制、动力和结构等多个分系统的设计，对空天飞机的研制具有重要意义。

轨迹优化是一类最优控制问题，属于复杂的大规模非线性优化范畴。轨迹优化设计方法主要有基于极大值原理的间接法和基于非线性规划理论的直接法。间接法的优点是计算精度高，并且能够满足最优必要性条件，但由于繁杂的数学推导和难于求解的两点边值问题，目前已很少使用间接法求解轨迹优化问题。随着计算机技术的发展，目前的研究更倾向于使用直接法求解。在采用直接法求解轨迹优化问题时，先将最优控制问题转化为非线性规划问题(Nonlinear Programming, NLP)，对控制变量和状态变量进行离散，对状态方程、端点约束以及路径约束进行节点转化，再选择优化算法求解决策向量和目标函数。

对于优化算法来说，目前主要采用两种优化算法：1)经典优化算法，例如序列二次规划算法(Sequential Quadratic Programming， SQP)[1]、梯度下降法、罚函数方法[2]等；2)智能优化算法，例如遗传算法[3-4](Genetic Algorithm，GA)、进化策略(Evolution Strategy， ES)、粒子群优化(Particle Swarm Optimization， PSO)等。经典优化算法的优点是收敛速度快，但是需要给出合理的初始值，才能收敛到最优解。对于空天飞机再入轨迹问题，有数百甚至数千个设计变量。面对这么多变量，以人工方式给出一个合理的初始值几乎不可能。智能优化算法的优点是无需给出初值，稳健性好，但收敛速度慢，优化结果不够精确。为了克服经典优化算法和智能优化算法的缺点，近年来有学者提出了混合优化方法[5]，即首先用智能优化算法获得一个初步优化结果，然后以该结果为初始点，利用经典优化算法，快速搜寻最优解。如Vavrina采用GA与SQP相结合的方法求解低推力轨迹优化问题[6]。但是，对于空天飞机再入轨迹优化问题，我们发现，当时间配点数较多时，混合优化方法也不能可靠地求解轨迹优化问题，其稳健性不能令人满意。

本文借鉴计算流体力学中预处理思路(变精度网格方法[7])，提出了一种变精度序列优化方法，目的是为空天飞机再入轨迹优化提供一种稳健的、收敛速度快且精度高的方法。

1 再入轨迹优化问题的描述

1.1 运动方程的建立

空天飞机再入大气层时的三自由度运动学和动力学方程组为[8-9]

(1)

式中：m为空天飞机质量；R为空天飞机质心到地心的距离；φ为经度；θ为纬度；v为速度；γ为航迹角；ψ为航迹方向角；β为倾侧角；D为阻力；L为升力；g为重力加速度。

D,L，g,H分别为

(2)

(3)

(4)

(5)

式中：g0=9.8 m/s2为地球表面的重力加速度；ρ为大气密度；H为空天飞机离地面高度；R0=6 357 km为地球半径；S为空天飞机参考面积；CD和CL分别为阻力系数与升力系数。

在大气飞行过程中的热流密度为[10]

(6)

式中：RN为飞行器头部前缘半径。

动压为

(7)

过载为

(8)

1.2 归一化处理

对于再入轨迹优化问题，不同变量之间的数值量级相差较大，例如R是106级别，而经度和纬度化为弧度后是1～10级别，二者数值上相差太大，会导致数值奇异，使优化计算失败。为了增加优化计算的稳健性，需要对原始模型进行归一化处理。令

(9)

升力和阻力按如下方法归一化处理：

(10)

(11)

1.3 空天飞机再入轨迹优化问题

(1)设计变量

对于升力式再入航天器的轨迹优化问题，由于航天器再入后依靠升力在大气层内做较长时间的滑翔，通过改变姿态来控制气动力的变化，以完成轨迹的调整和机动，因此通常将迎角α和倾侧角β作为控制变量。再入航天器的飞行轨迹可由航天器质心到地心的距离R、经度φ、纬度θ、速度v、航迹角γ、航迹方向角ψ等参数来描述，这组参数即为状态变量。控制变量和状态变量均是时间的函数，而状态变量随时间的变化是由微分形式的飞行动力学方程支配的。轨迹优化问题的设计变量由每个时间微段的控制变量和状态变量组成。

(2)优化目标

优化目标的选择主要由航天器的设计及任务要求确定，通常有如下几种优化目标：1)以飞行航程最大作为性能指标；2)以再入过程中总吸热量最小作为性能指标[12]；3)以最小能量消耗作为性能指标[13]；4)以到达指定目标点飞行时间最短为优化的性能指标。

(3)约束条件

再入轨迹优化中的约束包括终端约束、过载、动压及热流峰值约束和控制量约束[12]等。

(4)空天飞机再入段描述

再入段是引导空天飞机从巡航高度转至末端区域能量管理段的无燃料消耗飞行阶段，可视为常质量的无动力滑翔，再入段内无论高度还是速度都表现出递减的趋势，飞行包线呈收敛的走廊状。再入过程的初始段保持大迎角飞行姿态，且大气密度较小，空气动力作用不明显，但此时空天飞机具有很高的速度，因此首先受到热流约束的限制。再入飞行应最大限度地利用机体表面的防热材料，使空天飞机处于热流峰值的合理范围，并利用地球大气降低飞行速度，使其尽快通过高热流区域。再入过程的后半段处于稠密大气层中，空气动力的作用明显增强，动压和过载的共同约束逐渐加强。

2 变精度序列优化方法

“变精度”是指用粗离散(离散配点数较少)和细离散(离散配点数较多)2种方式对飞行时间进行离散。变精度序列优化方法的思路是：首先在粗离散情况下，应用混合优化方法进行轨迹优化，获得各时间点设计变量的初步优化结果；然后对初步优化结果应用插值处理方法，获得细离散情况下各设计变量的初始值；最后，在细离散情况下，应用经典优化算法，获得精度更高的优化结果。

智能优化算法采用一种最近新发展的优化方法——子集模拟优化算法[14]。有关研究表明：与GA相比，子集模拟优化算法具有更快的收敛速度[15]。经典优化算法采用SNOPT程序[16]，该程序采用了一种改进的SQP，适用于求解大规模非线性稀疏优化问题。

变精度序列优化方法的主要步骤如下：

(1)预优化

将飞行时间进行粗离散，应用混合优化算法求解轨迹优化问题。在粗离散过程中，离散时间段过多，则优化效率比较低；而离散点太少的话，得到的原始优化结果则不能足够准确地反应整个飞行过程，将这个原始优化结果作为后续SNOPT运算的初始值时SNOPT程序可能不会收敛。对于本文的问题，建议将飞行时间划分成20个时间微段。具体来说，首先采用智能优化算法迭代计算若干步(比如迭代100步、200步等)，计算结果不一定收敛，可以得到一个原始优化结果；然后以该原始优化结果作为SNOPT优化过程的初始值，进行优化计算直至收敛，获得初步优化结果。

(2)插值计算

根据轨迹设计精度的要求，进行细离散。对预优化的结果进行一元全区间插值处理，获得细离散情况下设计变量的值，然后将这个值设为精细优化的初始值。

(3)精细优化

在细离散情况下，应用SNOPT程序获得精度更高的计算结果。其中优化计算的初值为上一步插值计算中确定的初始值。变精度序列优化算法的流程如图1所示。

图1 变精度序列优化算法流程Fig.1 Variable-fidelity sequential optimization method flowchart

上述优化策略的优点是：既解决了空天飞机轨迹优化过程迭代初始值不易给定的问题，又可避免现有混合优化方法在离散配点数较大情况下不够稳健的缺点。

3 算例与验证

3.1 算例

空天飞机再入轨迹优化算例取自文献[17-18]，并在原始算例基础上增加了热流密度约束、动压约束和过载约束。

该飞行器的气动模型为

空天飞机参考面积和质量分别为S=250.237 m2、m=92 162 kg。空天飞机前缘半径为RN=1 m。不同高度下大气相关参数参考文献[19]。

该优化问题的初始状态为R=6 436 100 m，φ=0°,θ=0°,v=7 808 m/s,γ=-1°,ψ=90°。

末端状态为R=6 381 200 m,v=762.5 m/s,γ=-5°。

状态变量的约束为[17]

需要注意的是，这里的变量约束范围并未考虑工程实际情况，而是默认为在理想状态下，控制系统可以达到要求。

热流密度、动压、过载约束分别为

优化目标通过控制迎角α和倾侧角β使航程Ra最远。事实上，对于本算例来说，只需使末端纬度达到最大，就能使航程最远。故为方便起见，以末端纬度最大作为优化目标J：

J=maxθ(tf)

式中：θ(tf)为末端时刻tf的纬度。设计变量包括控制变量(迎角和倾侧角)和状态变量(参见1.3节)。

3.2 计算过程和结果分析

本算例中离散方法采用辛普森离散法。优化计算过程如下：

1)预优化。将时间离散成20段，采用子集模拟优化算法，迭代计算100步，得到原始优化结果(10s内就能计算完毕)。子集模拟优化算法的参数设定为：单层样本量取为100；条件概率取为0.51；最大层数取为100；终止精度取为10-6。

将子集模拟优化算法得到的原始优化结果作为SNOPT优化的初始点，进行优化计算，得到初步优化结果。

2)插值计算。将整个飞行时间离散为更多的时间微段。本算例考察了多种细离散情况，将整个飞行时间分别离散为30、50、60，80、100、150、200个时间微段。根据预优化后得到的结果，利用插值计算，获得各种细离散情况下设计变量的初始值。

3)精细优化。对于各种细离散情况，进行精细优化，优化结果见表1。计算机配置为3.10GHz,内存为2GB。从表1中可以看出，随着时间离散微段的增加，迭代步数总的趋势是逐渐增加的，精细优化时间也逐渐增加，得到的优化结果略有不同，目标值逐渐小幅上升。当时间离散为150个微段以上后，优化结果几乎不变，说明该离散段数足以满足精度要求。

表1 不同离散情况下精细优化结果

精细优化(细离散为200个时间微段)的优化结果轨迹如图2～图7所示。

图2 高度、航程时间关系Fig.2 Trajectory time histories for height and range

图3 速度、热流密度时间关系Fig.3 Trajectory time histories for velocity and heat flux

图4 经度、纬度时间关系Fig.4 Trajectory time histories for longitude and latitude

图5 航迹角、航迹方向角时间关系Fig.5 Trajectory time histories for flight path angle and azimuth

图6 迎角、倾侧角时间关系Fig.6 Trajectory time histories for angle of attack and bank angle

图7 动压、过载时间关系Fig.7 Trajectory time histories for dynamic pressure and overload

从图2～图5可以看出随着再入飞行过程的进行，高度、速度、航迹方向角均呈现下降趋势；航程、经度、纬度呈上升趋势；航迹角变化平稳，逐渐下降，且始终保持在小角度范围内，易于满足飞行过程的平衡要求。从图6可以看出下滑初期迎角较大，随着再入飞行过程的继续，迎角急剧减小，最后稳定在大约17°左右。这是由于在下滑初期，大气密度很低，操纵面气动效率较低，需要采用大迎角飞行。随着高度的降低，大气密度增大，当气动力能够克服重力后，轨迹进入一种近似滑翔的状态，迎角稳定在相对较小的角度；从倾侧角的变化曲线可以看出倾侧角变化相对平滑，易于控制，且满足控制量约束条件，下滑初始段，大倾侧角(绝对值最大达到78°左右)可防止再入时的轨迹跳跃，有利于尽快建立动压进入近似平衡滑翔状态，最后随着迎角降低，航向稳定性增大，横航向耦合下降，慢慢恢复至小倾侧角状态。从图3可知，在开始阶段，由于速度过大，热流密度急剧增加，一直增加到了约束最大值。虽然此后速度在减小，但是随着高度的降低，大气密度却在增加，所以导致热流密度一直保持在较大值。随着飞行过程的延续，在速度和大气密度的共同作用下，热流密度才逐渐减小。从图7中动压曲线可以看出，返回初期在很大的高度范围内动压较低，高度的降低导致大气密度呈指数增加，使动压呈指数增大，后来随着密度和速度的相互作用，动压呈现出一定的振荡特征，但始终在约束范围内；过载也呈现出递增趋势，由于在下滑初始阶段，迎角和速度比较大，大气密度也急剧增加，导致气动力急剧增大，过载急剧增加；后来在迎角、速度和大气密度的共同作用下，过载增加逐渐变缓，也呈现出一定的振荡现象，但仍在约束值3.0以内。

4 结束语

为了提高空天飞机再入轨迹优化的稳健性和精度，提出了一种变精度序列优化方法。变精度序列优化方法不仅可以用于空天飞机轨迹优化，其他大规模非线性优化问题也可以参考这种新的优化策略。最后用典型算例验证了这种方法的有效性，结果表明：

1)通过预优化，可为精细优化提供一个合理的初始值，解决了空天飞机优化过程迭代初始值难以给定的问题。对于不同配点数的细离散情况，精细优化计算过程均可收敛，获得优化结果，说明该方法具有很好的稳健性。

2)在精细优化中，当离散成200个时间微段时(离散后有1 609个设计变量，有1 803个约束方程)，仅仅迭代了223步就能收敛，优化计算效率高。

3)通过细离散和精细优化，可获得满足精度要求的优化结果。

References)

[1] 谷良贤, 赵吉松. 基于网格细化技术的地球-火星转移轨道优化[J]. 中国空间科学技术, 2013, 33(6): 17-25.

GU L X, ZHAO J S. Earth-to-Mars transfer orbit optimization based on mesh refinement[J]. Chinese Space Science and Technology, 2013, 33(6):17-25(in Chinese).

[2] 陈功, 傅瑜, 郭继峰.飞行器轨迹优化方法综述[J].飞行力学, 2011,29(4): 1-5.

CHEN G, FU Y, GUO J F. Survey of aircraft trajectory optimization methods[J]. Flight Dynamics, 2011, 29(4): 1-5(in Chinese).

[3] 胡海龙，南英，闻新.基于遗传算法的航天器再入动态终迹圈研究[J].中国空间科学技术, 2014, 34(1):35-41.

HU H L, NAN Y, WEN X. Research on dynamic reentry footprint for the spacecraft based on the genetic algorithm[J]. Chinese Space Science and Technology, 2014, 34(1):35-41(in Chinese).

[4] 闵学龙, 潘腾, 郭海林. 载人航天器深空飞行返回再入轨迹优化[J].中国空间科学技术, 2009, 29(4):8-12.

MIN X L, PAN T, GUO H L. Reentry trajectory optimization for manned deep space exploration[J]. Chinese Space Science and Technology, 2009, 29(4): 8-12(in Chinese).

[5] ZHANG D N, LIU Y. RLV reentry trajectory optimization through hybridization of an improved GA and a SQP algorithm[C]∥AIAA Guidance, Navigation, and Control Conference, 2011.

[6] VAVRINA M A，HOWELL K C. Global low-thrust trajectory optimization through hybridization of a genetic algorithm and a direct method[C]∥AIAA/AAS Astrodynamics Specialist Conference and Exhibit. Honolulu:[s.n.],2008: 1-27.

[7] BAGAI A, LEISHMAN J G. Adaptive grid sequencing and interpolation schemes for helicopter rotor wake analyses[J]. AIAA Journal, 1998,36(9): 1593-1602.

[8] RAHIMI A, KUMAR K D, ALIGHANBARI H. Particle swarm optimization applied to spacecraft reentry trajectory[J]. Journal of Guidance, Control, and Dynamics, 2013, 36(1): 307-310.

[9] ZHANG K N, CHEN W C. Reentry vehicle constrained trajectory optimization[C]∥The 17th AIAA International Space Planes and Hypersonic Systems and Technologies Conference, 2011.

[10] 吴德隆,王小军.航天器气动力辅助变轨动力学与最优控制[M]. 北京：中国宇航出版社,2006.

[11] 雍恩米, 唐国金, 陈磊.基于Gauss伪谱方法的高超声速飞行器再入轨迹快速优化[J].宇航学报, 2008, 29(6):1766-1772.

YONG E M, TANG G J, CHEN L. Rapid trajectory optimization for hypersonic reentry vehicle[J]. Journal of Astronautics, 2008, 29(6):1766-1772(in Chinese).

[12] 汤亮, 杨建民, 陈风雨. 多约束条件下的升力滑翔式再入轨迹优化[J]. 导弹与航天运载技术, 2013(1):1-5.

TANG L, YANG J M, CHEN F Y. Trajectory optimization with multi-constraints for a lifting reentry vehicle[J]. Missiles and Space Vehicles, 2013(1):1-5(in Chinese).

[13] 陈刚,万自明,徐敏,等.飞行器轨迹优化应用遗传算法的参数化与约束处理方法研究[J].系统仿真学报,2005,17(11):2737-2740.

CHEN G,WAN Z M, XU M, et al. Parameterization method and constraints transformation method for RLV reentry trajectory optimization using genetic algorithm[J]. Journal of System Simulation, 2005, 17(11):2737-2740(in Chinese).

[14] LI H S, AU S K. Design optimization using subset simulation algorithm[J]. Structural Safety, 2010, 32(6): 384-392.

[15] LI H S. Subset simulation for unconstrained global optimization[J]. Applied Mathematical Modelling, 2011, 35(10): 5108-5120.

[16] GILL P E, MURRAY W, SAUNDERS M A. User′s guide for SNOPT Version 7: Software for large-scale nonlinear programming[R].San Diego:University of California,2010.

[17] ZONDERVAN K P, BAUER T P, BETTS J T, et al. Solving the optimal control problem using a nonlinear programming technique Part 3: Optimal shuttle reentry trajectories: AIAA-1984-2039[R]. Reston: AIAA, 1984.

[18] BETTS J T. Practical methods for optimal control and estimation using nonlinear programming[M].Philadelphia: SIAM Press, 2010.

[19] 杨炳尉.标准大气参数的公式表示[J].宇航学报, 1983(1):83-86.

YANG B W. Formulization of standard atmospheric parameters[J]. Journal of Astronautics,1983(1):83-86(in Chinese).

(编辑：车晓玲、范真真)

A variable-fidelity sequential optimization method for reentry trajectory of space plane

ZHU Junjie,YU Xiongqing*,LUO Dongming, WANG Yu

KeyLaboratoryofFundamentalScienceforNationalDefense-AdvancedDesignTechnologyofFlightVehicle,NanjingUniversityofAeronauticsandAstronautics,Nanjing210016,China

A variable-fidelity sequential optimization method for reentry trajectory design of a space plane was proposed to increase robustness and accuracy of the optimization. The optimization process consisted of three steps. Firstly,there was a pre-optimization, in which for a coarse discrete case, the optimal solution of design variables for each discretized point was obtained by using a hybrid optimization method. Secondly,initial values of the discretized points in a refined discrete case were obtained by using linear interpolation from the pre-optimization solution. Finally, the trajectory optimization in the refined discrete case was conducted, and more accurate optimization results were obtained. A demonstration example was used to validate the proposed method. The results indicate that the reentry trajectory optimization problem can be solved robustly with satisfied accuracy and less computational burden.

space plane; reentry trajectory; optimization; hybrid optimization; interpolation

10.16708/j.cnki.1000-758X.2016.0032

2015-09-01；

2015-11-11；录用日期：2016-02-24；

时间：2016-04-29 10:49:39

http:∥www.cnki.net/kcms/detail/11.1859.V.20160429.1049.003.html

朱俊杰(1989-)，男，硕士研究生，1248400207@qq.com，主要研究方向为飞行器总体设计

*通讯作者：余雄庆(1965-)，男，教授，yxq@nuaa.edu.cn，主要研究方向为飞行器多学科优化设计和总体设计理论

朱俊杰,余雄庆,罗东明,等.空天飞机再入轨迹的变精度序列优化方法[J].中国空间科学技术,2016, 36(3):

50-56.ZHUJJ,YUXQ,LUODM,etal.Avariable-fidelitysequentialoptimizationmethodforreentrytrajectoryofspaceplane[J].ChineseSpaceScienceandTechnology, 2016, 36(3)：50-56(inChinese).

V221.6

A

http:∥zgkj.cast.cn