基于OPTIMUS的临近空间飞行器再入轨迹优化设计

张蕊 孙国庆 杨泽山

（北京航空航天大学宇航学院，北京 100191）

1 引言

临近空间飞行器再入时，由于需要较长一段时间在大气层中飞行，将受到热流密度、动压、过载的严格约束，以及需要进行主动姿态控制，给飞行器材料和内部设施热绝缘带来巨大的挑战，所以安全可靠的再入成为重要问题，再入轨迹的优化也成为一项关键技术。

现有的轨迹设计方法有工程规划法，依据工程经验将轨迹分为数段，然后分段求解［1］。轨迹优化设计主要有基于极大值原理的间接法［2-3］和基于非线性规划理论的直接法，直接法中的伪谱法［4-5］应用最为广泛，但无论是直接法还是间接法，均采用传统的寻优方法，根据工程经验设置参数初值，然后通过试凑法对这些参数进行调整。由于目标函数对待优化参数的微小变化十分敏感，这样在全局寻优时给轨迹优化带来很大问题，很难对整个的解空间进行搜索并找到全局最优解，即使能找到全局最优解，也需花费大量的时间寻优，而且优化结果的好坏，在很大程度上取决于某些参数的初始给定值，因而不便于工程应用。除了传统的直接法和间接法之外，还有许多其他优化方法应用也较广泛，如动态规划方法［6］、快速搜索 随机树 法［7］、遗传 算法［8］等，这些算法在OPTIMUS软件平台中都有体现。

OPTIMUS软件平台是一种过程集成与多学科优化的软件平台，主要功能包括过程自动化、设计空间开发、先进的概率和统计方法、数值优化算法等。OPTIMUS已经被广泛应用于航空航天领域的结构优化、涡轮气动优化设计［9］、弹道仿真控制系统参数优化［10］等方面。在再入轨迹优化问题中，基于OPTIMUS平台，可以充分应用其集成的数值优化算法，快速确定初始猜想值，并选取最佳优化方案进行求解。

本文针对再入轨迹优化中初值选取和计算效率的问题，利用OPTIMUS 寻优平台，找出轨迹优化的初值点，并确定最佳的优化算法，再依据算法编写优化程序，这种计算机集成的辅助优化设计可以取代大量的低效手动设计。仿真结果表明，基于OPTIMUS的轨迹优化方案计算快速、寻优准确，且计算结果精度高。

2 再入轨迹优化问题

飞行器再入过程中，由于各种状态变量的绝对值相差较大，这对于再入运动方程的数值计算非常不利，为提高计算精度、增强算法收敛性与快速性，需将实际再入运动方程进行归一化处理。因此，在考虑地球旋转的情况下，飞行器再入无量纲运动方程组为

上述再入运动方程组是关于无量纲时间τ＝t／的微分方程组，其中，R0为地球半径，g0是海平面重力加速度。式（1）～（6）方程组中：r为飞行器质心-地心的距离与地球半径R0的比值；θ和φ分别为再入过程中的经度纬度；V为飞行器相对地球的速度与归一化变量（圆周速度）的比值；γ为航迹倾角；ψ为航向角；Ω为地球自转角速度，σ为滚转角；L为飞行器的总升力加速度；D为阻力加速度。L、D计算公式为

式中：ρ为大气密度；m为飞行器质量；Sref为飞行器的参考面积；CL和CD分别为飞行器的升力系数和阻力系数，其与飞行器的马赫数Ma和攻角α有关，也即CL＝CL（Ma，α），CD＝CD（Ma，α）。

再入轨迹优化问题，即寻找最优控制律σ（t），使飞行器在再入过程中，满足要求的初始状态、目标状态和运动学方程，并使给定的性能泛函具有极值。本文中设计的性能泛函为在给定纵程下使再入横程最大。

为简化问题，再入轨迹优化时采用以下假设条件：重力加速度g为常值g0；满足准平衡滑翔条件γ≈0，γ≈0；由于高度相对地球半径很小，令r为常值且r＝1；不考虑地球自转，将哥氏加速度与牵连加速度影响因素忽略。简化后的再入运动方程为

过程约束考虑气动加热、过载及动压3种不等式约束。热流率约束为

动压约束q主要取决于机体结构和运载器表面防热材料的强度，qmax为允许的最大动压值。

过载约束n主要取决于飞行器结构强度和器载设备的过载范围，nmax为允许的最大过载值。

