基于同伦方法的地月系L2点小推力转移轨道优化

潘 迅，泮斌峰

(西北工业大学航天学院，西安 710072)

1 引言

在地月系三体模型中，L2点在地月连线上，且始终位于月球背面，不仅可以作为月球探测和深空探测的中转站，还可以在L2点建立月球观测站，实时监测月球背面情况，为月球背面着陆提供相关信息[1]。我国的嫦娥五号再入返回飞行试验器的服务舱曾于2014年11月28日到达地月系L2点，并在2015年1月4日离开L2转移至环月轨道[2-3]。

嫦娥五号再入返回试验器利用化学发动机经过多脉冲机动完成轨道转移，转移过程较为简单，所需时间较短。然而，随着推进技术的发展，电推进、等离子推进等持续小推力推进方式由于具有比冲大、控制精度高等优点，能极大地提高探测器的有效载荷的比重[4]，因此将在未来的航天任务，特别是深空探测任务中，扮演越来越重要的角色。

由于小推力发动机的推力幅值小，持续时间长，完成转移所需的轨道圈数多，在转移过程中初始状态量对转移轨道有较大影响，即对初始值很敏感，因此优化过程存在较大困难。此外，相比于传统的二体模型，限制性三体问题的非线性更强，导致该小推力转移轨道优化问题的求解难度增加。小推力转移轨道优化方法主要可以分为直接法、间接法和混合法[4]。直接法通过将状态量和控制量进行离散，将轨道优化问题转化为非线性规划问题，直接对离散参数进行寻解[5-6]；间接法基于最优控制理论的经典变分法推导得到最优控制的一阶必要条件，将其转换为两点边值问题，通过数值方法对该两点边值问题进行求解，从而得到转移轨道[7-8]；混合法只利用间接法推导得到的最优控制律，但不包含其他最优性条件，通过参数优化方法进行求解[9]。其中，间接法由于利用了最优性条件，可以保证所得解的局部最优性，且计算精度高。然而间接法将轨道优化问题转换成两点边值问题后，虽然理论上可以利用打靶法对该问题直接求解，但由于作为优化变量的协态变量没有物理意义，且算法的收敛域较小，因此，对于推力值小、转移圈数多的问题，求解过程存在很大困难。

为了增加优化算法的收敛域，降低求解难度，同伦方法被广泛应用于轨道优化。同伦方法，是指从较简单问题出发，通过改变同伦参数，求解对应的一系列子问题，最终得到原问题的解。其在小推力在转移轨道中的应用主要有两个方向：一是针对燃料最优问题，通过对性能指标进行同伦，先求解连续光滑控制问题，然后通过同伦迭代，最终得到不连续的Bang-Bang控制的燃料最优转移轨道[10-11]；二是针对推力很小的轨道优化问题，先求解较简单的大推力值对应的轨道优化问题，然后通过对推力幅值进行同伦，逐渐减小推力幅值，最终得到所需的很小推力幅值的转移轨道[12-13]。

目前对限制性三体问题中的小推力转移轨道优化问题研究较少，Caillau等[12]对从航天器从GEO轨道到地月系L1点和月球轨道的转移轨道进行相关研究，本文将其扩展到了航天器从地球GEO轨道到地月系L2点的时间最优转移轨道优化问题。首先基于最优控制理论，将地月系三体问题模型下的小推力转移轨道问题转换为两点边值问题，然后结合同伦方法对该两点边值问题进行求解。

2 动力学模型

考虑航天器同时受到地球引力和月球引力，忽略地球非球形摄动、太阳引力、太阳光压等其他摄动，并假设月球围绕地球的运行轨道为圆轨道，则可用圆限制性三体问题来描述航天器的运动。仅考虑航天器在月球轨道面内的平面运动，并且忽略转移过程中航天器的质量变化，则航天器在旋转坐标系下的动力学模型为式(1)～(2)：

(1)

(2)

