初值约束与两点边值约束轨道动力学方程的快速数值计算方法1)

张 哲 代洪华 冯浩阳 汪雪川 岳晓奎

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

(航天飞行动力学技术国家级重点实验室,西安 710072)

引言

轨道快速计算是航天工程领域的基础问题,广泛存在于地球绕行轨道设计及深空探测等各种航天工程任务中,具有重要意义.轨道动力学问题在形式上一般描述为常微分方程初值问题或边值问题,其中最具代表性的有轨道递推问题[1-2]、摄动Lambert问题[3-4]和三体模型下的转移轨道设计问题[5-6].轨道问题在早期一般通过微积分进行解析求解,牛顿与伯努利计算了万有引力作用下的二体轨道解析解[7],欧拉与拉格朗日求出了限制性三体问题的5 个拉格朗日点[8].但当考虑更为一般的三体系统以至N 体系统时,解析求解将不再适用[9].庞加莱发展了渐进法并将其用于求解限制性三体问题,但由于三体问题具有混沌效应,其长期轨道仍然难以通过渐进法精确预测.此外,早期对轨道动力学问题的求解都是基于理想的均匀圆球或椭球模型,而在实际航天任务中则需要考虑地球非球形等扰动项,这些干扰力使得航天器轨道难以求得精确解析解.由于存在上述困难,在实际工程任务中,航天器轨道均采用数值方法进行求解.

轨道动力学问题最为经典的数值解法是有限差分类方法,其中代表性的包括Euler 法、Runge-Kutta法等[10].此类算法将待求解问题划分为多个差分节点,逐点计算问题近似解.由于物理含义直观且递推方程较为简单,有限差分法受到了广泛关注,并在仿真软件与科学计算函数中得到大量应用[11-12].然而,有限差分法的计算精度严重依赖于小步长,方法性能受到限制,在诸如航天器追逃博弈以及失效航天器快速在轨维修等任务中,将难以满足任务实时性以及长期解算过程高效高精度的要求[13-14].

为克服有限差分法的原理性缺陷,学者们发展了一系列能实现大步长计算的迭代算法[15-21].1963年Clenshaw 与Norton[15]首次将Chebyshev 多项式与Picard 迭代法相结合,用于求解常微分方程的半解析解,该算法能够避免有限差分法逐步计算的缺点,但其系数计算方案较为复杂,且仅适用于求解一维方程.其后,日本Fukushima[16]对Chebyshev 多项式与Picard 迭代法的结合方式予以改进,简化了系数计算方案,并将求解对象拓展到高维常微分方程,但在这一改进方案中待解函数及其导数被分别估计,未能有效利用Chebyshev 多项式性质,引入的待定系数较多,计算量较大.此后,美国德州农工大学Bai[17]提出修正Chebyshev-Picard 迭代法(modified Chebyshev-Picard iteration method,MCPI),并将其用于求解轨道动力学问题,该算法利用Chebyshev多项式性质减少了待定系数,但由于在迭代公式推导中利用了反演公式,得到的迭代公式较为复杂,降低了算法效率.为进一步提高MCPI 算法计算效率,Woollands 等[18]对该算法进行改进,将导函数估计值与真实值之差作为迭代修正项,提高了算法收敛阶,但引入了雅可比矩阵,增加了算法计算量.Wang 等[19]将变分迭代法与时域配点法相结合,提出了局部变分迭代法(local variational iteration method,LVIM),该算法将计算残差用于误差修正,结合时域配点法简化了迭代公式推导过程,得到的计算公式较为简洁,在轨道计算问题中收敛速度较快,但该算法由于对广义拉格朗日乘子进行Taylor展开,同样引入了雅可比矩阵,在求解复杂高维问题时产生的计算量过大.2020 年,Wang 等[20]将牛顿法与Picard 迭代法相结合,提出了多步Newton-Picard法(multistep Newton-Picard method,MNP),该算法将导数方程在待定函数真解处进行Taylor 展开,并依据展开式建立了一系列衍生算法,但Taylor 展开引入的偏微分算子使得算法迭代公式较为复杂,且算法收敛阶较低.冯浩阳等[21]将局部变分迭代法、拟线性化法及叠加法相结合,提出了能用于求解边值问题的拟线性化-局部变分迭代法,拓展了局部变分迭代法的适用范围,但仍未能解决算法需要反复计算雅可比矩阵的弊端.

