基于序列凸规划的运载火箭轨迹在线规划方法

程晓明，尚 腾，徐 帆，王晋麟，姚寅伟

(北京航天自动控制研究所，北京 100854)

0 引 言

随着世界各航天强国对外层空间开发与利用的不断深化，运载火箭需要在处理故障与自主飞行方面具有更强的自主控制能力[1]。因此，轨迹在线规划技术已成为运载火箭未来进入空间的一项关键技术，该技术赋予火箭自主规划与故障应变能力，能够极大地提升火箭执行任务的灵活性与完成新型任务的能力，也可节省飞行成本[2]。

由于大气层内火箭轨迹规划模型的复杂性，在面对突发任务、故障情况时，如何快速、可靠地求解出最优飞行轨迹，是目前火箭控制中的难点之一[3]。近几十年来，国内外学者在轨迹规划方法上的研究主要分为基于最优控制原理的间接法和基于数值优化的直接法。如文献[4-8]基于间接法的思想，研究了气动影响下的最优轨迹规划方法，通过求解包含横截条件的两点边值问题的思路，对最优轨迹进行直接求解。这种方法形式简单，但较为依赖初始协态猜想，收敛速度较慢，导致其实用性大打折扣。

随着计算机技术的发展及其运算能力的提升，基于数值方法的直接法得到较快的发展[9]。如文献[10-15]利用直接多步打靶法、伪谱法等对火箭的轨迹规划进行了大量的研究，验证了直接法应用于火箭轨迹规划上的有效性。这些方法的优化结果精度较高，但由于离散变量多，计算量大，导致算法的实时性不足。另外，这些算法在理论上缺少收敛性证明。

为了进一步保证数值算法的收敛性与实时性，将轨迹规划问题转化为一个凸规划问题，能利用凸规划求解算法(如内点法)所具有的多项式时间收敛特性，以同时满足算法收敛性与实时性要求。近年来，文献[16-17]利用凸优化理论研究了火星着陆在线轨迹规划问题，文献[18-19]针对大气层外入轨问题的轨迹规划问题，研究了基于凸优化的轨迹规划方法。然而，对于地球大气层内的运载火箭轨迹规划问题，其气动模型的复杂性使其不能简易的转化为凸规划问题而得到求解，因此，需要进一步研究大气层内火箭轨迹规划问题的快速求解方法。

序列凸规划[20]是一项求解非凸规划问题的局部优化方法，通过序列逼近的方式，将原非凸规划问题转化为序列凸规划问题，并进行序列求解。文献[21]通过序列线性化的方式，将序列凸优化应用于航天器的轨迹优化中。文献[22-24]基于序列凸规划的思想，研究了火箭的自主轨迹规划方法。但在目前的这种方式中，模型的凸化均是基于对非线性模型的线性化方式，这种方式会造成任务状态变化较大时轨迹在线规划模型的失准，且依赖初始猜想轨迹，对任务的适应性较弱。

本文基于序列凸规划思想，研究了一种模型补偿序列凸规划方法，该方法不需要初始猜想，其核心思想是序列地利用前一次规划得到的最优轨迹，补偿下一次迭代时轨迹规划问题中动力学的非线性项(轴向力加速度、法向力加速度和重力加速度)与过程约束，直到前后两次迭代的最优解收敛。通过在序列凸规划问题中增加柯西约束，保证得到的序列最优解为一组柯西序列，从而保证了该方法的收敛性。最后针对火箭上升段与着陆段的轨迹规划问题对该算法进行仿真校验。

1 火箭轨迹规划问题描述

1.1 火箭运动方程

在发射点惯性坐标系下，建立火箭的运动方程为：

(1)

式中：r=[x,y,z]T表示位置，v=[vx,vx,vx]T表示速度，m为飞行器质量，g=[gx,gx,gx]T表示重力加速度，T=[Tx,Tx,Tx]T表示飞行器推力。Isp为飞行器的比冲，g0为海平面的重力加速度大小。A和N分别为飞行器气动力中的轴向力与法向力。

1.2 约束与性能指标

1.2.1控制约束与状态约束

假设飞行器的比冲不随飞行状态改变，但推力可以通过节流阀门进行调节，根据节流阀的最小与最大节流程度，可得推力的约束范围为：

(2)

