可回收火箭是指能够在执行完一次发射任务后，安全、稳定地重返地球，且具备可维护性、可重复利用性的新型航天任务发射运载器。近年来，商业火箭发射公司如Blue Origin、SpaceX等，相继成功研制了可垂直回收火箭，掀起了世界各航天大国可回收商业运载火箭的研制热潮。SpaceX研制的具备可重复利用的Falcon 9系列运载火箭，自2015年首次成功回收后的52次发射中，完成了44次一级火箭的陆地和海上回收^[1]。

火箭垂直回收的技术，最早可以追溯到阿波罗时代的月球软着陆研究，并在20世纪末和21世纪初的火星精确软着陆任务中得到进一步的发展。由于月球和火星表面地形复杂，且受到地面通信延迟等因素的影响，实现软着陆要求探测器具备一定的在线制导和障碍躲避等能力^[2]。早期的探测器月球软着陆制导算法，对着陆精度要求不高，主要目的是保证探测器尽可能接近零速度地降落在月球表面，例如月球垂线法^[3]、重力转弯制导等^[4]。为提高探测器的着陆精度，不断有新的制导方法被提出，例如基于标称轨迹的多项式制导方法等，一定程度上提高了着陆的精度和算法的实时性，但算法本身作了很多的假设，不具备应对多种未知不确定性误差的高鲁棒性^[2]。随着人类对太空探索的不断深入，火星探测等计划的提出和实施，使得精确软着陆技术逐渐成熟，制导算法的实时性和鲁棒性进一步得到提高。美国宇航局(National Aeronautics and Space Administration，NASA)的“好奇号”探测器，首次采用了大气制导技术，并结合标称轨迹跟踪，实现了2 km级别的着陆误差^[5-7]。为保证着陆器的安全着陆和落点精度，通常需要根据着陆器所携带的相机、雷达等，实时地根据地表特征进行障碍躲避、位置识别等，实现探测器的自主导航^[8-9]。NASA目前公布的下一代火星着陆任务中，提出了结合相对地形导航的实时在线制导技术，探测器能够在下降过程中对地形进行观测，并根据地形信息实时生成制导策略，从而尽可能地精确到达目标区域(百米量级)^[10-11]。相对于火星探测器软着陆，由于地球稠密大气等因素的存在，火箭垂直回收过程的动力学环境更为复杂，且火箭下降时间较短，制导方法面临着更为严峻的考验。传统的制导方法为了获得实时制导律，作了很多简化和假设，通常难以保证复杂动力学环境下的着陆精度，而基于更加精确模型下的数值和近似求解方法，则无法满足火箭着陆过程对制导方法实时性的要求。

凸优化最早是应用于求解火星探测器的有动力精确软着陆问题^[12-14]。文[15-16]提出了将火箭着陆的轨迹优化问题，转化为二阶锥问题(second order cone problem，SOCP)进行求解，进一步又提出了一种无损凸化的技术，推广至求解一般性轨迹优化问题。相对于传统求解轨迹优化问题的数值迭代算法，凸优化能够在保证解的等价性和最优性的前提下，使问题在多项式时间内得到求解，从而使得在线求解和实时制导成为可能^[17]。针对火箭垂直回收问题，国内Liu^[18]和张志国等^[12]分别建立了凸优化求解二维、三维问题的模型。NASA的喷气推进实验室(Jet Propulsion Laboratory，JPL)和欧空局(European Space Agency，ESA)等机构，也都开展了凸优化用于垂直起降火箭制导算法的相关研究，将问题的维度进一步拓展至六自由度^[19-20]。利用凸优化求解火箭着陆问题时，重点解决了制导算法的问题，但通常对模型做很多假设，例如忽略大气阻力作用、近似重力加速度为常值、要求推力连续可调等，使得制导算法能够得到快速收敛。然而在实际飞行过程中，该制导算法的鲁棒性及适应边界，尚没有明确界定。

