基于hp自适应伪谱法的固体运载火箭轨迹优化

洪 蓓 辛万青

北京宇航系统工程研究所，北京100076

多级固体运载火箭(Multistage Solid Launch Vehicle，MSLV)具有响应快速、机动性强、成本低、可靠性高的特点，目前已成为各国研制的热点。运载能力是运载火箭设计的重要指标，而主动段轨迹优化设计对提高MSLV 运载能力具有重要的意义。

MSLV 主动段轨迹优化实质上是一类终端状态固定、终端时刻自由，带有状态约束、控制约束的多阶段多约束非线性最优控制问题。最优控制问题的求解方法一般可分为间接法和直接法两种［1］。间接法根据极小值原理，将最优控制问题转化为两点边值问题，其解的精度高，但收敛半径小，且协态变量由于没有物理意义而很难准确估计。直接法不需对协态变量初值进行猜测，采用离散化方法将最优控制问题转化为非线性规划问题(NLP)，其解的收敛半径大，但解的最优性不能保证。

作为直接法的一种，全局伪谱法近年来得到了广泛的应用，其通过较少的节点便可获得较高精度的解。但在求解大型复杂问题及非光滑问题时，其运用效果并不理想。hp 自适应伪谱法融合了全局伪谱法与有限元法的优点，采用双层策略对最优控制问题进行求解，相比于全局伪谱法其占用更少的计算时间却可以获得更高精度的解［2］。

MSLV 采用三级固体助推+液体末修的推进方式，整个飞行过程分为5个阶段:一级、二级、三级动力飞行段、无动力滑行段和末级飞行段。对hp 自适应伪谱法在MSLV 主动段轨迹优化中的应用进行了研究，仿真结果表明各级之间过渡平滑，终端入轨条件和过程约束得到很好满足，hp 自适应伪谱法具有很高的计算精度和较快的收敛速度，是一种行之有效的新型优化算法。

1 主动段轨迹优化模型

1.1 动力学模型

在地心惯性坐标系OE－ XIYIZI下建立MSLV动力学模型，简单且利于计算。地心惯性坐标系定义［3］为:原点位于地心OE，OEXI轴于赤道平面内指向春分点;OEZI轴垂直于赤道平面，与地球自转角速度方向一致;OEYI轴由右手定则确定。MSLV 三自由度质点动力学方程如式(1)所示，该方程基于下列前提或假设:

1)地球是一个绕自身轴旋转的均匀圆球;

2)推力方向总是与飞行器体轴方向重合;

3)MSLV 在飞行过程中以零侧滑角飞行。

由上式可知，该动力学方程共包含7个状态变量，其中r=［rxryrz］T，v =［vxvyvz］T分别为位置和速度矢量，m 为运载火箭质量;此外，μ =GM为引力常数，T 为发动机推力，u=［uxuyuz］T为发动机推力方向，η 为发动机秒流量，D 为气动阻力，其计算公式如下:

式(2)中，Cd为阻力系数，Sref为运载火箭参考面积，ρ 为大气密度，其计算采用杨炳尉在“标准大气参数的公式表示”一文中给出的拟合公式［4］，vref为相对运动速度，即:

ωe为地球相对惯性空间的旋转角速度。

1.2 优化模型

(1)性能指标

在发动机性能参数一定的情况下，MSLV 主动段轨迹优化的目的是为了提高运载能力，由于MSLV 前三级均为耗尽关机，因此选取末级助推消耗推进剂质量最小为优化指标，其可以等价描述为末级助推结束后MSLV 剩余质量最大，即状态变量中质量的终端时刻值为最大，故指标泛函选择如下:

式中，tf为MSLV 主动段飞行结束的终端时刻。

(2)控制变量

选取发动机推力方向u =［uxuyuz］T及无动力滑行时间thx为待优化的控制变量。

(3)约束条件控制变量约束:

终端入轨约束:

