基于双层凸优化模型的月面运载器大范围转移最优轨迹规划方法

乔衍迪，张泽旭

（哈尔滨工业大学 深空探测基础研究中心，哈尔滨 150001）

引言

月球作为地球的天然卫星，被视为探索前沿航天科技的实验基地和未来火星探测任务的跳板，“嫦娥”工程作为中国深空探测的基础工程，在过去的几十年间取得了举世瞩目的巨大成功[1]，随着深空发射、自主变轨、自主导航、制导与控制技术的成熟，月球南极探测和月球基地的建设从论证阶段被列入航天发展的实施计划，各个国家均提出了自己的月球探测任务或月球基地建设规划[2]。未来聚焦在月球南极的科学探测任务，面临着昼夜温差大、地形复杂等重重挑战，在结合多项数据后，多数科研团队将月球基地的选址确认在南极点附近的沙克尔顿（Shackleton）及其附近的几个陨石坑[3]。月球资源勘探开采、月球基地建设等一切月基活动的关键基础是长时间稳定且充足的月基原位能源供给。月球水冰的大规模开采实现以电解水为原理的燃料能源系统构建是实现月球基地初期探测任务的主要途径，燃料系统以氢氧燃料电池为主。为了节约水冰运输与开采成本，各国都在研究可重复运载器，作为月球表面可靠运载工具[4]。本文提出了一种基于可重复运载器的月面大范围资源运载方案。

实现月球自主高精度垂直起降任务，可以将轨迹规划问题表述为具有非线性动力学、大量状态和控制约束的连续时间最优控制问题，通过直接法和间接法来解决，其中直接法是通过离散化方法将最优问题进行迭代求解，文献[5]将多阶段轨迹优化问题利用动态规划的思想，利用Legendre-Gauss-Radau配点法对多阶段轨迹规划问题离散化，通过求解非线性规划问题得到全局近似最优解。

近年来，随着凸优化理论的发展，提出了许多基于序列凸优化和无损凸优化方法的轨迹规划方法，文献[6]提出实现了6-DOF自主一般全尺度的轨迹优化问题。文献[7]将轨迹规划问题扩展为6-DOF的运载器动力学与运动学，并比较了不同离散化方法对于算法性能的影响；文献[8]在凸优化的基础上，利用虚拟控制函数和信赖域法，强化了轨迹规划算法的鲁棒性；文献[9]进一步将自适应信赖域半径法降低了凸优化算法的收敛性与收敛速度；文献[10]引入复合状态触发约束的概念，将其表示为连续凸化函数，在复杂约束条件下实现快速收敛。文献[11]利用神经网络优化了序列凸优化算法的收敛效率和实时性，提高了扩展性，在飞行器有效地求解了避障过程中复杂约束下的轨迹优化问题。

文献[12]针对飞行器垂直起降过程中大姿态翻转问题，提出了1种考虑气动力矩影响的平面运动模型，并利用凸优化方法在线规划出了飞行轨迹。问题是只考虑了二维轨迹问题。文献[13]将多阶段轨迹优化问题作为完整的系统进行优化，并在分阶段利用单独的模型替代原始模型，利用傅里叶级数降低求解难度，通过预设计得到最优发射再入轨迹。文献[14]针对月面飞跃问题，利用黄金分割法搜索上升着陆时间的最优解，在固定上升着陆时间的条件下，利用凸优化方法对分段求解拼接全飞行轨迹，不足是只得到了分段燃料最优解。

在求解月面大范围转移轨迹优化问题时，其一，利用凸优化工具求解，往往通过固定终端时刻来进行线性搜索，得到最优解，且模型非线性程度越高，规划时间越久，既无法保证算法的实时性，也无法保证整条轨迹的最优性；其二，考虑水冰采集任务中，运载器的飞行时间长，横向行程远大于纵向行程，着陆精度要求高，需要合理规划主推进器点火时长与时间点，以节省燃料消耗；其三，多数轨迹规划问题将任务简化为一个二维剖面，得到的推力分配方案是非现实的。基于上述分析，本文借鉴双层凸优化模型的思想，将大范围垂直起降问题转化为内外层凸优化问题，分段求解，在保证分段燃料最优性的同时，实现全局最优轨迹规划。

1 月球南极水冰采集任务概述

1.1 问题提出

