一种用于连续推力航天器燃料最优轨迹优化的协状态初始化方法

诸葛运成,任 飞*,李安庆,王 叶

齐鲁工业大学(山东省科学院) a.机械工程学院;b.生物工程学院,山东 济南 250353

航天器的轨道转移技术在各种空间任务中扮演着重要角色。在轨道转移过程中,对航天器的轨迹进行优化设计可以有效减少燃料的消耗,提高空间操作的效率。连续推力作用下航天器轨迹优化问题是近几年的研究热点。目前,航天器的轨迹优化方法主要可分为两类：直接法和间接法。直接法是指通过数值离散将原问题转化为一个参数优化问题,再运用非线性规划[1]之类的优化方法进行求解。通常,直接法的计算量巨大,而且也很难得到最优解。间接法是指通过运用庞特里亚金极大值原理[2]将最优控制问题转化为两点边值问题并进行求解[3-4]。理论上间接法是优于直接法,因为间接法中,控制策略是由极大值原理优化而来,未知参数数量也远远比直接法少。不过间接法求解的两点边值问题收敛域极窄,对协状态变量的初值非常敏感。极大值原理揭示燃料最优问题的推力幅值要么为最大上限,要么为零,即所谓的bang-bang控制。它使得微分方程的右端项不连续,使用传统的数值积分方法很难保证精度,尤其是对于时间较长,发动机开关频繁的情况。Bertrand等[5]提出了一种平滑技术来求解bang-bang控制问题,这种技术称为同伦法。Jiang等[6]和Tang等[7]后续利用同伦法研究了一些连续小推力航天器的轨道转移问题。同伦法的主要思路为：先求解较为简单的平滑的能量最优问题,后通过迭代的方式将问题过渡至不连续的燃料最优问题。

同伦法大大提高了间接法的实用性。然而,在求解能量最优问题时,虽然能量最优问题相对容易求解,但其初值敏感的问题依然存在,如何对协状态变量进行初始化是一个非常棘手的问题。目前常用的办法是将能量最优问题看作是一个优化问题并利用智能优化算法求解[8-11],然后将用智能优化算法求得的解作为同伦法的初始猜测值。然而,智能优化算法一般是进化类算法,其解具有随机性,难以保证同伦法最终的收敛性和最优性。另外,使用智能优化算法时,需要对状态方程和协状态方程反复进行数值积分来计算适应度,计算效率很低。

庞特里亚金极大值原理揭示燃料最优控制问题的最优控制为bang-bang控制,但对于推力在周向、径向和法向的3个分力并非是bang-bang控制,3个分力是随时间连续变化的函数。现实中设在3个方向的发动机推力目前极难做到连续变化控制,最好的方法是主推发动机是bang-bang控制,通过调整航天器的姿态来满足3个方向推力的最优控制。而对于3个方向推力幅值随时间的变化较为剧烈的情况,显然很难通过调整航天器的姿态来实现推力的最优控制。因此本文将3个推力进行各自独立的bang-bang控制,但是这使问题更加复杂,间接法的收敛域更小,初始猜测值的选取非常困难。Ren等[12]尝试将序列二次规划算法嵌入到同伦法的迭代中,一定程度上提高了收敛性,但是计算量剧增。

Casalino[13]针对“近圆近共面”轨道转移问题做了详细的研究。受到Casalino的启发,考虑到整个轨道转移过程中轨道半长轴的变化是最主要的因素,本文提出一种能量最优问题协状态初始化方法。首先,对航天器动力学方程做了合理的简化,与简化动力学方程相应的协状态可以解析地计算出来,从而避免数值积分;然后,将利用简化动力学方程得到的协状态作为能量最优问题的协状态初始猜测值,求解能量最优问题;最后,利用同伦法最终求得燃料最优问题的解。

1 轨道动力学模型及最优控制

1.1 改进春分点根数的轨道动力学模型

绕地飞行轨道动力学模型通常基于经典轨道根数或改进春分点轨道根数建立[14]。当轨道偏心率或倾角为0时,经典轨道根数的动力学方程会发生奇异。因此,采用改进春分点轨道根数建立动力学模型。航天器的动力学方程可表示为：

其中,Fr、Fu和Fh分别为航天器的径向、周向和法向最大推力,g0为重力加速度,Isp为发动机比冲,m为航天器质量。推力大小由uj∈[-1,1] (j=r,u,h)控制。x=[p,f,g,h,k,L]T为航天器的状态变量,p、f、g、h、k和L代表6个改进春分点形式的轨道要素。改进春分点轨道要素和经典轨道要素之间的关系为：p=a(1-e2),f=ecos (ω+θ),g=esin (ω+θ),h=tan (i/2)cosΩ,k=tan (i/2)sinΩ,L=ω+Ω+θ。a、e、Ω、i、ω和θ分别代表轨道半长轴、偏心率、升交点赤经、倾角、近地点幅角和真近点角。