本文将误差反馈机制与Picard 迭代法相结合,提出一种新的反馈迭代算法,避免对原函数进行泰勒展开,消除迭代公式中的雅可比矩阵,以保证计算效率;利用配点法与Chebyshev 多项式性质,将迭代公式推导中涉及的复杂符号运算转化为线性代数运算,简化推导过程,降低迭代公式的复杂度;结合拟线性化法与叠加法,拓展算法适用范围,使之能够求解两点边值问题;将变参数计算方案引入局部配点反馈迭代法,使算法能够在运行过程中自适应选择计算参数,进一步提高算法计算效率与计算精度,从而为在轨航天器的长期轨道快速预测提供一种高性能数值计算方法.

1 局部配点反馈迭代法

1.1 修正Chebyshev-Picard 迭代法

修正Chebyshev-Picard 迭代法是用于迭代求解常微分方程初值问题的数值算法[17],该算法求解常微分方程的具体过程简述如下.

对于具有如下形式的一阶常微分方程

假设未知函数及其导数可以用正交多项式的线性组合来近似,即

式中,Ti(t)为一组正交多项式组的第i 项,βi及Fi为Ti(t)的系数.

Picard 迭代法计算公式为

将式(2)与式(3)代入式(4),并令式(4)两边正交多项式对应项系数在一系列给定时刻相等,进一步整理可得未知函数y(t)的迭代计算公式为

1.2 局部配点反馈迭代法

在式(4)中,迭代公式通过对导数向量积分得到函数值,利用新的函数值更新导数向量并再次进行积分迭代,直到达到目标精度.本节将误差反馈机制引入迭代过程,结合时域配点法与局部离散化思想,建立一种新的局部配点反馈迭代法.

Picard 迭代公式(4)与下式等价

对式(6)进行整理,以导数估计值 y˙n(t) 与利用状态估计值计算得到的导数值 g(yn(t),t) 之差作为修正项,将式(6)改写为如下迭代公式

称式(7)对应的迭代方案为反馈迭代法.

对于[t0,tf]范围内的任意时刻tm,可以使用下式将其转化到区间[-1,1]范围内

通过式(8) 将求解区间进行转化后,使用Chebyshev 正交多项式的线性组合作为未知函数y(t)的估计解,即近似认为

式中,cos(karccost) 为第k 项Chebyshev 多项式在t 时刻的值,lk为第k 项Chebyshev 多项式的系数.

为求解式(9)中N 个未知多项式的系数lk(k=1,2,···,N),对式(9)采用时域配点法[22],即令正交多项式线性组合所构成的函数估计解在N 个给定时刻t1,t2,···,tN与函数真值相等,从而构造由N 个线性方程构成的线性代数方程组

将式(10)记为如下矩阵形式

式中,y 为N 个时刻 t1,t2,···,tN对应的函数真值构成的真值向量,C 为Chebyshev 多项式矩阵,其中第i 行j 列元素Cij=cos(jarccosti),L 为系数向量L=[l1,l2,···,lN]T.

在给定时刻,未知函数的真值与Chebyshev 多项式矩阵C 可以预先求得,因此可以利用下式求解系数向量L

将式(12)得到的系数向量代入式(11),再将式(11)代入式(7),得到配点反馈迭代法的迭代公式

由于Chebyshev 多项式形式已知,对其进行微积分运算可得

由式(14)及式(15)可以得到Chebyshev 多项式矩阵的微分形式及积分形式.据此可将式(5)推导中涉及的符号运算全部转化为代数运算,最终得到配点反馈迭代法的迭代计算式

式中,P为积分形式的Chebyshev 多项式矩阵,计算公式为P=

在计算过程中对式(16)反复迭代,直至计算结果精度达到目标要求.此时,得到的系数向量L 中每一项即为使用式(9)作为未知函数估计解时,满足精度要求的Chebyshev 多项式对应项系数.