本文针对火箭垂直回收着陆问题，从凸优化求解算法出发，建立凸优化求解火箭着陆制导问题模型，并利用内点法(interior point method, IPM)定制求解器，以保证算法求解的实时性。为测试算法在复杂动力学环境下的鲁棒性，搭建了制导、导航与控制(guidance, navigation, and control, GNC)一体化闭环仿真平台。测试了在大气随机扰动、发动机节流特性、导航系统随机偏差以及系统整体时延等多重不确定性因素作用下算法的鲁棒性，提出了发动机推力偏差和节流性能受限情况下的推力执行策略，以及处理时延的状态估计制导策略。

1 火箭着陆段的动力学模型

火箭的动力着陆段，是指火箭在垂直回收过程中，完成主动减速、调姿等操作后，开启发动机完成精确软着陆的最后阶段。本节将分析火箭着陆段的动力学模型及受到的约束，并在此基础上建立燃料最优的最优控制问题模型。

1.1 火箭动力学及约束

本文主要考虑火箭中心质点的运动问题。建立以目标落点为原点的站心直角坐标系，y轴与重力方向相反，垂直当地水平面向上，x轴和z轴与当地水平面平行且分别指向正北、东方向，三轴呈右手坐标系，如图 1所示。由于火箭动力着陆段的起始高度通常为百米到公里量级，因此重力的作用在本阶段简化为常值处理，且不考虑由地球自转引起的非惯性力。该阶段的动力学方程可以表示为

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\dot r}} = \mathit{\boldsymbol{v}},}\\ {\mathit{\boldsymbol{\dot v}} = \mathit{\boldsymbol{g}} + \frac{\mathit{\boldsymbol{T}}}{m} + {\mathit{\boldsymbol{a}}_{\rm{d}}},}\\ {\dot m = - \frac{{\left\| \mathit{\boldsymbol{T}} \right\|}}{{{I_{{\rm{sp}}}}{g_0}}}.} \end{array}} \right. $

(1)

其中：r和v分别表示火箭运动的位置和速度；g表示重力加速度；T表示火箭发动机推力；m表示火箭的质量；a_d为气动加速度；I_sp表示火箭发动机比冲；g₀取9.801 m/s²，表示地球海平面重力加速度常数。

火箭着陆段运动轨迹受到地面避障约束角β的限制，如图 1所示。假设火箭发动机推力方向始终沿着箭体轴心向下，并通过调整姿态的变化来改变推力的方向。因此，为保证着陆段火箭垂直降落，箭体允许发生的最大倾斜角度，即推力方向的角度幅值约束如图 1所示，用γ表示。此外，火箭着陆段的约束还包括软着陆约束，即到达目标落点时的速度为0。

1.2 最优控制问题模型

基于式(1)动力学方程约束，考虑燃料最优的性能指标，火箭着陆段轨迹优化问题可以描述为

$ \min J = - m\left( {{t_f}} \right). $

(2)

状态约束：

$ {\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{r}}\left( {{t_0}} \right) = {\mathit{\boldsymbol{r}}_0},}\\ {\mathit{\boldsymbol{v}}\left( {{t_0}} \right) = {\mathit{\boldsymbol{v}}_0},}\\ {m\left( {{t_0}} \right) = {m_0}.} \end{array}} \right.} $

(3)

$ {\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{r}}\left( {{t_f}} \right) = 0,}\\ {\mathit{\boldsymbol{v}}\left( {{t_f}} \right) = 0,}\\ {m\left( {{t_f}} \right) \ge {m_{{\rm{dry}}}}.} \end{array}} \right.} $

(4)

$ {\arctan \frac{{\sqrt {r_x^2(t) + r_z^2(t)} }}{{{r_y}(t)}} \le \beta .} $

(5)

其中m_dry表示火箭干重。

控制量约束：

$ {{T_{\min }} \le \left\| {\mathit{\boldsymbol{T}}(t)} \right\| \le {T_{\max }}.} $

(6)

$ {\arctan \frac{{\sqrt {T_x^2(t) + T_z^2(t)} }}{{{T_y}(t)}} \le \gamma .} $

(7)

其中T_min和T_max分别表示推力幅值的上下限约束。