给定纵程下最大横程问题的描述如下（见图1）：设飞行器的初始再入点为O点，其对应的经纬度为θO和φO；要求到达的目标点为T，其对应的经纬度为θT和φT；P点为经过O点和T点的大圆弧上的点，且满足O点到P点的距离为给定的纵程值SP，其对应的经纬度为θP和φP；F点即为给定纵程值SP所对应的最大横程点，其对应的经纬度为θF和φF。再入轨迹优化问题即为寻求最优滚转角控制量σ（t），使飞行器在满足初始状态（O点）和给定的纵程SP的情况下，到达最大横程点F处，且满足再入过程约束和再入运动方程。

图1 最大横程示意图Fig．1 Maximum crossrange diagram

本文主要是基于最优控制理论，设计滚转角控制量，通过选择合适的求解方法求解最优控制问题，使再入飞行器在给定纵程SP下，获得最大横程SF。优化目标是横程SF最大，优化变量为滚转角σ（｜σ｜＜π／2），性能泛函即为

终点约束条件主要有2个：①达到预设的末端能量；②保证连接∠OPF为直角。

式中：rF为终端高度；VF为终端速度；eF为终端能量，哈密尔顿函数H写为

式中：pθ，pφ，pV，pψ为共轭状态变量，共轭状态方程为

根据最优条件，最优滚转角控制律要满足

根据运动的积分以及整理，具体参见文献［11］，可得到滚转角σ最优控制律为

在滚转角表达式（22）中，可以看出分子分母是常数c1，c2和c3的线性表达式，由此可以得出，滚转角的值只与其中2个独立变量有关，最终优化目标也即通过搜索找到2个独立变量的值。假设：c2＝sinλ，c3＝cosλ。终点约束条件2可以改写为

至此，给定纵程下最大横程优化问题可表述为：在满足终点约束条件的前提下，利用寻优算法求解最优控制参数c1和λ的值。c1作用域为［-1，1］，λ作用域为［-π，π］。对实际飞行器轨迹优化而言，可以取

其中K1，K2为权值，可取为K1＝1，K2＝1。需要寻找两个参数c1和λ，使Δ（c1，λ）最小。利用OPTIMUS集成的优化平台，可以方便快速地搜索出最优参数c1和λ。

3 基于OPTIMUS软件平台的再入轨迹优化

在OPTIMUS软件平台中，集成了多种优化算法，分为两大类：全局优化算法和局部优化算法。全局优化算法主要包括：差分进化、自适应进化、模拟退火、遗传算法［12］等；局部优化算法主要包括：模式搜索法［13］、序列二次规划、非线性二次规划、广义简约梯度等。全局优化算法往往收敛速度较慢，而局部优化算法收敛速度较快，但如果初始点选取不够准确，局部优化算法将有可能使搜索进入局部极值，导致结果不准确。本文结合全局优化算法和局部优化算法各自的优缺点，提出了基于全局优化算法和局部优化算法相结合的混合优化的概念，全局优化算法采用遗传算法，局部优化算法采用模式搜索法，将遗传算法得到的结果作为模式搜索法的初始点，并进一步精确搜索得到最优解。采用混合优化算法，不仅缩短了仿真时间、降低了错误率，而且提高了优化精度，满足了全局最优性。下面主要介绍在OPTIMUS软件平台下优化设计的详细步骤。

3．1 建立工作流

为了应用OPTIMUS软件平台，通过可视化的建模方式，将仿真程序和工作流集成在OPTIMUS平台的框架中，首先建立M 文件global．variable．dat，其中包含有待优化的控制参数c1和λ的值。接着将事先编好的数学模型仿真程序导入OPTIMUS平台，搭建系统工作流。在OPTIMUS 中选择控件，并将所有控件用连接线进行连接，搭建OPTIMUS工作流。工作流控件中包含设计变量、输入输出文本文件、仿真程序模块、输出向量、输出变量以及连接线。搭建工作流的目的，在于以图像化的方式描述仿真中计算的逻辑和数据流程，仿真工作流在输入文件中更新输入参数、进行仿真计算，从结果文件中读入计算结果。

在手动输入完成后，OPTIMUS 就开始自动化的计算。先对c1和λ取值，修改global．variable．dat中相应的数值。之后运行M 文件搭建的动力学仿真模块，待程序运行结束后读取Output．dat中的数据，根据事先设置的算法，自动计算出c1和λ更新值，回到第一步进行数据刷新，然后重复这个工作流程，直至得到符合条件的解。整个工作不需要人工干预，由计算机独立完成长时间的运算，提高了工作效率。

