变比冲连续小推力轨迹优化方法综述*

迟哲敏，李俊峰，蒋方华，宝音贺西

(清华大学 航天航空学院·北京·100084)

0 引 言

目前，各国正在积极开展深空探测任务，对具有科学价值的天体的探测已经成为深空探测领域的研究热题。我国计划在2020年发射第一个火星探测器“天问”一号，同时也筹划着未来的小行星探测任务。远距离探测活动的开展离不开推进系统的选择。1998年NASA成功发射的“深空”一号首次验证了电推进技术在深空探测领域的应用。由于电推进的比冲是化学推进的10倍左右，所以电推进在减少燃料消耗、增加有效载荷质量方面更具有优势。电推进相较于化学推进可以提供“连续”的小推力，因此发动机的寿命较长，可以工作15000小时以上，适合完成长距离、长时间的飞行任务。

常见的采用电推进方式的发动机有三种，首先是美国艾德阿斯特拉火箭公司在研的磁等离子体火箭(VASIMR)[1]，它在恒定的大功率下工作，通过控制比冲来调节推力等参量的大小，计划采用核能作为其动力来源，因此有利用VASIMR完成30天到达火星的构想。另外两种电推进发动机——霍尔效应推进器和离子电推进发动机均已经应用于深空探测任务中。霍尔效应推进器利用基于电磁角力的霍尔效应原理来提供飞行动力[2]。2003年欧空局发射的第一个月球探测器“智慧”一号(SMART-1)就是将霍尔效应推进器作为主推进器，成功完成了月球探测任务[3]。而离子电推进发动机是在高温下将工质电离，再令离子在电场中加速喷出从而产生推力，目前普遍选择氙气作为工质。近几年，离子电推进发动机越来越广泛地应用于深空探测任务中，NASA先后研制了NSTAR(NASA Solar Technology Application Readiness)和NEXT(NASA Evolutionary Xenon Thruster)离子发动机，NEXT发动机在推进效能和燃料节省方面均优于NSTAR发动机，表1展示的是NSTAR和NEXT发动机的部分性能参数[4]。“深空”一号装载一个NSTAR 发动机[5]，发动机工作的功率由太阳能帆板提供。“黎明”号携带3个NSTAR发动机先后探测了灶神星和谷神星[6]，成为了第一个环绕2颗地球外天体的探测器。NEXT发动机也在2010年测试成功，原定被应用于2024年探测彗星67P的“凯撒”号任务中，然而任务竞标失败。它的商业化版NEXT-C计划于2021年作为“双小行星重定向测试”任务(Double Asteroid Redirection Test, DART)[7]的主推进器。JAXA的小行星采样任务Hayabusa和Hayabusa2均采用μ10离子电推进发动机[8]，Hayabusa在2005年11月完成了对小行星Itokawa的采样任务，而Hayabusa2也于2019年11月在小行星Ryugu上进行采样，Hayabusa2所采用的μ10推进器是Hayabusa任务推进器的升级版[9]，它们均是变推力变比冲发动机。英国的T6离子发动机[10]正应用于在飞的Bepicolombo任务中，Bepicolombo借助一次地球、两次金星和六次水星引力辅助计划于2025年到达水星。

表1 NSTAR推进器和NEXT推进器的性能特征对比

上述电推进发动机的比冲均是变化的，因此变比冲电推进模型是更符合工程实际的推进系统模型。而对于变比冲连续小推力轨迹优化的研究是具有工程应用意义的。由于连续小推力自身的特点和太阳能提供功率的限制，以及考虑多目标或多次引力辅助探测问题而引起的多内点约束，变比冲电推进探测任务的设计仍面临着困难和挑战。传统的直接法和间接法也可以适用于解决变比冲电推进轨迹优化问题，同时也需要针对具体模型和具体问题进行方法上的改进。本文对目前常见的变比冲电推进模型以及常用的变比冲连续小推力轨迹优化方法的研究现状进行了相关介绍。

1 变比冲电推进模型

目前，在深空轨迹计算中，比较常用的变比冲电推进简化模型有三种，分别对应着不同真实发动机系统的工作原理和方式。