由于着陆段火箭运动的速度通常较小，式(1)中的阻力项在求解时忽略，闭环仿真中作为扰动项加以分析考虑。通过观察上述最优控制问题，式(5)和式(7)均为二阶锥形式的约束，符合SOCP问题的约束类型，而非线性动力学方程的等式约束和控制量幅值的不等式约束均为非凸约束，需通过凸化处理后，才能用凸优化方法求解。

2 原问题的凸化处理

采用凸优化方法求解一般优化问题的基本思路，即通过凸化手段，将原问题中的非凸约束，采用变量替换、松弛、线性化等凸化手段进行处理，并将问题离散后采用IPM求解器进行求解。对于原问题中非凸因素的处理，通常有逐次凸化和无损凸化2种凸化方法^[21]。其中，逐次凸化主要用于处理动力学方程等一些非线性约束，通过逐次线性化等方式，将非线性约束通过线性化和逐次逼近的方法处理，具有全局收敛性；而无损凸化，主要采用松弛、变量替换等方式，处理推力幅值的非凸约束，不影响问题的最优性。本节将分别采用逐次凸化和无损凸化方法，处理非线性动力学方程和推力幅值约束。

2.1 动力学方程处理

式(1)的动力学方程表示为时间连续的等式非线性约束，且时间为自由变量。在凸化处理时，通常需要将时间离散，转化为离散点上的等式线性约束。由于火箭着陆问题的复杂性，飞行时间难以预先设定，因此末端时刻为待求量。对于这种问题的处理，文献中通常有2种处理方法，一种是通过固定末端时刻，将问题凸化处理求解后，再将末端时刻作为参数量进行一维搜索^[22]；另一种方法，则是在动力学方程凸化处理前，先将非线性微分方程，转化为离散的差分代数方程约束，保留时间项作为优化变量，并通过线性化处理该非线性代数方程约束^[23]。对比2种处理手段，后者省去了将凸优化求解程序外加一维循环的流程，因此具有更好的收敛特性，本文将采取该方法处理动力学方程的约束。

取变量x表示动力学方程的状态量，即火箭的位置、速度矢量和质量，可将式(1)的动力学方程写如下等价的紧凑形式：

$ \mathit{\boldsymbol{\dot x}}(t) = \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}}(t),\mathit{\boldsymbol{u}}(t)). $

(8)

其中u表示归一化的推力矢量控制变量，定义为

$ \mathit{\boldsymbol{u}}(t) = \frac{{\mathit{\boldsymbol{T}}(t)}}{{{T_{\max }}}}. $

(9)

取时间域上的N个离散点，相邻离散点间的时间间隔可以表示为

$ \Delta t = \frac{{{t_f}}}{{N - 1}}. $

(10)

将式(8)用梯形公式离散后可以表示为

$ \mathit{\boldsymbol{x}}[k + 1] = \mathit{\boldsymbol{x}}[k] + \frac{1}{2}\Delta t(\mathit{\boldsymbol{f}}[k] + \mathit{\boldsymbol{f}}[k + 1]). $

(11)

其中k表示离散点序号，取值范围为[1, N-1]。

在式(11)中，等式右端第二项，包含非线性的动力学方程表达式，以及系数变量Δt，仍为非线性约束，因此可以通过线性化，取状态量x、控制量u和时间项Δt的展开式的一阶近似，得出如下紧凑形式的表达式：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}_k}\mathit{\boldsymbol{x}}[k] + {\mathit{\boldsymbol{A}}_{k + 1}}\mathit{\boldsymbol{x}}[k + 1] + {\mathit{\boldsymbol{B}}_k}\mathit{\boldsymbol{u}}[k] + }\\ {{\mathit{\boldsymbol{B}}_{\mathit{\boldsymbol{k}} + 1}}\mathit{\boldsymbol{u}}[k + 1] + \mathit{\boldsymbol{C}} \cdot \Delta t + \mathit{\boldsymbol{D}} = 0.} \end{array} $

(12)

其中：