通常，飞行器的当前状态作为初始端点约束

s(t0)=s0∈R7

(3)

(4)

令飞行器的净重为mdry，初始质量m0,飞行器的燃料约束描述如下：

M={m(t)∈R+:mdry≤m(t)≤m0, ∀t∈[0,tf]}

(5)

1.2.2过程约束

火箭在大气层内飞行时，考虑三种常见的过程约束：动压约束、轴向过载约束和弯矩约束。

(6)

(7)

|qα|≤Qα max

(8)

式中：qmax，amax和Qα max分别为动压约束、轴向过载约束和弯矩约束的最大允许值。α为飞行器的攻角，由飞行器体轴与速度轴的夹角构成。

1.2.3性能指标

考虑最省燃料这一性能指标：

(9)

由于飞行器的质量是单调递减的，因此可将最省燃料的性能指标(9)转化为梅耶(Mayer)形式：

minJ=-m(tf)

(10)

因此，本文火箭轨迹规划问题的主要目标便是求解最优的推力矢量T(t)∈R3，使其满足各约束的同时，令性能指标(式(10))最小。

2 基于序列凸规划的轨迹规划方法

2.1 模型补偿序列凸规划方法

模型补偿序列凸规划方法的核心是如何序列地补偿大气层内轨迹规划问题中的非线性项，使其逼近真实的气动加速度模型与重力加速度模型，因此，为了在后续方法设计中对这些非线性项进行较为直观的处理，将火箭的动力学方程重新描述为：

(11)

并将控制量增广为：η=[υ,uT]T，控制约束可以重新写为：

(12)

0≤Tmine-Z≤υ≤Tmaxe-Z

(13)

U(t)={η(t)=[υ(t),uΤ(t)]Τ∈R4:

0≤Tmine-Z[k-1](t)≤υ(t)≤Tmaxe-Z[k-1](t)

(14)

由于轴向力加速度a、法向力加速度n、重力加速度g以及过程约束的非线性，导致该问题无法得到快速求解，因此本文研究了模型序列补偿的方式，以迭代的方式对这些非线性项进行逐步逼近，最终收敛到真实的模型，并提出了针对大气层内火箭轨迹规划问题的模型补偿序列凸规划方法。

2.1.1模型序列补偿方法

模型序列补偿方法的实现策略是，在第一次迭代中，将动力学中的非线性项建模为时间的线性函数，并将非凸最优控制问题转化为一个凸规划问题进行求解；然后在后续迭代中，利用前一次迭代产生的最优解对这些非线性项模型进行序列补偿，直到其收敛到真实模型。

为了便于描述，将动力学中的非线性项(包括重力加速度、气动力)以一个新的变量表示：

ψ(t)=g(t)+a(t)+n(t)

(15)

在序列迭代补偿方案设计中，在序列补偿的第一次轨迹优化时，根据轨迹规划任务的始末状态，将自动给出这些非线性项的初始线性猜想，其次基于该模型形成的凸规划问题，利用凸优化方法进行快速轨迹求解，然后将前一次得到的结果对下一次优化时动力学模型中的ψ(t)进行补偿。表示如下：

ψ(k+1)(t,s(k))=g(k)(t,s(k))+a(k)(t,s(k))+

n(k)(t,s(k))

(16)

通过上述补偿方式，可以将火箭轨迹规划问题转化为序列求解凸规划问题，直到模型ψ(t)收敛。

2.1.2过程约束的序列凸化

在模型补偿序列凸规划方法中，对过程约束的凸化采用了利用前一次迭代产生的最优解对后续迭代中的过程约束进行序列凸化的策略。主要包括：

1)动压约束

通过分析动压函数对高度和速度的海森阵可知，动压约束(6)不是一个凸约束。

为了满足本文所提出的算法在求解最优控制问题时对凸模型的需求，利用前一次迭代得到的最优解，对动压约束进行如下近似：

(17)

2)轴向加速度约束

根据前文所述的变量替换方法，轴向加速度约束可以基于前一次迭代得到的重力加速度，写为：

(18)

因此，轴向加速度约束可以简单的转化为凸约束。

3)弯矩约束

(19)

(20)

此时，通过利用前一次迭代的最优轨迹结果将过程约束转化为凸约束

2.1.3凸模型下的大气层内火箭轨迹规划问题