首先是比冲在区间内变化的推进模型。在这一模型下，将比冲作为控制量，控制其在区间范围内变化，优点是若在需要减小燃料消耗的时候，令比冲取范围中的较大或是最大值，若需要在一个位置加速，那么又会取比冲范围中相应的较小值，根据不同情况可以灵活调节[11]。最符合这一模型的是磁等离子体电推进发动机，因为考虑采用核能供电，所以发动机输入功率是恒定的，比冲是可控的。Seywald等[12]就是构建了VASIMR发动机模型来求解共面圆轨道之间的燃料最优轨迹，Park等[13]也是应用VASIMR模型求解到达越地小行星的交会轨迹，为地球防御任务的实施提供参考；Kechichian[14]也应用了恒定功率变比冲模型，在他的研究中，将春分点轨道根数代替位置速度来表示用于优化的状态量。然而也存在另一种功率模型假设：在太阳能电推进模式下，功率是星日距离的函数，但比冲仍是在区间内变化并可被当作控制量[15-17]。Casalino等[18]在此变比冲模型的基础上构建了双比冲模型，即比冲只能取范围的最大和最小值，再对比了两种模型下的燃料消耗；他们还将变比冲模型用于求解从近地轨道向地球同步轨道转移的多圈问题[19]。Mengali等[20]在变功率变比冲模型的基础上还考虑了推进效率的影响，推进效率是比冲的多项式函数，根据拟合程度有线性、二次、三次等函数关系，陈杨[21]在博士论文中将其中的线性关系用于行星际间的轨迹优化。Taheri等[22]在研究中考虑了多模式的推进系统模型，在他们的模型构建中，分别将发动机输入功率和比冲作为控制量来求解燃料最优轨迹。

另一种是比冲随发动机输入功率变化的数学模型。这种情况下，发动机是由太阳能供电，因此比冲和推力也是星日距离的函数。此模型是最贴合光栅离子电推进发动机工作原理的，发动机的输入功率作用于两部分：等离子体的产生和离子的加速。等离子体生成器是在等离子体生成功率的作用下工作的，在燃烧室中，燃料中的中性原子在电子轰击下电离，生成的离子在一组光栅中静电加速，这组光栅的控制电压为高电压，相对应的功率取决于高电压和静电加速的电流。一旦给定一个高电压，性能参数：推力和比冲，确定为是发动机输入功率的函数。图1是NSTAR光栅离子发动机的原理图[23]。

图1 NSTAR离子发动机[23]

发动机试验得出的结果是一系列相关的推力随功率以及比冲随功率变化的离散点(参见图2)，因此相应的函数关系是通过曲线拟合得到，最常见的是多项式拟合。需要解释的是，一般通过多项式拟合的是推力和质量流率对输入功率的函数关系，而比冲是通过推力和质量流率计算得到。William等[24]在1997年已经对将应用于“深空”一号探测任务的NSTAR离子电推进发动机的参数进行线性拟合和评估，Petropoulos等[25-26]采用此线性拟合模型进行小推力引力辅助轨迹设计。William等[27]在火星探测轨迹的求解中应用了NSTAR的质量流率是功率的二次多项式拟合函数模型。更精确一些的NSTAR和NEXT发动机的四次多项式拟合曲线[28]也用于求解行星际转移的燃料最优轨迹。David等[29]分别根据用于NASA的五个深空探测任务的NEXT发动机的参数离散点给出相应的“大推力”和“高比冲”四次多项式拟合模型。Saripalli等[30]分别拟合了不同电压下的、“大推力”和“高比冲”情况下的推力-功率以及比冲-功率四次曲线，并且将它们用于任务计算来对比各种模型得到的轨迹的最优性。一系列霍尔效应发动机的参量拟合关系也用于轨迹优化的计算中[31-34]。Casalino等[35]在电推进三次曲线拟合模型下研究了功率作为控制量的多模式推进系统的轨迹优化问题。国内的孙冲等[36]应用了太阳帆电推进的二次多项式拟合模型进行虚拟引力场的轨道设计。除此之外，Karthik等[37]对微小卫星上的电推进系统在考虑了电压电流等各种因素的影响下进行了精确建模。

(a)推力变化图

最后是功率离散点模型或称为功率分档模型。这种模型对应两种情况：首先是一个发动机具有多个工况点，推力和比冲随着功率分档变化；第二种情况是一个推进系统由多个恒定推力恒定比冲的发动机组成，发动机工作的功率由太阳能帆板提供，所以随着功率减小，发动机依次关闭，也就导致推力随着功率分档变化[38]。Brophy等[39]给出了NSTAR发动机所有工作周期内的6个工况点数据，Chi等[38]应用此6个离散点的发动机模型求解火星探测的时间最优与燃料最优轨迹。Polk等[40]提供了用于“深空”一号任务的NSTAR的16个工况点试验数据，并且将它们用于任务轨迹的验证和评估。Quarta等[41-42]应用NEXT发动机的41个工况离散点[43]求解行星际转移的燃料最优轨迹，他们对代入离散点数据(功率、推力和比冲)的哈密顿函数进行排序来优化轨迹。在不考虑功率和比冲的情况下，Axelrod等[44]研究了带有推力离散档位的发动机模型下的轨迹优化方法，该方法既适用于多推力档位的推进系统又适用于多个单向发动机组成的推进系统。