Cuthrell 和Biegler[23-24]的相关研究结果指出,将整个函数待解区域作为配点区间,直接得到全局估计解的方法仅适用于求解光滑问题,而对于非光滑或难以直接全局近似估计的问题,应当对全局区域进行离散,并在离散后得到的每个较小子区间内使用正交多项式进行分段估计.考虑到将区间离散化后在每个较小子区间内对未知函数进行估计能够取得更好的估计效果,在应用配点反馈迭代法时,首先将整个待解区间离散化为多个子区间,从第一个子区间开始对未知函数进行迭代计算,之后将第一个子区间终点函数值作为下一子区间初值,再次使用配点反馈迭代法估计该区间内的未知函数数值解,重复上述过程,直至求得未知函数在全部待解区间内的数值解.这一结合局部离散思想的配点反馈迭代法称为局部配点反馈迭代法(local collocation feedback iteration method,LCFI),算法计算过程示意图见图1.

图1 局部配点反馈迭代法示意图Fig.1 Illustration of LCFI

本节提出的局部配点反馈迭代法与Picard 迭代法在数学上虽然等价,但式(7)通过结合误差反馈,将导数残差作为修正项引入迭代公式,实际计算效率具有明显优势.此外,在Bai[17]提出的修正Chebyshev-Picard 迭代算法中,由于其利用多项式正交性计算系数,推导过程十分繁琐,且迭代公式较为复杂,增加了编程工作量,降低了算法的计算效率.本文通过结合配点法,直接计算正交多项式系数向量,在实际计算中效率更高.汪雪川[25]通过相关数值实验证明了局部变分迭代法的计算效率高于修正Chebyshev-Picard 迭代算法,在本文下一节的数值仿真中可以看到,本文提出的局部配点反馈迭代法计算效率远高于局部变分迭代算法.上述事实证明相比于修正Chebyshev-Picard 迭代算法及变分迭代法,本文算法在计算效率上具有显著的优越性.同时,通过比较本文算法迭代公式与文献[25]中局部变分迭代算法迭代公式不难发现,本文算法消除了雅可比矩阵,这一特点使得算法能够更容易的被应用于求解复杂高维问题.

2 拟线性化-局部配点反馈迭代法

局部配点反馈迭代法能够计算常微分方程初值问题,但不能直接求解Lambert 问题等边值问题,本节介绍一种将拟线性化法及叠加法与局部配点反馈迭代法相结合的方案,使算法能够处理两点边值约束下的轨道计算问题.

2.1 拟线性化法

拟线性化法是一种将非线性微分方程近似线性化的数学方法[21].

对具有如下形式的二阶非线性微分方程边值问题

式(17)可以改写为

令y0为未知函数y 的初始估计解,yn和yn+1分别表示第n 次和第n+1 次迭代结果,将第n+1 次迭代结果在处展开,略去二阶及以上偏导数,可以得到将非线性二阶微分方程(17)拟线性化后的迭代计算公式

式(19)对应的边界条件为 yn+1(a)=ya,yn+1(b)=yb.

利用式(19)迭代求解未知函数y,求解过程中,对于第n 次迭代得到的yn,在代入式(19)进行第n+1次迭代时,若将yn视为已知量,则式(19)可视为仅关于yn+1的二阶微分方程,求解该方程即可得到yn+1.重复上述迭代过程,当yn+1=yn时,迭代结束,此时yn+1即为微分方程(17)的解.

2.2 叠加法

叠加法能够将微分方程边值问题转化为初值问题[26].对于如下形式的微分方程边值问题

假设解的形式为

将式(21)代入式(20),可将边值问题(20)转化为求解使如下两个二阶常微分方程成立的未知函数y1与y2

满足边界条件的y1与y2可由如下两组微分方程初值问题分别求得

满足边界条件的μ 由下式计算得到

在计算时,首先求解式(24)和式(25)两个微分方程初值问题得到y1(t)与y2(t),之后将b 点的两个函数值y1(b)与y2(b)代入式(26)计算μ,最后,将μ及y1(x)与y2(x)代入式(21),得到未知函数y(x)的解.

对于Lambert 问题等两点边值条件约束下的轨道动力学问题,首先利用拟线性化法将其转化为二阶线性微分方程边值问题,再利用叠加法将二阶线性微分方程边值问题转化为一组初值问题,使用局部配点反馈迭代法对初值问题分别求解,再依据式(21)即可得到两点边值问题的解.

3 数值仿真结果

3.1 轨道递推

轨道递推问题的动力学方程为