3．2 试验设计与响应面模型的搭建

完成建模后要进行试验设计，通过较少的试验次数，比较准确地反映输入参数和输出参数的关系，并且通过响应面模型使结果可视化。OPTIMUS提供了多种试验设计（Design of Experiment，DOE）方法，包括全参数、随机试验设计、部分参数等。试验设计的目的在于分析参数相关性，为响应面的建立提供样本点，保证在设计变量变化时（特别是边界上）不发生干涉、参数不匹配及命令不执行等情况。此外还应对网格尺度、计算步长、差分格式及数等进行匹配调试，以便在计算时间、精度及稳定性等方面得到平衡。

优化变量c1可行域为［-1，1］，λ可行域为［-π，π］，图2为拉丁超立方样本点分布，图3为参数相关性。

图2 拉丁超立方样本点分布Fig．2 Latin super cubic sample point distribution

图3 拉丁超立方参数相关性Fig．3 Latin super cubic parameters correlation

通过实验设计分析得出：该飞行器轨迹优化的目标泛函对选择的参数敏感性非常高，参数响应面存在很多局部的急剧变化，因此选择一种合适的优化算法显得尤其重要。试验点的选择应针对不同模型有侧重地进行考虑，虽然尽可能地将试验点个数减少，但前提是体现整个参数可行域内的分布情况。

选择相应的设计实验方法可以得到响应面模型，来模拟一系列输入与输出参数之间的响应关系。由于该飞行器优化模型非线性极强，对于参数的响应面波动十分剧烈，所以需要相当精确的拟合方法对其进行拟合，以减少后续分析处理中由于模型误差引入的不确定性。以插值法中的RBF 神经网络方法为例，形成的插值响应面如图4所示。

图4 RBF神经网络方法插值响应面Fig．4 RBF neural network method interpolation response surface

响应面模型搭建的出发点在于利用一系列的拟合算法，在若干试验点的基础上，对待优化的精确响应面模型进行离散样本点拟合，用数理分析统计的方法代替实际飞行器模型。虽然舍弃了小部分精度，但换来的是后续的高效运算，极大地减少了数据存储和运算量，基于该拟合模型，就可以将输入参数和输出参数关联起来，从而具有较快的响应速度和较好的响应精度，降低了运算成本和响应时间。

完成了上述的工作流的建立、以及试验设计和响应面模型的搭建之后，就可以选择OPTIMUS软件平台内部自带的遗传算法和模式搜索法对参数c1和λ进行寻优计算，并最终得到参数的最优值为c1＝0．784 7，λ＝0．524 6。

4 编程实现轨迹优化与仿真结果

将OPTIMUS软件平台优化出的参数c1和λ的值，代入到公式（22）的滚转角最优控制律中，并编写相应的MATLAB 再入轨迹仿真程序，即得到满足性能泛函和过程约束的再入轨迹。

MATLAB 仿真用到的再入飞行器参数为：参考面积149．388m2，飞行器质量37 362．9kg，最大热流率794 005 W／m2，最大动压15 425N／m2，最大过载5gn。攻角α变化范围为［0°，45°］；马赫数Ma变化范围为［3，25］；高度变化范围为［0km，120km］。初始条件为：高度123．104km，经度和纬度分别为-122．498°和-33．262 9°，速度为7 625．15m／s，航迹倾角为-1．249 4°，航迹偏角为46．061 7°。终端条件为：高度30．427km，经度和纬度分别-62．346 5°，35．270 6°，速度为908．151 6m／s，纵程为9 871．4km。图5～8是仿真结果，对任意的给定纵程，都有左右两边的两个最大横程。图5是高度-速度空间中的再入轨迹，从图中可以看到满足热流率、动压和过载约束；图6是对应的最优滚转角；图7是航迹倾角；图8是星下点轨迹（给定纵程下对应的最大横程有左右2条）。

图5 再入轨迹图Fig．5 Entry trajectory

图6 滚转角控制曲线Fig．6 Bank angle profile

图7 航迹倾角变化曲线Fig．7 Flight path angle profile

图8 星下点轨迹Fig．8 Trajectory of ground track

5 结束语