$ {{\mathit{\boldsymbol{A}}_k} = \mathit{\boldsymbol{I}} + {{\left. {\frac{1}{2}\Delta {t^\prime }\left( {{\mathit{\boldsymbol{f}}^\prime }[k] + {\mathit{\boldsymbol{f}}^\prime }[k + 1]} \right) \cdot \frac{{\partial f}}{{\partial x}}} \right|}_{{\mathit{\boldsymbol{x}}_k},{\mathit{\boldsymbol{u}}_k}}},} $

(13)

$ {{\mathit{\boldsymbol{A}}_{k + 1}} = {{\left. {\frac{1}{2}\Delta {t^\prime }\left( {{\mathit{\boldsymbol{f}}^\prime }[k] + {\mathit{\boldsymbol{f}}^\prime }[k + 1]} \right) \cdot \frac{{\partial \mathit{\boldsymbol{f}}}}{{\partial \mathit{\boldsymbol{x}}}}} \right|}_{{\mathit{\boldsymbol{x}}_{k + 1}},{\mathit{\boldsymbol{u}}_{k + 1}}}} - \mathit{\boldsymbol{I}},} $

(14)

$ {{\mathit{\boldsymbol{B}}_k} = {{\left. {\frac{1}{2}\Delta {t^\prime }\left( {{\mathit{\boldsymbol{f}}^\prime }[k] + {\mathit{\boldsymbol{f}}^\prime }[k + 1]} \right) \cdot \frac{{\partial \mathit{\boldsymbol{f}}}}{{\partial u}}} \right|}_{{\mathit{\boldsymbol{x}}_k},{\mathit{\boldsymbol{u}}_k}}},} $

(15)

$ {{\mathit{\boldsymbol{B}}_{k + 1}} = {{\left. {\frac{1}{2}\Delta {t^\prime }\left( {{\mathit{\boldsymbol{f}}^\prime }[k] + {\mathit{\boldsymbol{f}}^\prime }[k + 1]} \right) \cdot \frac{{\partial \mathit{\boldsymbol{f}}}}{{\partial u}}} \right|}_{{\mathit{\boldsymbol{x}}_{k + 1}},{\mathit{\boldsymbol{u}}_{k + 1}}}},} $

(16)

$ {\mathit{\boldsymbol{C}} = \frac{1}{2}\left( {{\mathit{\boldsymbol{f}}^\prime }[k] + {\mathit{\boldsymbol{f}}^\prime }[k + 1]} \right),} $

(17)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{D}} = {\mathit{\boldsymbol{x}}^\prime }[k] - {\mathit{\boldsymbol{x}}^\prime }[k + 1] - {\mathit{\boldsymbol{A}}_k}{\mathit{\boldsymbol{x}}^\prime }[k] - }\\ {{\mathit{\boldsymbol{A}}_{k + 1}}{\mathit{\boldsymbol{x}}^\prime }[k + 1] - {\mathit{\boldsymbol{B}}_k}{\mathit{\boldsymbol{u}}^\prime }[k] - }\\ {{\mathit{\boldsymbol{B}}_{k + 1}}{\mathit{\boldsymbol{u}}^\prime }[k + 1] - \mathit{\boldsymbol{C}}\Delta {t^\prime }.} \end{array} $

(18)

其中上标“^′”用来表示前一次迭代求解得到的已知量。

从式(13)—(18)的形式可以看出，式(12)中所有未知变量前的系数均为常数矩阵，方程为等式线性约束，满足凸优化问题的等式约束要求。此外，为保逐次凸化处理后问题的收敛性和等价性，通常需要将方程的输入量施加信赖域约束，即：

$ {\left\| {\delta \mathit{\boldsymbol{u}}[k]} \right\| \le {\eta _u}[k],} $

(19)

$ {|\delta \Delta t| \le {\eta _{\Delta t}}.} $

(20)

其中η_u和η_Δt需要通过惩罚项增广至目标函数中。

由于逐次凸化过程中，问题被做了大量的近似处理，且状态量和控制量均受到严格的约束，在逐次逼近迭代求解过程中，极易出现前几次迭代无可行解的问题。为避免此类问题，Szmuk等^[23]提出了一种松弛技术，即在控制力之外，添加一项虚拟加速度，并在目标函数中对该项施加较大的惩罚项系数，可以提高算法的收敛性。

2.2 推力幅值约束处理