式中,r=[x,y,z]T,表示航天器的位置矢量,μ=3 986 004.418 × 108m3/s2,表示地球引力常数,r=‖r‖,表示地心与航天器质心间的距离,t0为轨道递推初始时刻.aJ为摄动项,本算例中考虑J2项摄动.

轨道递推为初值问题,在给定初值条件的情况下,可以使用局部配点反馈迭代法直接求解.本文通过连续100 次仿真实验比较算法性能.为了与同类方法的计算性能进行对比,本算例使用文献[21]中轨道递推问题初始位置,相关计算参数见表1,100 次仿真平均单次计算时间比较结果见表2.

表1 递推轨道计算参数Table 1 Calculation parameters of LCFI in orbit propagation

表2 轨道递推计算效率对比Table 2 Comparation of efficiency on orbit propagation

由表2 可见,在相同配点条件及计算精度下,针对轨道递推问题,局部配点反馈迭代算法(LCFI)的计算速度为局部变分迭代法(LVIM)的1.5 倍以上(1.78 倍),在相同精度要求下,LCFI 算法的计算效率更高.进一步的仿真结果表明,当配点个数上升到32 个时,在表1 计算精度条件下,LCFI 算法计算步长可以达到12 000 s,而基于有限差分原理的ode113 算法与ode45 算法计算步长均小于1200 s(最大步长分别为889.9 s 与1 147.3 s),本文算法在上述问题中能将计算步长提高到有限差分类算法的10 倍以上.

LCFI 与LVIM 计算结果及误差见图2～ 图4,对于7 日内的轨道递推结果,两种算法所得目标轨道最大绝对误差小于0.4 m,最大相对误差小于4.9 ×10-8,这样的计算精度在工程实际中是可以接受的.同时,从图3 绝对误差变化曲线与图4 相对误差变化曲线可以看出,算法在迭代过程中能够对误差进行周期性自动修正,从而降低误差发散速度.

图2 LCFI 与LVIM 所得递推轨道对比Fig.2 Comparison between orbits propagated via LCFI and LVIM

图3 LCFI 与LVIM 所得递推轨道绝对误差Fig.3 Absolute error of the orbits propagated via LCFI and LVIM

图4 LCFI 与LVIM 所得递推轨道相对误差Fig.4 Relative error of the orbits propagated via LCFI and LVIM

3.2 摄动Lambert 轨道

Lambert 轨道转移问题是经典的轨道力学问题,也是航天器轨道动力学快速计算方法研究中常用的求解对象[21,25].本节对摄动Lambert 转移轨道及转移初速度进行求解,并与拟线性化-局部变分迭代算法求解结果及计算效率进行对比,验证LCFI 算法的精确性及快速性.

摄动Lambert 问题的动力学方程与式(27)所示二体轨道动力学系统的递推问题相同,将式(27)右侧记为g(r),拟线性化处理后的Lambert 问题动力学方程及其边界条件表达式为

依据叠加法,将 rn+1写为如下形式

其中 r1与 r2利用如下两个微分方程初值问题求得

式(29)中μ 由下式求得

为与同类方法的计算性能进行对比,本算例同样使用文献[21]中Lambert 转移问题对应的始末位置条件,坐标系采用地球固连坐标系,经过连续100 次仿真,单次平均用时见表3,初始条件及相关计算参数见表4,计算所得轨道结果及不同算法所得轨道计算误差见图5～ 图7.

图5 LCFI 与LVIM 所得Lambert 转移轨道结果对比Fig.5 Comparison between Lambert orbits obtained via LCFI and LVIM

图6 LCFI 与LVIM 所得Lambert 转移轨道绝对误差Fig.6 Absolute error of Lambert orbits obtained via LCFI and LVIM

图7 LCFI 与LVIM 所得Lambert 转移轨道相对误差Fig.7 Relative error of Lambert orbits obtained via LCFI and LVIM

表3 摄动Lambert 问题计算效率对比Table 3 Comparation of efficiency on Lambert problem

表4 摄动Lambert 转移轨道计算参数Table 4 Calculation parameters of LCFI in perturbed Lambert problem

表3 表明,在相同计算参数及计算精度条件下,针对摄动Lambert 轨道转移问题,局部配点反馈迭代法(LCFI)计算速度为文献[21]中拟线性化-局部变分迭代算法(QL-LVIM) 的2 倍以上(2.68 倍).

