0 引 言
小天体（小行星、彗星）是深空探测的前沿目标。开展小天体探测，对各类小天体进行全面的研究，能够促进对行星演化和太阳系的形成的了解。近年来，小天体探测方法主要集中在着陆和采样两个方面。截至目前，已有5次探测任务完成了在小天体上的着陆，包括“近地小行星交会”（Near Earth Asteroid Rendezvous，NEAR）任务、“隼鸟”（Hayabusa）任务、“隼鸟”2任务、“罗塞塔”（Rosetta）任务和“奥西里斯-雷克斯”（OSIRIS-REx）任务。美国的NEAR探测器通过预先制定的变轨序列进行4次制动，并在着陆过程中采用惯性定向姿态，成功在爱神星（Eros）上着陆[1]。在欧空局的Rosetta任务中，“菲莱”（Philae）着陆器依靠与轨道器分离时的速度，以自由落体的方式降落在彗星上，并在着陆过程中保持纵轴惯性定向[2-3]。日本宇宙航空研究开发机构研制的Hayabusa探测器[4-5]和Hayabusa 2探测器[6]采用地形跟踪和目标标志器跟踪的导航方式，使用闭环自主制导和控制技术，实现了在小天体上的接触着陆。美国国家航空航天局的OSIRIS-REx小行星探测器根据事先制定的TAG（Touch-and-Go）序列，通过3次制动，将探测器从1 km轨道依次导引至检查点（Checkpoint）和匹配点（Matchpoint），并成功在小天体上接触着陆[7]。在以上3种着陆策略中，闭环制导的着陆精度最高，Hayabusa 2探测器可以达到3 m的着陆精度[6]，而OSIRIS-REx任务的着陆精度低于25 m[7]。此外，2025年5月，我国发射了“天问”二号探测器，对小行星2016HO3和彗星311P进行探测，其主任务是实现对小行星的测绘、下降、着陆和采样，其着陆精度小于1 m。综上所述，高控制精度的闭环制导方式是实现在小天体上定点着陆和采样的重要途径。
考虑到小天体的不规则形状和高速自旋，着陆制导问题需要综合考虑多种约束，以确保探测器的安全性和满足观测目的。因此，需将这样一个复杂的制导问题转化为满足各种约束条件的轨迹优化问题。轨迹优化问题的求解方法可分为直接法和间接法两类。间接法计算速度快，但收敛效果通常不佳[8-9]。直接法是将原问题转化为非线性规划问题（Non-Linear Programming，NLP）求最优解，通常计算速度较慢，且对初始估计不敏感[10-11]。在求解NLP问题的方法中，凸优化方法因能在多项式时间内求解[12]，被越来越多地应用于着陆轨迹规划问题[13-16]。凸优化方法要求代价函数和不等式约束均为凸函数，等式约束均为仿射函数，因此需先将一般非凸NLP问题转化为凸问题[17]，然后通过离散化将连续时间凸优化问题转化为离散时间问题，并通过求解器（如内点法、序列非线性规划法等）进行解算[18-20]。在离散方法中，伪谱法提供了一种可获得平滑结果的离散连续时间问题的方法[21]。该方法利用非均匀高斯积分网格和拉格朗日差分多项式对变量进行近似，并使用NLP问题的稀疏结构进行计算，所需计算机存储量少，计算精度高[22]。
本文结合凸优化方法与伪谱法，求解小天体着陆轨迹规划问题。首先，建立具有动力学、约束和性能指标函数的小天体着陆轨迹规划的最优控制问题，并给出将一般NLP问题转化为凸优化问题的方法，进而构建小天体着陆轨迹规划的凸优化问题。基于Radau伪谱方法，建立分段状态约束下的小天体着陆轨迹规划伪谱序列凸优化问题，并给出求解过程。最后，通过仿真验证该算法在小天体自主着陆中的可行性和有效性。
1 分段状态约束小行星着陆轨迹规划最优控制问题
1.1 动力学
首先，给出建立动力学模型的两个坐标系：小行星固连坐标系${\mathcal{F}_{{\rm{L}}}}$ 和探测器本体坐标系${\mathcal{F}_{{\rm{B}}}}$ 。${\mathcal{F}_{{\rm{L}}}}$ 为原点${o_{{\rm{L}}}}$ 位于小天体的质心，${o_{{\rm{L}}}}{z_{{\rm{L}}}}$ 轴沿小行星的自转轴，${o_{{\rm{L}}}}{x_{{\rm{L}}}}$ 沿小天体最大（最小）惯量主轴，${o_{{\rm{L}}}}{y_{{\rm{L}}}}$ 轴构成右手坐标系。${\mathcal{F}_{{\rm{B}}}}$ 为原点${o_{{\rm{B}}}}$ 位于探测器的质心，${o_{{\rm{B}}}}{z_{{\rm{B}}}}$ 轴沿探测器的最大惯量轴，${o_{{\rm{B}}}}{x_{{\rm{B}}}}$ 沿探测器的最小惯量主轴，${o_{{\rm{B}}}}{y_{{\rm{B}}}}$ 轴构成右手坐标系。
在${\mathcal{F}_{{\rm{L}}}}$ 系下以矢量描述探测器的质心动力学，给出探测器的位置矢量${{{\boldsymbol{r}}}_{{\rm{L}}}}(t)$ 和速度矢量${{{\boldsymbol{v}}}_{{\rm{L}}}}(t)$ ，则质心动力学为：
式中，$ {{\boldsymbol{\omega}}_e} $ 为${\mathcal{F}_{{\rm{L}}}}$ 系下的小天体自转角速度，$ {T_{{\rm{L}}}}(t) $ 和$ \nabla U(r) $ 为${\mathcal{F}_{{\rm{L}}}}$ 系下探测器的推力矢量和所受到的小行星引力。
在${\mathcal{F}_{{\rm{B}}}}$ 系下以矢量描述探测器的姿态动力学，并采用修正罗德里格参数（Modified Rodrigues Parameters，MRPs）${{{\boldsymbol{\sigma}}} _{{\rm{B}}{\rm{L}}}}(t)$ 来描述${\mathcal{F}_{{\rm{B}}}}$ 系的空间指向，则姿态动力学为：
式中，$J$ 为探测器的惯性张量，${{\boldsymbol{\omega }}_{{\rm{B}}}}(t)$ 和${{\boldsymbol{M}}_{{\rm{B}}}}(t)$ 为${\mathcal{F}_{{\rm{B}}}}$ 系下探测器的转动角速度和控制力矩；${{\boldsymbol{\omega}} _{{\rm{B}}{\rm{L}}}}(t)$ 为${\mathcal{F}_{{\rm{L}}}}$ 系相对于${\mathcal{F}_{{\rm{B}}}}$ 系的转动角速度，考虑到小行星的自转，${{\boldsymbol{\omega}} _{{\rm{B}}{\rm{L}}}}(t)$ 和$ {{\boldsymbol{\omega}} _e} $ 有如下关系：
式中，${\boldsymbol{C}}_{{\rm{L}}}^{{\rm{B}}}({{\boldsymbol{\sigma}} _{{\rm{B}}{\rm{L}}}})$ 为${\mathcal{F}_{{\rm{B}}}}$ 系和${\mathcal{F}_{{\rm{L}}}}$ 系的方向余弦阵，${\boldsymbol{C}}_{{\rm{L}}}^{{\rm{B}}}({{\boldsymbol{\sigma}} _{{\rm{B}}{\rm{L}}}})$ 和式（3）中${\boldsymbol{R}}({{\boldsymbol{\sigma}} _{{\rm{B}}{\rm{L}}}})$ 表达式为：
式中，${{\boldsymbol{\sigma}} _{{\rm{B}}{\rm{L}}}}(t) = {\left[ {{{\boldsymbol{\sigma}} _1}, {{\boldsymbol{\sigma}} _2}, {{\boldsymbol{\sigma}} _3}} \right]^{\mathrm{T}}}$ 。
探测器的推力方向由姿态控制，定义${\mathcal{F}_{{\rm{B}}}}$ 系下的推力矢量为$ {{\boldsymbol{T}}_{{\rm{B}}}}(t) $ ，则$ {{\boldsymbol{T}}_{{\rm{L}}}}(t) $ 和$ {{\boldsymbol{T}}_{{\rm{B}}}}(t) $ 的关系为：
式中，$ {\boldsymbol{C}}_{{\rm{B}}}^{{\rm{L}}}({{\boldsymbol{\sigma}} _{{\rm{L}}{\rm{B}}}}) = {\boldsymbol{C}}_{{\rm{L}}}^{{\rm{B}}}{({{\boldsymbol{\sigma}} _{{\rm{B}}{\rm{L}}}})^{ - 1}} = {\boldsymbol{C}}_{{\rm{L}}}^{{\rm{B}}}{({{\boldsymbol{\sigma}} _{{\rm{B}}{\rm{L}}}})^{\mathrm{T}}} $ 。
此外，推力器推力的输出会使探测器的质量发生变化，假设燃料消耗的速率与推力矢量的大小成正比，则质量变化可表示为：
式中，${I_{{\mathrm{sp}}}}$ 为探测器推力器的真空比冲，${g_0}$ 为地球的标准重力常数。
1.2 状态约束
由于小天体一般外形不规则且自旋速度较快，在着陆过程中需确保探测器不与小天体表面发生碰撞，并实现对着陆点的观测。小天体的着陆过程可以分为两段：下降段和最终着陆段，在不同的着陆阶段，需采用不同的状态约束来满足安全和观测的需求。
在下降段，探测器距离着陆点较远，且通常不在小天体着陆区上方。此阶段的目的是降低探测器的高度，并转移到着陆区上方，同时消除探测器的水平方向速度，使其与小天体同步。因此，可以引入一个包围小行星的椭球，限定探测器在椭球外部飞行，以防止探测器与小天体发生碰撞。椭球防碰撞状态约束可以描述为：
式中，${{\boldsymbol{R}}_e} = {{diag}}(1/{a^2}, 1/ {b^2} , 1/ {c^2} )$ ，$\left( {a, b, c} \right)$ 为椭球半长轴，${t_s}$ 为下降段结束时间。
在最终着陆段，探测器已经从较远的位置转移到着陆点上空并会在椭球内飞行，所以引入一个以着陆点为顶点和锥角为$2\theta $ 的圆锥约束，限定探测器在此圆锥内飞行，则该圆锥防碰撞状态约束为：
式中，${{\boldsymbol{r}}_{{\rm{L}},ls}}$ 为着陆点的位置，${{\boldsymbol{n}}_{{\rm{L}},ls}}$ 为指向小行星外部的着陆点法向量，${t_f}$ 为着陆结束的时间。
此外，在最终着陆段，需要对着陆点进行跟踪以实现定点着陆，即要求着陆点一直在机载着陆相机的视场内，该视场角约束可描述为：
式中，${{\boldsymbol{\rho}} _{{\rm{B}}}}$ 为着陆相机在探测器上的安装位置矢量，${{\boldsymbol{d}}_{{\rm{B}}}}$ 为相机光轴矢量。
1.3 控制约束和边界条件
探测器的机动能力由推力器的输出控制力和飞轮输出的控制力矩决定，推力和控制力矩的大小受到限制，所以控制约束为：
式中，${T_{{\mathrm{max}}}}$ 为推力器的最大推力，${M_{{\mathrm{max}}}}$ 为飞轮的最大控制力矩。
对于边界条件，探测器从空间某观测点(悬停点、转移点）开始着陆过程，该观测点的探测器状态已知，探测器的初始状态是该观测点的状态；而着陆过程随探测器接触小天体上某已知着陆点结束，探测器的终端状态就是该着陆点的状态。边界条件可以表示为：
式中，$\left( {{{{\boldsymbol{r}}}_{{\rm{L}},op}},\; {{{\boldsymbol{v}}}_{{\rm{L}},op}},\; {{\boldsymbol{\omega}}_{{\rm{B}},op}}, {{\boldsymbol{\sigma}} _{{\rm{B}}{\rm{L}},op}}} \right)$ 为着陆起始观测点的状态值，$\left( {{{{\boldsymbol{r}}}_{{\rm{L}},ls}},\; {{{\boldsymbol{v}}}_{{\rm{L}},ls}},\; {{\boldsymbol{\omega}}_{{\rm{B}},ls}}, {{\boldsymbol{\sigma}} _{{\rm{B}}{\rm{L}},ls}}} \right)$ 为着陆点的状态值。
1.4 小天体着陆轨迹规划最优控制问题
以推进剂最优为性能指标函数，根据式（9），推进剂消耗会使探测器质量发生变化，推进剂最优为整个着陆期间探测器的质量变化最小，即有：
将式(9）代入式(16），则推进剂最优的性能指标函数为：
式中，$\alpha = {1}/({{{I_{{sp}}} \cdot {g_0}}})$ 。
则分段状态约束小天体最优着陆轨迹规划最优控制问题可以表述为：满足式（1）~（4）的质心动力学和姿态动力学、式（10）~（12）的状态约束和式（13）的控制约束，以及式（14）~（15）的边界条件的以探测器推进剂消耗最少为性能指标式（17）的最优控制问题。
2 小天体着陆轨迹规划凸优化问题
2.1 凸优化
凸优化问题的定义如下：
式中，函数${f_0}, \cdot \cdot \cdot ,{f_m}:{{\bf{R}}^n} \to {\bf{R}}$ ，${f_i}$ 满足：
式中，$x,y \in {{\bf{R}}^n}$ ，$\alpha ,\beta \in {\bf{R}}$ ，$\alpha + \beta = 1$ ，$\alpha \geqslant 0$ ，$\beta \geqslant 0$ 。式（19）是凸函数的定义，即凸优化中要求性能指标函数和所有的约束函数都是凸函数。
2.2 小天体着陆轨迹规划凸优化问题
由于第1.4节中描述的小天体着陆轨迹规划最优控制问题中，其动力学方程式（1）~（4）、椭球约束式（10）和视场角约束式（12）均为非凸函数，因此在凸优化框架下解决该最优控制问题，需先将该非凸问题转化成为凸问题。
定义${\boldsymbol{x}}\left( t \right) = {\left[ {{{{\boldsymbol{r}}}^{\mathrm{T}}_{{\rm{L}}}}{{(t)}}, {{{\boldsymbol{v}}}^{\mathrm{T}}_{{\rm{L}}}}{{(t)}}, {{\boldsymbol{\omega}}^{\mathrm{T}}_{{\rm{B}}}}{{(t)}}, {{\boldsymbol{\sigma}} ^{\mathrm{T}}_{{\rm{B}}{\rm{L}}}}{{(t)}} } \right]^{\mathrm{T}}}$ 为探测器的状态向量，定义探测器的控制向量为${\boldsymbol{u}}\left( t \right) = {\left[ {{{\boldsymbol{T}}^{\mathrm{T}}_{{\rm{B}}}}{{(t)}}, {{\boldsymbol{M}}^{\mathrm{T}}_{{\rm{B}}}}{{(t)}}} \right]^{\mathrm{T}}}$ ，则探测器动力学式（1）~（4）可转化成${\boldsymbol{x}}\left( t \right)$ 和${\boldsymbol{u}}\left( t \right)$ 描述的方式，则有：
式中，$f\left[ {{\boldsymbol{x}}\left( t \right), {\boldsymbol{u}}\left( t \right)} \right] \in {{\bf{R}}^{12}}$ 为非线性非凸的函数，可以采用线性化的方式将其转化成凸函数。
给出参考轨迹$\left\{ {\hat {\boldsymbol{x}}\left( t \right),\;\hat {\boldsymbol{u}}\left( t \right),\;\hat m\left( t \right)} \right\}$ ，在参考轨迹附近对式（20）右侧进行一阶泰勒展开，则有：
对于式（10）中的椭球约束，同样可以采用线性化的方式对不等式左边进行凸化，则凸化的椭球防碰撞状态约束为：
式中，${{\boldsymbol{C}}_r} = \left[ {{{\boldsymbol{E}}_3},\;{{\boldsymbol{0}}_{3 \times 9}}} \right]$ ，${{\boldsymbol{P}}_e} = {{\boldsymbol{C}}_e}^{\mathrm{T}}{{\boldsymbol{R}}_e}{{\boldsymbol{C}}_e}$ 。
对于式（11）的圆锥防碰撞状态约束，将其转换成用状态$x\left( t \right)$ 描述的凸二阶锥函数形式：
式中，$x\left( {{t_f}} \right)$ 为探测器的终端状态。
对于式（12）中的视场角约束，将其转化成为以下形式：
式中，${{\boldsymbol{l}}_{{\rm{B}}}}\left[ {{\boldsymbol{x}}\left( t \right)} \right] = {\boldsymbol{C}}_{{\rm{L}}}^{{\rm{B}}}\left[ {{{\boldsymbol{C}}_{\boldsymbol{\sigma}} }x\left( t \right)} \right]{{\boldsymbol{C}}_r}\left[ {{\boldsymbol{x}}\left( t \right) - {\boldsymbol{x}}\left( {{t_f}} \right)} \right]$ ，${{\boldsymbol{C}}_{\boldsymbol{\sigma}} } = \big[ {{{\boldsymbol{0}}_{3 \times 9}}, {{\boldsymbol{E}}_3}} \big]$ 。由于${l_{{\rm{B}}}}\left[ {x\left( t \right)} \right]$ 为非凸函数，所以式（24）也是非凸的。可以对${l_{{\rm{B}}}}\left[ {x\left( t \right)} \right]$ 线性化使其变成凸函数，则式（24）可以转化成二阶锥形式的凸函数。凸化的视场角约束为：
定义探测器的初始状态和终端状态为：
用控制向量${\boldsymbol{u}}\left( t \right)$ 描述控制约束式（13）和性能指标函数式（17）：
式中，${{\boldsymbol{C}}_T} = \left[ {{{\boldsymbol{E}}_3}, {\boldsymbol{{0}}_3}} \right]$ ，${{\boldsymbol{C}}_M} = \left[ { {{\boldsymbol{0}}_3}, {{\boldsymbol{E}}_3}} \right]$ 。
通过对小天体着陆轨迹规划最优控制问题凸化后，小天体着陆轨迹规划凸优化问题陈述为：满足式（21）的探测器动力学，式（22）、式（23）、式（25）的状态约束和式（29）的控制约束，以及式（27）~（28）的边界条件的以探测器燃料推进剂最少为性能指标式（30）的凸优化问题。
3 分段状态约束小天体着陆轨迹规划伪谱序列凸优化求解
由于构建的小天体着陆轨迹规划凸优化问题为连续时间问题，需将其转化为有限维离散时间问题才能进行求解。传统的离散方法，如有限差分法，其误差通常仅为代数速率收敛，导致求解速度较慢；而伪谱法具有指数收敛性，尤其采用满足协态映射定理的Legendre-Gauss-Radau（LGR）配点的Radau伪谱法，其精度和收敛速度均优于采用其他配点的伪谱法。因此，利用Radau伪谱法对连续时间问题进行离散，只需更少的时间节点即可达到所需精度，从而减少了需要处理的状态量，加快求解速度。
小天体着陆轨迹规划伪谱凸优化问题，正是结合Radau伪谱法，将连续的小天体着陆轨迹规划凸优化问题离散化，转化为有离散、限维的二阶锥规划（Second-Order Cone Programming，SOCP）问题。该问题可以使用内点法进行求解，且其解为全局最优解。然而，在该问题转化为SOCP的过程中，采用线性化方法对原小天体着陆最优控制问题进行凸化处理，虽然该方法大幅降低了非凸问题的求解难度，但对初始参考轨迹的选择却极为苛刻。为此，采用序列凸优化的方法对该伪谱凸优化问题进行序列迭代求解，以确保所获得的解能够收敛至原问题的解，并且具有较好的局部最优性。
3.1 Radau伪谱法
伪谱法是一种利用多项式插值来逼近状态和控制变量，并用正交配点对微分方程进行离散化处理，将连续时间的问题转化为离散的NLP问题的最优控制求解方法。Radua伪谱法作为一种非对称伪谱方法，允许对问题的初始条件有一个自然而直接的定义，且与其他方法相比无须积分计算终端状态以施加终端约束。
Radua伪谱法中的配点是Legendre-Radau多项式的根，其定义为：
式中，$ {\tilde L_n}\left( \tau \right) = \dfrac{1}{{{2^n}n!}}\dfrac{{{d^n}}}{{d{\tau ^n}}}{\left( {{\tau ^2} - 1} \right)^n} $ 。
伪谱时间区域$ \tau \in \left[ { - 1, 1} \right] $ 可由仿射变换给出，其表达式为：
则最优控制问题中的微分方程用$\tau $ 表示为：
最优控制问题中的约束方程用$\tau $ 表示为：
最优控制问题中的边界用$\tau $ 表示为：
所以可以用$ \tau \in \left[ { - 1, 1} \right] $ 这个离散区域对${\boldsymbol{x}}\left( \tau \right)$ 进行连续近似重构，有：
式中，${{\boldsymbol{X}}_i}$ 为第i个配点的状态。
因为配点是与时间无关的，对式（36）进行关于$\tau $ 微分，可得：
结合式（37）和式（33）合，则式（33）描述的微分方程可替换为：
式中，${D_{ki}} = {\dot P_i}\left( \tau \right)$ ，$k = 1,2,\cdots,n$ ，$i = 0,1,\cdots,n$ 。定义${\boldsymbol{D}} = \left[ {{D_{ki}}} \right] \in {{\bf{R}}^{n \times \left( {n + 1} \right)}}$ ，D被称为Radau伪谱状态微分矩阵，${D_{ki}}$ 可以用重心坐标拉格朗日（Barycentric Lagrange）插值有效地计算。
除微分运算外，还需进行积分运算。该运算符出现在如式（30）的性能指标函数中，包含拉格朗日项，需适当离散化处理。利用高斯积分公式，对于Radau伪谱法，可以用离散和来代替连续积分，其表达式为：
式中，权重$ {{\boldsymbol{\omega}}_k} $ 可以用式（40）计算：
其中，$flip $ 表示中心反对称操作符。
则一个连续时间的最优控制问题通过伪谱法转化成为一个NLP问题，表达式为：
式（41）和式（42）分别为用配点${{\boldsymbol{X}}_i}$ 表示的约束方程和边界条件。
3.2 分段状态约束小天体着陆轨迹规划伪谱序列凸优化问题
序列凸优化是通过求解一系列相关的凸子问题，间接地解决非凸优化问题。在本文探讨的小天体着陆轨迹规划问题中，每个相关的凸子问题被定义为离散的小天体着陆轨迹规划的SOCP问题。伪谱序列凸优化的过程如下。
首先，引入三重表示法来区分迭代过程、迭代中的各个段以及段中的各个点，并定义每个配点的状态和控制为：
式中，i表示迭代次数，p为最终迭代次数，j表示飞行阶段，k为每个飞行阶段中的配点数。
将小天体着陆轨迹规划凸优化问题的时间在Radua伪谱法的配点上进行离散，每次迭代中每个飞行阶段的飞行时间定义为：
令第（i−1）的迭代结果为第i次迭代中对动力学和非凸状态约束的线性化的参考轨迹，则对于式（21）中的动力学，用Radua伪谱法离散后，有：
式中，${D^{i,j}}$ 为第i次迭代中第j段的Radau伪谱状态微分矩阵，$ V_k^{i,j}$ 为第i次迭代中第j段的第k个点的虚拟控制，其作用是用来消除对动力学线性化中的人工不可行性，保证每次迭代中每个小天体着陆轨迹规划的凸子问题都有解。
用伪谱法离散后的椭球防碰撞约束可以表示为：
用伪谱法离散后的圆锥防碰撞约束可以表示为：
用伪谱法离散后的视场角约束可以表示为：
用伪谱法离散后的边界条件可以表示为：
用伪谱法离散后的控制约束可以表示为：
每个飞行阶段的连接条件为：
对于式（30）描述的性能指标函数，用Radua伪谱法离散后，可表示为：
式中，性能指标函数中增加了虚拟控制$ {\left\| { V_k^{i,j}} \right\|_1} $ 项，$ {W_v} $ 为虚拟控制项的权重。可以验证，增加虚拟控制项的优化问题的最优解，同样是原问题（即无虚拟控制的性能指标的优化问题）的最优解。
对小天体着陆轨迹规划凸优化问题，在离散化处理后，分段状态约束下的伪谱序列凸优化问题可表述为：需满足式（45）的等式约束、式（46）~（48）和式（52）的不等式约束，以及式（50）~（51）的边界条件和式（53）~（54）的连接条件，并以式（55）作为性能指标的序列凸优化问题。
4 仿真与结果分析
本节将通过仿真实验来验证伪谱凸优化方法在求解小天体着陆轨迹规划问题中的可行性。该仿真采用Matlab软件凸规划软件包CVX及SDPT3 4.0求解器进行计算。首先，给出探测器的物理参数如表1 所示。
表 1 探测器参数Table 1 Probe parameters |
| 参数 | 数值 |
| 初始惯性张量/(kg·m2) | |
| 初始质量/(kg) | 1 400 |
| 推力约束Tmax/(N) | 20 |
| 推力矩约束Mmax/(N·m) | 0.2 |
| 斜坡约束锥角 | 15 |
| 相机视场角约束 | 25 |
| 相机安装位置 | [0.9, 0, −1.0]T |
| 相机光轴矢量 | [0, 0, −1.0]T |
| 推力器真空比冲 | 225 |
| 地球标准重力常数 | 9.806 65 |
对于每次迭代的质量变化，可用式（56）表达：
着陆目标小行星选取为围绕太阳公转的小行星4769（Castalia）。在仿真过程中，探测器的初始点状态向量${\boldsymbol{X}}_0^{i,1}$ 、预着陆点的状态向量${\boldsymbol{X}}_n^{i,2}$ 及预着陆点的法向量${\boldsymbol{n}}_{{\rm{L}},ls}^{i,2}$ 如表2 所示。
表 2 探测器边界条件Table 2 Probe boundary conditions |
| 状态 | 数值 |
| [−225.202 6, −925.325 3, −679.098 0, 1.753 6, −1.767 4, 0.118 1, 0, 0, 0, 0.471 1, 0.328 4, −0.100 8]T | |
| [34.935 5, −292.355 7, −255.997 8, 0, 0, 0, 0, 0, 0, 0.548 9, 0.542 4, −0.422 1 ]T | |
| [-0.463 8, −0.741 2, −0.485 ]T |
给出探测器的飞行时间$ {t_s} = 360 $ s和$ {t_f} = 600 $ s，每个飞行阶段的配点数分别为36个和24个，仿真结果如图1 ~5 所示。在仿真过程中，设置最大迭代次数为10，求解精度（即相邻两次迭代中状态偏差的无穷范数）为10−5。当迭代次数大于10或相邻两次迭代的状态误差小于10−5时，迭代结束，并输出最后一次迭代的结果为最优轨迹。
图 2 探测器着陆过程控制量变化Fig.2 Changes in control variables during the landing process of the probe |
从仿真结果可以看出，采用伪谱序列凸优化算法规划的小天体着陆轨迹中，所有的状态约束和控制约束均得到满足。此外，每个飞行阶段的开始和结束节点分布密集，而中间部分节点较为稀疏，这与采用非均匀高斯积分网格的特点相一致，从而使得生成的轨迹更为平滑。这说明伪谱凸优化方法在求解小天体着陆轨迹规划问题上具有有效性和可行性。
另外在图1 ~ 3 中，同时给出了采用一般NLP方法及基于Radua伪谱法的h个时间区间、p个控制参数的hp-伪谱法进行计算的小天体着陆规划结果。在伪谱序列凸优化算法中，迭代次数为6，中央处理器计算时间为45.785 352 s；而在hp-伪谱法中，迭代次数为3，中央处理器计算时间为115.945 200 s，伪谱序列凸优化的迭代次数较多，但总体计算时间较短，这符合该方法计算速度快的特点。此外，从图1 ~3 可以看出，采用伪谱序列凸优化方法规划的探测器轨迹更为平缓，且探测器重量变化更小，即探测器的推进剂消耗更低。这表明伪谱序列凸优化算法在计算速率和规划轨迹的优越性方面均表现更佳，从而为未来探测器在小天体的实时自主降落提供了可能性。
5 结 语
本文提出了一种基于伪谱序列凸优化的小天体着陆轨迹规划方法。首先，建立了基于MRPs的六自由度探测器着陆动力学模型，给出了分段的状态约束及控制约束。结合边界条件，构建了分段状态约束下的小天体着陆轨迹规划推进剂最优控制问题。为在凸优化框架内求解此最优控制问题，阐述了原问题的凸化过程，形成了小天体着陆轨迹规划的凸优化模型。然后，详细推导了采用Radau伪谱法将连续时间凸优化问题离散化的表达式，给出了具有分段状态约束的小天体着陆轨迹规划伪谱凸优化子问题的序列求解步骤。通过仿真验证了该算法的有效性，即该算法能够在多项式时间内解决复杂的小天体着陆轨迹规划问题，并生成平滑的轨迹。