火箭发动机推力幅值约束如式(6)所示，属于非凸约束。引入松弛变量Γ，将原约束改写为如下形式：

$ {\left\| {\mathit{\boldsymbol{T}}[k]} \right\| \le \mathit{\Gamma }[k],} $

(21)

$ {{T_{\min }} \le \Gamma [k] \le {T_{\max }}.} $

(22)

同时，式(1)中质量变化率中的推力项也由新的松弛变量代替：

$ \dot m = - \frac{\mathit{\Gamma }}{{{I_{{\rm{sp}}}}{g_0}}}. $

(23)

通过Pontryagin极大值原理可以证明，上述松弛变换前后的两个问题，最终的收敛解具有等价性，即松弛变换后的问题，一定会收敛到使得不等式(21)活跃的最优解上，因此被称为无损松弛^{[16, 21]}。

2.3 凸优化问题模型

基于上述对问题的凸化处理，火箭着陆段连续时间轨迹优化问题，可以写成如下有限维形式的SOCP问题：

$ \begin{array}{c} \begin{array}{*{20}{c}} {\min J = - m[N] + {\omega _u} \cdot \left\| {{\eta _u}} \right\| + }\\ {{\omega _{\Delta t}} \cdot \left| {{\eta _{\Delta t}}} \right| + {\omega _{av}} \cdot \left\| {{\mathit{\boldsymbol{a}}_v}} \right\|,}\\ {{\rm{s}}.{\rm{t}}.} \end{array}\\ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{A}}_k}\mathit{\boldsymbol{x}}[k] + {\mathit{\boldsymbol{A}}_{k + 1}}\mathit{\boldsymbol{x}}[k + 1] + {\mathit{\boldsymbol{B}}_k}\mathit{\boldsymbol{u}}[k] + }\\ {{\mathit{\boldsymbol{B}}_{k + 1}}\mathit{\boldsymbol{u}}[k + 1] + C \cdot \Delta t + D = {\bf{0}},}\\ {\mathit{\boldsymbol{r}}[1] = {\mathit{\boldsymbol{r}}_0},\mathit{\boldsymbol{v}}[1] = {\mathit{\boldsymbol{v}}_0},m[1] = {m_0},}\\ {\mathit{\boldsymbol{r}}[N] = {\bf{0}},{\bf{v}}[N] = {\bf{0}},m[N] \ge {m_{{\rm{dry}}}},}\\ {\sqrt {r_x^2[k] + r_z^2[k]} - {r_y}[k]\tan \beta \le 0,}\\ {\sqrt {T_x^2[k] + T_z^2[k]} - {T_y}[k]\tan \gamma \le 0,}\\ {\left\| {\mathit{\boldsymbol{T}}[k]} \right\| \le \mathit{\Gamma }[k],{T_{\min }} \le \mathit{\Gamma }[k] \le {T_{\max }},}\\ {\left\| {\delta \mathit{\boldsymbol{u}}[k]} \right\| \le {\eta _u}[k],|\delta \Delta t| \le {\eta _{\Delta t}}.} \end{array}} \right. \end{array} $

(24)

其中ω_u、ω_Δt和ω_av分别表示控制量信赖域、时间信赖域和虚拟加速度的惩罚项系数。

3 凸化处理问题的求解

上述SOCP形式的问题，可以采用现有成熟的求解器，例如面向MATLAB接口的通用求解器CVX^[24]求解。此类通用求解器往往不满足实时性和在线计算的需求，因此仅适用于算法的验证。为了进一步提高凸优化问题的求解效率，满足在线计算和实时制导的需求，通常需要针对特定的问题，利用IPM方法定制求解器^[25]。本文采用了一种基于IPM方法面向嵌入式计算的SOCP求解器ECOS(embedded conic solver)^[26]，将问题进行定制求解，使算法具备实时性。

3.1 问题定制求解过程

