刘柯邑，岳晓奎，张 莹

（西北工业大学 航天学院，陕西 西安 710072）

小推力转移轨道设计与优化问题逐步成为研究的热点[1-2]，与传统的火箭发动机相比，小推力发动机具有比冲高、质量轻的特点，可有效减少燃料消耗，提高任务的有效载荷，降低发射成本。但由于其推力小（通常为mN量级）、作用时间长导致小推力轨道的动力学非线性严重，这给小推力轨道的设计带来了难题。小推力转移轨道设计问题可归结为求解有约束的连续最优控制问题。学者们针对该问题所采用的方法包括间接法和直接法[3-4]。

间接法基于变分法将最优控制问题最终转换为一个两点边值问题，但这种方法收敛半径小，对协状态的初始猜测值很敏感，对于复杂约束难以处理。相对间接法，直接法在收敛的鲁棒性和解决实际问题的适应性上具有优势。直接法是将最优控制问题离散成参数优化问题，然后利用非线性规划方法求解。高斯伪谱方法[5-6]是最近研究较多的直接配点优化方法，相对于其他直接配点法[7]的优势是可以用较少的节点获得较高的精度。麻省理工学院的Benson从理论上证明了高斯伪谱法的KKT条件与最优控制理论中的一阶必要条件是一致的；并且，由于其无需猜测协状态的初始值，这就大大降低了求解最优控制问题的难度。

文中研究了基于高斯伪谱法的燃料最省小推力转移轨道设计问题。首先建立了星际小推力轨道优化模型，对转移轨道参数进行了无量纲化处理以提高数值计算精度。然后利用高斯伪谱方法将转移轨道优化问题转化为多约束的非线性规划问题，并应用SNOPT优化软件包求解了最优转移轨道的推力控制角变化规律。

1 小推力轨道优化问题

只考虑探测器在星际巡航段的转移轨道设计问题，采用笛卡尔坐标系下的小推力探测器动力学模型：

其中，r=[x，y，z]T为探测器在日心黄道惯性系中的位置矢量，v=[vx，vy，vz]T为发动机推力在各个方向上的加速度矢量，m为飞行器和推进剂的总质量； 为太阳引力系数； 为推力幅值；α=[αx，αy，αz]T为发动机推力在各个方向的加速度矢量；g0为海平面重力加速度；Isp为发动机比冲。

假定小推力发动机的输入功率P随着探测器日心距离的增大而减小，它与探测器日心距离的平方成反比，即

其中，P0为探测器日心距离为1AU时的输入功率，r采用天文单位。发动机推力T为：

其中，η为推进器的工作效率。

在推进器连续工作条件下，探测器从地球发射与目标星体的燃料最省转移轨道的性能指标可以写成：

探测器的发射时间和到达时间均不固定，对于任何的发射时间t0，探测器受到地球星历约束，探测器的位置和速度需要满足下面的等式约束

其中，rsc，vsc是发射时刻探测器在日心黄道坐标系中的位置和速度；rE，vE是地球的位置和速度；v∞是地球逃逸双曲线超速，其大小已知，方向不确定。由公式（4），（5）构成了探测器轨道的初始边界条件：

当探测器与目标星体交会时，探测器的位置和速度需要满足下面的等式约束：

其中，tf是探测器与目标星体的交会时间；rt（tf），vt（tf）是探测器交会时目标星体的位置和速度。由公式（7）、（8）构成了探测器轨道的末端边界条件：

在转移过程中，推力方向矢量必须满足如下条件：

则星际探测任务的燃料最省转移轨道可以描述为：寻求最优的发射时间t0，连续的轨道状态 xsc（t），连续的推力方向矢量α（t）和交会时间tf，使探测器轨道在满足动力学约束（1），边界约束（6）和（9），推力方向矢量约束（10）的条件下，使性能指标（3）达到最小。

寻优参数可以表示为：

为提高计算精度，避免利用高斯伪谱方法求解过程中出现病态矩阵，需要对轨道参数进行无量纲化处理，使各状态参数的取值范围都在一个相近的数值范围内。为此，文中取量纲质量单位[M]=m0；量纲距离单位[L]=AU；量纲时间单位[T]=。下面将采用高斯伪谱法对上述小推力轨道优化问题进行处理。