式中：r=[x,y]T和v=[vx,vy]T分别为航天器的位置矢量和速度矢量，m0为航天器初始质量，Tmax为航天器可以提供的最大推力，u为转移过程中的实际推力大小，α=[αx,αy]T为单位推力方向。则可定义航天器有效势能函数U如式(3)，g(r)和h(v)分别表示为式(4)、(5)：

(3)

(4)

h(v)=[2vy,-2vx]T

(5)

需要注意的是，在动力学模型式(1)～(2)中，各变量均为单位化后的值，其中长度归一化单位为1 DU=3.8440×105km，时间归一化单位为1 TU=3.75676967×105s，质量归一化单位为1 MU=M0kg(M0为航天器实际质量)。归一化后，地球在旋转坐标系的位置为(-μ,0)，月球的位置为(1-μ,0)。

本文针对航天器从地球GEO轨道到L2平动点的时间最优小推力转移轨道进行研究。在初始时刻t0，航天器在GEO轨道上的位置自由，其状态量可以表示为式(6)：

x0=[r0cosθ0-μ,r0sinθ0,-v0sinθ0,v0cosθ0]T

(6)

时间最优问题的性能指标可以表示为式(7)：

(7)

时间最优小推力转移轨道优化问题的本质，就是在满足初始状态约束式(6)和终端状态约束xf的条件下，通过优化变量u和α，使性能指标最小。根据庞德里亚金极小值原理，引入拉格朗日乘子，构造哈密尔顿函数为式(8)：

(8)

式中：λr=[λx,λy]T和λv=[λvx,λvy]T也称为协态变量，其微分方程为式(9)～(10)：

(9)

(10)

为使哈密尔顿函数为极小值，从式(8)可知，推力大小和推力方向的最优控制律应为式(11)～(12)：

u*=1

(11)

(12)

对于航天器初始时刻的状态量x0，由于存在优化变量θ0，需要对横截条件进行推导。根据最优控制理论，可以得到对应的横截条件为式(13)：

(x0+μ)λy0-y0λx0+vx0λvy0-vy0λvx0=0

(13)

对于时间最优问题，转移时间自由，对应的横截条件为式(14)：

H(tf)=0

(14)

至此，将最优控制问题转换为两点边值问题，其初始时刻的状态值由式(6)确定，终端时刻的约束包括状态约束xf和横截条件式(13)～(14)。该两点边值问题可利用打靶法进行求解，打靶函数表示为式(15)：

F(y)=O

(15)

式中y=[θ0,λx0,λy0,λvx0,λvy0,tf]T为优化变量。

3 推力幅值的同伦方法

将时间最优小推力转移轨道优化问题转换为两点边值问题后，虽然理论上可以利用打靶法对式(15)进行求解得到优化变量y的值，再对初值进行积分得到转移轨道。然而推力越小，转移时间越长，转移圈数越多，打靶过程对初始猜测值就越敏感，打靶的收敛域就越小。此外，相比于传统的只考虑中心引力场的二体模型，限制性三体模型的非线性更强，从而增加了求解难度。因此，为了克服小推力转移轨道优化问题的初始值敏感问题，本文利用同伦方法进行求解。

3.1 构造同伦函数

基于所需求解的原问题，构造容易求解的问题作为同伦初始问题，而同伦过程中具有全局收敛性，因此同伦方法可以有效增大打靶过程收敛域，减小求解难度。同伦函数的一般构造形式可以表示为式(16)：

H(y,κ)=κF(y)+(1-κ)G(y)

(16)

式中：y∈n为优化变量，κ为同伦参数，F:n→n和G:n→n均为光滑的映射函数，分别表示所需求解的原问题和构造的同伦初始问题，H:n×→n为构造得到的同伦函数。根据初始问题G的不同，一般可以将同伦方法划分为牛顿同伦(Newton Homotopy)，定点同伦(Fixed-point Homotopy)和尺度不变仿射同伦(Scale-invariant Affine Homotopy)[14]。针对本文的时间最优小推力转移轨道优化问题，在推力幅值中引入同伦参数κ，构造的同伦函数的动力学模型为式(17)～(18)：