式(1)中,矩阵Mr、Mu、Mh和D为：

1.2 庞特里亚金极大值原理和同伦法

航天器轨迹优化的优化指标通常有两种：能量最优指标和燃料最优指标。能量最优指标是对推力幅值平方或加速度幅值平方的积分,燃料最优指标是对推力幅值绝对值或加速度幅值绝对值的积分。两种指标分别对应了能量最优问题和燃料最优问题。相对于燃料最优问题,能量最优问题是比较容易求解的。这是因为能量最优问题的控制是连续的,而燃料最优问题的控制是不连续的bang-bang控制。利用同伦法建立能量最优问题和燃料最优问题的联系,可构造如下优化指标[5-7]：

其中,t0和tf为初始和终端时刻。当ε=1时,优化指标对应能量最优问题,当ε=0时,优化指标对应燃料最优问题。使ε由1减小至0(或接近0),能量最优问题则过渡至燃料最优问题。

基于庞特里亚金极大值原理,构造哈密顿函数：

其中,λx=[λp,λf,λg,λh,λk,λL]T和λm为状态变量和质量所对应的协状态变量。将Hr、Hu和Hh分别命名为径向哈密顿函数、周向哈密顿函数和法向哈密顿函数。以Hu为例：

最优控制应该使哈密顿函数取最小值,那么当-1≤uu<0时：

当0≤uu≤1时：

所以：

相同的方法,可以得出径向和法向的最优控制。

哈密顿函数(式(5))对各状态变量求偏导,可以得到协状态方程：

终端物理约束描述为：

Ψ[x(tf),m(tf),tf]=0。

(11)

通常,航天器的终端质量并不做约束,那么λm(tf)为0。

对于轨道转移问题,航天器的初始状态x(t0)=[p(t0),f(t0),g(t0),h(t0),k(t0),L(t0)]T以及初始质量m0=m(t0)是已知的,目标状态(终端状态)xf=[pf,ff,gf,hf,kf,Lf]T是已指定好的。如果给定一组协状态初始值λ(t0)=[λp(t0),λf(t0),λg(t0),λh(t0),λk(t0),λL(t0),λm(t0)]T,对式(1)和式(10)进行积分,那么可以得到终端状态和终端协状态,同时通过式(9)得到最优控制。终端状态和终端协状态需要满足终端物理约束Ψ以及λm(tf)=0,那么此时轨迹优化问题转变成了一个两点边值问题,其靶方程表示为：

其中,z表示待求解的初始协状态z=λ(t0)。求解靶方程可以使用Levenberg-Marquardt算法(LMA)。LMA能提供非线性最小化(局部最小)的数值解,可用于求解非线性方程组和非线性优化问题。LMA在高斯-牛顿算法的基础上引入阻尼参数,借由执行时修改阻尼参数达到结合高斯-牛顿算法和梯度下降法的优点,并改善两者之不足。具体算法细节可参阅文献[15]。

从式(7)和式(8)可以看出,只要ε的值不为0,最优控制力的幅值虽然对时间不可导,但总是连续的。经验表明,ε=1时(能量最优问题)的靶方程相比ε=0(燃料最优问题)的情况要易求得多。因而可求得ε=1时的解,再逐渐减少ε的值,并以前一步得到的解作为初始猜测值进行迭代求解直至得到ε的值很接近0时的解,即燃料最优问题的解。ε的减小可以按照这样的原则：先选用一个较大的步长减小,随着ε减小,减小其减小的步长。采用的ε的减小规则为：

先利用LMA求得ε0=1所对应的靶方程的解,即能量最优解;后利用LMA求解ε1至ε60所对应的60个靶方程的解,求解时,以εn-1对应的解作为初始猜测值去求解εn对应的解;最终得到ε60=10-4所对应的靶方程的解,即燃料最优解。LMA求解每个靶方程时,都需要初始猜测值,或者说需要对待求解量进行初始化。因此,在整个计算流程的开始,还需要有能量最优问题的λ(t0)的初始猜测值。实际上,此初始猜测值的选取是非常困难的,且初始猜测值的质量好坏直接影响最终能否得到收敛的燃料最优解。