月球基地是目前最为可行的地月空间基础设施开发的工程起点和最佳实践场。为了满足月球基地健康可持续运维能力，月球基地应具备勘探艾特肯盆地（Aitken Basin）任意陨石坑的能力，因此综合考虑重复往返各陨石坑水冰采集区的运输成本，选择Shackleton陨石坑的坑缘区域[15]，作为月球基地的长期运维地点，如图1所示。

图1 月球基地与水冰采集区域示意图（图中红色星星表示月球基地选址，黄色星星表示水冰采集区域选址）Fig.1 Illustration of lunar base and in-situ water excavation locations(red star indicates lunar base location;yellow stars indicate in-situ water excavation locations)

本文中将水冰采集区域放在距离月球基地较远的舒梅克（Shoemaker）陨石坑，德格拉什（de Gerlache）陨石坑，斯维德鲁普（Sverdrup）陨石坑等，如图1黄色星星所示，为了满足较长距离的可重复运输，运载器在设计时应满足此类任务的上限。

1.2 任务场景概述

如图2所示，运载器大范围运载任务可以描述成：小型液体火箭运载器分别装载完有效载荷和燃料后（LO2/LH2），氢氧发动机点火，保障运载器垂直上升，到达一定高度后，主发动机关机，运载器姿控发动机工作，帮助运载器重新定向，使运载器的头部指向月球基地，之后姿控发动机关闭，主发动机第2次点火，推动运载器飞向驻人月球基地，横向速度提升到一定数值后，主发动机关机，运载器无动力飞行，且在此过程中，姿控发动机工作，在运载器的尾部指向月球基地后，姿控发动机关机，在距离月球基地的一定范围时，主发动机重新点火，在水平速度消减至0 m/s后，关机，接着运载器姿控发动机工作，调整运载器尾部指向目标着陆点，主发动机点火，实现垂直软着陆。

图2 月面运载器大范围转移过程示意图Fig.2 Scheme of lunar large-scale transportation trajectory designation

水冰运输任务中，位于陨石坑底的水冰采集区与月球基地之间的海拔落差不超过10 km，但是运载器的飞行范围跨度可达上百km，运载器的主推进器角度可调节范围有限，无法及时完成运载器转身调姿任务，需要借助姿控发动机，进行翻转操作；其次，考虑到运载器燃料原位制备的可行性，采用液氢液氧发动机，优点是发动机比冲高，可达500 s，推力范围大，可达上百tnf，可重复点火，但是液体发动机燃料消耗速率相应也很快，无法完成连续推力控制；最后，考虑到任务的安全性，运载器的发射与着陆姿态均应垂直于月面。

基于上述分析，将运载器姿态调整与动力飞行拆分，分别进行优化设计，运载器飞行轨迹划分为垂直上升段、姿态调整段、动力飞行段和垂直下降段。

1.3 可重复运载器设计

为了满足月面大范围运载任务，包括GNC系统的质量在内，有必要定义一下实际的可重复运载器模型。由于此类运载器的没有公开且可靠的数据，参考经典的发动机信息。因此，参考文献[16]，将液氢液氧作为运载器发动机的燃料，该运载器为单级结构，包含助推器和箭身小型喷气式推进器。火箭发动机可包括液体燃料火箭发动机，例如LO2/LH2发动机。可包括多个小型推进器，用于控制火箭的上升和下降轨迹。

本节的目的是初步确定可重复使用运载器规格，以支持任务需求。包括运载器发动机的推力、推进剂的体积、流量、混合比，有效载荷的体积、质量，以实现推力矢量控制。同时，燃料与有效载荷的参数将直接决定运载器主体的直径，长度和其它构型。

根据图中飞行剖面，确定任务所需的最大推力，并确定燃料的质量与体积。主要考虑2个主要的动力段，即动力上升段与垂直下降段，忽略运载器姿态翻转所需的少量燃料。

本文从运载器到达初始滑行速度或动力上升段结束时的速度Vpa后的无动力上升段开始计算，通过式（1）确定无动力上升段的初始速度Vpa：

其中：hf为最大上升高度；hpa为动力上升段的高度；Vpa为动力上升段结束时运载器的速度的纵向分量；gm为月球的重力加速度，gm=1.62 m/s2。

为了保证运载器能够达到最大高度，运载器可以综合考虑高度与速度要求，动力上升段hpa结束得越早，对段终的速度Vpa要求就越高。通常，无动力上升段的高度应该大于动力上升段的高度。

在动力上升段，运载器的主推进器第一次点火，当确定了动力上升段结束时的速度Vpa后，再来确定动力上升段的时间，即发动机点火时长。