2 变比冲小推力轨迹优化方法

2.1 轨迹优化方法

变比冲电推进轨迹优化问题的本质是连续小推力轨迹优化问题。这里的“连续”和“小”是相对于脉冲化学推进而言的。常见的求解连续小推力最优轨迹的数值方法分为两种：一是将问题进行转换，比如直接法，将原先连续控制问题转化为参数优化问题；第二种是针对问题的特点给出方法，连续小推力轨迹优化问题等效于泛函极值问题，所以在此可应用基于变分法和最优控制原理的间接法。下面分别介绍一下直接法和间接法的基本原理和特点。

2.1.1 直接法

直接法实际上是将连续控制问题转化为参数优化问题。对于包含状态量、控制量、过程约束和性能指标的小推力轨迹优化问题，直接法是将状态量和(或)控制量离散，并将状态约束条件和性能指标写为数值微分或是数值积分的代数形式，所以原问题转化为求解多个离散后序列的非线性规划问题(Nonlinear Programming, NLP)[45]。常见的求解非线性规划问题的方法是序列二次规划算法(Sequence Quadratic Program, SQP)，常用的SNOPT软件就是基于此算法来求解大规模的非线性规划问题[46]。

通过上述描述，直接法包含一个离散的过程，这就涉及到两个问题，首先是离散的过程中并不能确保状态量和控制量的匹配，也不能保证得到的解可以满足最优控制的一阶必要条件，所以直接法不能保证解的最优性。再者为了提高求解的精度，需要增加离散点数量，会导致计算量增大，相应地也就会降低计算效率。但直接法也有很多优点，比如离散过程简单，解的收敛性高等。

直接法包括几种常用的方法，比如直接配点法和直接打靶法。它们的区别体现在前者将动态约束条件写为数值微分(有限差分或拉格朗日插值)的形式，后者是通过数值积分的形式来实现。伪谱法就是直接配点法的一种[47-49]，它是利用基于全局近似的拉格朗日插值的数值微分方法重构动力学方程。伪谱法可以改善上一段所叙述的应用直接法存在的问题，首先伪谱法自身构造的一阶必要条件与最优控制的一阶必要条件等价[50]，所以可以提高解的最优性；其次是在相同数量离散点的情况下，伪谱法相比其他直接法的求解精度更高[51]。Rao等[52-53]基于伪谱法和自适应离散网格法创建了MATLAB的GPOPS(Gauss Pseudospectral Optimization Software)软件包，目前该软件在轨迹优化中广泛应用。伪谱法还有一个优势在于，由于一阶必要条件的等价性，所以伪谱法的解可以为间接法求解提供协态初值猜测[54]。但是利用伪谱法求解得到的Bang-Bang控制曲线存在数值振荡等问题，所以结果的精度还是比间接法弱一些，Bai等[55]和郭铁丁等[50,56]对此提出了相应的改进方法。

另外，基于直接法原理还开发了其他常用的程序包，比如MALTO、EMTG等，在后文的叙述中会对它们做简单介绍。

2.1.2 间接法

应用间接法的求解思路是构造对应于状态变量的协态变量，在性能指标的基础上建立哈密顿函数，再推导欧拉-拉格朗日方程(协态微分方程)。同时积分状态和协态微分方程，在边界点处打靶，打靶方程包括状态约束以及与协态变量相关的一阶必要条件。通过庞得里亚金极大或是极小值原理推导得到最优控制律。因此，“间接”是通过构造辅助优化的协态变量来体现的。

间接法是将最优控制问题转化为两点边值问题(Two-Point Boundary-Value Problem, TPBVP)或多点边值问题(Multi-Point Boundary-Value Problem, MPBVP)，MPBVP对应于多内点约束问题，比如多目标交会或是多次引力辅助问题。其本质都是求解非线性方程组，在此可以应用牛顿迭代法[57-58]等。Jiang等[59]采用Fortran里面的Minpack-1程序包求解TPBVP和MPBVP，Minpack-1是基于结合了牛顿法和梯度法的Powell混合法的非线性求解器[60]。Zhang等[61]通过构造解析的雅各比矩阵来估算打靶过程的梯度信息，以此提高牛顿迭代的收敛效率。