通过对控制量约束、动力学模型与过程约束进行序列凸化后，非凸的大气层内火箭轨迹规划问题可以建模为凸最优控制问题(21)。

η∈U(t),s∈Χ(t)，[s,η]∈P(t)，s(t0)=s0

U(t)={η(t)=[(υ(t),uΤ(t)]Τ∈R4:

0≤Tmine-Z[k](t)≤υ(t)≤Tmaxe-Z[k](t)

Χ(t)={s(t)∈R7:smin≤s(k)≤smax,

Z(t)≥ln(mdry)}

(21)

2.2 基于模型补偿序列凸规划对火箭轨迹规划算法描述

采用等间隔方法，对凸最优控制问题(21)进行离散。通过将时间区间进行离散、并将最优控制问题中的各种约束限制在每个离散点上的方式，将无穷维的连续线性凸最优控制问题转化为有限维的参数优化问题。由于该最优控制问题中的约束均为线性约束或二阶锥约束，因此离散后将得到一个二阶锥规划问题，从而可以使用有效的凸优化求解方法对其进行实时求解。

对于给定的时间区间[0,tf]，取等间隔时间区间为Δt，那么每个离散点对应的时刻为：

tk=nΔt,n=0,…,N

(22)

式中：N=tf/Δt为离散点的个数。控制变量采用零阶保持，即在时间区间t∈[tn,tn+1)内，控制量为ηn=η(tn)。

基于上述状态变量与控制变量的离散和已知的时间区间，式(21)所描述的连续凸最优控制问题可以转化为一个离散凸最优控制问题：

(23)

式中：上标(k+1)代表第(k+1)次迭代，即当前迭代，下标n代表第n个离散点。由式(23)确定的最优控制问题中的动力学方程可以看作是一个线性时不变方程。仅有ψ(k+1)会随着时间改变，但其各离散点的值均已知。性能指标函数为关于终端状态的凸函数，可以将其离散为关于最后一个离散点处状态的函数，该函数的凸性不变。因此，可以认为式(23)这一离散凸最优控制问题可以逼近原连续凸最优控制问题(21)。

式(23)中的离散状态空间矩阵为：

(24)

其中，η=[uΤ,υ]Τ;矩阵Ac,Bc和Cc分别为

(25)

(26)

(27)

在完成凸规划问题的离散建模后，该问题可利用原始-对偶内点法对其进行快速求解。因此，在进行火箭轨迹规划算法设计时，基于模型序列补偿方法，对离散凸规划问题进行重复求解，当max(s(k+1)-s(k))<ϑ时，即可认为算法收敛，其中，ϑ为设定的状态收敛阈值，此时算法直到收敛到一个定点最优解Y*=φ(Y*)。

定理1.当给定的初始离散设计变量在可行域内时，即Y(0)∈M⊆RNY，模型补偿序列凸规划算法在有限次迭代后能够收敛到定点最优解Y*∈M。

证.由于本文提出的模型补偿序列凸规划算法增加了柯西约束：

(28)

在迭代解的更新中，序列凸规划算法自映射迭代算子：φ:M⊆RNY→M⊆RNY表示离散可行设计变量Y(k)的更新，即输入由前一次迭代k得到的全体设计变量Y(k)，该算子能够在迭代k+1中得到状态量Y(k+1)：

Y(k+1)=φ(Y(k))

(29)

当增加了式(28)描述的柯西约束后，可知离散凸优化模型补偿迭代序列Y(k)=φ(Y(k-1))=φ(k)(Y(0)),k=1,2,…,M为一个柯西序列。当式(29)给出的柯西序列有界时，该柯西序列可收敛到定点最优解Y*=φ(Y*)∈M。

定理2.假设由式(23)表示的离散凸最优控制问题通过模型补偿序列凸规划算法，收敛到一个离散最优解Y*，那么由这个解逼近得到的连续解([r*,v*,m*],T*)是原非凸最优控制问题的一个局部最优解。

证.首先假设由式(23)表示的离散凸最优控制问题，在模型补偿序列凸规划算法下，收敛到一个离散的最优解Y*=(s*,η*)。根据本文对连续最优控制问题的离散化描述，离散方法能够逼近连续凸最优控制问题(21)，因此连续凸最优控制问题(21)的最优解(s*,η*)可以由Y*逼近得到，并在此假设(s*,η*)与Y*等价。

对于由式(23)描述的最优控制问题，矩阵Ac,Bc和Cc均为时不变矩阵，该问题的可控性克莱姆矩阵为：

(30)

式中：Φ(t,τ)为状态转移矩阵,Φ(t,τ)=eAc(t-τ)。

由于状态转移矩阵Φ(t,τ)在任意时间区间处均为满秩，即rank(eAc(t-τ))=7，且rank(Bc)=4，可知，null(eAc(t-τ)Bc)=0，因此可知可控性克莱姆矩阵是W(t)非奇异的，最优控制问题(23)中由[Ac,Bc]确定的连续线性动力学是可控的。

在可控性基础上，根据文献[16]可知，变量替换后控制量约束的松弛凸化是无损的，因此可知在每次迭代Ψ(k+1)都已知的情况下，式(21)表示的连续凸最优控制问题的最优解(s*,η*)，也是控制量约束“松弛”处理前的线性非凸最优控制问题的最优解(s*,u*)。

由离散最优解Y*是一个定点最优解，可知(s*,u*)也为定点最优解。这意味着动力学非线性项Ψ(k+1)通过序列凸化，最后能够收敛，因此(s*,u*)也是原始非凸最优控制问题的局部最优解。此时，通过变量反替换可知，与(s*,u*)对应的([r*,v*,m*],T*)是原非凸最优控制问题的一个局部最优解。

总结上述基于模型序列补偿凸规划的火箭轨迹规划算法，下面将算法流程描述如下。

第二步：基于线性化后的重力/气动加速度猜想值与过程约束凸化方法，利用无损凸化与变量替换技术，将轨迹规划问题转化为凸规划问题，并进行求解后储存为轨迹k，轨迹状态为s(k)。

第三步：计算轨迹k上状态s(k)对应的重力/气动力加速度值g(s(k)),a(s(k))和n(s(k))，将其补偿到最优控制问题(23)的建模中。

第四步：求解式(23)，得到轨迹k+1。

第五步：如果轨迹k+1与轨迹k中对应状态的差小于允许的阈值，则停止程序；否则，令k=k+1，返回第三步。

3 数值仿真与分析

为了检验本文所提出模型补偿序列凸规划方法的有效性，以火箭大气层内轨迹在线规划任务为背景，对该方法进行仿真校验。数值仿真平台为PC机，性能为Inter Core i7-6700 3.4 GHz，仿真软件采用了MATLAB2014a，其中的凸优化求解软件采用了CVX工具箱[25]。

本节主要通过大气层内轨迹规划的两个方面进行数值仿真实验设计：1)火箭上升段轨迹自主规划仿真；2)火箭着陆段轨迹自主规划仿真。