(17)

(18)

式中：TL≫Tmax为足够大的推力值，在该推力值下，转移轨道优化问题较容易直接求解，Tmax为所需求解问题的实际推力幅值。对应于式(16)中的一般同伦形式，F为所需求解的推力值为Tmax的原问题，G为容易求解的推力值为TL的初始同伦问题。在时间最优的性能指标下，得到的哈密尔顿函数为式(19)：

(19)

协态变量的微分方程保持式(9)～(10)的形式不变，同时最优推力大小u*和最优推力方向α*和式(11)～(12)保持一致。边界约束条件和横截条件也与原问题保持不变。至此，同伦问题转换成一系列随着κ值变化的两点边值问题，通过求解这一系列的同伦子问题，最终可得到原问题的解。

3.2 同伦曲线跟踪方法

同伦方法除了构造合理的同伦函数，还需要选取合适的同伦曲线跟踪方法。同伦曲线的跟踪方法主要可以分为离散同伦和连续同伦。离散同伦可以描述为，将同伦参数从0到1划分为m个离散节点，即0=κ1<κ2<…<κm-1<κm=1，先求解κ1对应的初始问题的解，然后将得到的解y1作为κ2对应的同伦子问题的初始猜测值，求解得到该问题的解y2，然后再将y2作为κ3对应的子问题的初始猜测值，依次计算，最终得到κm对应的原问题的解ym[7]。离散同伦简单直观，但若相邻两个κ节点之间距离较远，则以前一个节点的解作为当前节点子问题的猜测值时收敛性不能保证；更重要的是，当同伦曲线存在拐点时，离散同伦将在拐点处失效。

本文利用连续同伦中的伪弧长法跟踪同伦路径[15]。不同于离散同伦根据同伦参数κ来选取同伦节点，伪弧长法沿着同伦曲线的切线方程，根据伪弧长的步长Δs进行迭代计算。伪弧长法是一种预测-校正方法。首先在当前节点计算表征梯度信息的雅克比矩阵，将步长Δs分配到各个优化变量上，得到下一个节点的信息，然后再利用牛顿法等校正方法对下一个节点的值进行修正，从而得到该节点所对应的解[15]。如图1所示，假设在第i个节点，(κi,yi)为非线性方程组H(κ,y)=O的已知解，定义该点处的单位切向量为式(20)：

(20)

图1 伪弧长法示意图Fig.1 Geometric interpretation of pseudo arclength continuation method

(21)

H(κ,y)=O

(22)

由于κ值的变化方向由同伦曲线的切线方法确定，在同伦过程中κ可以增加，也可以减小，因此可以用于跟踪存在拐点的同伦路径。

4 数值仿真验证

本文的仿真算例为地月系平面限制性三体问题下航天器从GEO轨道到L2的时间最优小推力转移轨道，航天器的初始质量为1500 kg，发动机最大推力值为Tmax=1 N，对应的初始推重比为6.8027×10-5。航天器初始位置由式(6)确定，终端状态约束为航天器终端时刻状态量与L2点重合，GEO轨道半径为r0=0.109689855932071 DU，航天器在GEO轨道上的初始速度为v0=3.000969693845573 DU/TU。

首先选取较大的推力值为TL=10 N，此时利用打靶法，可得初始问题的解为y0=[-7.4762,-3.7986,6.9355,-0.3092,-0.0844,1.4413]T。对应的转移轨道如图2所示，其中蓝色实心圆点表示地球位置，红色实心圆点表示月球位置，航天器到达L2点所需的转移轨道圈数为1圈。

图2 TL=10 N时从GEO到L2点的转移轨道Fig.2 Transfer trajectory from GEO to L2 libration point with TL=10 N