间接法可以有效克服应用直接法产生的困难，比如采用间接法得到的解可以满足一阶必要条件，至少能够保证是局部最优的，再者因为无需进行多参数离散求解，计算效率会相应地提高。但是因为引入协态变量且协态变量自身物理意义不明确，所以无法解析给定，如果对其猜测不准确，那么将一直无法收敛到最优解。对于多内点约束问题，随着协态初值数量的增多，收敛难度增大。因此应用间接法的最大困难是由于引入协态变量导致的较低的收敛效率。目前采用协态初值估计[59,62-63]、同伦方法[15,59,64-66]以及开关检测[59,61,67]等一系列策略来提高间接法求解的计算精度和收敛效率。

最后简单介绍一下混合法，它是直接法和间接法的结合，用直接法离散协态变量，通过间接法推导最优控制律，最后利用非线性规划算法对问题进行求解。混合法克服了直接法和间接法的一些弊端，比如解的最优性难以保证，协态初值猜测困难等，但是因为计算参数较多，求解效率仍不够高。Kluever和Gao等[68-71]提出并改进了混合法用于求解行星际间转移的轨迹优化问题。

2.2 比冲区间变化的轨迹优化方法

对于比冲区间变化的模型下的轨迹优化问题，在所查找到的文献中基本都应用间接法来处理[12-15,17-20,72]。应用间接法在处理区间变比冲与恒定比冲轨迹优化问题的不同之处在于，在此将比冲作为优化过程中的控制量，因此需要用庞得里亚金极大值原理推导关于比冲的最优控制律。在使用间接法的基础上辅助有其他策略，可以提高求解效率。比如对于VASIMR恒定功率变比冲模型，Seywald等[12]在求解共面圆轨道之间的最优转移轨迹时，先假设飞行器的最优推力方向始终与速度方向同向或反向；Kechichian等[14]在春分点轨道根数下推导了采用间接法的最优控制律；Casalino等[19]将由经典轨道根数表示的Edelbaum方法扩展到变比冲最优轨迹的求解中。对于太阳能电推进下的变功率变比冲的情况，Mengali等[20]和陈杨[21]在推导最优控制律的判断条件上考虑了推进效率的影响，推进效率是比冲的函数；Chi等[15-16]结合了间接法、同伦法和协态归一化方法，提高了对于TPBVP和MPBVP的求解效率，并且他们改进了对数同伦函数以适用于比冲区间变化的情况。对于这种比冲区间变化的情况，将比冲作为控制量，所以轨迹优化过程中至少有三个控制量存在，Taheri等[17,22,72]提出了适用于多控制量优化的复合平滑控制技术，它是以双曲正切函数近似的平滑思想[66]为基础，应用该策略可以辅助间接法更有效地求解多控制量以及多模式情况下的优化问题。

对于多次引力辅助轨迹优化问题，Chi等[16]用间接法推导了多次引力辅助的一阶必要条件，通过随机协态初值猜测进行轨迹优化。然而对于三次及以上引力辅助问题，随机猜测难以收敛，因此采用单段拼接的方法[73]获得全局可行解。并且他们在单次引力辅助最优轨迹的计算中，利用前文提到的基于伪谱法的GPOPS软件包[52-53]进行了算例对比。Chi等[74]之后采用协态变量转换的方法[65]重新推导了一阶必要条件，因此单段的最优解可以为多段协态初值提供参考，有效地提高了收敛效率。

2.3 比冲随距离变化的轨迹优化方法