公式（2）为加入重力修正的Tsiolkowsky方程，提供了推进剂质量与动力上升段发动机开机时长tpa的关系，假设该方程只作用于运载器的纵轴

同时，根据式（3），可以得到动力上升段高度与发动机开机时长的关系

根据式（1）与（2）判断动力上升段结束时的速度与高度关系，在确定运载器的燃料载荷后，可以通过式（2）与式（3）得到运载器的总重和动力上升段的发动机点火时长的曲线，在后续轨迹规划过程中做出合理性判断依据。

2 月面运载器大范围垂直起降轨迹建模

2.1 坐标系与运动学方程

如前所述，构建局部参考坐标系，如图3所示，假设月球引力场均匀，且考虑月球自转的影响。

图3 坐标系示意图Fig.3 Guidance reference frame

运载器主发动机初始推力方向与运载器纵轴指向一致；

1）发射点坐标系—OX0Y0Z0

如图3所示，坐标原点为运载器在水冰采集区的发射起点，Z0轴为坐标原点与月心的连线且垂直向上，X0指向运载器的着陆点方向，Y0轴与X0轴、Z0轴构成右手正交坐标系。

运载器的运动学方程为

2.2 轨迹优化问题建模

运载器大范围垂直起降轨迹最优考虑的是燃料的消耗问题，姿态调整段、动力飞行段等飞行阶段的初始条件和终端条件不作约束，只对运载器发射点的位置、速度以及目标位置和速度做出约束

其中：r0和v0为运载器在发射点的位置和速度；rf和vf为运载器在着陆点的位置和速度；m0为运载器在初始质量，包含有效载荷和燃料以及运载器干重。

对于运载器的主推进器推力约束可以表示为

其中，Tmin和Tmax为主推进器的推力范围。

各飞行阶段均将燃料最优作为轨迹规划目标，则目标函数为

则动力学方程、控制约束、状态约束与目标函数，共同构成了月面运载器大范围垂直起降的轨迹优化问题方程。

2.3 问题的凸优化

为了运用凸优化算法对轨迹优化问题求解，需要将观察到的非凸控制约束推力进行凸化，即引入松弛变量 Γ，令

则原始轨迹优化问题中的推力约束可以写作

目标函数为

式（11）所描述的优化问题是对原始优化问题的松弛，原始优化问题的解显然也满足式（4）。此外，由文献[17]可知，对于松弛后的优化问题，当取最优解时，对于松弛后的轨迹优化问题，若存在最优解，则

成立，同时满足边界条件（10）。即原始最优问题与凸化后的最优问题可以看作是等价的。

为了数值求解最优轨迹规划问题，定义一组新的变量σ(t)，将非线性规划问题转换为具有凸性目标函数的连续时间最优控制问题，将状态约束与控制约束转换为二阶锥约束或线性约束

则运载器动力学方程转换为

令

则运载器推力约束可以写作

等式右侧为非凸域，利用泰勒展开实现二阶锥和线性逼近，构造可行凸域，即z0(t)=ln(m0-αTmaxt)

其中，为运载器燃料最大消耗量，z(t)的约束为

则最终运载器大范围垂直起降燃料最优轨迹规划问题可以表述为

动力学方程转换为（9），控制约束为（11）、（12），燃料消耗约束为（13），状态变量约束条件不变。

为了求解最优轨迹规划问题，需要将问题离散化，即将[t0,tf]拆分为等距的时间间隔，并在时间节点处施加线性或二阶锥约束，将无限维优化问题转化为有限维二阶锥优化问题求解，离散化的最优轨迹规划问题可以表述为

动力学方程

控制约束

其中，

状态约束

边界条件不变。

3 基于双层凸优化模型的梯度下降算法

在求解上述离散化最优轨迹规划问题时，凸优化方法无法在姿态调整段与动力飞行段处理推力分配，导致在燃料有限的情况下，轨迹优化问题无可行解；其次，在求解最优燃料问题时，凸优化算法在设置燃料最优指标下，推力分配不符合bang-bang特性，导致运载器上升高度过高，造成燃料浪费或燃料消耗大于搭载值。为了解决上述问题，引入双层凸优化模型的思想，将大范围轨迹优化问题看作双层凸优化问题，并利用梯度下降法进行最优燃料轨迹规划求解。

3.1 双层凸优化模型构造

将多段轨迹优化问题转换为双层凸优化问题，即在内层光滑凸函数的约束极小值集上实现外层光滑凸函数最小化的问题，并利用梯度下降法得到优化问题的数值解。对于内层凸函数，利用内点法进行最优问题求解，再得到可行解后，加入惩罚因子，利用梯度下降法求得外层凸函数的最优解。