另外,当ε接近0时,控制律趋于bang-bang控制,微分方程的右端虽然理论上仍然连续,但是在控制切换点处变化非常剧烈,常用的自适应步长数值积分方法通常难以应对这种情况。因此,本文中的数值积分方法采用定步长四阶龙格库塔法,积分步长为0.1 TU。

2 能量最优问题的协状态初始化方法

轨道转移过程中,如果初始轨道和目标轨道的偏心率较小且轨道倾角差较小,那么轨道半长轴的变化是影响燃料消耗的主要因素。将f、g、h、k简化为0,p简化为一个常数P=0.5[p(t0)+pf]。这样可以得到简化的动力学方程：

将式(18)代入式(14)得：

将式(19)整理成如下形式：

式(20)中,矩阵N为：

其中,

对于轨道转移问题,航天器的初始状态x(t0)和目标状态xf是已知的,那么由t0时刻(一般t0为0)至tf时刻,航天器的状态改变量为：

3 数值仿真

考虑一个连续小推力长时间多圈飞行的异面轨道转移问题,此轨道转移问题在我们以往的研究中已有所涉及[16]。本文将以此轨道转移问题为例,尝试将航天器3个方向的推力控制分开进行独立优化。传统的主推力开关及其方向最优控制方法可参阅文献[16]。航天器的初始质量为1 000 kg,发动机比冲为400 s。径向,周向和法向最大推力分别为2、4和4 N。轨道转移时间为570 TU(约5 d 8 h),即要求航天器在570 TU时到达目标位置。航天器的初始状态,目标状态以及计算结果(终端状态)如表1所示。能量最优解的燃料消耗为264.12 kg,燃料最优解的燃料消耗为199.50 kg,绕地球飞行33圈。通过比较终端状态和目标状态可以发现解的精度很高,利用简化动力学模型生成的初始猜测值可以得到收敛且可行的能量最优解和燃料最优解。

表1 航天器初始状态与终端状态

能量最优飞行轨迹及其最优控制如图1和图2所示,燃料最优飞行轨迹及其最优控制如图3和图4所示。通过使用同伦法,连续变化的能量最优控制(图2)成功转变成为燃料最优bang-bang控制(图4)。图4中可以看出,径向推力Fr始终为0。实际上,轨道转移过程中径向推力不是必要的,仅使用周向和法向推力即可完成轨道转移,周向推力改变轨道形状,法向推力改变轨道面。如果径向推力出现,必然会造成额外的燃料消耗,与我们节省燃料的目的相悖。本算例中故意设置了2 N的径向推力,而得出的最优控制中并未使用这2 N的径向推力,验证了最终燃料最优解的合理性以及算法的有效性。图3可看出,周向推力(淡蓝色正方形)主要出现在近地点附近区域。本算例中初始轨道半长轴比目标轨道半长轴小,即航天器由低轨向高轨转移。近地点处航天器的速度快,于近地点处施加推力有利于航天器加速向高轨转移。法向推力(红色圆形)主要出现在初始和目标轨道面的交线附近。按照轨道动力学理论,轨道面的改变实际上是轨道面以航天器法向推力施加点和地心的连线为轴发生转动。在航天器异面轨道转移问题中,使初始轨道面绕着初始和目标轨道面的交线转动至目标轨道面无疑是最高效的,因此法向推力出现在初始和目标轨道面的交线附近。综上,本算例中航天器的周向和法向推力控制是合理的,进一步验证了算法的有效性。

图1 航天器能量最优飞行轨迹

图2 能量最优控制

图3 航天器燃料最优飞行轨迹

图4 燃料最优控制

相较于将能量最优问题看作优化问题并利用智能优化算法搜索协状态初值,提出的基于简化动力学模型的初始化方法可以解析地生成初始猜测值,避免了反复数值积分计算适应度的过程,因此大大提高了计算效率。另外,本算例考虑了实际工程中的需求,分开对3个方向的推力控制进行优化,这大大增加了问题的复杂程度,间接法初值敏感问题变得更加明显。利用提出的协状态初始化方法可以顺利求解这种复杂问题。针对本算例,也尝试了随机初始化和利用智能优化算法进行初始化,但最终未能得到收敛的燃料最优解。

4 结 论

针对连续小推力多圈飞行异面轨道转移这类复杂轨迹优化问题进行了研究。考虑到实际工程中的需求,对3个方向的推力控制分开进行优化。采用间接优化方法,并利用同伦法将能量最优解过渡至燃料最优解。提出了一种基于简化动力学模型的初始化方法,可以快速生成高质量的能量最优协状态初始猜测值,有效解决了间接法的初值敏感问题。仿真结果验证了这种初始化方法的有效性和高效性。