SOCP问题求解的一般过程如图 2所示。通过上述对原优化问题的凸化、离散，使得原始连续无限维的优化问题，转化为有限维离散点的SOCP问题。在给定该问题的常量、初值等外部数据后，通用IPM求解器一般会根据式(24)解释性语言的问题描述，在求解器内部转化为对应维有限个SOCP子问题，然后再根据IPM方法，将该簇SOCP求解，并生成满足通用求解器输出接口的解。而求解SOCP问题完成对原问题的求解，通常需要经过数次迭代，则需要将上次求解的结果作为最新初值，重复上一过程，直到满足优化问题的收敛条件。

分析上述通用求解器求解一般SOCP问题的过程，存在大量的重复性计算，求解效率较低，且往往存在对特定软件工具包的依赖性，不具备在线制导的要求。本文针对ECOS的格式化接口，将式(24)描述的问题直接描述为一系列子问题实例。对于每个离散点上的子问题实例约束中出现的系数矩阵，采用列压缩算法，提取反映稀疏矩阵的全部信息数组，极大地节约了存储空间^[26]。此外，该数组信息对特定的SOCP问题，只有记录非零元素值的数组会随着迭代而变化，其余2个数组均为固定的常数数组。因此，在反复的迭代过程中，仅需要修改一个数组的信息即可，极大地提高了迭代计算效率。

3.2 仿真算例

以坐标轴原点作为火箭着陆段的目标落点，考虑火箭以起始位置r₀=[500 m, 1 600 m, 0 m]^T，初速度v₀=[-5 m/s, 26.7 m/s, 0 m/s]^T为着陆段起始状态。部分仿真参数如表 1所示。假设火箭发动机推力具备深度可调节能力，上限为235 kN，下限为0 kN，比冲为218 s，火箭初始质量为24 t，干重10 t。本算例考虑从该起始位置首次产生整个着陆段的制导律，为末端时刻自由的燃料最优问题。

本仿真算例中，飞行时间收敛的精度设置为1×10^-6，控制量迭代精度为1×10^-6，离散点数取30个。用于求解效率对比的通用求解器采用了CVX中的SDPT3求解器。仿真平台为配置4核3.4 GHz i7 4770处理器、8 G内存的PC机。仿真结果及计算效率对比如图 3及表 2所示。

从图 3的位置速度求解结果对比中可以看出，2种求解器得到了一致的优化结果。从表 2中的计算效率对比来看，2种求解器所用迭代次数相同，但CPU有效计算时间方面，通用求解器达到了百毫秒的量级，由于复杂的前处理和数据迭代，总计算时长超过了22 s。而定制求解器的总计算时间为63.720 ms，单次平均迭代时间为10.620 ms，具备实时在线制导的条件。

4 GNC一体化闭环数值仿真

本文建立了火箭着陆段的燃料最优控制问题模型，并转化为SOCP问题通过定制求解器完成了快速求解，使火箭着陆段制导问题具备了实时在线求解的条件。然而，上述过程仅针对火箭着陆段制导问题完成一次求解。火箭着陆过程中，需要不断地根据当前位置实时地进行在线制导，反复重复前述凸优化求解过程。此外，在前述开环仿真算例中，未考虑环境大气的干扰、推力执行误差和节流能力等因素，算法在闭环控制系统中的鲁棒性尚未可知。本节将搭建GNC一体化数值闭环仿真平台，以测试前述制导算法在实际的应用过程中对大气扰动、推力误差、导航系统误差、制导过程延时等非理想因素作用下的鲁棒性。

4.1 GNC一体化闭环数值仿真框架

在火箭着陆闭环仿真过程中，以某一状态触发精确软着陆制导程序后，即进入GNC一体化闭环数值仿真程序。在一个控制周期内，基于当前状态，制导模块采用前述凸优化在线制导算法，生成燃料最优推力曲线，控制模块则根据当前周期的最优推力曲线，输出推力执行指令。在该闭环仿真框架内，火箭模型模块根据控制模块输出的推力执行指令，进行数值积分，得到一个周期后火箭的状态。导航模块则负责采集一个周期后的火箭状态，并传递给制导模块，从而进入下一个制导周期，直到火箭的位置速度满足精确软着陆的要求，或火箭飞行高度达到0。本文所搭建的GNC一体化闭环数值仿真框架如图 4所示。