在计算结果精度方面,LCFI 算法计算得到的轨道转移初速度为[-1 687.309 900 000,-0.024 275 234,2 922.608 000 000] m/s,LVIM 算法计算得到的轨道转移初速度为[-1 687.309 900 000,-0.024 305 875,2 922.607 900 000] m/s,两种算法求得的转移初速度绝对误差二范数小于1.05 × 10-4m/s,相对误差二范数小于3.10 × 10-8,速度计算误差远小于实际测量与控制精度,这样的计算误差在实际应用中是可以忽略的.在计算误差变化趋势方面,由图6 及图7 可见,针对本算例中的摄动Lambert 轨道,两种算法对位置计算结果的最大绝对误差不超过0.83 m,最大相对误差不超过2.1×10-8,且整个计算过程中误差经过初期增长后受到抑制,在后续计算过程中呈下降趋势.

3.3 圆型限制性三体模型约束下转移轨道

当两个主要天体围绕其系统质心做圆周运动,同时存在一个质量可忽略的第三天体在两个主要天体运动的公共平面内运动,这样的问题就构成了平面圆形限制性三体问题[27].月球探测航天器在地-月转移过程中,运动范围始终位于地球影响球内,而地-月系本身围绕系统质心旋转,因此在初步设计地-月间循环轨道时,可以将航天器的运动简化为地-月-航天器圆形限制性三体问题模型(circular restricted three body problem,CR3BP),该模型也是描述星际探测航天器轨道转移问题的经典非线性模型[28-33],本节通过解算该模型验证LCFI 算法的有效性.

对于地球、月球、航天器构成的三体系统,令地-月系质心为坐标系原点,x 轴正方向由地球质心指向月球质心,地-月系轨道平面为(x,y)平面,z 轴与地月系统轨道面垂直,建立空间右手直角坐标系.设地球质量为m1,月球质量为m2,地球到坐标系原点距离为l1,月球到坐标系原点距离为l2,将地月间距离及地月系质量均无量纲化为1,即令月球的无量纲质量为μm,地球的无量纲质量为1-μm,地球到原点的距离为μ,月球到原点的距离为1-μ,航天器的无量纲形式运动方程为[33]

式(33)表征的圆形限制性三体问题动力学模型中,在其L1,L2以及L3拉格朗日点附近存在对应于系统某一周期解的Halo 轨道.由于Halo 轨道附近存在着类似管道的不变流形,航天器在流向Halo 轨道的稳定流形中运动时几乎不需要消耗能量,能够显著节约燃料[34],因此Halo 轨道及其不变流形常常被用于设计低成本登月轨道[9,35-38].

在由近地轨道或绕月轨道向不变流形转移的过程中,需要对航天器运动状态进行多次调整,这一过程中的转移轨道需要精确计算[21].本节算例通过计算探月航天器从185 km 高度地球停泊轨道向流向L1点的稳定流形机动时的转移轨道,验证算法在圆形限制性三体问题动力学模型中轨道计算的的高效性.为了与同类方法的计算性能进行对比,本算例使用文献[21]中该算例计算参数,具体计算参数见表5,其中转移时间以地月系统运动周期作为系统的无量纲时间2π,相关转移时间的计算以此为基准,连续100 次计算时间比较结果见表6.计算所得轨道结果及不同算法所得轨道计算误差见图8～ 图10.

表5 圆型限制性三体模型约束下转移轨道计算参数Table 5 Calculation parameters of LCFI in transfer orbit calculation problem with the constraint of circular restricted three-body model

表6 圆形限制性三体模型约束下转移轨道计算效率对比Table 6 Comparation of efficiency on transfer orbit calculation problem with the constraint of three-body model

图8 LCFI 与LVIM 轨道计算结果对比Fig.8 Comparison between orbits obtained via LCFI and LVIM

图9 LCFI 与LVIM 对平面圆形限制性三体模型约束下转移轨道计算绝对误差Fig.9 Absolute error of transfer orbits under the constraint of circular restricted three body model obtained via LCFI and LVIM

图10 LCFI 与LVIM 对平面圆形限制性三体模型约束下转移轨道计算相对误差Fig.10 Relative error of transfer orbits under the constraint of circular restricted three body model obtained via LCFI and LVIM