下面介绍双层凸优化问题的构造，首先，将外层函数f1(x)的可行解集表示为

其中，g(x)为内层函数的不等式约束，则内层函数约束可以定义为f2(x)

双层凸优化问题可利用梯度投影法求解，即迭代解算双层凸优化问题的基函数，其中D为闭合凸子集。

其中：η为惩罚因子；f1(x)为内层光滑凸函数，即原始分段轨迹优化问题；f2(x)为外层光滑凸函数。

为了满足内层函数的解算的收敛性，需要保证

为了得到外层函数的最优解，在求解内层目标函数的基础上，需要适当选择函数的可行解集D或适当调整式（28）中的惩罚因子 η，其中，D可以是非负性约束，或球形约束，或某个R维状态空间。

在解算内层优化问题时，直接求解带有边界约束的凸优化问题，会有效提高算法的计算效率，因此只在外层优化问题加入惩罚因子，将优化问题写作

在求解过程中，双层凸优化问题采用显式求解方法，在 ηk的迭代过程中不计算最小正则化函数，不存在迭代收敛问题，降低了计算成本。其次，利用梯度下降法解算外层凸优化问题时，可以不求解最优化内层函数，因此不会增加任何“隐形”计算成本。

综上，本文利用的双层凸优化算法思想，利用内点法求解内层复杂约束下的二阶锥规划问题，内点法的单次迭代特点，可以很快得到可行解，保证算法的较高的求解效率；利用梯度下降法求解外层凸优化问题，通过惩罚因子将外层约束优化问题转换为无约束问题，利用梯度下降法多次迭代尝试得到完整轨迹的燃料最优解。内点法与连续时间最速下降法共同得到一组最优轨迹的可行解，其中惩罚因子的更新方法可以对照序列二次规划法。

3.2 算法流程

本文提出一种基于双层凸优化模型的梯度下降算法，内层循环利用内点法求解分段轨迹优化问题的可行解，外层循环利用梯度下降法迭代更新完整轨迹规划问题的最优解。

双层凸优化问题中，内层凸优化问题的目标函数与外层凸优化问题一致，如式（15）所示，其它约束为式（17）～（20），边界条件不变。但是外层凸优化问题的状态变量必须为内层凸优化问题的可行解的集合，可表述为

下面对基于双层凸优化模型的最优轨迹规划算法总结如下：

步骤1：通过内点法求解从初始状态x0到指定着陆点的分段最优燃料标称轨迹，令k=0。

步骤2：判断是否满足外层凸优化问题的约束；若不满足，则重新划定外层约束，求解分段最优燃料标称轨迹；若满足，则给定一组参数，，η 和 θ，并规定＞0，η∈(0,1)，θ∈(0,1)

步骤3：令m=T到0，计算

并得到满足条件的最小负整数m，令

sk=ηmk s

步骤4：更新

步骤5：判断是否满足f1(xk+1)＞f1(xk)，若满足，则k=k+1，再次计算分段最优燃料标称轨迹，直到f1(xk+1)≤f1(xk)。

算法流程如图4所示。

图4 基于双层凸优化模型的轨迹优化算法流程图Fig.4 Flow diagram of the bilevel model optimal trajectory optimization algorithm

4 仿真试验

下面给出仿真试验的条件：

将水冰采集区设计在Shoemaker坑的永久阴影区，具体位置为88.051 7 °S，45.488 1 °E，海拔为–4 200 m，着陆点为Shackleton陨石坑的坑沿位置，具体位置为89.678 2 °S，143.915 6 °W，海拔为1 660 m。运载器高2.65 m，直径0.8 m，可运输水冰400 kg，装载燃料100 kg，发动机比冲450 s，推力为3.2 t。

图5～8为基于双层凸优化算法得到的轨迹规划曲线，根据任务描述，发动机点火，以最大推力水平工作4.2 s，如图5所示，主要作用为将运载器尽快带离发射点，之后主推进器关机，利用姿控发动机单独对运载器进行姿态调整，调整时间为8.7 s，主推进器再次开机，推力水平为2.61 tnf，时长为5.9 s，此时横向速度与位置开始变化，如图6～7所示，当X向速度上升至175 m/s，Y向速度上升至265 m/s时，主推进器关机，运载器进入滑行状态，至纵向位置上升至9 899 m时，主推进器再次开机，推力水平为2.21 tnf，时长为5.7 s，将X向速度与Y速度降至0，且运载器到达目标着陆点上空，此时利用姿控发动机再次调整运载器姿态，调整时间为8.2 s，当运载器机头垂直向上时，主推进器开机，推力水平为3.2 tnf，时长为3.0 s，此时纵向速度降为0，实现运载器安全着陆。图8为单次仿真完整的轨迹规划示意图。