直接法和间接法均曾用于求解比冲随距离变化的发动机模型下的轨迹优化问题。其中的一些研究是应用直接法软件MALTO(Mission Analysis Low Thrust Optimization)[75]来计算变比冲最优轨迹[28,34,76]，MALTO是根据Sims等[77]提出的将轨迹按控制点划分成多段，每一段都假设在中途位置施加一个脉冲的离散方法而创建的，之后再通过SNOPT软件包求解此离散变量的非线性规划问题，此外，MALTO非常适用于求解小推力多目标交会以及多次引力辅助问题，控制点可以根据交会天体和甩摆天体来确定。Conway等[27]提出了用解析梯度代替直接配点法中采用的有限差分近似，在求解连续小推力的轨迹优化问题中，无论在收敛率还是求解效率方面，解析梯度法更具有优势。Saripalli等[30]应用Englander等[32,78]开发的EMTG(Evolutionary Mission Trajectory Generator)软件重新计算NEXT发动机多种拟合曲线下的Dawn任务和NEARER任务轨迹，EMTG是基于直接法分为外环和内环优化程序集合，外环采用非支配排序遗传算法离散变量，内环通过单调域跳跃(MBH)和非线性规划算法得到收敛解。Petropoulos等[25-26]提出了用基于形状的标称轨迹算法筛选出合适的轨迹，为直接法的精确求解提供初值。Karthik等[37]在对电推进发动机精确建模的基础上，用MATLAB的工具包DIRETTO对参量进行离散，再用工具包IPOPT求解此非线性规划问题。

Casalino等[35]和Mengali等[33]应用间接法求解比冲随距离变化模型下的燃料最优轨迹。Williams等[27]采用混合法，即直接的遗传算法优化为间接法求解提供初值参考，在此应用软件包SETOP(Solar Electric Propulsion Trajectory Optimization Program)进行求解，SETOP是间接法软件VARITOP[79]的改进版。Taheri等[80]将他提出的复合平滑控制技术[22]用于求解光栅离子发动机模型下的轨迹优化问题。另外，孙冲等[36]在虚拟引力场下通过平均算法推导了轨道设计的解析解，该解析方法可以为直接法或间接法提供初值。

对于复杂的多次引力辅助轨迹优化问题, 目前的很多研究应用形状曲线近似比冲随距离变化的连续小推力多次引力辅助轨迹的方法。Petropoulos等[26,81]应用基于形状的指数正弦曲线逼近最优轨迹，该方法不仅可以对轨迹进行快速的全局搜索，也可以为直接法或间接法轨迹优化提供初值。普渡大学开发的一款轨迹优化软件GALLOP(Gravity-Assist Low-thrust Local Optimization Program)[82]就是根据基于形状的全局搜索和局部直接法优化相结合的思想而编制的。

2.4 功率分档轨迹优化方法

最后对于功率分档发动机模型，很难找到相关的应用直接法的研究。Quarta等[41-42]在NEXT发动机41个工况离散点模型下采用间接法求解时间最优轨迹，41个工况点对应各自的功率、比冲和推力值，因此轨迹优化的控制量为41个工况点以及推力方向，在满足功率约束的前提下，对哈密顿函数按照离散点进行排序并寻找到使哈密顿函数最大的工况点。Chi等[38]应用间接法结合同伦方法以及功率档位检测方法优化行星际转移的时间最优和燃料最优轨迹，这套方法既适用于单个发动机具有多个工况档位的情况，也适合于多个恒定比冲恒定推力发动机带有功率约束的情况。Axelrod等[44]针对离散推力模型提出了基于间接法求解的闭环控制技术，并分别讨论了应用离散控制律和次优控制律的求解效果。

3 结 论

本文对常用的变比冲电推进简化模型和变比冲连续小推力轨迹优化方法进行了简介和研究现状综述。变比冲电推进简化模型有比冲在区间内变化、比冲随距离变化和功率分档变化三种，国内外的很多学者针对此三种模型进行了轨迹优化研究。常见的轨迹优化方法分为直接法和间接法，直接法收敛效率相对较高，但是求得的解无法保证局部最优性；应用间接法可以满足局部最优的一阶必要条件，但是引入协态变量会导致收敛效率较低。由于方法的不足和模型的特殊性，很多研究针对具体问题提出了相应的处理办法。例如，对于比冲在区间内变化的模型，采用同伦法以克服初值猜测的敏感性，但由于考虑比冲作为控制量，所以需要对常用的同伦函数进行改进；对于功率分档发动机模型，提出功率档位检测技术来解决被积项离散不连续的问题，从而可以提高积分精度和收敛效率。此外，对于任务设计中的多次引力辅助问题，采用相应的直接法程序包进行求解，或在应用间接法的过程中增加提高计算效率的策略。

目前没有针对变比冲电推进模型的初始轨道设计方法，所以仍采用适用于连续小推力的脉冲搜索和形状曲线近似等方法。并且，对于考虑更为实际的变比冲发动机模型和更为复杂的多目标探测任务，连续小推力轨迹优化问题的求解仍是一个难题，高效的求解方法和策略仍是需要研究的目标。