由表6 可见,在相同计算参数条件下,针对平面圆形限制性三体模型约束下的轨道转移问题,LCFI 的计算速度为LVIM 的6 倍以上(6.92 倍),即在相同精度要求下,LCFI 算法的计算效率更高.从图9 及图10 误差变化趋势来看,误差在计算初期增长后迅速受到抑制,最大绝对误差小于279.5 m,最大相对误差小于5.2 × 10-6.

4 变参数计算方案及仿真实例

本文提出的局部配点反馈迭代法相比于传统数值积分方法具有大步长计算特性,为发挥算法大步长收敛优势,在满足工程计算精度要求的基础上尽可能减少计算量,本文借鉴ph 网格细化法(ph mesh refinement method)[39-40]计算思想及针对该方法的相关改进策略,将变参数计算方案引入本文算法,从而使得局部配点反馈迭代法能够自适应选择具有更高计算效率的参数组合.

4.1 变参数计算方案

参数变化算法流程图见图11,具体的自适应参数调节过程分为两部分.

图11 变参数局部变分迭代算法流程图Fig.11 Flow chart of adaptive local collocation feedback iteration method

(1)在一个子区间内,迭代计算精度大于迭代终止误差限,但迭代次数超过预定的迭代计算终止次数Nmax时,将当前区间的计算步长减小为原来的a 倍(a＜1),并对当前计算区间内各配点时刻对应的函数值重新进行计算.

(2)在一个子区间内,迭代计算精度达到迭代终止误差限,且迭代次数小于预设定的迭代计算终止次数Nmax时,首先对子区间内M 个节点状态进行插值,利用插值函数计算子区间配置M+1 个节点时各节点状态估计值.其次对子区间内M 个节点的导数估计结果进行插值,并利用插值函数计算子区间配置M+1 个节点时各节点的导数估计值,对导数估计值积分后得到子区间内M+1 个节点函数状态估计值.根据相关研究结果,两次估计结果差的模值与当前对M 个节点的估计结果的误差基本相等[39],因此可以使用该方法对当前计算结果的误差进行估计,上述相对误差计算过程的具体计算公式如下

当计算结果的相对误差向量er中值最大的元素高于规定的误差上限e1时,此时计算结果误差过大,需要增加计算精度以满足计算精度要求.当计算结果的相对误差向量er中最大元素低于规定的误差下限e2时,需要适当降低计算精度以提高计算速度.依据ph 网格细化法,需用节点个数计算公式为

式中,er＞e1时P=lg(er/e1),er＜e2时P=lg(er/e2),M2为计算所得的需用节点个数,M1为当前的节点个数.

当计算结果的相对误差er大于规定的误差上限e1且需用节点个数大于节点数上限Mmax时,令当前计算区间的步长降为当前步长的b 倍(b＜1),节点数等于最小节点个数,并重新计算当前区间内各个节点函数值.当计算结果的相对误差er大于规定的误差上限e1且需用节点个数小于节点数下限Mmax时,令节点数等于需用节点个数,并重新计算当前区间结果.当计算结果的相对误差er低于规定的误差下限e2且需用节点个数小于节点数下限Mmin时,令当前计算区间的步长增加到当前步长的c 倍(c＞1),节点数等于最小节点个数,并重新计算当前区间结果.当计算结果的相对误差er低于规定的误差下限e2且需用节点个数大于节点数下限Mmin时,令节点数等于需用节点个数,并令当前计算区间的步长增加到当前步长的d 倍(d＞1),重新计算当前区间结果.当计算结果的相对误差er高于规定的误差下限e2且低于规定的误差上限e1时,保留当前区间计算结果与计算参数设置情况,并进行下一区间的计算.

4.2 仿真实例

近年来,微小卫星在航天工程领域受到了广泛关注.缺乏推进系统的微小卫星在轨运行状况发射前需要依赖轨道动力学模型经过数值计算得到.目前,微小卫星大多运行于近地轨道,在该轨道上航天器受到的最主要干扰力为地球非球形摄动,其大小比其他摄动力高出3 个量级以上[41],本文基于40 阶EGM2008 球谐重力场模型,在相同初始计算参数及迭代计算精度下,计算两组不同初始条件下卫星轨道经过12 960 400 s (15 天)后的目标位置.初始计算参数见表7 及表8,变参数LCFI 相关计算参数见表9,计算用时见表10,不同算法计算结果见图12 及图13,计算过程中参数变化情况见图14～ 图17.