2 高斯伪谱方法

高斯伪谱法的基本求解思路为：将未知的轨道状态和控制时间历程在一系列高斯点上离散化，而后用这些离散的状态与控制分别构造Lagrange插值多项式去逼近真实的轨道状态与控制，再通过对状态量求导来代替动力学微分方程。

由于高斯点分布在（-1，1）区间，若要把动力学方程在高斯点上进行离散，首先需要把所研究问题的时间区间t∈[t0，tf]转换到τ∈[-1，1]，这个转化可以通过下式完成：

高斯伪谱法取 N 个高斯点 τ1，τ2，…，τN和初始端点 τ0=-1上的离散状态构造Lagrange插值多项式去近似状态的时间历程，即

式中，x（τ）为真实的状态时间历程，X（τ）为由 Lagrange插值多项式近似得到的状态时间历程；Li（τ）为Lagrange插值基函数，i=0，1，…，N。

类似地，控制量的时间历程也可以用离散的控制构造Lagrange插值多项式来近似，即

式中，u（τ）为真实的控制时间历程，U（τ）为由 Lagrange插值多项式近似得到的控制时间历程；（τ）为Lagrange插值基函数，i=1，…，N。

需要注意的是，构造状态用了N个高斯点与初始端点共N+1个点，而构造控制只用了N个高斯点，因而方程（13）的脚标是从0开始，而方程（15）的脚标却从1开始。

式（13）对时间求导，得

定义微分近似矩阵D，其中元素为：

式中，D 为 N×（N+1）二维矩阵，k=1，…，N，i=0，1，…，N。

由式（17）、（18）将轨道动力学约束转化为一系列代数约束：

式中，Xk=X（τk），Uk=U（τk）；k=1，2，…，N， f为轨道动力学方程右端函数。

终端状态由高斯求积公式得到：

式中，X0=X（-1），ωk为高斯积分权重，有

式（19）将轨道动力学约束转化为只与节点处离散状态和控制变量相关的代数约束，同样，性能指标和约束条件也可以处理为只与节点处状态和控制变量相关的形式。

如此，星际小推力轨道优化问题可以转化为一个多约束的参数优化问题，寻优参数为发射时间、交会时间以及各节点上的状态和控制变量：

性能指标可以表示为：

始末端边界条件可以表示为：

推力矢量约束可以表示为：

离散后的小推力轨道最优控制问题为典型的非线性规划问题，本文采用TOMLAB/PROPT优化软件包进行求解。

3 数值仿真与分析

为了验证本文所提设计方法的正确性和有效性，对小推力探测器的星际转移轨道进行设计与优化。探测器主要参数如表1所示。

表1 探测器主要参数表Tab.1 Parameters of space detector

算例一：地球-火星交会任务

采用直接打靶法和高斯伪谱法对地球-火星交会任务的燃料最省小推力转移轨道分别进行设计。

探测任务要求探测器在2009年1月1日到2009年12月31日之间从地球出发，出发时刻认为探测器的日心位置和速度与地球相同。地球-火星交会的初始转移时间设为365天。采用直接打靶法和高斯伪谱法优化的轨道参数如表2所示。

由表2可以看出，采用直接打靶法和高斯伪谱法的优化结果基本一致，这表明采用高斯伪谱法求解行星际小推力转移轨道问题是有效地。采用直接打靶法求解该问题时，所用计算时间为1 712.343 8 s；采用高斯伪谱法时计算时间为56.140 6 s，这主要是由于直接打靶法在迭代过程中需要对状态动力学方程进行积分造成的，同时也验证了高斯伪谱法在求解小推力转移轨道问题上的优势。

表2 地球-火星轨道优化结果对比Tab.2 Orbit optimization result of earth-mars