以初始解y0为同伦初始点，利用伪弧长法跟踪同伦轨迹，可以得到从κ=0到κ=1的同伦路径，其中优化变量中的转移时间tf的同伦路径如图3所示，初始时刻的协态变量λ0=[λx0,λy0,λvx0,λvy0]T如图4所示，表示初始时刻航天器在GEO轨道上的位置的变量θ0的同伦曲线如图5所示。从图中可以看出，该同伦曲线并不是随着κ单调递增的，而是存在多个拐点。在同伦初始阶段，κ值沿着同伦路径增加，但同伦曲线在κ=0.8588时遇到第一个拐点，改变了同伦路径的方向，使κ值沿着同伦路径减小，直到在κ=0.8408时遇到第二个拐点，使κ值重新沿着同伦路径增加，如此往复，直到同伦路径到达κ=1，最终得到Tmax=1 N的转移轨道。

图3 转移时间tf的同伦曲线Fig.3 Homotopy curve of the transfer time tf(κ)

图4 初始时刻协态变量λ0的同伦曲线Fig.4 Homotopy curves of costate λ0(κ)

图5 表示航天器在GEO上的位置的θ0的同伦曲线Fig.5 Homotopy curve of θ0(κ)

值得注意的是，同伦路径上的每一个点都对应该点κ值下的一个解，即在对应推力值下，满足打靶方程的一条时间最优转移轨道。这是因为间接法的最优控制律和横截条件是基于一阶最优性条件推到得到的，可以保证满足这些约束条件的解至少是极值点，即局部最优解，但不能保证是全局最优解。而在轨道转移优化问题中，特别是转移圈数较多的问题，一般都存在多个局部最优解。从同伦曲线中可以看出，在第二个拐点κ=0.8408之后，每个κ值对应的推力下都存在多个解。以κ=0.90为例，其对应的推力值为1.9 N，在这条同伦路径中共存在13个解，对应的转移时间范围为5.0480～ 6.1190 TU，通过对性能指标即转移时间进行对比，即可在这些局部最优解中选取得到最优解。对于原问题中的Tmax=1 N的转移轨道，虽然在这条同伦路径中只存在1个局部解，但可以预见的是，如果继续跟踪同伦路劲，随着同伦过程的进行，将会得到更多的局部解。根据tf的同伦路径的变化规律，当前得到的这个解，对应的转移时间将始终小于其他解，即可认为当前得到的解是最优解。

从图5中可知，与图2对应的10 N时的转移轨道相比，Tmax=1 N的转移圈数增加了9圈。该解对应的转移轨道如图6所示，所需转移时间为8.7113 DU，即37.8777天。转移过程中航天器发动机的推力幅值始终保持最大值1 N，推力方向α*的变化情况如图7所示。

图6 Tmax=1 N时的从GEO到L2点的转移轨道Fig.6 Transfer trajectory from GEO to L2 libration point with Tmax=1 N

图7 Tmax=1 N的转移轨道推力方向α*随时间的变化曲线Fig.7 Variation of thrust direction α* along the transfer trajectory with Tmax=1 N

5 结论

1) 在转移轨道优化问题中，针对推力大小进行同伦时，同伦曲线存在拐点，离散同伦不能得到完整的同伦曲线，因此需要选取连续同伦轨迹跟踪方法，利用伪弧长法可以跟踪得到连续同伦路径，得到包括拐点在内的完整同伦信息，有利于对同伦过程及同伦结果进行分析；

2) 对于多圈的小推力转移轨道优化问题，可能存在多个满足最优性一阶必要条件的局部最优解，需要对性能指标进行对比才能确定最优解。若直接对小推力转移轨道进行求解，并得到其中一个解，不能从该解本身判断其是否为全局最优解，但如果利用连续同伦法跟踪完整同伦路径得到小推力转移轨道，则可以依据该解在同伦路径的位置及同伦路径变化规律来判断该解是否为全局最优解。