高度、动压和过载等过程约束分别为:

状态变量连接点约束:

式(5)～(9)中，tmin，tmax分别为无动力滑行时间的上下限，Hf，Vf，Θf分别为入轨点高度、速度和当地弹道倾角，af，ef，if分别为轨道的半长轴、偏心率和轨道倾角，可根据终端位置rf和速度vf求得，Re为地球半径，qmax为动压上限值，nmax为法向过载上限值，ms为分离时抛掉的结构质量。

2 hp 自适应伪谱法

2.1 概述

近年来，在求解最优控制问题时，直接法得到了越来越广泛的应用，文献［5］从有限元法的思想出发，将直接法分为h 方法、p 方法和hp 方法3 种。

h 方法先将优化问题进行单元剖分，每个单元上使用低阶插值基函数进行离散，在优化过程中保证基函数阶次不变而通过逐步加密网格剖分对单元长度h 进行细分(即h－细化)。常用的有限元法一般采用h 方法。采用h 方法离散后的NLP 一般都是稀疏的，计算效率高，但其为多项式级的收敛速度，在求解大型复杂问题时将带来维数灾难问题。

p 方法正好相反，其将优化问题作为一个单元(或少量单元)处理，每个单元上采用高阶插值基函数进行离散，优化过程中保持单元剖分不变，通过增大基函数的阶次p 以提高精度(即p－细化)。全局伪谱法就是一种典型的p 方法，其具有简单的结构、较高的精度及指数性的收敛速度，但对于非光滑问题，收敛速度非常慢，即使采用高阶的插值基函数也可能无法得到满足精度的解。

hp 方法允许单元长度h 和基函数阶次p 在优化过程中同时进行自适应调节。作为hp 方法的一种，hp 自适应伪谱法将有限元法与全局伪谱法的优点进行融合，运用双层策略决定单元数和基函数的阶次以满足精度要求。该方法包含两部分内容:hp－LG 离散方法和hp－网格细化方法。

2.2 hp－LG 离散方法

(1)时域变换

假设整个最优控制问题在t∈［t0，tf］上被分成K个单元，即t0＜t1＜… ＜tK=tf。对于任意单元k，通过下式将时间区间由t∈［tk－1，tk］转换到τ∈［－1，1］。

转换后的性能指标为:

系统微分方程变为:

边界条件:

不等式约束:

此外，还存在内点约束为:

(2)状态变量与控制变量的近似

设x(k)(τ)和u(k)(τ)分别表示单元k 上的状态变量和控制变量。在单元k 上以Nk个LG 点及τ0= －1为配点，构造Nk+1 阶Lagrange 插值多项式并以此为基函数近似状态变量为:

其中

类似的，控制变量也可以近似如下:

(3)动力学微分方程约束转换为代数约束

定义微分近似矩阵D(k)为:

系统微分方程动力学约束可进一步转换成一系列代数形式:

(4)离散条件下的终端状态约束

由于在状态变量的近似中没有考虑终端时刻τf=1，因此将终端状态约束离散并近似为:

(5)离散条件下的性能指标

将Bolza 性能指标中的Lagrange 积分项用Gauss 积分近似，即

这样就将一个最优控制问题转化为参数优化问题，选择合适的非线性规划算法即可进行求解。

2.3 hp－网格细化方法

hp－网格细化方法的核心在于如何判定单元k上到底是进行h－细化还是p －细化。对单元k 上的约束满足程度进行考核，在单元k 上任意选取L个节点，且满足:

其中，i=1，…，n;j=1，…，s;l=1，…，L。

取e(k)的算术平均值为:

定义

定义误差指标ρ，下面分3 种情况进行讨论:

3 仿真计算及结果分析

3.1 仿真条件

目标轨道选择为700km 高度的太阳同步圆轨道，终端入轨约束条件如表1 所示。

表1 终端入轨约束

3.2 仿真结果分析

(1)仿真精度分析

