A Multi-stage Robust Optimization Guidance Method for Endoatmospheric Powered Descent of Reusable Rockets
-
摘要: 针对大气层内可回收火箭的动力下降问题, 提出一种多阶段的鲁棒优化(Robust optimization, RO)方法. 由于大气层内存在未知风场, 如何在火箭下降段考虑这种不确定性具有十分重要的意义. 首先, 建立一个关于高度的不确定风场模型, 在该风场下给出火箭动力下降的鲁棒最优控制问题. 为了求解该问题, 使用一种对不等式约束采取一阶近似并将一阶项作为安全裕量加入约束的鲁棒优化方法, 得到一个可以求解的单阶段鲁棒优化算法. 其次, 定量给出安全裕量的上界, 基于该上界提出一种多阶段鲁棒优化算法, 避免单阶段鲁棒优化算法中安全裕量可能过大导致无法求解的问题. 最后, 通过仿真对比各个算法在多个实际风场下的性能, 结果表明所提出的多阶段鲁棒优化方法同时具有较高的落点精度和对于不同风场的鲁棒性.Abstract: This paper proposes a multi-stage robust optimization (RO) method for the endoatmospheric powered descent problem of reusable rockets. Due to the unknown wind field in the atmosphere, it is of great significance to consider this uncertainty during the rocket descent phase. Firstly, a model of uncertain wind field with respect to altitude is established, and a robust optimal control problem for rocket powered descent is formulated under this wind field. To solve this problem, a tractable single-stage RO algorithm is developed by approximating the inequality constraints using a first-order expansion and incorporating the first-order term as a safety margin. Secondly, an upper bound on the safety margin is quantitatively derived. Based on this upper bound, a multi-stage RO algorithm is proposed, which avoids the infeasibility problem caused by the excessively large safety margin in the single-stage RO algorithm. Finally, simulation results are presented to compare the performance of each algorithm under various actual wind fields. The results demonstrate that the proposed multi-stage RO method achieves both high landing accuracy and robustness against different wind fields.
-
近些年来, 由于可以显著降低太空探索的成本, 关于可回收火箭的研究受到广泛的关注. 20世纪50 ~ 60年代, 美国在航天飞机上首次实现航天运载器的部分重复使用, 随后进行大量的可重复使用运载器的研究[1]. 21世纪以来, SpaceX、Blue Origin、Maste Space System等商业公司同样开展大量关于可回收火箭的研究[2-3]. 动力下降作为可回收火箭中的一个重要问题, 要求利用稠密大气的空气阻力和火箭发动机的推力来实现火箭的定点软着陆, 是火箭着陆的最后一环. 由于环境存在许多干扰, 在动力下降段中常常存在大量的不确定性, 尤其是未知的风场, 这给动力下降段的制导控制带来挑战.
目前针对火箭动力下降问题的制导方法主要可以分为三类: 解析制导方法, 基于轨迹优化的制导方法, 基于学习的制导方法. 解析制导方法是指制导指令由火箭状态量解析表示的一系列方法, 包括多项式制导[4]、重力转弯制导[5]和近似最优解析制导[6] (零控位置/速度误差制导)等. 这些制导方法具有形式简单、易于实现等优点, 已经在阿波罗登月、嫦娥登月等实际场景中得到应用. 但是为了获得解析的制导指令, 往往使用较为简单的问题模型, 难以处理大气层内的动力下降问题. 基于学习的制导方法[7-8]则是近年来随着深度学习兴起而发展起来的方法, 通过训练可以显著提高方法对初始状态和模型不确定性的适应性. 但是基于学习的方法具有难以设计奖励函数、可解释性差等缺点, 因此目前还难以在实际场景中进行应用.
基于轨迹优化的制导方法[9-10]将控制和状态作为待优化变量, 在一定的性能指标下通过求解优化问题来获得制导指令, 其主要分为间接法和直接法. 间接法是将动力下降问题建模为最优控制问题, 通过极大值原理将最优控制问题转化为两点边值问题进行求解[11-12]. 但是对于较为复杂的问题模型来说, 求解转化后的两点边值问题同样较为困难, 因此间接法更多只能用在一些较为特殊的、化简后的模型. 直接法则是指通过对原始最优控制问题进行直接描述和寻优的方法. 凸优化是一种在动力下降问题中得到广泛应用的直接法, 基于凸优化的轨迹规划将原始的非凸问题转化或近似转化为凸问题, 再通过离散化将无穷维的连续系统转化为有限维的离散系统, 最后通过成熟的凸优化求解算法进行快速求解. 文献[13-14]首次在动力下降问题中提出无损凸化的概念, 通过引入松弛变量将控制幅值约束转化为凸约束并推广到更一般的情形, 从理论上严格证明转化前后最优解的等价性. 文献[15]则将上述结果推广到同时具有控制和状态约束的情况, 同样证明转化前后最优解的等价性. 但是上述无损凸化的方法只对特定的约束形式成立, 难以推广到动力着陆段的一般情形. 文献[16]提出序列线性化的方法, 针对非凸的、一般的非线性不等式, 每次迭代使用上一次迭代的轨迹并在其附近进行一阶近似展开, 求解近似的凸优化问题, 使用信赖域技术, 通过多步迭代使近似最优解收敛. 随后同样有一系列工作使用凸优化作为轨迹优化的工具[17-19]. 目前, 基于凸优化的轨迹规划方法已经被广泛应用于可回收火箭的动力着陆问题中. 但是上述方法针对的均是确定性问题, 未考虑不确定性的影响.
对于轨迹规划中不确定性的研究则相对较少. 一种考虑不确定性的方法是基于敏感度矩阵(系统状态关于未知参数的雅可比矩阵)的方法[20-23]. 文献[20]将敏感度矩阵作为增广状态并在原始的最优控制问题中增加关于增广状态的初始和终端约束, 通过求解增广的最优控制问题来进行鲁棒轨迹规划. 文献[21]则设计一个基于敏感度矩阵和风场观测器的临近最优制导律, 并且为处理真实风场的不确定性, 加入扰动观测补偿器以提升鲁棒性. 文献[22]考虑存在参数不确定性的轨迹跟踪问题, 通过最小化状态敏感度和输入敏感度进行鲁棒轨迹规划. 文献[23]将敏感度矩阵作为增广状态, 对不等式约束进行一阶近似并使用敏感度矩阵计算一阶项, 将原问题转化为一个易求解的优化问题. 此外, 基于多项式混沌展开(Polynomial chaos expansion, PCE)的方法[24]则是轨迹规划中处理不确定性的一种典型方法, 其主要思想为, 将存在不确定性的随机过程展开为一系列正交的随机多项式基底函数的加权和, 从而建立不确定输入与不确定输出的映射关系. 文献[25]则将PCE与凸优化的方法相结合, 在仅考虑初值不确定性的情况下, 通过序列线性化的方法降低求解问题的计算代价. 文献[26]通过引入谱分解技术来自适应更新基函数, 以处理长时间的优化场景可能导致PCE方法发散的问题. 文献[27]通过特征选择和序列采样提出序列增强PCE方法, 相比于原始的PCE方法具有更高的计算效率和精度.
可以看出, 现有对存在不确定性的动力下降问题或者轨迹规划的研究要么对不确定性(或者风场)的建模较为简单[23-24], 要么使用随机的PCE方法, 通过期望与方差等统计特性来衡量性能指标. 这种方法需要把原始的随机问题转化为更高维的确定性问题, 大大增加求解的复杂度. 此外, PCE方法的最优指标是在期望和方差等统计学意义下的平均最优, 难以应对极端场景. 本文则使用一种确定性的鲁棒优化(Robust optimization, RO)方法[23], 该方法转化后的问题维度远低于PCE方法, 并且考虑在参数最差时的最优性, 可以处理更加极端的情况. 本文的主要工作包括: 首先, 建立更加贴近实际场景的风场模型并在确定性框架下给出动力下降的RO问题; 其次, 使用一种对不等式约束采取一阶近似的RO方法, 得到一个可以易于求解的单阶段RO算法; 然后, 通过理论分析单阶段RO算法在终端约束范围较小时可能无解的缺陷, 进一步提出一种改进的多阶段RO算法; 最后, 通过仿真说明多阶段RO算法对未知不确定风场同时具有较高的最终状态精度和较强的鲁棒性.
1. 火箭着陆段问题建模
1.1 火箭动力学模型
大气层内的可回收火箭的动力下降问题是指火箭在发动机推力、重力和空气阻力的影响下从高空状态被引导到指定的着陆点. 在本文中, $ x $轴、$ y $轴和$ z $轴分别指向东、北和天空. 为了简便, 文中的下标$ (\cdot)_x $、$ (\cdot)_y $和$ (\cdot)_z $分别代表对应量在$ x $轴、$ y $轴和$ z $轴上的分量. 考虑如下大气层内火箭的动力学模型
$$ \begin{cases} \dot{{\boldsymbol{r}}} = {\boldsymbol{v}}\\ \dot{{\boldsymbol{v}}} = {\boldsymbol{g}} + \dfrac{{\boldsymbol{u}}}{m} + \dfrac{{\boldsymbol{D}}( r_z,{\boldsymbol{v}},{\boldsymbol{v}}_w )}{m}\\ \dot{m} = -\dfrac{1}{\kappa} \Vert{\boldsymbol{u}} \Vert_2 \end{cases} $$ (1) 其中, $ {\boldsymbol{r}} = \left[r_x,r_y,r_z\right]^{\rm{T}}\in {{\bf{R}}}^3 $, $ {\boldsymbol{v}} = \left[v_x,v_y,v_z\right]^{\rm{T}}\in {{\bf{R}}}^3 $分别是位置向量和速度向量, $ {\boldsymbol{g}} $代表重力加速度, $ m\in {{\bf{R}}} $是火箭质量, $ \kappa\in {{\bf{R}}} $代表火箭恒定的燃料消耗系数, $ {\boldsymbol{u}} = \left[u_x,u_y,u_z\right]^{\rm{T}}\in {{\bf{R}}}^3 $是系统的控制输入, 代表火箭引擎提供的推力, 其满足如下幅值约束
$$ \left| u_i\right| \leq \bar{U}_i,\; i = x, y, z $$ $ \bar{U}_i $代表各个方向对应的推力幅值. 此外, 式(1)中
$$ {\boldsymbol{D}}(r_z,{\boldsymbol{v}},{\boldsymbol{v}}_w ) = -\frac{1}{2}\rho\Vert {\boldsymbol{v}}-{\boldsymbol{v}}_w \Vert_2 S_{\mathrm{ref} } C_D \left({\boldsymbol{v}} -{\boldsymbol{v}}_w\right) $$ (2) 其中, $ {\boldsymbol{D}}(r_z,{\boldsymbol{v}},{\boldsymbol{v}}_w ) $代表空气阻力, $ {\boldsymbol{v}}_w (r_z)\in {{\bf{R}}}^3 $代表风场向量, $ \rho(r_z)\in {{\bf{R}}} $是大气密度, $ S_{\mathrm{ref} }\in {{\bf{R}}},\; C_D\in {{\bf{R}}} $分别代表火箭的参考面积和空气阻力系数.
设初始时刻为$ t_0 = 0 $, 该问题的初始条件为
$$ {\boldsymbol{r}}(t_0) = {\boldsymbol{r}}_0, \; {\boldsymbol{v}}(t_0) = {\boldsymbol{v}}_0,\; m(t_0) = m_0 $$ 终端时刻$ t_f $自由, 在终端时刻需要满足终端条件
$$ {\boldsymbol{r}}(t_f) = {\boldsymbol{r}}_f,\; {\boldsymbol{v}}(t_f) = {\boldsymbol{v}}_f $$ (3) $ {\boldsymbol{r}}_f,\; {\boldsymbol{v}}_f $代表期望的终端位置和速度. 性能指标为终端时刻和推力二次型积分之间的加权和
$$ J = t_f+ \int_{t_0}^{t_f} \gamma\; {\boldsymbol{u}}^{\rm{T}} {\boldsymbol{u}} \; {\rm{d}}t $$ (4) 其中, $ \gamma $代表加权权重.
记$ {\boldsymbol{X}} = [{\boldsymbol{r}}^{\rm{T}},{\boldsymbol{v}}^{\rm{T}},m]^{\rm{T}} $, $ {\boldsymbol{X}}_0 = [{\boldsymbol{r}}_0^{\rm{T}},{\boldsymbol{v}}_0^{\rm{T}},m_0 ]^{\rm{T}} $, 那么动力学模型和初值可以写为如下紧凑形式
$$ \begin{align} &\dot{{\boldsymbol{X}}} = {\boldsymbol{f}}({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{v}}_w) \end{align} $$ (5) $$ \begin{align} &{\boldsymbol{X}}(t_0) = {\boldsymbol{X}}_0 \end{align} $$ (6) 现假设有一个标称风场$ \hat{{\boldsymbol{v}}}_w (r_z) $已知, 那么在$ \hat{{\boldsymbol{v}}}_w (r_z) $下的火箭动力下降问题为:
标称优化(Nominal optimization, NO)问题1. 给定初始状态$ {\boldsymbol{X}}_0 $以及目标终端状态$ {\boldsymbol{X}}_f $, 选择合适的控制律$ {\boldsymbol{u}} $, 在$ {\boldsymbol{v}}_w = \hat{{\boldsymbol{v}}}_w (r_z) $的情况下, 满足式(3)、式(5)和式(6)的同时最小化性能指标(4).
NO问题使用标称风场信息进行计算, 因此可以直接使用标准优化算法进行求解. 基于NO问题可以得出如下NO算法.
NO算法. 输入$ {\boldsymbol{X}}_0 $, $ {\boldsymbol{r}}_f $, $ {\boldsymbol{v}}_f $, 求解NO问题, 输出最优控制和相应的状态轨迹与终端时刻.
1.2 不确定风场模型
第1.1节中的NO问题使用标称风场, 但在实际场景中, 风场具有很强的不确定性. 为了在模型中考虑不确定性的影响, 在风场模型中引入不确定参数.
假设风场在$ y $轴方向存在不确定性, 不确定性风场为$ v_{w_y} = v_{w_y} (r_z,{\boldsymbol{s}}) $, $ {\boldsymbol{s}} = \left[s_1,\cdots,s_d \right]^{\rm{T}} \in {{\bf{R}}}^d $是未知参数, $ d $代表未知参数的维度. 具体而言
$$ v_{w_y} (r_z,{\boldsymbol{s}}) = s_1 r_z^{*d-1} +s_2 r_z^{*d-2} + \cdots +s_d $$ 其中, $ r_z^* = (r_z-\mu) / \sigma $代表标准化后的高度, $ \mu,\; \sigma $均为已知参数, 分别代表拟合数据中高度的均值和方差. 可以证明, 对于从零开始、等间距采样的数据点, $ r_z^* $满足
$$ \left|r_z^*\right|\leq \sqrt3 $$ (7) 未知参数$ {\boldsymbol{s}} $的取值范围满足
$$ {\boldsymbol{s}}\in {\boldsymbol{S}}_\tau := \left\{\hat{{\boldsymbol{s}}}+ \tau {\boldsymbol{A}} {\boldsymbol{\delta}}:\Vert {\boldsymbol{\delta}}\Vert_\infty \leq 1 \right\} $$ (8) 其中, $ \hat{{\boldsymbol{s}}} = \left[\hat{s}_1,\cdots,\hat{s}_d\right]^{\rm{T}} $代表参数的标称值, $ \tau>0 $代表参数变化的幅值, $ {\boldsymbol{\delta}}\in {{\bf{R}}}^d $代表参数集合变量, $ {\boldsymbol{A}} $为根据参数变化特点选择的矩阵, 很多情况下, 可以将$ {\boldsymbol{A}} $选为正定对角矩阵, 此时, $ {\boldsymbol{S}}_\tau $是一个矩形立方体集合. 上述模型可以通过事先测量的风场数据得到.
此外, 由于实际场景中$ z $轴方向的风速可以忽略不计, 因此假设$ z $轴方向不存在风. 对于$ x $轴方向的风场, 为了简单起见, 本文假设$ x $轴方向的风场是确定已知的, 满足
$$ v_{w_x}(r_z) = kr_z+b $$ 其中, $ k,\; b $均为已知参数. 最后, 记$ {\boldsymbol{v}}_w (r_z,{\boldsymbol{s}}) $代表上述定义的、包含不确定参数$ {\boldsymbol{s}} $的风场.
1.3 基于鲁棒优化的问题建模
在第1.2节给出的不确定风场$ {\boldsymbol{v}}_w (r_z,{\boldsymbol{s}}) $下, 考虑如下鲁棒优化问题:
鲁棒优化(RO)问题2.
$$ \begin{split} & \min\limits_{{\boldsymbol{u}}} \; J = t_f+\int_{t_0}^{t_f} \gamma \; {\boldsymbol{u}}^{\rm{T}} {\boldsymbol{u}} \; {\rm{d}}t\\ &\text{s.t.}\\ &\dot{{\boldsymbol{X}}} = {\boldsymbol{f}}({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{v}}_w)\\ &{\boldsymbol{X}}(t_0) = {\boldsymbol{X}}_0\\ & \left\{\begin{aligned} &\left| r_i(t_f) - r_{f_i} \right| \leq L_{ri}\\ &\left| v_i(t_f) - v_{f_i} \right| \leq L_{vi} \end{aligned}\right.\; \forall {\boldsymbol{s}} \in {\boldsymbol{S}}_\tau , i = x,y,z \end{split} $$ (9) 其中, $ L_{ri},\; L_{vi} $代表事先确定的终端约束范围.
注1. RO问题2在满足对任意参数都不违背约束条件的可行解中寻求最小化优化指标的解, 这使得RO问题2的最优解即使在最差的参数条件下也能拥有较好的效果.
注2. 注意到RO问题2的终端约束和NO问题1的终端约束不同, 这是因为在风场存在不确定性的情况下, 使火箭严格满足给定的终端条件式(3)是难以实现的, 而要求火箭尽可能落在给定终端状态的一定范围之内更加具有可行性.
注3. 在式(9)中选择不同的约束范围$ L_{ri} $, $ L_{vi} $, 是因为第1.2节定义的风场不确定性出现在$ y $轴方向, 所以相比于$ y $轴方向而言, 火箭由风场不确定性导致的、在$ x $轴和$ z $轴方向的偏差较小, 自然也应该取不同的约束范围.
一般而言, 直接获得RO问题2的最优解是十分困难的, 特别是当系统的动态方程较为复杂时. 接下来将简要介绍一种对不等式约束采用一阶近似的RO方法, 该方法可以将RO问题2近似转化为一个易求解的最优控制问题, 并且转化后的问题维度低于PCE方法. 该方法更详细的推导过程见文献[23]. 考虑如下一般形式的非线性优化问题
$$ \begin{split} & \min\limits_{{\boldsymbol{X}},{\boldsymbol{u}}} \; \phi({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{s}}) \\ &\text{s.t.}\\ &{\boldsymbol{F}}({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{s}}) = {\bf{0}}\\ &{\boldsymbol{G}}({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{s}})\leq {\bf{0}},\; \forall {\boldsymbol{s}} \in {\boldsymbol{S}}_{\tau} \end{split} $$ (10) 在上述优化问题中, 通过等式约束确定一个$ {\boldsymbol{X}} $与$ {\boldsymbol{u}} $, $ {\boldsymbol{s}} $之间的隐式关系式, 不妨设为
$$ {\boldsymbol{X}} = {\boldsymbol{X}}({\boldsymbol{u}},{\boldsymbol{s}}) $$ 将$ {\boldsymbol{F}}({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{s}}) = {\bf{0}} $两端同时关于$ {\boldsymbol{s}} $求微分, 可得$ {\boldsymbol{X}}_{{\boldsymbol{s}}}({\boldsymbol{u}},{\boldsymbol{s}}) $满足
$$ {\boldsymbol{F}}_{{\boldsymbol{X}}}({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{s}}){\boldsymbol{X}}_{{\boldsymbol{s}}}+ {\boldsymbol{F}}_{{\boldsymbol{s}}}({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{s}}) = {\bf{0}} $$ (11) 其中, 变量下标代表关于该变量求偏导数.
另一方面, 在隐式关系$ {\boldsymbol{X}} = {\boldsymbol{X}}({\boldsymbol{u}},{\boldsymbol{s}}) $下, 令
$$ \hat{{\boldsymbol{G}}}({\boldsymbol{u}},{\boldsymbol{s}}) = {\boldsymbol{G}}({\boldsymbol{X}}\left({\boldsymbol{u}},{\boldsymbol{s}}),{\boldsymbol{u}},{\boldsymbol{s}}\right) $$ (12) 设$ \hat{g}_i ({\boldsymbol{u}},{\boldsymbol{s}}) $为$ \hat{{\boldsymbol{G}}}({\boldsymbol{u}},{\boldsymbol{s}}) $的第$ i $个元素, $ i = 1,\cdots, m $. 那么不等式约束$ {\boldsymbol{G}}({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{s}})\leq {\bf{0}} $等价于
$$ \max\limits_{{\boldsymbol{s}} \in {\boldsymbol{S}}_\tau } \hat{g}_i({\boldsymbol{u}},{\boldsymbol{s}})\leq 0,\; i = 1,\cdots,m $$ 其中, $ {\boldsymbol{S}}_\tau $由式(8)定义. 将$ \hat{g}_i({\boldsymbol{u}},{\boldsymbol{s}}) $一阶泰勒展开可得
$$ \hat{g}_i({\boldsymbol{u}},{\boldsymbol{s}}) \approx \hat{g}_i({\boldsymbol{u}},\hat{{\boldsymbol{s}}}) + \tau \langle \nabla_{{\boldsymbol{s}}} \hat{g}_i({\boldsymbol{u}},\hat{{\boldsymbol{s}}}),{\boldsymbol{A}}{\boldsymbol{\delta}} \rangle $$ 其中, $ \nabla_{{\boldsymbol{s}}} $代表关于$ {\boldsymbol{s}} $的梯度算子, $ \langle\cdot,\cdot\rangle $代表内积. 因此
$$ \begin{split} & \max\limits_{{\boldsymbol{s}} \in {\boldsymbol{S}}_\tau } \hat{g}_i({\boldsymbol{u}},{\boldsymbol{s}}) \approx \hat{g}_i({\boldsymbol{u}},\hat{{\boldsymbol{s}}})\;+\\ & \quad \tau\max\limits_{\Vert {\boldsymbol{\delta}}\Vert_\infty\leq 1} \langle {\boldsymbol{A}}^{\rm{T}}\nabla_{{\boldsymbol{s}}} \hat{g}_i({\boldsymbol{u}},\hat{{\boldsymbol{s}}}),{\boldsymbol{\delta}} \rangle = \\ &\quad \hat{g}_i({\boldsymbol{u}},\hat{{\boldsymbol{s}}})+\tau \Vert {\boldsymbol{A}}^{\rm{T}}\nabla_{{\boldsymbol{s}}} \hat{g}_i({\boldsymbol{u}},\hat{{\boldsymbol{s}}}) \Vert_1 \end{split} $$ (13) 其中, 等号成立是因为赫尔德不等式(Holder inequality). 为了求取$\nabla_{{\boldsymbol{s}}} \hat{g}_i({\boldsymbol{u}},\hat{{\boldsymbol{s}}})$, 将式(12)关于$ {\boldsymbol{s}} $求偏导可得
$$ \begin{split} \left(\nabla_{{\boldsymbol{s}}} \hat{g}_i({\boldsymbol{u}},\hat{{\boldsymbol{s}}})\right)^{\rm{T}} = &\; {\boldsymbol{e}}_i^{\rm{T}} \hat{{\boldsymbol{G}}}_{{\boldsymbol{s}}}({\boldsymbol{u}},\hat{{\boldsymbol{s}}} ) = \\ & \left.{\boldsymbol{e}}_i^{\rm{T}} \left({\boldsymbol{G}}_{{\boldsymbol{X}}} {\boldsymbol{X}}_{{\boldsymbol{s}}}+ {\boldsymbol{G}}_{{\boldsymbol{s}}}\right)\right|_{{\boldsymbol{s}} = \hat{{\boldsymbol{s}}}} \end{split} $$ (14) 其中, $ {\boldsymbol{e}}_i\in {{\bf{R}}}^m $是$ m $阶单位矩阵的第$ i $列, $ m $是$ {\boldsymbol{G}} $的维数. 将式(14)代入式(13)中, 再结合式(11)和原始问题式(10), 可以得到原始问题在一阶近似意义下的鲁棒优化问题
$$ \begin{split} & \min\limits_{{\boldsymbol{X}}, {\boldsymbol{X}}_{{\boldsymbol{s}}}, {\boldsymbol{u}}} \; \phi({\boldsymbol{X}},{\boldsymbol{u}},\hat{{\boldsymbol{s}}})\\ &\text{s.t.}\\ &{\boldsymbol{F}}({\boldsymbol{X}},{\boldsymbol{u}},\hat{{\boldsymbol{s}}}) = {\bf{0}}\\ &{\boldsymbol{F}}_{{\boldsymbol{X}}}{\boldsymbol{X}}_{{\boldsymbol{s}}}+ {\boldsymbol{F}}_{{\boldsymbol{s}}} = {\bf{0}}\\ &g_i({\boldsymbol{X}}, {\boldsymbol{u}}, \hat{{\boldsymbol{s}}})+\tau \Vert {\boldsymbol{e}}_i^{\rm{T}} \left({\boldsymbol{G}}_{{\boldsymbol{X}}} {\boldsymbol{X}}_{{\boldsymbol{s}}}+ {\boldsymbol{G}}_{{\boldsymbol{s}}}\right){\boldsymbol{A}} \Vert_1\leq 0\\ &i = 1,\cdots,m \end{split} $$ (15) 其中, $ g_i $代表$ {\boldsymbol{G}} $的第$ i $个元素, $ {\boldsymbol{F}}_{{\boldsymbol{X}}} $, $ {\boldsymbol{F}}_{{\boldsymbol{s}}} $, $ {\boldsymbol{G}}_{{\boldsymbol{X}}} $和$ {\boldsymbol{G}}_{{\boldsymbol{s}}} $均使用参数标称值$ \hat{{\boldsymbol{s}}} $计算, $ {\boldsymbol{X}}_{{\boldsymbol{s}}} $是新的待优化变量, 代表通过等式约束$ {\boldsymbol{F}}({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{s}}) = {\bf{0}} $和式(11)所导出的$ {\boldsymbol{X}} $关于$ {\boldsymbol{s}} $的偏导数.
注4. 式(15)中不等式约束的第二项可以看作是式(10)的不等式约束对应于参数变化的“安全裕量”, 它代表$ g_i $对于参数改变的敏感度. 虽然在计算中使用的是参数标称值, 但是通过在式(15)中添加安全裕量来尽可能保证实际的参数不会导致优化结果违背不等式约束.
文献[23]指出, 虽然上述鲁棒优化方法不能保证不等式约束$ {\boldsymbol{G}}({\boldsymbol{X}}({\boldsymbol{u}},{\boldsymbol{s}}),{\boldsymbol{u}},{\boldsymbol{s}})\le {\bf{0}} $对$ \forall {\boldsymbol{s}}\in {\boldsymbol{S}}_\tau $都一定成立, 但是在一定的连续性假设下, 它能保证超出约束范围的值具有一个关于$ \tau^2 $的上界.
为了叙述简便, 令
$$ {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}: = {\boldsymbol{v}} - {\boldsymbol{v}}_w = \left[ \Delta_{v_{x}}, \Delta_{ v_y}, \Delta_{ v_z}\right]^{\rm{T}} $$ $ {\boldsymbol{\Delta}}_{{\boldsymbol{v}}} $代表风和火箭的相对速度, $ \Delta_{v_{x}}, \Delta_{ v_y}, \Delta_{ v_z} $分别代表$ {\boldsymbol{\Delta}}_{{\boldsymbol{v}}} $在$ x,\; y,\; z $轴上的分量. 由于风场在$ z $轴方向风速为零, 因此$ \Delta_{v_z} = \; v_z $. 另外, 令
$$ \alpha(r_z) : = \frac{1}{2}\rho(r_z)S_{\text{ref}}C_D $$ (16) 那么式(2)中的空气阻力可以表示为
$$ {\boldsymbol{D}} = -\alpha(r_z) \Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}} \Vert_2 {\boldsymbol{\Delta}}_{{\boldsymbol{v}}} $$ 此外, 对于一个矩阵$ {\boldsymbol{P}} $, 引入符号$ {\boldsymbol{P}}(i:j,k:l) $ 代表由$ {\boldsymbol{P}} $的第$ i $行到第$ j $行、第$ k $列到第$ l $列组成的子矩阵, $ {\boldsymbol{P}}(i:j,\; :\; ) $代表由$ {\boldsymbol{P}} $的第$ i $行到第$ j $行构成的子矩阵. 将上述鲁棒优化方法应用到RO问题2中, 可以得到:
RO问题3.
$$ \begin{split} &\min\limits_{{\boldsymbol{u}}} \; J = t_f+\int_{t_0}^{t_f} \gamma\; {\boldsymbol{u}}^{\rm{T}} {\boldsymbol{u}} \; {\rm{d}}t\\ &\text{s.t.}\\ &\dot{{\boldsymbol{X}}} = {\boldsymbol{f}}({\boldsymbol{X}},{\boldsymbol{u}},{\boldsymbol{v}}_w(r_z,\hat{{\boldsymbol{s}}}))\\ & \dot{{\boldsymbol{X}}}_{{\boldsymbol{s}}} = {\boldsymbol{f}}_{{\boldsymbol{X}}} {\boldsymbol{X}}_{{\boldsymbol{s}}} +{\boldsymbol{f}}_{{\boldsymbol{s}}} \\ &{\boldsymbol{X}}(t_0) = {\boldsymbol{X}}_0,\; {\boldsymbol{X}}_{{\boldsymbol{s}}}(t_0) = {\bf{0}} \\ &\left\{\begin{aligned} &\left| r_i(t_f) - r_{f_i} \right| + \tau \Vert {\boldsymbol{r}}_{{\boldsymbol{s}}_i} (t_f) {\boldsymbol{A}}\Vert_1 \leq L_{ri} \\ &\left| v_i(t_f) - v_{f_i} \right| + \tau \Vert {\boldsymbol{v}}_{{\boldsymbol{s}}_i} (t_f) {\boldsymbol{A}}\Vert_1 \leq L_{vi} \end{aligned}\right.\\ &i = x,y,z \end{split} $$ (17) 其中, $ {\boldsymbol{X}}_{{\boldsymbol{s}}} = \left[{\boldsymbol{r}}_{{\boldsymbol{s}}}^{\rm{T}}, {\boldsymbol{v}}_{{\boldsymbol{s}}}^{\rm{T}}\right]^{\rm{T}} \in {{\bf{R}}}^{6\times d} $为新的状态变量, $ d $代表不确定参数的维度.
RO问题3中, $ {\boldsymbol{X}}_{{\boldsymbol{s}}} $的等式约束具体表达式为式(18) ~ (19), 其中, $ {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,\; :\; ) $各行表达式为式(20) ~ (22), 式(23)给出$ {\boldsymbol{f}}_{{\boldsymbol{s}}}(4:6,\; :\; ) $的表达式, 而$ \frac{\partial f_4}{\partial r_z} $, $ \frac{\partial f_5}{\partial r_z} $和$ \frac{\partial f_6}{\partial r_z} $的表达式在式(24) ~ (26)给出. 需要注意的是, 虽然按照本节给出的方法, $ {\boldsymbol{X}}_{{\boldsymbol{s}}} $的形式应该为$ {\boldsymbol{X}}_{{\boldsymbol{s}}} = \left[{\boldsymbol{r}}_{{\boldsymbol{s}}}^{\rm{T}},{\boldsymbol{v}}_{{\boldsymbol{s}}}^{\rm{T}},{\boldsymbol{m}}_{{\boldsymbol{s}}}^{\rm{T}}\right]^{\rm{T}} $, 但是容易得到在$ {\boldsymbol{X}}_{{\boldsymbol{s}}} $的等式约束中, $ {\boldsymbol{m}}_{{\boldsymbol{s}}} $所对应的动态方程全都为零. 又因为$ {\boldsymbol{m}}_{{\boldsymbol{s}}} $的初值也为零, 因此$ {\boldsymbol{m}}_{{\boldsymbol{s}}} $恒为零, 所以在考虑$ {\boldsymbol{X}}_{{\boldsymbol{s}}} $的等式约束时可以忽视$ {\boldsymbol{m}}_{{\boldsymbol{s}}} $对应的行与列, 同时对问题不产生影响. 因此在本问题中可以不失一般性地令${\boldsymbol{X}}_{{\boldsymbol{s}}} = [{\boldsymbol{r}}_{{\boldsymbol{s}}}^{\rm{T}},{\boldsymbol{v}}_{{\boldsymbol{s}}}^{\rm{T}}]^{\rm{T}}$
$$ \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}_{3\times 3}} = {\boldsymbol{v}}_{{\boldsymbol{s}}_{3\times 3}} $$ (18) $$ \begin{split} \dot{{\boldsymbol{v}}}_{{\boldsymbol{s}}_{3\times 3}} =\;& {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,\; :\; )_{3\times 6} \left[\begin{matrix} {\boldsymbol{r}}_{{\boldsymbol{s}}}\\ {\boldsymbol{v}}_{{\boldsymbol{s}}} \end{matrix}\right]_{6\times 3} +\\ &{\boldsymbol{f}}_{{\boldsymbol{s}}}(4:6,\; :\; )_{3\times 3} \end{split} $$ (19) $$ \begin{split} \frac{\partial f_4}{\partial {\boldsymbol{X}}} =\;& {\boldsymbol{f}}_{{\boldsymbol{X}}}(4,\; :\; ) = \\ &\left[ 0,0, \frac{\partial f_4}{\partial r_z},-\frac{\alpha}{m}\left(\frac{\Delta_{v_x}^2}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2} + \Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 \right),\right.\end{split} $$ $$ \left. -\frac{\alpha}{m} \frac{\Delta_{v_x} \Delta_{v_y} } {\Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2},-\frac{\alpha}{m} \frac{\Delta_{v_x} v_z} {\Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2} \right] $$ (20) $$ \begin{split} \frac{\partial f_5}{\partial {\boldsymbol{X}}} =&\; {\boldsymbol{f}}_{{\boldsymbol{X}}}(5,\; :\; ) = \\ &\left[ 0,0, \frac{\partial f_5}{\partial r_z}, -\frac{\alpha}{m}\frac{\Delta_{v_x}\Delta_{v_y}}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2},\right.\\ &\left.-\frac{\alpha}{m}\left(\frac{\Delta_{v_y}^2}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2} + \Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 \right), -\frac{\alpha}{m} \frac{\Delta_{v_y} v_z} {\Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2} \right] \end{split} $$ (21) $$ \begin{split} \frac{\partial f_6}{\partial {\boldsymbol{X}}} =&\; {\boldsymbol{f}}_{{\boldsymbol{X}}}(6,\; :\; ) =\\ &\left[ 0,0, \frac{\partial f_6}{\partial r_z}, -\frac{\alpha}{m}\frac{\Delta_{v_x}\Delta_{v_z}}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2},\right. \\ &\left.-\frac{\alpha}{m} \frac{\Delta_{v_y} \Delta_{v_z} } {\Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2}, -\frac{\alpha}{m}\left(\frac{v_z^2}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2} + \Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 \right) \right] \end{split} $$ (22) $$\begin{split} {\boldsymbol{f}}_{{\boldsymbol{s}}}(4:6,\; :\; ) =\; &\frac{\alpha}{m}\text{ diag} \left\{ \frac{\Delta_{v_x}\Delta_{v_y}}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2}, \frac{\Delta_{v_y}^2}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2}\;+\right.\\ &\left. \Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 ,\frac{\Delta_{v_y} v_z} {\Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2} \right\} \left[\begin{matrix} r_z^{*2} &r_z^* & 1 \\ r_z^{*2} &r_z^* & 1\\ r_z^{*2} &r_z^* & 1 \end{matrix}\right] \end{split} $$ (23) $$ \begin{split} \frac{\partial f_4}{\partial r_z} = \;& \frac{\partial }{\partial r_z}\left( \frac{u_x}{m}- \frac{\alpha}{m}\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2\Delta_{v_x}\right) = \\ & - \frac{1 }{m}\Bigg( \frac{\partial \alpha}{\partial r_z}\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2\Delta_{v_x} -\alpha \frac{\partial\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2}{\partial{\boldsymbol{v}}_w} \frac{\partial{\boldsymbol{v}}_w}{\partial r_z}\Delta_{v_x} -\\ &\alpha \Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 \frac{\partial \Delta_{v_x}}{\partial r_z}\Bigg) = - \frac{1 }{m} \Bigg( \frac{\partial \alpha}{\partial r_z}\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2\Delta_{v_x} -\\ &\alpha \frac{ k \Delta_{v_x}+\frac{2\hat{s}_1 r_z^*+\hat{s}_2}{\sigma}\Delta_{v_y}}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2}\Delta_{v_x} -\alpha k \Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 \Bigg)\\[-18pt] \end{split} $$ (24) $$ \begin{split} &\frac{\partial f_5}{\partial r_z} = - \frac{1 }{m} \Bigg( \frac{\partial \alpha}{\partial r_z}\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2\Delta_{v_y} -\\ &\;\;\alpha \frac{ k \Delta_{v_x}+\frac{2\hat{s}_1 r_z^*+\hat{s}_2}{\sigma}\Delta_{v_y}}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2}\Delta_{v_y} -\alpha\frac{2\hat{s}_1 r_z^*+\hat{s}_2}{\sigma} \Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 \Bigg) \end{split}$$ (25) $$ \begin{split}\frac{\partial f_6}{\partial r_z} = \;&- \frac{1 }{m} \Bigg( \frac{\partial \alpha}{\partial r_z}\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 v_z \;-\\ &\alpha \frac{ k \Delta_{v_x}+\frac{2\hat{s}_1 r_z^*+\hat{s}_2}{\sigma}\Delta_{v_y}}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2} v_z \Bigg)\end{split} $$ (26) RO问题3使用标称风场信息进行计算, 因此可以直接使用传统优化算法求解. 基于RO问题3可以得出如下单阶段RO算法.
单阶段RO算法. 输入$ {\boldsymbol{X}}_0 $, $ {\boldsymbol{r}}_f $, $ {\boldsymbol{v}}_f $, $ L_{ri} $, $ L_{vi} $, $ i = x,y,z $, 求解RO问题3, 输出最优控制、相应的状态轨迹与终端时刻.
2. 多阶段鲁棒优化算法
第1.3节得到了参数的变化导致的终端约束函数变化, 即注4所说的安全裕量, 也即式(17)中的$ \Vert {\boldsymbol{r}}_{{\boldsymbol{s}}_i} (t_f) {\boldsymbol{A}}\Vert_1 $, $ \Vert {\boldsymbol{v}}_{{\boldsymbol{s}}_i} (t_f) {\boldsymbol{A}}\Vert_1 $. 由式(13)的推导可知, 这种方法在每个时刻总会选取最差的参数(尽管实际场景中很难出现这种情况), 并且随着时间推移, 这种最差的参数对约束函数带来的变化也会逐渐累积, 最终会导致安全裕量比较大. 这就导致RO问题3出现矛盾: 一方面, 终端约束范围$ L_{ri} $和$ L_{vi} $, $ i = x,y,z $越小越好, 因为这代表着能使火箭更加精准地以给定速度落到给定地点; 另一方面, 如果$ L_{ri} $和$ L_{vi} $比安全裕量还小, 那么式(17)中的不等式约束将无法满足, 导致优化问题无解, 因此获取更多安全裕量的信息十分重要. 事实上, 本节将会定量地给出安全裕量的一个上界. 为此, 首先给出一个假设和一个引理.
假设1. 参数的维度$ d = 3 $, 此外, 由式(16)定义的$ \alpha(r_z) $, $ \alpha(r_z) $关于高度的导数$ \dfrac{\text{d}\alpha}{{\text{d}r}_z} $, $ \Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 $和质量均有界, 满足
$$ \begin{aligned} & \left|\alpha(r_z)\right|\leq \alpha_0,\; \left| \frac{\text{d} \alpha}{\text{d} r_z}\right|\leq \alpha_0'\\ & \Vert {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 \leq \bar{\Delta}_v,\; m\geq m^*\end{aligned} $$ 在假设1下, 可以证明有如下引理成立.
引理1. 若假设1成立, 在RO问题3中, $ {\boldsymbol{X}}_{{\boldsymbol{s}}} $等式约束内的$ {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,\; :\; ) $和$ {\boldsymbol{f}}_{{\boldsymbol{s}}}(4:6,\; :\; ) $满足
$$ \begin{split} &\Vert {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,1:3)\Vert_2 < M_1 \\ &\Vert {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,4:6)\Vert_2 \leq M_2\\ &\Vert {\boldsymbol{f}}_{{\boldsymbol{s}}}(4:6,\; :\; )\Vert_2 \leq M_3 \end{split} $$ 其中
$$ \begin{split} & M_1 = \frac{\bar{\Delta}_v}{m^*} \left( \alpha_0' \bar{\Delta}_v+ \alpha_0\left(\left|k\right|+\frac{2\sqrt3\left|\hat{s}_1\right|+\left|\hat{s}_2\right|}{\sigma} \right) \right) \\ & M_2 = 2 \alpha_0 \frac{\bar{\Delta}_v}{m^*} \\ & M_3 = \left(8+2\sqrt{3}\right)\alpha_0\frac{\bar{\Delta}_v}{m^*} \end{split} $$ $ k $, $ \hat{s}_1 $, $ \hat{s}_2 $, $ \sigma $均为第1.2节定义的已知参数.
证明. 对于$ \Vert {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,1:3)\Vert_2 $, 由于该矩阵的2范数代表矩阵的特征值的绝对值的最大值, 因此$\Vert {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,1:3)\Vert_2 = \left|\frac{\partial f_6}{\partial r_z}\right|$.
由式(26)和假设1, 有
$$ \begin{split} &\left|\dfrac{\partial f_6}{\partial r_z}\right| < \; \frac{1}{m} \Bigg(\alpha_0' \bar{\Delta}_v^2\;+\\ &\;\;\;\alpha_0 \left(\left|k\right|+ \frac{2\left |\hat{s}_1\right|\left|r_z^*\right|+\left|\hat{s}_2\right|}{\sigma}\right)\left|v_z \right|\Bigg)\leq\\& \;\;\;\frac{\bar{\Delta}_v}{m^*} \left(\alpha_0' \bar{\Delta}_v+\alpha_0\left( \left|k\right|+ \frac{2\sqrt3\left |\hat{s}_1\right|+\left|\hat{s}_2\right|}{\sigma}\right) \right) = M_1\end{split} $$ 其中, 第二个不等式是因为式(7). 因此
$$ \begin{array}{l} \Vert {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,1:3)\Vert_2 < M_1 \end{array} $$ 对于$ {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,4:6) $, 可以得到
$$ \begin{split} {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,4:6) = &\; -\frac{\alpha}{m}\left( \frac{{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}^{\rm{T}}}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2}+ \Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2 {\bf{1}}_{3\times 3}\right) = \\ & -\frac{\alpha}{m}\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2\left( \frac{{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}^{\rm{T}}}{\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2^2}+{\bf{1}}_{3\times 3}\right)\\[-18pt] \end{split} $$ (27) 其中, $ {\boldsymbol{\Delta}}_{{\boldsymbol{v}}} = \left[\Delta_{v_{x}},\Delta_{v_y}, v_z\right]^{\rm{T}} $. 注意到$ {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}^{\rm{T}}/\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2^2 $和$ {\boldsymbol{\Delta}}_{{\boldsymbol{v}}}^{\rm{T}}{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}/\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2^2 $有相同的非零特征值, 而后者为标量且值为1. 因此式(27)中最后一个等号右端括号内矩阵的特征值为2和1, 进而可以得到
$$ \Vert {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,4:6)\Vert_2 = \; 2\frac{\alpha}{m}\Vert{\boldsymbol{\Delta}}_{{\boldsymbol{v}}}\Vert_2\leq M_2 $$ 对于$ {\boldsymbol{f}}_{{\boldsymbol{s}}}(4:6,\; :\; ) $, 在式(23)中, 如果记
$$ {\boldsymbol{R}}_z^* : = \begin{bmatrix} r_z^{*2} &r_z^* & 1 \\ r_z^{*2} &r_z^* & 1\\ r_z^{*2} &r_z^* & 1 \end{bmatrix} $$ 那么可以得到
$$ {\boldsymbol{R}}_z^* = \begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix} \begin{bmatrix} r_z^{*2},r_z^*,1 \end{bmatrix} $$ 注意到$ {\boldsymbol{R}}_z^* $和 $ \left[r_z^{*2},r_z^*,1\right]\times\left[1,1,1\right]^{\rm{T}} $具有相同的非零特征值, 而该矩阵唯一的特征值为$ 1+r_z^*\;+ r_z^{*2} $, 因此
$$ \Vert{\boldsymbol{R}}_z^*\Vert_2 = \left| 1+r_z^*+r_z^{*2}\right|\leq 1+\sqrt3 +3 = 4+\sqrt3 $$ 易知式(23)中的对角矩阵的2范数小于等于$ 2\bar{\Delta}_v $, 所以可以得到
$$ \Vert {\boldsymbol{f}}_{{\boldsymbol{s}}}(4:6,\; :\; ) \Vert_2 \leq \frac{\alpha_0}{m^*}2\bar{\Delta}_v \left( 4+\sqrt3\right) = M_3 $$ □ 引理1说明式(17)中的$ {\boldsymbol{f}}_{{\boldsymbol{X}}} $, $ {\boldsymbol{f}}_{{\boldsymbol{s}}} $都是有界的, 因此可以通过$ \dot{{\boldsymbol{X}}}_{{\boldsymbol{s}}} = {\boldsymbol{f}}_{{\boldsymbol{X}}} {\boldsymbol{X}}_{{\boldsymbol{s}}} +{\boldsymbol{f}}_{{\boldsymbol{s}}} $求出$ {\boldsymbol{r}}_{{\boldsymbol{s}}_i} $, $ {\boldsymbol{v}}_{{\boldsymbol{s}}_i} $, $ i = x, y,z $的上界, 进而求出式(17)中安全裕量的界. 具体而言, 给出如下定理.
定理1. 若$ {\boldsymbol{X}}_{{\boldsymbol{s}}} $满足
$$ \dot{{\boldsymbol{X}}}_{{\boldsymbol{s}}} = {\boldsymbol{f}}_{{\boldsymbol{X}}} {\boldsymbol{X}}_{{\boldsymbol{s}}} +{\boldsymbol{f}}_{{\boldsymbol{s}}}, \; {\boldsymbol{X}}_{{\boldsymbol{s}}}(t_0) = {\boldsymbol{0}} $$ 那么对$ \forall t_f \geq 0 $, 式(17)中的$ \Vert {\boldsymbol{r}}_{{\boldsymbol{s}}_i} (t_f){\boldsymbol{A}}\Vert_1 $和$ \Vert {\boldsymbol{v}}_{{\boldsymbol{s}}_i} (t_f){\boldsymbol{A}}\Vert_1 $有界, 满足
$$ \Vert {\boldsymbol{r}}_{{\boldsymbol{s}}_i} {\boldsymbol{A}}\Vert_1 < 3\sqrt3 \Vert{\boldsymbol{A}}\Vert_1 M_3 \sqrt{\frac{t_f}{ C}}\left(\text{e}^{Ct_f} -1 \right)^{\frac{1}{2}} $$ (28) $$ \Vert {\boldsymbol{v}}_{{\boldsymbol{s}}_i}{\boldsymbol{A}} \Vert_1 < 3\sqrt {3\; t_f} \Vert{\boldsymbol{A}}\Vert_1 M_3\text{e}^{\frac{1}{2}Ct_f} $$ (29) 其中, $ i = x,y,z $, $ C = 1+ \left(t_f+1\right) M_1 + 2M_2 $.
证明. 定理中的结果需要证明$ {\boldsymbol{r}}_{{\boldsymbol{s}}} $和$ {\boldsymbol{v}}_{{\boldsymbol{s}}} $所有行向量的1范数有界, 因此只需证明$ \Vert{\boldsymbol{r}}_{{\boldsymbol{s}}} \Vert_\infty $和$ \Vert{\boldsymbol{v}}_{{\boldsymbol{s}}} \Vert_\infty $有界. 为了证明$ \Vert{\boldsymbol{r}}_{{\boldsymbol{s}}} \Vert_\infty $和$ \Vert{\boldsymbol{v}}_{{\boldsymbol{s}}} \Vert_\infty $有界, 首先证明它们的Frobenius范数(F范数)有界, 再通过矩阵范数之间的关系(对任意$ n $阶方阵$ {\boldsymbol{A}} $, 有$ \Vert {\boldsymbol{A}}\Vert_\infty\leq \sqrt n\Vert {\boldsymbol{A}} \Vert_2 \leq \sqrt n\Vert {\boldsymbol{A}} \Vert_{\text{F}} $)得到$ \Vert {\boldsymbol{r}}_{{\boldsymbol{s}}} \Vert_\infty $和$ \Vert{\boldsymbol{v}}_{{\boldsymbol{s}}} \Vert_\infty $的界. 为了简便, 在本证明中, 令$ \left|\; \cdot\; \right| $代表向量或者矩阵的2范数, 对于标量则代表取绝对值, 并记$ \; {\boldsymbol{r}}_{{\boldsymbol{s}}^j} = {\boldsymbol{r}}_{{\boldsymbol{s}}} (\; :\; ,j), {\boldsymbol{v}}_{{\boldsymbol{s}}^j} = {\boldsymbol{v}}_{{\boldsymbol{s}}} (\; :\; ,j),\; j = 1,2,3 $代表$ {\boldsymbol{r}}_{{\boldsymbol{s}}} $和$ {\boldsymbol{v}}_{{\boldsymbol{s}}} $的各个列向量, 注意和$ \; {\boldsymbol{r}}_{{\boldsymbol{s}}_i} $, $ {\boldsymbol{v}}_{{\boldsymbol{s}}_i} $区分, $ i = x,y,z $. 它们代表$ {\boldsymbol{r}}_{{\boldsymbol{s}}} $和$ {\boldsymbol{v}}_{{\boldsymbol{s}}} $的各个行向量. 将式(19)拆为三个列向量对应的微分方程可得
$$ \begin{split} \dot{{\boldsymbol{v}}}_{{\boldsymbol{s}}^j} = &\; {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,1:3){\boldsymbol{r}}_{{\boldsymbol{s}}^j}+ {\boldsymbol{f}}_{{\boldsymbol{X}}}(4:6,4:6){\boldsymbol{v}}_{{\boldsymbol{s}}^j}\;+\\ & {\boldsymbol{f}}_{{\boldsymbol{s}}}(4:6,\; :\; ) \end{split} $$ 两端同时取2范数并由引理1可得
$$ \left| \dot{{\boldsymbol{v}}}_{{\boldsymbol{s}}^j} \right| < M_1\left|{\boldsymbol{r}}_{{\boldsymbol{s}}^j}\right|+ M_2\left|{\boldsymbol{v}}_{{\boldsymbol{s}}^j}\right|+ M_3 $$ 或者等价地
$$ \left|\ddot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right| < M_1\left|{\boldsymbol{r}}_{{\boldsymbol{s}}^j}\right|+ M_2\left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|+ M_3 $$ 对$ \left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|^2 $关于时间求导可得
$$ \frac{{\rm{d}}}{{\rm{d}}t}\left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|^2 = \frac{{\rm{d}}}{{\rm{d}}t}\langle \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}, \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j} \rangle = 2\langle \ddot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}, \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j} \rangle $$ $ \langle\cdot, \cdot\rangle $代表内积. 对上式两端在$ \left[0,t\right] $上积分, 注意到$ \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}(0) = {\bf{0}} $, 可得
$$ \begin{split} \left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|^2 = & \; 2 \int_0^t \langle \ddot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}, \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j} \rangle {\rm{d}}l\leq 2 \int_0^t \left| \ddot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|\left| \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j} \right|{\rm{d}}l<\\ &\; 2 \int_0^t \left( M_1\left|{\boldsymbol{r}}_{{\boldsymbol{s}}^j}\right|+ M_2\left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|+ M_3 \right) \left| \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j} \right|{\rm{d}}l\leq \\ &\; M_1 \int_0^t \left(\left|{\boldsymbol{r}}_{{\boldsymbol{s}}^j}\right|^2 + \left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|^2\right){\rm{d}}l\;+\\ & 2M_2 \int_0^t \left| \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|^2 {\rm{d}}l + M_3^2t+\int_0^t \left| \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|^2 {\rm{d}}l = \\ &\; M_1\int_0^t \left|{\boldsymbol{r}}_{{\boldsymbol{s}}^j}\right|^2{\rm{d}}l\;+ \\ & \left(1+ M_1+2M_2 \right)\int_0^t \left| \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|^2 {\rm{d}}l+M_3^2t\\[-17pt] \end{split} $$ (30) 对于上式中的$ \left|{\boldsymbol{r}}_{{\boldsymbol{s}}^j}\right|^2 $, 有
$$ \left|{\boldsymbol{r}}_{{\boldsymbol{s}}^j}\right|^2 = \; \left| \int_0^l \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}(\xi){\rm{d}}\xi \right|^2 \leq\; \int_0^l \left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}(\xi)\right|^2{\rm{d}}\xi $$ 因此
$$ \begin{split} \int_0^t \left|{\boldsymbol{r}}_{{\boldsymbol{s}}^j}(l)\right|^2{\rm{d}}l \leq&\; \int_0^t \int_0^l \left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}(\xi)\right|^2{\rm{d}}\xi {\rm{d}}l\leq\\ &\; \int_0^t \int_0^t \left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}(\xi)\right|^2{\rm{d}}\xi {\rm{d}}l\leq\\ &\; t \int_0^t \left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}(\xi)\right|^2{\rm{d}}\xi\leq\\ &\; t_f \int_0^t \left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}(\xi)\right|^2{\rm{d}}\xi \end{split} $$ 将上式代入式(30)中, 得到
$$ \begin{split} \left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|^2 < &\; \left(1+ \left(t_f+1\right) M_1+2M_2 \right)\int_0^t \left| \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|^2 {\rm{d}}l\;+\\ &M_3^2t \leq C \int_0^t \left| \dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}\right|^2 {\rm{d}}l+M_3^2t_f \end{split} $$ 由Gronwall不等式可得
$$ \left|{\dot{{\boldsymbol{r}}}}_{{\boldsymbol{s}}^j}(t)\right|^2 < M_3^2 t_f \text{e}^{Ct} $$ (31) 也就是
$$ \left|{\boldsymbol{v}}_{{\boldsymbol{s}}^j}(t)\right| < M_3 \sqrt{t_f } \text{e}^{\frac{1}{2}Ct} $$ 回顾$ {\boldsymbol{v}}_{{\boldsymbol{s}}^j} $代表$ {\boldsymbol{v}}_{{\boldsymbol{s}}} $的列向量, 由F范数的定义可知
$$ \Vert {\boldsymbol{v}}_{{\boldsymbol{s}}}(t)\Vert_\text{F} \leq \sum\limits_{j = 1}^3 \left|{\boldsymbol{v}}_{{\boldsymbol{s}}^j}(t)\right|\leq 3 M_3 \sqrt{t_f } \text{e}^{\frac{1}{2}Ct} $$ (32) 又因为矩阵的2范数总是小于矩阵的F范数, 所以
$$ \begin{split} \Vert {\boldsymbol{v}}_{{\boldsymbol{s}}}(t)\Vert_\infty \leq &\; \sqrt 3 \Vert {\boldsymbol{v}}_{{\boldsymbol{s}}}(t)\Vert_2 \leq\\ & \sqrt 3 \Vert {\boldsymbol{v}}_{{\boldsymbol{s}}}(t)\Vert_\text{F}\leq\\ & 3\sqrt {3\; t_f} M_3\text{e}^{\frac{1}{2}Ct} \end{split} $$ (33) 由式(33)即可得知定理1的式(29)成立. 对于$ {\boldsymbol{r}}_{{\boldsymbol{s}}} $, 由式(31)可得
$$ \begin{split}& \left| {\boldsymbol{r}}_{{\boldsymbol{s}}^j}(t)\right|^2 \leq \int_0^t \left|\dot{{\boldsymbol{r}}}_{{\boldsymbol{s}}^j}(l)\right|^2 {\rm{d}}l\leq \int_0^t M_3^2 t_f \text{e}^{Cl}{\rm{d}}l = \\ &\quad\qquad\qquad\frac{M_3^2t_f}{C}\left(\text{e}^{Ct} -1 \right) \\ &\left| {\boldsymbol{r}}_{{\boldsymbol{s}}^j}(t)\right| \leq M_3 \sqrt{\frac{t_f}{ C}}\left(\text{e}^{Ct} -1 \right)^{\frac{1}{2}} \end{split} $$ 类似于式(32)和式(33)的推导, 有
$$ \begin{split} \Vert {\boldsymbol{r}}_{{\boldsymbol{s}}}(t)\Vert_{\rm{F}} \leq &\; \sum\limits_{j = 1}^3 \left|{\boldsymbol{r}}_{{\boldsymbol{s}}^j}(t)\right|\leq 3 M_3 \sqrt{\frac{t_f}{ C}}\left(\text{e}^{Ct} -1 \right)^{\frac{1}{2}}\\ \Vert {\boldsymbol{r}}_{{\boldsymbol{s}}}(t)\Vert_\infty \leq &\; \sqrt 3 \Vert {\boldsymbol{r}}_{{\boldsymbol{s}}}(t)\Vert_\text{F}\leq\\ & 3\sqrt3 M_3 \sqrt{\frac{t_f}{ C}}\left(\text{e}^{Ct} -1 \right)^{\frac{1}{2}} \\[-17pt]\end{split} $$ (34) 由式(34)即证定理1的式(28)成立.
□ 从定理1可以看出, 式(28) ~ (29)中给出的界是随时间指数增长的, 对于一些需要较长时间才能完成的优化场景, 显然会导致必须选取较大的终端约束范围优化问题才能有解, 但是大的终端约束范围又会降低落点的精度. 对于优化时间较长的场景, 为了通过缩短优化时长使问题具有可行性, 本文提出一种多阶段RO算法. 其主要思想是通过将整个优化进程拆分为多个串联的子阶段, 在每个子阶段中分别求解RO子问题, 缩短优化时长, 进而避免优化时间较长致使RO问题3无解.
具体而言, 在如图1所示的$ K $阶段RO算法流程图中, 选取每个阶段的预估终端时刻$ t^1 , \; t^2 , \cdots , t^K $, 使得每个阶段的预估时长为标称轨迹时长的$ 1/K $. 将标称轨迹在$ t^1 $时刻的状态$ {\boldsymbol{r}}_{\text{NO}}(t^1) $、$ {\boldsymbol{v}}_{\text{NO}}(t^1) $作为第一阶段的目标终值, 将初始状态$ {\boldsymbol{X}}_0 $作为第一阶段的初值, 使用单阶段RO算法进行优化, 将最优解$ {\boldsymbol{u}}^1_{\text{RO}} $输入到火箭中, 得到火箭在第一阶段末的位置$ {\boldsymbol{X}}^1_{\text{RO}}(t^1_{f_{\text{RO}}} ) $. 在第二阶段中, 则以$ {\boldsymbol{X}}_{\text{RO}}^1(t^1_{f_{\text{RO}}}) $为第二阶段的初值, 以$ {\boldsymbol{r}}_{\text{NO}}(t^2) $, $ {\boldsymbol{v}}_{\text{NO}}(t^2) $作为第二阶段的目标终值, 同样使用单阶段RO算法求解, 如此进行下去, 直到第$ K $阶段结束. 需要注意的是, 最后一个阶段的目标终值满足
$$ {\boldsymbol{r}}_{\text{NO}}(t^K) = {\boldsymbol{r}}_f,\; {\boldsymbol{v}}_{\text{NO}}(t^K) = {\boldsymbol{v}}_f $$ 通过该多阶段RO算法将会得到一组控制序列$ {\boldsymbol{u}}_{\text{RO}}^1 $, $ {\boldsymbol{u}}_{\text{RO}}^2 , \cdots , {\boldsymbol{u}}_{\text{RO}}^K $. 这样, 通过将原始的优化问题分为多个串联的子优化问题进行求解, 使每个阶段的优化时长得以缩短, 进而使算法可以求解具有更小终端约束的优化问题. 多阶段RO算法的具体实施流程如算法1所示.
算法1. 多阶段RO算法
输入. 初始状态$ {\boldsymbol{X}}_0 $, 初始时刻$ t_0 = 0 $, 期望终端状态$ {\boldsymbol{r}}_f $, $ {\boldsymbol{v}}_f $, 每个子阶段的终端约束范围, $ L_{ri} $, $ L_{vi} $, $ i = x,y,z $, 阶段总数$ K $.
步骤1. 将$ {\boldsymbol{X}}_0 $, $ {\boldsymbol{r}}_f $, $ {\boldsymbol{v}}_f $输入NO算法(见第1.1节)进行求解, 得到一条标称风场下的标称轨迹 $ {\boldsymbol{X}}_{\text{NO}} $和其优化时长$ t_{f_{\text{NO}}} $. 计算各个阶段的预估终端时刻$ t^1,t^2,\cdots,t^K $, 使得$ t^{k+1}-t^k = t_{f_{\text{NO}}}/K $, $ k = 1,2,\cdots,K-1 $. 令$ k = 1 $.
步骤2. 将$ {\boldsymbol{X}}_0 $, $ {\boldsymbol{r}}_{\text{NO}}(t^k) $, $ {\boldsymbol{v}}_{\text{NO}}(t^k) $, $ L_{ri} $, $ L_{vi} $ 输入单阶段RO算法(见第1.3节)中进行求解, 得到第$ k $阶段的控制$ {\boldsymbol{u}}_{\text{RO}}^k $和终端时刻$ t_{f_{\text{RO}}}^k $.
步骤3. 将$ {\boldsymbol{u}}_{\text{RO}}^k $输入至火箭中, 得到第$ k $阶段末实际状态$ {\boldsymbol{X}}^k_{\text{RO}}(t^k_{f_{\text{RO}}} ) $. 令$ {\boldsymbol{X}}_0 = {\boldsymbol{X}}^k_{\text{RO}}(t^k_{f_{\text{RO}}} ) $.
步骤4. 判断: $ k = K? $否: $ k = k+1 $, 返回步骤2; 是: 结束算法.
注5. 关于阶段数的选择, 在实际场景中应综合考虑多方面因素. 一方面, 阶段数量越多, 每个阶段的优化时长就越短, 根据本节的分析, 就可以求解更小终端约束范围的优化问题, 进而使最终状态更加精确. 另一方面, 阶段数量进一步增加也会带来一些问题. 首先, 阶段数增多会使每个阶段的优化时长缩短, 缩短优化时长是为了缩小约束范围$ L_{ri} $, $ L_{vi} $, $ i = x,y,z $. 但是从定理1可以看出, 当优化时长$ t_f $减小并趋于零, 式(17)中的安全裕量同样会减小并趋于零, 若$ L_{ri} $, $ L_{vi} $同样减小并趋于零, 此时RO问题3将近似等价于NO问题, 这表明当优化时间过短时, RO算法相比于NO算法便不再具有优势. 其次, 当阶段数增加时, 整个优化进程需要求解的优化问题也会增加, 这对计算性能有着更大的需求. 因此, 在实际场景中, 设计算法时应综合考虑火箭初始状态、能够容许的终端约束范围大小、计算能力等条件来选择阶段数. 本节算法将每个阶段的预估优化时长取为相等, 正是考虑到每个阶段的优化时长不能过长、也不能过短的因素.
3. 仿真结果与分析
首先说明仿真所使用的风场模型. 本文通过实际测量出的6个时刻风场数据拟合出6个风速与高度之间的函数关系, 将这些拟合函数各个参数对应的最大、最小值的平均值作为参数标称值. $ x $轴方向的风场参数选为拟合出的标称参数. 对于$ y $轴方向的不确定风场, 将各个拟合参数对应的最大、最小值分别作为该参数取值范围的上、下界, 拟合出的$ y $轴风场与高度之间的函数关系如图2所示. 如果令风场序号0代表标称风场, 风场序号1到6代表实际风场, 那么各个实际风场的参数与标称参数在无穷范数意义下的距离如图3所示. 可以看出在6个实际风场中, 风场1、2、6与标称参数差距较大, 风场4、5与标称参数差距较小. 因此可以认为, 相比于标称风场而言, 风场1、2、6是较差的风场, 风场4、5是较好的风场.
图 2 $ y$轴的六个实际风场和一个标称风场Fig. 2 Six actual wind fields and one nominal wind field on the y-axis
图 3 实际风场参数与标称风场参数的距离Fig. 3 Distances between parameters of actual wind fields and parameters of nominal wind field本节的仿真参数如表1所示, 图4则给出在表1的仿真参数条件下随机生成的一组风场, 可以看出本节的参数条件较好地覆盖了6个实际风场.
表 1 仿真参数值Table 1 Simulation parameter values参数 取值 初始状态 ${\boldsymbol{r}}_0=\left[2~968~{\rm{m}},263~{\rm{m}},4~326~{\rm{m}}\right]^{\rm{T}} $, ${\boldsymbol{v}}_0 = \left[-295~{\rm{m/s}},27~{\rm{m/s}},-296~{\rm{m/s}}\right]^{\rm{T}}$, $m_0=48~185$ kg 期望落点状态 ${\boldsymbol{r}}_f =\left[0~{\rm{m}},0~{\rm{m}},0~{\rm{m}}\right]^{\rm{T}}$, $ {\boldsymbol{v}}_f =\left[0~{\rm{m/s}},0~{\rm{m/s}},0~{\rm{m/s}}\right]^{\rm{T}}$ 火箭燃料消耗常数 $\kappa=2~975$ 火箭参考面积 $S_{{{\rm{ref}}}}=8.814\ {\rm{m}}^2 $ 空气动力学阻力系数 $C_D=4.5$ 标称风场参数 $\hat{{\boldsymbol{s}}}=\left[-0.073,4.460,5.230\right]^{\rm{T}},\ \mu=2~234,\ \sigma=1~635$, $ k=2.16\times{10}^{-3}, b=0.35$ ${\boldsymbol{S}}_\tau $中的参数 $\tau=1,\ {\boldsymbol{D}}=\text{ diag}\left\{0.83,\ 1.20,\ 2.66\right\}$ 终端约束范围 单阶段RO算法${\boldsymbol{L}}=\left[L_{rx},L_{ry},L_{rz},L_{vx},L_{vy},L_{vz}\right]=\left[7,35,7,35,35,35\right]$
多阶段RO算法${\boldsymbol{L}}=\left[4,20,4,4,4,4 \right]$多阶段RO算法的阶段数 2 
此外, 本节使用Matlab软件包GPOPS来求解最优控制问题. GPOPS是一个基于伪谱法和序列二次规划方法来求解最优控制问题的算法, 具有求解速度快、精度高等特点.
注6. 在表1中, 单阶段RO算法的终端约束范围大于多阶段RO算法, 这是因为当单阶段RO算法和多阶段RO算法取相同的终端约束范围参数时, 单阶段RO算法将由于约束范围太小而无法求解.
3.1 单阶段RO算法和多阶段RO算法的对比
在本节, 将通过对比单阶段RO算法和多阶段RO算法来显示多阶段RO算法在最终落点状态精度上的优越性. 单阶段RO算法和多阶段RO算法在不同风场下最终位置和期望位置的误差如图5所示, 最终速度和期望速度的误差如图6所示, 它们在各个实际风场下的轨迹如图7所示, 控制–时间关系则如图8、图9所示. 从图5和图6可以看出, 多阶段RO算法无论是位置还是速度, 最终状态的精度均高于单阶段RO算法. 并且随着风场的变化, 多阶段RO算法的最终位置误差波动也明显小于单价段RO算法. 这显示多阶段RO算法相比于单阶段RO算法在最终落点状态精度上的优越性.
图 5 单阶段RO算法和多阶段RO算法在不同风场下的最终位置误差Fig. 5 The terminal position errors of single-stage RO algorithm and multi-stage RO algorithm under different wind fields
图 6 单阶段RO算法和多阶段RO算法在不同风场下的最终速度误差Fig. 6 The terminal velocity errors of single-stage RO algorithm and multi-stage RO algorithm under different wind fields
图 7 单阶段RO算法和多阶段RO算法在各个实际风场下的轨迹Fig. 7 The trajectories of single-stage RO algorithm and multi-stage RO algorithm under each actual wind field此外, 单阶段RO算法和多阶段RO算法的运行时间如表2所示. 从表2可以看出, 由于多阶段RO算法需要求解更多优化问题, 所以具有更高的计算代价, 这也表明在设计算法时阶段数不宜过多.
表 2 单阶段RO算法和多阶段RO算法的运行时间Table 2 Running time of single-stage RO algorithm and multi-stage RO algorithm算法 运行时间(s) 单阶段RO算法 6.61 多阶段RO算法 13.14 3.2 RO算法的约束满足情况以及与NO算法的对比
在本节, 将通过单阶段RO算法与多阶段RO算法在不同风场下的终端约束满足情况以及单阶段RO算法与NO算法的对比, 来显示RO算法的有效性. 单阶段RO算法在各个风场下最终位置、最终速度和相应的约束范围如图10、图11所示, 其中, 虚线代表该方向对应的约束范围, 圆圈代表对应的位置或速度. 多阶段RO算法在各个风场下最终位置、最终速度和相应的约束范围如图12、图13所示, 其中, 虚线代表该方向对应的约束范围, 叉号代表对应的位置或速度. 从图10 ~ 图13可以看出, 无论是单阶段RO算法还是多阶段RO算法, 它们在不同的风场下均未违背对应的终端状态约束.
图 10 单阶段RO算法在不同风场下的最终位置与约束范围Fig. 10 The terminal position and its constraint range of single-stage RO algorithm under different wind fields
图 11 单阶段RO算法在不同风场下的最终速度与约束范围Fig. 11 The terminal velocity and its constraint range of single-stage RO algorithm under different wind fields
图 12 多阶段RO算法在不同风场下的最终位置与约束范围Fig. 12 The terminal position and its constraint range of multi-stage RO algorithm under different wind fields
图 13 多阶段RO算法在不同风场下的最终速度与约束范围Fig. 13 The terminal velocity and its constraint range of multi-stage RO algorithm under different wind fields图14则对比在不同风场下单阶段RO算法与NO算法的最终位置误差. 从图14可以看出, 除了标称风场(风场0), 单阶段RO算法在其他风场下的最终位置精度都高于NO算法, 并且在较差的风场条件(风场1、2、6)下, 两算法精度差距较大. 这表明NO算法在面对未知的风场时, 受风场变化的影响较大, 而RO算法则对变化的风场具有更好的鲁棒性. 上述讨论说明了RO算法的有效性.
图 14 单阶段RO算法与NO算法在不同风场下的最终位置误差Fig. 14 The terminal position errors of single-stage RO algorithm and NO algorithm under different wind fields4. 总结
本文考虑处于不确定风场中可回收火箭动力下降的鲁棒制导问题, 主要工作包括:
1)建立一个包含不确定性的风场模型用来刻画未知风场, 在该模型的基础上给出火箭动力下降的鲁棒优化问题. 进一步, 使用一种对不等式约束采取一阶近似并将一阶项作为安全裕量加入约束中的鲁棒优化方法, 得到一个易于求解的单阶段鲁棒优化算法.
2)从理论上定量地给出安全裕量的上界, 借助该上界提出一种多阶段鲁棒优化方法, 该方法相比于单阶段鲁棒优化算法, 可以求解具有更小的终端约束范围的优化问题.
3)通过仿真对比分析单阶段、多阶段鲁棒优化算法和标称风场下的优化算法, 结果表明提出的多阶段鲁棒优化方法在落点精度较高的同时, 对不同的风场具有较强的鲁棒性.
-
图 2 $ y$轴的六个实际风场和一个标称风场
Fig. 2 Six actual wind fields and one nominal wind field on the y-axis
图 3 实际风场参数与标称风场参数的距离
Fig. 3 Distances between parameters of actual wind fields and parameters of nominal wind field
图 5 单阶段RO算法和多阶段RO算法在不同风场下的最终位置误差
Fig. 5 The terminal position errors of single-stage RO algorithm and multi-stage RO algorithm under different wind fields
图 6 单阶段RO算法和多阶段RO算法在不同风场下的最终速度误差
Fig. 6 The terminal velocity errors of single-stage RO algorithm and multi-stage RO algorithm under different wind fields
图 7 单阶段RO算法和多阶段RO算法在各个实际风场下的轨迹
Fig. 7 The trajectories of single-stage RO algorithm and multi-stage RO algorithm under each actual wind field
图 10 单阶段RO算法在不同风场下的最终位置与约束范围
Fig. 10 The terminal position and its constraint range of single-stage RO algorithm under different wind fields
图 11 单阶段RO算法在不同风场下的最终速度与约束范围
Fig. 11 The terminal velocity and its constraint range of single-stage RO algorithm under different wind fields
图 12 多阶段RO算法在不同风场下的最终位置与约束范围
Fig. 12 The terminal position and its constraint range of multi-stage RO algorithm under different wind fields
图 13 多阶段RO算法在不同风场下的最终速度与约束范围
Fig. 13 The terminal velocity and its constraint range of multi-stage RO algorithm under different wind fields
图 14 单阶段RO算法与NO算法在不同风场下的最终位置误差
Fig. 14 The terminal position errors of single-stage RO algorithm and NO algorithm under different wind fields
表 1 仿真参数值
Table 1 Simulation parameter values
参数 取值 初始状态 ${\boldsymbol{r}}_0=\left[2~968~{\rm{m}},263~{\rm{m}},4~326~{\rm{m}}\right]^{\rm{T}} $, ${\boldsymbol{v}}_0 = \left[-295~{\rm{m/s}},27~{\rm{m/s}},-296~{\rm{m/s}}\right]^{\rm{T}}$, $m_0=48~185$ kg 期望落点状态 ${\boldsymbol{r}}_f =\left[0~{\rm{m}},0~{\rm{m}},0~{\rm{m}}\right]^{\rm{T}}$, $ {\boldsymbol{v}}_f =\left[0~{\rm{m/s}},0~{\rm{m/s}},0~{\rm{m/s}}\right]^{\rm{T}}$ 火箭燃料消耗常数 $\kappa=2~975$ 火箭参考面积 $S_{{{\rm{ref}}}}=8.814\ {\rm{m}}^2 $ 空气动力学阻力系数 $C_D=4.5$ 标称风场参数 $\hat{{\boldsymbol{s}}}=\left[-0.073,4.460,5.230\right]^{\rm{T}},\ \mu=2~234,\ \sigma=1~635$, $ k=2.16\times{10}^{-3}, b=0.35$ ${\boldsymbol{S}}_\tau $中的参数 $\tau=1,\ {\boldsymbol{D}}=\text{ diag}\left\{0.83,\ 1.20,\ 2.66\right\}$ 终端约束范围 单阶段RO算法${\boldsymbol{L}}=\left[L_{rx},L_{ry},L_{rz},L_{vx},L_{vy},L_{vz}\right]=\left[7,35,7,35,35,35\right]$
多阶段RO算法${\boldsymbol{L}}=\left[4,20,4,4,4,4 \right]$多阶段RO算法的阶段数 2 
下载: 导出CSV
表 2 单阶段RO算法和多阶段RO算法的运行时间
Table 2 Running time of single-stage RO algorithm and multi-stage RO algorithm
算法 运行时间(s) 单阶段RO算法 6.61 多阶段RO算法 13.14
下载: 导出CSV
-
[1] McCurdy D, Roche J. Structural sizing of a horizontal take-off launch vehicle with an air collection and enrichment system. In: Proceedings of the 40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. Fort Lauderdale, USA: AIAA, 2004. Article No. 3394 [2] Reynolds T, Malyuta D, Mesbahi M, Acikmese B, Carson J M. Funnel synthesis for the 6-DOF powered descent guidance problem. In: Proceedings of the AIAA Scitech 2021 Forum. Virtual Event: AIAA, 2021. Article No. 0504 [3] Scharf D P, Acikmese B, Dueri D, Benito J, Casoliva J. Implementation and experimental demonstration of onboard powered-descent guidance. Journal of Guidance, Control, and Dynamics, 2017, 40(2): 213-229 doi: 10.2514/1.G000399 [4] Klumpp A R. Apollo lunar descent guidance. Automatica, 1974, 10(2): 133-146 doi: 10.1016/0005-1098(74)90019-3 [5] McInnes C R. Direct adaptive control for gravity-turn descent. Journal of Guidance, Control, and Dynamics, 1999, 22(2): 373-375 doi: 10.2514/2.4392 [6] Zhou L Y, Xia Y Q. Improved ZEM/ZEV feedback guidance for Mars powered descent phase. Advances in Space Research, 2014, 54(11): 2446-2455 doi: 10.1016/j.asr.2014.08.011 [7] Liu Y H, Wang H L, Fan J X, Wu J F, Wu T C. Control-oriented UAV highly feasible trajectory planning: A deep learning method. Aerospace Science and Technology, 2021, 110: Article No.106435 doi: 10.1016/j.ast.2020.106435 [8] You S X, Wan C H, Dai R, Rea J R. Learning-based onboard guidance for fuel-optimal powered descent. Journal of Guidance, Control, and Dynamics, 2021, 44(3): 601-613 doi: 10.2514/1.G004928 [9] 周宏宇, 王小刚, 单永志, 赵亚丽, 崔乃刚. 基于改进粒子群算法的飞行器协同轨迹规划. 自动化学报, 2022, 48(11): 2670-2676Zhou Hong-Yu, Wang Xiao-Gang, Shan Yong-Zhi, Zhao Ya-li, Cui Nai-Gang. Synergistic path planning for multiple vehicles based on an improved particle swarm optimization method. Acta Automatica Sinica, 2022, 48(11): 2670-2676 [10] Wei H L, Shi Y. MPC-based motion planning and control enables smarter and safer autonomous marine vehicles: Perspectives and a tutorial survey. IEEE/CAA Journal of Automatica Sinica, 2023, 10(1): 8-24 doi: 10.1109/JAS.2022.106016 [11] Lawden D F. Optimal Trajectories for Space Navigation. London: Butter Worths, 1963. 42−51 [12] Lu P. Propellant-optimal powered descent guidance. Journal of Guidance, Control, and Dynamics, 2018, 41(4): 813-826 doi: 10.2514/1.G003243 [13] Acikmese B, Ploen S R. Convex programming approach to powered descent guidance for Mars landing. Journal of Guidance, Control, and Dynamics, 2007, 30(5): 1353-1366 doi: 10.2514/1.27553 [14] Acikmese B, Blackmore L. Lossless convexification of a class of optimal control problems with non-convex control constraints. Automatica, 2011, 47(2): 341-347 doi: 10.1016/j.automatica.2010.10.037 [15] Harris M W, Acikmese B. Lossless convexification of non-convex optimal control problems for state constrained linear systems. Automatica, 2014, 50(9): 2304-2311 doi: 10.1016/j.automatica.2014.06.008 [16] Liu X F, Lu P. Solving nonconvex optimal control problems by convex optimization. Journal of Guidance, Control, and Dynamics, 2014, 37(3): 750-765 doi: 10.2514/1.62110 [17] Liu X, Li S, Xin M. Mars entry trajectory planning with range discretization and successive convexification. Journal of Guidance, Control, and Dynamics, 2022, 45(4): 755-763 doi: 10.2514/1.G006237 [18] Szmuk M, Reynolds T P, Acikmese B. Successive convexification for real-time six-degree-of-freedom powered descent guidance with state-triggered constraints. Journal of Guidance, Control, and Dynamics, 2020, 43(8): 1399-1413 doi: 10.2514/1.G004549 [19] 王祝, 徐广通, 龙腾. 基于定制内点法的多无人机协同轨迹规划. 自动化学报, 2023, 49(11): 2374-2385Wang Zhu, Xu Guang-Tong, Long Teng. Customized interior-point method for cooperative trajectory planning of unmanned aerial vehicles. Acta Automatica Sinica, 2023, 49(11): 2374-2385 [20] Boscariol P, Richiedei D. Robust point-to-point trajectory planning for nonlinear underactuated systems: Theory and experimental assessment. Robotics and Computer-Integrated Manufacturing, 2018, 50: 256-265 doi: 10.1016/j.rcim.2017.10.001 [21] Chen X L, Zhang R, Xue W C, Li H F. Robust neighboring optimal guidance for endoatmospheric powered descent under uncertain wind fields. In: Proceedings of the 34th Chinese Control and Decision Conference (CCDC). Hefei, China: IEEE, 2022. 1491−1496 [22] Brault P, Delamare Q, Giordano P R. Robust trajectory planning with parametric uncertainties. In: Proceedings of the IEEE International Conference on Robotics and Automation (ICRA). Xi＇an, China: IEEE, 2021. 11095−11101 [23] Zhang Y. General robust-optimization formulation for nonlinear programming. Journal of Optimization Theory and Applications, 2007, 132: 111-124 doi: 10.1007/s10957-006-9082-z [24] Prabhakar A, Fisher J, Bhattacharya R. Polynomial chaos-Based analysis of probabilistic uncertainty in hypersonic flight dynamics. Journal of Guidance Control and Dynamics, 2010, 33(1): 222-234 doi: 10.2514/1.41551 [25] Wang F G, Yang S X, Xiong F F, Lin Q Z, Song J M. Robust trajectory optimization using polynomial chaos and convex optimization. Aerospace Science and Technology, 2019, 92: 314-325 doi: 10.1016/j.ast.2019.06.011 [26] Jiang X Q, Li S, Furfaro R, Wang Z B, Ji Y D. High-dimensional uncertainty quantification for Mars atmospheric entry using adaptive generalized polynomial chaos. Aerospace Science and Technology, 2020, 107: Article No.106240 doi: 10.1016/j.ast.2020.106240 [27] Zhang W, Wang Q, Zeng F Z, Yan C. An adaptive sequential enhanced PCE approach and its application in aerodynamic uncertainty quantification. Aerospace Science and Technology, 2021, 117: Article No.106911 doi: 10.1016/j.ast.2021.106911 期刊类型引用(1)
1. 高多志，龚有敏，郭延宁，李传江，马广富. 考虑导航误差的火星探测器动力下降段鲁棒轨迹优化. 中国科学:物理学 力学 天文学. 2025(02): 174-186 . 
百度学术其他类型引用(0)
-
下载:
计量
- 文章访问数: 465
- HTML全文浏览量: 212
- PDF下载量: 152
- 被引次数: 1