对于小推力轨道设计问题，通常采用俯仰和偏航两个控制角来定义推力的指向。定义俯仰控制角α为推力矢量在轨道平面内的投影与当地水平之间的夹角；偏航控制角β为推力矢量与密切轨道平面之间的夹角，则推力的单位方向矢量就可以表示为：a=[sinαcos β cosαcos β sin β]。 直接打靶法和高斯伪谱法的控制角随时间的变化规律如图1所示。

由于火星轨道面与地球轨道面之间的夹角很小，偏航控制角在整个飞行过程中均较小，最大时约为20°。探测器的推力方向随飞行轨迹的变化如图2所示。

图1 直接打靶法和高斯伪谱法控制轨线对比Fig.1 Compare of orbital control trajectory

图2 探测器飞行轨迹与推力方向Fig.2 Trajectory and thrust orientation

算例二：地球-金星交会任务

为了进一步考察高斯伪谱法求解行小推力轨道优化问题的有效性，采用高斯伪谱法对地球-金星交会任务的燃料最省小推力转移轨道进行设计。探测任务要求探测器在2021年1月1日到2021年12月31日之间从地球出发，出发时刻认为探测器的日心位置和速度与地球相同，即探测器的地球逃逸速度为零。初始假设探测器飞行 天到达金星，此时要求探测器的日心速度等于金星的日心速度，与金星进行交会。采用高斯伪谱法优化的轨道参数如表3所示。

表3 地球-金星轨道优化结果Tab.3 Orbit optimization result of earth-mars

表2和表3中转移时间和剩余质量的优化结果对比，说明了发动机开机时间长短对太阳能电推进器的燃料消耗的影响不大；由于发动机推力大小与探测器日心距离的平方成反比，相比于地球-火星段轨道，探测器地球-金星轨道段的推力比较大，所以剩余质量相差不大。

探测器在地球-金星交会任务中的飞行轨迹如图3所示。

图3 探测器飞行轨迹Fig.3 Trajectory of space detector

图4 地球-金星轨道控制轨线Fig.4 Control orbital trajectory of earth-venus

从飞行轨迹图3可以看出，在转移轨道初始阶段，发动机推力方向主要指向轨道内以减小轨道的半长轴；轨道中段，推力方向主要指向轨道面外以增加轨道倾角；在探测器接近金星的最后阶段，推力方向主要指向轨道外以降低探测器的速度，避免撞击金星。探测器在地球-金星交会任务中的发动机推力矢量控制角随时间的变化如图4所示。由于金星轨道面与地球轨道面之间的夹角相对较大，在整个飞行过程中，偏航控制角最大时约为。

4 结 论

文中研究了基于高斯伪谱法的小推力转移轨道设计与优化问题。通过对燃料消耗最省火星和金星交会任务的轨道设计与分析表明：采用高斯伪谱法进行求解时，仿真初值选取相对自由，并且收敛速度快，采用此优化算法，可保证所得解最优的基础上，有效提高问题的计算效率，具有良好的应用前景。

[1]Rayman M D，Williams S N.Design of the first interplanetary solar electric propulsion mission[J].Journal of Spacecraft and Rockets，2012，39（4）:589-595.

[2]Kugelberge J，Bodin P，Persson S，et al.Accommodating electric propulsion on SMART-1[J].Acta Astronautica，2004，55（2）:121-130.

[3]Ranieri C L，Ocampo C A.Indirect optimization of spiral trajectories[J].Journal of Guidance， Control and Dynamics，2006，29（6）:1360-1366.

[4]Betts J T，Erb S O.Optimal low thrust trajectories to the moon[J].SIAM Journal of Applied Dynamical Systems，2003，2 （2）:144-170.

[5]Benson D A.A gauss pseudospectral transcription for optimal control[D].Massachusetts:Department of Aeronautics and Astronautics， Massachusetts Institute of Technology，2005.

[6]Benson A，Thorvaldsen T，Rao V.Direct trajectory optimization and costae estimation via an orthogonal collocation method[J].Journal of Guidance， Control and Dynamics，2006，29 （6）:1435-1440.

[7]Betts J T.Very low-thrust trajectory optimization using a direct SQP method[J].Journal of Computation&Applied Mathmatics，2009，120（1）:27-40.