在各分系统的测试中，火箭着陆段的起始状态均设置为r₀=[0 m, 1 760 m, 0 m]^T，v₀=[0 m/s, -29.3 m/s, 0 m/s]^T，其余火箭参数和制导算法参数设置同表 1。在本文的仿真算例中，考虑火箭落地后位置误差小于1 m、速度误差小于1 m/s即为满足精确软着陆的需求。如无特殊说明，本文中闭环控制周期均假设为200 ms(控制频率为5 Hz)。

4.2 GNC分系统测试

本节将针对上述GNC一体化仿真平台的各个分系统进行测试，分析大气扰动，制导、控制与导航模块各种干扰因素下制导算法的鲁棒性。

4.2.1 环境扰动适应性

首先，考虑算法在大气扰动下的鲁棒性。为了收敛性和计算效率，在线制导算法并未考虑大气扰动因素。在闭环仿真时，将该因素考虑到火箭的运动模型中，会导致实际的飞行状态与制导算法规划的路径不同，因此需要制导模块根据当前状态，实时的更新制导律。考虑火箭飞行过程中受大气阻力，以及侧向10 m/s的切变风扰动(风力约6级)，测试结果如图 5所示。

由于受到大气阻力和侧向风影响，如按照初次规划的飞行制导方案执行，火箭实际飞行轨迹则会发生较大的水平位移，且垂直方向的速度提前减速到0。通过200 ms制导和控制周期的实时在线修正，火箭最终的着陆误差为3.252×10^-2 m和7.941×10^-2 m/s，满足精确软着陆的需求。

4.2.2 控制系统误差适应性

为测试在线制导算法对控制系统误差的适应性，考虑闭环控制框架下控制系统的推力的执行方案和执行误差对仿真结果的影响。仍假设以200 ms为制导和控制周期，测试对比了3种不同的推力执行方案：

1) 控制周期内推力常值：考虑闭环控制周期内推力固定为一常值，不可连续变化。该常值取值策略为控制周期内最优推力曲线的插值平均数。

2) 推力分档：考虑推力不可连续调节的节流特性，且只能按照0.2、0.4、0.6、0.8、1.0倍最大幅值等几个固定档位执行。推力分档策略为，根据方案(1)中的常值对应取最近的推力档位执行。

3) 推力分档且存在最大±5%的幅值偏差：在考虑推力分档的方案(2)基础上，施加±5%以内的幅值随机偏差。

3种情况下推力执行方案以及对应的闭环仿真结果如表 3所示。

从上述结果中可以看出，对于推力控制周期内常值，以及推力分档执行且幅值存在±5%的偏差时，火箭均能实现精确软着陆，着陆误差为10^-2量级(m或m/s)，可见该算法对控制系统的误差适应性较好。

4.2.3 导航系统误差适应性

假设火箭飞行过程中，导航系统存在一定的测量偏差。综合考虑绝对误差和相对误差，导航模块获取火箭的实时状态后，传递给制导模块前施加如下形式的误差：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{r}}_{{\rm{test }}}} = {\mathit{\boldsymbol{r}}_{{\rm{true }}}} \pm \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \min \left( {10,\left| {{\mathit{\boldsymbol{r}}_{{\rm{true }}}}} \right| \cdot \varepsilon \% } \right) \times {\mathop{\rm rand}\nolimits} (0,1),\\ {\mathit{\boldsymbol{v}}_{{\rm{test }}}} = {\mathit{\boldsymbol{v}}_{{\rm{true }}}} \pm \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \min \left( {1,\left| {{\mathit{\boldsymbol{v}}_{{\rm{true }}}}} \right| \cdot \varepsilon \% } \right) \times {\mathop{\rm rand}\nolimits} (0,1). \end{array} \right. $

(25)

其中：ε表示相对误差的最大范围，在仿真中分别取不同值以测试算法对导航系统误差的适应边界；r_test、v_test和r_true、v_true分别表示导航模块输出和输入的位置速度矢量。测试结果如表 4所示。