hp 自适应伪谱法得到的终端入轨参数与目标值相比，偏差很小，详见表2 所示，说明该算法具有非常高的计算精度。从图1 和图2 中可以看出，各阶段之间衔接过渡平滑，高度和速度的终端约束得到了很好的满足，三级与末级之间加入无动力滑行段，起到了提升高度和降低速度的作用。在无动力滑行段，由于发动机推力为零，推力方向不再是控制量，在满足控制约束条件下可以任意取值，为方便起见，在优化过程中将uz置为1，ux，uy置为0，因此控制变量在无动力滑行段出现了跳变。

表2 优化轨道根数与目标轨道根数的偏差

(2)配点情况分析

Gauss 伪谱法的配点为LG 点，该配点呈现出两边密中间疏的特点［6］，而hp 自适应伪谱法在优化过程中能够对单元数和基函数的阶次进行自动调节，配点分布与Gauss 伪谱法有很大的不同。在无动力滑行段，MSLV 为纯惯性飞行，运动相对简单，配点数也较少，在其余的4个飞行阶段，运动相对复杂，因此配点数也较多。从配点分布情况来看，配点不再呈现出两边密中间疏的特点，而是在梯度变化小的地方配点数较少，变化大的地方配点数相对较多，说明该方法具有良好的自适应调节配点数的能力。

(3)最优性分析

图4 ～图5 为状态变量对应的协态变量曲线，hp 自适应伪谱法可以对协态变量特别是协态初值进行准确的估计。将协态变量初值与状态变量初值一起代入状态变量与协态变量组成的微分方程组，积分求解即可得到主动段最优轨迹，将该轨迹与优化得到的轨迹进行比较可知两者基本一致。图6 为Hamilton 函数变化曲线，在每一个阶段，Hamilton 函数都保持常值，这也从反面验证了hp自适应伪谱法与最优控制理论中的一阶必要条件是一致的。

图1 高度随时间变化曲线

图2 速度随时间变化曲线

图3 控制变量随时间变化曲线

图4 协态变量λx，λy，λz随时间变化曲线

4 结论

本文以MSLV 运载能力最大为优化目标，建立了基于hp 自适应伪谱法的固体运载火箭主动段轨迹优化模型，并进行了仿真计算。仿真结果表明:

1)优化后的轨迹各级之间连接过渡平滑，严格遵循过程约束条件，且终端入轨约束满足精度高;

图5 协态变量λvx，λvy，λvz随时间变化曲线

图6 Hamilton函数随时间变化曲线

2)hp 自适应伪谱法的配点分布能根据实际情况进行自适应调节，相比全局伪谱法两端密中间疏的情况，该分布更为合理;

3)hp 自适应伪谱法对协态变量的初值能够进行准确的估计，可将该算法用于伪谱最优反馈控制，实现对轨迹的实时跟踪;

4)采用hp 自适应伪谱法可降低对变量初始值的敏感性，具有快速、全局的优点，具有良好的通用性和可扩展性。

［1］Betts J T. Survey of Numerical Methods for Trajectory Optimization［J］. Journal of Guidance，Control and Dynamics，1998，21(2):193-207.

［2］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.

［3］龙乐豪，等.总体设计(上)［M］. 北京:宇航出版社，1989.

［4］杨炳尉. 标准大气参数的公式表示［J］. 宇航学报，1983，(1):83-86.(Yang bingwei.Formulization of Standard Atmospheric Parameters［J］. Journal of Astronautics，1983，(1):83-86.)

［5］Darby C L，Hager W W，Rao A V.Direct Trajectory Optimization Using a Variable Low-Order Adaptive Pseudospectral Method［J］. Journal of Spacecraft and Rockets，2011，48(3):433-445.

［6］David Benson. A Gauss Pseudospectral Transcription for Optimal Control［D］. Cambridge，Massachusetts:Massachusetts Institute of Technology，2005.