3.1 火箭上升段轨迹自主规划仿真与分析

火箭上升段轨迹自主规划仿真采用的火箭模型以及基本任务参数设置如表1所示。

表1 上升段火箭模型与任务参数设置

图1中，标记“·”的曲线代表最后一次迭代得到的最优推力大小曲线。图2中，标记“·”的曲线代表最后一次迭代得到的最优俯仰角和偏航角指令曲线。

图1 上升段推力大小变化曲线

图2 上升段姿态角变化曲线

经过仿真校验，基于该推力曲线，模型补偿序列凸规划方法得到的最优轨迹能够精确满足表1中给定的火箭上升段终端状态约束与控制量约束。

图3中，标记“·”的曲线代表最后一次迭代得到的最优解对应的气动加速度变化曲线。从图3可以看出，气动加速度曲线呈现递减的趋势，这主要是随着飞行高度的增高，大气密度快速减小导致的。从气动加速度的迭代过程图可以看出，模型补偿序列凸优化方法通过对气动力的不断补偿进行序列轨迹规划，对气动力补偿的收敛意味着该序列方法的收敛。

图3 上升段气动力加速度在发惯系三个方向分量的变化曲线

在轨迹规划收敛时间方面，该任务仿真时通过7次迭代后收敛，共耗时1.5 s，平均每次迭代耗时0.21 s；本次数值仿真对应50个离散点(以1 s为等时间间距进行离散)，当实际任务需要进一步减少轨迹生成时间时，可以在轨迹精度允许范围内，增大离散点间隔，从而通过减少离散点个数的方式加快轨迹生成速度。