从表 4可以看出，当随机误差小于5%时，火箭着陆的误差能够达到10^-2 m(或m/s)量级；而当误差增加到10%及15%时，着陆的位置误差和速度误差均明显增大，最高达到了3.573×10^-1 m的位置误差和4.018×10^-1 m/s的落地速度。当导航系统误差较大时，由于火箭无法准确判定当前的位置和速度，因此精确软着陆也就难以实现。

4.2.4 制导系统延时适应性

从节3的仿真结果来看，制导算法每次计算更新制导方案均需要一定的时间。为提高制导方案的计算效率，每次更新制导方案时，将上次生成的制导方案通过时间插值作为下次计算制导方案的迭代初值，可将迭代次数降低至1~2次，计算时间约20~30 ms。然而，考虑到箭载计算机的计算能力有限，存在无法在控制周期(例如本文假设的200 ms)内完成一次制导方案生成的可能性，本文提出了一种基于状态估计器处理系统延时的制导方案，该方案的流程示意图如图 6所示。

相比图 4的闭环方案，该制导方案增加了一个基于状态估计器的外部制导回路。考虑到系统不同程度的延时，可设置不同的制导周期T_g。在一个制导周期内，在线制导算法以状态估计器生成的状态估计作为输入量进行制导方案的生成与更新，此后，闭环仿真系统仍以原来的控制周期按照当前制导方案执行。本文采用的状态估计方法为，依据导航系统测得的当前状态，以10步四阶Runge-Kutta(RK4)定步长积分一个制导周期T_g，估计一个制导周期后的状态，传递给导航模块。由于RK4定步长积分仅需要进行少量的代数运算，因此该方法可避免由于大量复杂计算带来的效率损失，且在系统延时较短时具有较高的精度。为测试该方案的可行性和准确性，本文分别设置3个不同的制导周期进行仿真，仿真结果如表 5所示。

可以看出，在系统延时较短的情况下，设置制导周期T_g为0.5 s，能够保证火箭着陆较高的精度；而当T_g设置为2 s时，着陆速度误差随之增大，达到了7.239×10^-1 m/s。

4.3 GNC一体化闭环数值仿真

综合考虑上述各分系统，对本文所提出的在线制导算法进行完整系统的GNC一体化闭环数值仿真测试，各分系统的仿真设置为：制导周期为1 s；考虑推力在控制周期内为常值；考虑大气阻力及5 m/s的侧向切变风；考虑±3 m(或±3%)的测量误差。

火箭起始状态设置为r₀=[0 m, 1 760 m, 0 m]^T，v₀=[0 m/s, -29.3 m/s, 0 m/s]^T，其余火箭参数和制导算法参数设置同表 1，闭环控制周期设置为200 ms。仿真结果对比如表 6所示。

火箭着陆过程的飞行轨迹及位置速度状态对比如图 7所示。可以看出，由于空气阻力的存在，火箭实际的飞行时长比初次预计飞行时间长约1.4 s，但由于实际的推力执行方案中，推力幅值普遍与初次规划的推力幅值较低，因此实际的燃料消耗少于初次规划预计的燃料消耗。从仿真结果来看，火箭最终着陆的位置误差和速度误差均满足预定要求，证明了本文所提出的在线制导算法和一体化闭环仿真系统，在环境大气扰动、控制系统不同情况下的误差、导航系统偏差及系统延时等各种复杂扰动的综合影响下，仍可以根据实时在线制导，完成火箭着陆段的精确软着陆，具备较高的算法鲁棒性和准确性。

5 结论

本文针对可回收火箭动力着陆段的制导问题，建立了最优控制问题模型，并通过无损凸化和逐次凸化方法，将该问题转化为二阶锥优化问题，通过内点法定制求解器快速求解。仿真结果表明：通过定制化求解，该制导算法能够在保证得到最优解的同时，实现了相比通用凸优化求解器更高的求解效率，达到毫秒量级的计算时间，具备了实时在线制导的条件。基于该制导算法，搭建了制导、导航与控制一体化数值仿真平台，以测试该算法的鲁棒性。闭环仿真结果表明：该算法能够适应大气扰动、推力执行系统误差、导航系统误差和系统延时等多重复杂干扰的影响，具备较高的鲁棒性。