图5 推力示意图Fig.5 Results of thrust magnitude

图6 三轴速度示意图Fig.6 Results of velocity magnitude

图7 位置示意图Fig.7 Results of trajectory magnitude

图8 运载器飞行轨迹示意图Fig.8 Three dimensional trajectory history

随机生成动力上升段结束时刻产生的500个随机位置误差，其中X、Y向的位置误差标准差为±500 m，Z向位置的标准差±200 m，符合正态分布，利用蒙特卡洛试验验证算法的鲁棒性。

如图9和图10所示，运载器水平方向位置误差在±0.6 m以内，速度误差在±0.02 m/s以内实现了运载器的高精度着陆精度要求；如图11与图12所示，由于最终着陆阶段时间较短，纵向位置误差在±0.2 m以内，误差多大于目标高度，这对于软着陆任务是有利的，避免运载器直接接触月面而坠毁，纵向速度误差±6 × 10–3nm/s以内，保障了运载器的安全着陆；根据剩余质量约束，仿真的成功率为97.4%，其中有13次仿真试验剩余燃料小于0，其次燃料剩余最多分布在410 kg附近，根据液体燃料即时加注的特点，燃料的利用率达到80%以上。

图9 最终横向位置误差柱状图Fig.9 Histogram of final horizon position errors

图10 最终横向速度误差柱状图Fig.10 Histogram of final horizon velocity errors

图11 最终纵向位置误差柱状图Fig.11 Histogram of final vertical direction errors

图12 最终纵向速度误差柱状图Fig.12 Histogram of final Z direction velocity errors

接着，随机生成动力上升段结束时刻产生的500个随机速度误差，其中X、Y向的位置误差标准差为±70 m/s，Z向速度误差的标准差±50 m/s，符合正态分布，利用蒙特卡洛试验验证算法的鲁棒性。

如图13和图14所示，运载器水平方向位置误差在±0.4 m以内，速度误差在±0.4 m/s以内实现了运载器的高精度着陆精度要求，且由于Y向行程较长，位置和速度误差分布优于X向误差，尤其是速度误差，64%的误差在±0.1 m/s以内；如图15与图16所示，由于最终着陆阶段时间较短，纵向位置误差在±0.4 m以内，纵向速度误差±1.5 × 10–3m/s以内，虽然精度有所下降，但依然保障了运载器的安全着陆；根据剩余质量约束，仿真的成功率为96.7%，其中有16次仿真试验剩余燃料小于0，其次燃料剩余最多分布在407 kg附近，根据液体燃料即时加注的特点，燃料的利用率达到80%以上。

图13 最终横向速度误差柱状图Fig.13 Histogram of final horizon velocity errors

图14 最终横向位置误差柱状图Fig.14 Histogram of final horizon position errors

图15 最终纵向速度误差柱状图Fig.15 Histogram of final Z direction velocity errors

图16 最终纵向位置误差柱状图Fig.16 Histogram of final vertical direction error

如图17所示为500次仿真试验得到的运载器大范围垂直起降的轨迹图，为了更明显的显示初始位置不同，采用了对数坐标，可以看到得到的轨迹全部收敛于最终着陆点，证明算法具有良好的收敛性与稳定性。

图17 运载器飞行轨迹图Fig.17 The 3-Dimensional Trajectory of the vehicle

4 结论

本文利用双层凸优化思想构造了月面运载器大范围最优轨迹规划问题，将大范围飞行轨迹划分为动力上升段、大范围动力飞行段、垂直下降段的最优化问题，并将其作为内层凸优化问题，通过内点法进行求解，同时，结合全局燃料最优化目标，设计目标函数，作为外层凸优化目标函数，利用最速下降法的思想，进行迭代求解，得到了大范围垂直起降的燃料最优轨迹。仿真试验证明，本文提出的算法可以实现运载器大范围垂直起降轨迹优化问题，且在初始位置误差在700 m以上时规划出满足着陆约束，并实现水平位置误差小于0.6 m，高度误差小于0.2 m的运载器轨迹，且燃料利用率达80%以上，剩余燃料有97.4%的概率满足约束。