3.2 火箭垂直着陆段轨迹自主规划仿真与分析

对于火箭垂直着陆任务，需要在火箭着陆时刻保持火箭速度为0，且箭体竖直向上。另外，通常要求火箭在着陆阶段燃料消耗最省，以节省火箭的发射成本、提高火箭的有效载荷重量。

仿真采用的火箭模型以及基本任务参数设置如表2所示。

表2 着陆段火箭模型与基本任务参数设置

图4和图5中，标记“·”的曲线代表最后一次迭代得到的最优飞行轨迹和最优速度在发惯系三个方向上的分量，没有标记“·”的曲线则代表迭代过程中得到的飞行轨迹和速度曲线。

图4 垂直着陆段位置在发惯系三个方向分量的变化曲线

图5 垂直着陆段速度在发惯系三个方向分量的变化曲线

每幅图有多条曲线，代表序列迭代过程中每次迭代产生的最优解曲线。虽然每次迭代产生的解均能满足着陆约束，但中间迭代采用的气动模型还处于不断补偿修正中，并不能代表真实的模型，只有最后一次模型补偿收敛后得到的解，才能算是具有真实意义的最优解。

图6和图7中，标记“·”的曲线分别代表最后一次迭代得到的最优推力大小曲线和最优俯仰角、偏航角指令曲线。仿真结果表明，基于该推力大小和姿态角形成的推力矢量，模型补偿序列凸规划方法得到的最优轨迹能够精确满足火箭垂直着陆所需要的终端状态约束条件。

图6 垂直着陆段推力大小变化曲线

图7 垂直着陆段姿态角变化曲线

图8中，标记“·”的曲线代表最后一次迭代得到的最优解对应的气动加速度变化曲线。从图8可以看出，火箭在着陆过程中，气动加速度也会随着速度的降低而减小，最后时刻减小为0。从气动加速度的迭代过程图可以看出，模型补偿序列凸规划方法对气动加速度进行了序列的补偿并最终收敛。

图8 垂直着陆气动力加速度在发惯系三个方向分量迭代更新图

最后，在轨迹规划收敛时间方面，该任务仿真时通过8次迭代后收敛，共耗时2 s，平均每次迭代耗时0.25 s；每次迭代所需时间与离散点个数、约束数量相关，本次数值仿真对应60个离散点(以0.5 s为等时间间距进行离散)，当实际任务需要进一步减少轨迹生成时间时，可以在轨迹精度允许范围内，通过增大离散点间隔，减少离散点个数的方式进行处理。

综合以上仿真结果可知，本文提出的模型补偿序列凸规划方法可以通过序列补偿的方式处理大气层内非线性气动力对轨迹规划带来的不利影响。在利用该方法对轨迹规划问题序列求解时，火箭的气动力通过不断的补偿最终收敛到准确的气动模型，同时也完成了方法的收敛迭代过程，这一快速收敛性将有利于满足火箭轨迹自主规划任务对轨迹优化实时性的需求。

4 结 论

针对大气层内的火箭轨迹自主规划问题，本文研究了一种模型补偿序列凸规划方法。该方法对火箭轨迹规划问题的非线性项(轴向力加速度、法向力加速度和重力加速度)和过程约束进行了序列补偿方案设计，结合无损凸化和变量替换技术，将火箭轨迹规划问题转化为一个序列凸规划问题。通过在序列凸规划问题中增加柯西约束，保证得到的序列最优解为一组柯西序列，从而进一步保证了该方法对应求解算法的收敛性。从理论分析可知，模型补偿序列凸规划方法不需要初始猜想轨迹，具有设计简单，自主性强的特点。最后，分别通过对大气层内上升段与垂直着陆段的轨迹自主规划任务的仿真，对本文提出的模型补偿序列凸规划方法进行了校验，得到结论如下：

1)该方法能够较好地处理大气层内气动力影响下的火箭轨迹规划问题。

2)该方法具有较好的收敛性，能够满足火箭轨迹自主规划任务对轨迹规划问题求解的收敛性要求与实时性要求。

3)该方法对火箭的着陆段与大气层内上升段的轨迹规划任务都有较好的适应性，对在线轨迹规划在工程上的应用提供了理论基础。