本文基于OPTIMUS优化软件平台，提出了临近空间飞行器的再入轨迹优化方案。与传统的优化方法不同，先利用OPTIMUS平台优化求出初值，并选出最佳优化算法，再编写优化程序进行求解。OPTIMUS平台在轨迹优化中的应用，可以快速地搭建起整个数据流程，通过应用多种不同的优化算法，对同一个问题进行求解分析，可以更加全面地了解轨迹优化问题各个环节的细节，并获得优化问题的初始猜想，然后针对具体问题设计优化方案，所开发的算法程序具有计算高效、便捷准确的特点。

（References）

［1］赵汉元，飞行器再入动力学和制导［M］．长沙：国防科学技术大学出版社，1997

Zhao Hanyuan．Flight vehicle entry dynamics and guidance［M］．Changsha：University of National Defense Science and Technology Press，1997（in Chinese）

［2］Vinh N X．Optimal trajectories in atmospheric flight［M］．NewYork：Elsevier Scientific Publishing Company，1981

［3］Lu P，Xue S B．Rapid generation of accurate entry landing footprints［J］．Journal of Guidance，Control and Dynamics，2010，33（3）：756-767

［4］雍恩米，唐国金，陈磊．基于Gauss伪谱方法的高超声速飞行器再入轨迹快速优化［J］．宇航学报，2008，29（6）：1766-1772

Yong Eenmi，Tang Guojin，Chen Lei．Rapid trajectory optimization for hypersonic reentry vehicle via Guass pseudospectral method［J］．Journal of Astronautics，2008，29（6）：1766-1772（in Chinese）

［5］宗群，田栢苓，窦立谦．基于Gauss伪谱法的临近空间飞行器上升段轨迹优化［J］．宇航学报，2010，31（7）：1775-1781

Zong Qun，Tian Bailing，Dou Liqian．Ascent phase trajectory optimization for near space vehicle based on Gauss pseudospectral method［J］．Journal of Astronautics，2010，31（7）：1775-1781（in Chinese）

［6］John T B．Survey of numerical methods for trajectory optimization［J］．Journal of Guidance，Control and Dynamics，1998，21（3）：193-207

［7］周克强，高晓光，白奕．改进的快速扩展随机树在航迹规划中的应用［J］．系统工程与电子技术，2006，28（10）：1538-1540

Zhou Keqiang，Gao Xiaoguang，Bai Yi．Trajectory planning using improved rapidly exploring random tree［J］．System Engineering and Electronics，2006，28（10）：1538-1540（in Chinese）

［8］周浩，陈万春．高超声速飞行器滑行航迹优化［J］．北京航空航天大学学报，2006，32（5）：513-517

Zhou Hao，Chen Wanchun．Optimization of glide trajectory for a hypersonic vehicle［J］．Journal of Beijing University of Aeronautics and Astronautics，2006，32（5）：513-517（in Chinese）

［9］严俊峰，吴宝元，逯婉若．基于OPTIMUS的涡轮气动优化设计［J］．火箭推进，2008，34（2）：13-17

Yan Junfeng，Wu Baoyuan，Lu Wanruo．Turbo aerodynamics optimization design based on OPTIMUS flat［J］．Journal of Rocket Propulsion，2008，34（2）：13-17（in Chinese）

［10］周健，姜晓菡，周喻虹．基于OPTIMUS的控制系统参数优化计算［J］．弹箭与制导学报，2008，28（3）：13-18

Zhou Jian，Jiang Xiaohan，Zhou Yuhong．Optimization of parameters for control system based on OPTIMUS［J］．Journal of Projectiles，Rockets，Missiles and Guidance，2008，28（3）：13-18（in Chinese）

［11］Lu P．Entry trajectory optimization with analytical feedback bank angle law，AIAA-2008-7268［R］．Washington：AIAA，2008

［12］袁麟博，章卫国，李广文．一种基于遗传算法-模式搜索法的无人机路径规划［J］．弹箭与制导学报，2009，29（3）：279-282

Yuan Linbo，Zhang Weiguo，Li Guangwen．A path planning for UAV based on genetic-pattern searching algorithm［J］．Journal of Projectiles，Rockets，Missiles and Guidance，2009，29（3）：279-282（in Chinese）

［13］张煜东，吴乐南，王水花．基于遗传算法与模式搜索的混合优化算法［J］．南京信息工程大学学报，2012，4（1）：34-39

Zhang Yudong，Wu Lenan，Wang Shuihua．A hybrid optimization method based on genetic algorithm and pattern search［J］．Journal of Nanjing University of Information Science ＆ Technology，2012，4（1）：34-39（in Chinese）