表7 圆轨道递推初始计算参数Table 7 Calculation parameters of orbit propagation

表8 椭圆轨道递推初始计算参数Table 8 Calculation parameters of orbit propagation

表9 变参数LCFI 递推轨道计算参数Table 9 Calculation parameters of adaptive LCFI in orbit propagation

表10 不同算法计算时间Table 10 Time cost of different numerical method

图12 变参数LCFI 及定参数LCFI 圆轨道递推计算结果Fig.12 Circular orbit obtained via Adaptive LCFI and LCFI

图13 变参数LCFI 及定参数LCFI 椭圆轨道递推计算结果Fig.13 Elliptical orbit obtained via adaptive LCFI and LCFI

图14 变参数LCFI 圆轨道递推计算步长变化Fig.14 Step size of ALCFI in the propagation of a circular orbit

图15 变参数LCFI 圆轨道递推计算配点数变化Fig.15 Variation of the collocation point number in the propagation of a circular orbit

图16 椭圆轨道递推计算步长变化Fig.16 Step size of ALCFI in the propagation of an elliptical orbit

图17 椭圆轨道递推计算配点数变化Fig.17 Variation of the collocation point number in the propagation of an elliptical orbit

根据图14～图17 可见,针对考虑40 阶EGM2008地球球谐重力场模型的轨道递推问题,在相同初始计算参数及迭代计算精度条件下,变参数LCFI 算法能够自适应选择更优的计算步长及基函数阶次,从而降低计算量.

根据表10 可见,变参数LCFI 计算速度是定参数LCFI 的6 倍以上(6.57 倍与7.04 倍),定参数LCFI 由于初始计算参数设置的合理性不足,限制了算法自身大步长高速计算特性的发挥,该问题通过本节变参数计算方案得到了有效解决.

从图12 及图13 可以看出,由于变参数计算方案中引入了误差评估机制,在航天器高精度轨道预测问题中,相同迭代计算精度下,变参数局部配点反馈迭代法的(ALCFI)计算结果基本落在精确参考轨道附近,而采用定参数局部配点反馈迭代法的(LCFI)的轨道计算结果则产生了更大程度的偏离(图12 及图15 中精确参考轨道均使用Matlab 软件中ode45 程序计算得到,ode45 设置精确参考轨道计算相对误差2.3 × 10-14,绝对误差1 × 10-14m),表明变参数计算方案引入的误差评估与调节机制能够进一步提高局部配点反馈迭代法计算结果精度.

5 结论

针对在轨航天器轨道动力学系统快速解算需求,本文基于积分补偿迭代思想,提出了一种适用于轨道动力学快速计算的高性能数值积分方法,有效克服了传统数值差分类方法的缺陷.主要工作与结论总结如下.

(1) 结合Picard 迭代法、误差反馈思想和时域配点法,提出了一种能够高效求解微分方程初值问题的局部配点反馈迭代法.

(2) 将局部配点反馈迭代法和拟线性化法及叠加法相结合,求解了轨道转移两点边值问题.

(3) 通过解算递推轨道、摄动Lambert 轨道转移问题以及圆型限制性三体模型约束下的转移轨道问题验证了本文算法的有效性.计算结果显示在相同迭代精度及计算参数设置条件下,本文算法在上述轨道问题中计算效率比局部变分迭代算法提高50%以上.

(4) 将变参数计算方案引入局部配点反馈迭代法,使算法能够自适应调节计算参数.相关计算结果表明,本文引入的变参数计算方案能够在保证计算精度条件下发挥算法大步长快速收敛优势,并进一步提高算法计算精度.

局部配点反馈迭代法对轨道动力学初值问题及两点边值问题的计算效率相较于现有方法具有明显优势,在计算能力较弱的星载计算机中[42],该优势对计算时间的节约程度将更为显著,这对执行失效卫星抓捕等空间快速机动任务具有重要意义.在后续工作中,将探究局部配点反馈迭代法的并行算法设计等改进方案,以进一步提高算法性能.