精细指数积分法在卫星编队飞行动力学中的应用

邓子辰　李庆军

西北工业大学工程力学系, 西安 710072; † E-mail: dweifan@nwpu.edu.cn



精细指数积分法在卫星编队飞行动力学中的应用

邓子辰†李庆军

西北工业大学工程力学系, 西安 710072; † E-mail: dweifan@nwpu.edu.cn

编队飞行卫星间的距离远小于卫星的轨道半径, 其动力学方程表现为弱非线性。针对弱非线性方程的求解, 提出精细指数积分方法, 用精细积分法求解指数积分方法中的指数矩阵。用精细指数积分法和 Runge-Kutta方法, 在不同条件下求解弱非线性方程的算例, 验证了精细指数积分法的有效性。通过Lagrange方程,建立卫星编队飞行动力学方程的半线性形式, 用精细指数积分方法与 Runge-Kutta 方法求解方程。数值计算结果表明, 与同阶的 Runge-Kutta 求解弱非线性微分方程相比, 精细指数积分法具有更高的精度, 为卫星编队飞行动力学仿真提供了一种有效的数值算法。

指数积分方法; 精细积分法; 卫星编队飞行; Runge-Kutta方法

北京大学学报（自然科学版）第52卷第4期2016年7月

Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 52, No. 4 (July 2016)

卫星编队飞行指若干个小卫星按特定的编队飞行, 各卫星协同工作, 形成一个分布式的编队卫星群。与大卫星相比, 编队卫星群成本低、可靠性强、灵活性高、编队自由, 能完成大卫星所不能完成的特殊任务。因此, 编队卫星技术在导航、通信、观测、侦查等军事和民用方面都扮演着重要角色, 具有广阔的应用前景[1]。

编队的控制是编队卫星群运行的重要前提, 包括编队调整与编队保持, 编队控制的基础则是卫星间的相对运动动力学关系。1960 年, Clohessy 等[2]最早开展卫星飞行的相对运动动力学关系研究, 并针对空间两临近飞行器交会问题, 提出 Clohessy-Wiltshire方程(C-W方程)。C-W方程是线性化、常系数的微分方程组, 形式简单, 存在解析解[1], 在早期的卫星编队飞行和卫星空间交会的动力学与控制中得到广泛应用[3-4]。然而, 万有引力与卫星轨道半径之间的关系是非线性的, 因此描述卫星编队飞行的动力学方程是非线性方程[5]。相对于卫星的轨道半径, 编队卫星群的卫星间距很小, 若将万有引力引起的非线性项做泰勒展开, 高阶项部分数量级相对较小, 其线性部分占有主导地位, 因此, 卫星编队飞行的动力学方程呈现弱非线性的特点。在研究卫星编队飞行动力学时, 若采用 C-W 模型而忽略非线性项的影响, 势必带来较大的误差, 其后果是消耗更多的燃料以维持卫星编队, 影响功能的充分发挥。对卫星编队飞行动力学方程的积分, 应当充分利用其弱非线性的特点, 首先将其做泰勒展开后写成半线性的形式(即方程包含线性部分和非线性部分), 然后用精细积分法[6-7]求解线性部分,用合适的积分方法求解高阶项的非线性部分。目的是精确计算数量级较大的线性部分, 而由于高阶非线性部分的数量级较小, 误差也比较小。

针对线性常系数微分方程, 经典常微分方程理论的主要困难是指数矩阵的计算。1994年, 钟万勰等[6-7]针对指数矩阵的计算提出精细积分法(precise integration method, PIM), 用矩阵的加法代替矩阵乘法, 巧妙地避免了舍入误差。精细积分法有多方面的优势, 如可以采用大步长, 计算精度高, 稳定性好[8]。另外, Ashi 等[9]比较了 6 种常用的指数矩阵的求解方法, 其中矩阵分解(matrix decomposition)方法无论在大步长还是小步长都能精确计算指数矩阵, 但当矩阵独立特征向量个数小于维数时, 此方法的应用比较困难。精细积分法则不受此限制, 无论步长大小都能精确计算指数矩阵, 是一种有效的计算方法。

针对半线性微分方程, Hersch[10]于 1958 年首次提出指数积分的思想。Certaine[11]于 1960 年用常数变异公式, 首次给出指数 Adams-Moulton 方法。Hochbruck 等[12]于 2010 年综述了指数积分方法的发展, 并系统地总结了指数积分法的计算格式。由于指数积分方法的思想是精确地求解微分方程的线性部分, 而方程的刚性和高振荡性通常表现在其线性部分, 所以在指数积分方法提出的初期, 主要用于解决刚性问题或高振荡问题[13]。在指数积分方法的计算格式中, 指数矩阵往往难以精确计算, 尤其在维数较大的情况下。这限制了指数积分方法的进一步发展, 因此在指数积分方法提出后的很长时间都没有得到重视。21世纪初, 由于新的针对指数矩阵的计算方法的提出[9], 指数积分方法得到快速发展。目前指数积分方法在各个领域得到广泛应用, 如延时问题[14]、高振荡问题[15]、抛物型方程的求解[16]、李群算法[17]以及 Hamilton 系统的保辛算法的构造[18]等。针对此方法的研究重点由计算格式的构造逐渐转移到其计算性能的研究, 如稳定性[19-20]、收敛性[21]等。

本文在精细积分法与指数积分方法的基础上,提出精细指数积分方法, 并将其用于弱非线性方程的求解, 为弱非线性方程的求解提供一种有效的计算方法, 为卫星编队飞行动力学仿真提供一种精确的算法。

1　卫星编队飞行动力学方程

卫星编队飞行动力学主要研究伴随卫星(follower satellite)相对于参考卫星(reference satellite)的相对运动与受力之间的关系。如图 1 所示, 假设参考卫星在圆形轨道上运行, 轨道半径为 r = |r|, O2X2Y2Z2为轨道坐标系, 其中O2Y2轴的方向与地心到参考卫星的方向一致, O2X2方向为参考卫星运动的反方向, O2Z2方向由右手螺旋定则确定。在轨道坐标系中, 伴随卫星的位置矢量为r2= [x, y, z]T, 在不考虑摄动力和控制力的情况下, 伴随卫星相对于参考卫星运动的Lagrange函数为

其中,μ为地球引力常数, m为伴随卫星的质量, θ˙为参考卫星的定常轨道角速度, 等号右边第一项为伴随卫星的动能, 第二项为势能。由Lagrange方程可导出伴随卫星相对于参考卫星的动力学方程:

方程(2)的非线性项在于万有引力引起的方程右端的分母。卫星编队飞行的距离通常比参考卫星的轨道半径小得多, 即, 故方程(2)右端的分母部分表现为弱非线性。为了简化方程, 寻求解析解, 一般可将其做泰勒展开并保留线性项, 可得

由于方程的线性化会引入误差, 当误差积累到一定程度时需要对卫星施加控制, 从而消耗更多的燃料以维持编队, 甚至影响编队卫星群功能的实现。在求解卫星编队飞行的非线性方程时, 为了方便精细指数积分方法的使用, 将方程(2)写成半线性的形式:

方程(4)右端的第一项为方程(2)的线性部分; 第二项和第三项合起来为方程(2)的高阶非线性部分, 由于方程表现为弱非线性, 该部分的数值比线性部分小若干个数量级。

2　精细积分法

对于卫星编队飞行的线性化动力学方程(式(3)), 可以利用精细积分法, 计算得到伴随卫星相对于参考卫星的近似运动规律。精细积分法的详细推导过程[6-7]如下。

将方程(3)降阶成一阶微分方程组:

其中,

若能精确求解矩阵 T, 则可以精确求解常系数齐次线性微分方程(5)。精细积分法的核心正在于矩阵T的精确求解, 求解思路是利用指数函数的加法原理, 将时间步长τ分成2n等分, 令则有

由于tΔ是一个很小的时间段, 所以在 exp(LΔt)的泰勒展开式中直接取其前 5 项, 即可达到足够的精度, 即

其中, I为单位矩阵。与I相比, T1是一个很小的矩阵。为防止计算过程中“大数吃小数”, 要避免T1与I直接相加。由方程(7)和(8)可知:

若有矩阵序列:

则根据式(9)有如下递推公式:

通过求此矩阵序列得到 Tn+1, 最后与单位矩阵I相加, 便可计算矩阵 T 的值。在此过程中, 巧妙地运用指数函数的性质, 避免了严重的舍入误差, 所以这种算法有很高的精度, 一般取n=20可保证算法良好的性能。由于计算过程中在每个时间步长τ 内插入220个点, 所以即使步长比较大, 也不影响算法的精确性。精细积分法求解一般的常系数齐次线性微分方程, 能够达到计算机意义上的精度。

3　精细指数积分方法

精细积分法只能处理卫星编队飞行的线性化动力学方程。为了更加精确地求解伴随卫星相对于参考卫星的动力学规律, 还需考虑其非线性因素。将半线性微分方程(4)降阶, 且写成如下形式:

其中, L与方程(5)中的L相同, N(u, t)为方程的非线性部分, 当方程为弱非线性时有

对时间离散后, 方程(11)的精确解可写为如下形式:

其中, 右端的积分项称为Duhamel积分。用不同的方法求解, 可得到不同的指数积分方法, 目前应用较广的是指数 Runge-Kutta 方法。s 级的指数 Runge-Kutta方法格式如下:

其中, aij(Lτ)和bi(Lτ)均为待定矩阵, 且为 L 和τ 的函数, ci为待定常数。当N(u, t)=0时, 指数Runge-Kutta 方法退化为常系数齐次线性微分方程的精确解, 当 L=0 时, 此方法退化为经典的 Runge-Kutta方法。

根据 aij(Lτ)和 bi(Lτ)给出的方法不同, 指数Runge-Kutta方法可以分成不同的类型, 其中包括两种最常用的类型: 一种是积分因子法(integration factor, IF), 也称为 Lawson 方法[22], 这种方法的待定矩阵完全根据已有的 Runge-Kutta 方法确定, 构造和应用都比较简单; 另一种是指数时间微分方法(exponential time differencing, ETD), 此方法通过常数变异公式, 将非线性部分用代数插值多项式代替, 需要定义φ函数[12,19]。本文只对积分因子法做简单介绍。

积分因子法中的 aij(Lτ)和 bi(Lτ)按如下规律[19]给出:

员工是企业发展的关键，企业发展的动力是广大员工，企业员工的管理工作十分重要。当然，在知识经济的大环境下，对员工管理所存在的各种问题，人们逐渐有了全新的认识。当前社会，很多企业都十分关心企业客户的满意度。随着整个市场竞争活动的不断加剧，越来越多的企业认识到完善的经营理念并非只有外部因素。因此，企业经营者不仅需要关注外部环境的要求，还需要关注客户与企业内部的员工管理活动，探索适合当前企业经营的新道路。

其中, aij, bi, ci与经典Runge-Kutta方法的系数相同,即给定一种经典的Runge-Kutta方法就有一种积分因子法与之对应, 且当 N(u, t)=0 时积分因子法退化为相应的经典的Runge-Kutta方法。例如, 经典四级四阶的Runge-Kutta方法所对应的积分因子法如下:

4　卫星编队飞行数值仿真

在本节中, 先通过一个弱非线性的算例, 验证本文提出的积分因子法的有效性, 然后将此方法用于卫星编队飞行动力学方程的求解。

4.1积分因子法有效性的验证

考虑以下弱非线性微分方程:

方程有解析解:

当 a 取值较小时, 方程的非线性部分数量级比较小, 此处取 a=0.1。分别采用经典的四级四阶Runge-Kutta 方法(RK4_4)、八级六阶 Runge-Kutta方法(RK8_6)[23]、三级六阶隐式 Runge-Kutta 方法(RK3_6)[23]、与 RK4_4 对应的积分因子法(IF4_4)、与 RK3_6 对应的积分因子法(IF3_6), 在不同的时间步长下求解方程(16), 结果如图2所示。

图2为以上各种方法相对误差绝对值的最大值与积分时间步长之间的关系。由图 2 可知, RK4_4与IF4_4、RK3_6与IF3_6基本上平行, 说明积分因子法与对应的 Runge-Kutta 方法具有相同的阶数。在相同的步长时, IF4_4比RK4_4、IF3_6比RK3_6的计算误差都低4个数量级左右, 说明在弱非线性情况下, 指数 Runge-Kutta 方法比 Runge-Kutta 方法在精度上有比较明显的优势。将 IF4_4与 RK3_6 进行比较, 在较大步长时, IF4_4 方法比RK3_6 更精确; 在中等步长时, 两种方法精度相当;而在较小步长时, RK3_6更有优势。

为了分析方程非线性部分的数值大小与各种方法求解相对误差之间的关系, 取步长为0.1τ=, a在0.01~1 之间变化, 同样用上述5种方法求解方程(16), 结果如图3所示。

由图 3 可知, 对于普通的 Runge-Kutta 方法,方程非线性的强弱对求解结果的相对误差影响不大, 然而, 精细指数积分方法求解结果的相对误差随非线性部分的数值增大而增加。方程的非线性越弱, 精细指数积分方法求解结果越准确。与同阶的Runge-Kutta 方法相比, 精细指数积分方法在精度上有优势。

4.2卫星编队调整动力学仿真

在卫星编队飞行的非线性动力学方程中, 对于Runge-Kutta 方法, 方程(2)与方程(4)没有本质上的差别, 为了便于计算, 可直接积分方程(2); 对于指数Runge-Kutta方法, 先将方程写成半线性形式(即方程(4)), 再进行计算。由于非线性动力学方程没有解析解, 所以用RK3_6和RK8_6两种高阶算法的计算结果作为参考, 主要比较RK4_4与IF4_4的计算结果。另外, 用精细积分法(PIM)求解线性化的动力学方程(方程(3)), 研究方程的线性化带来的误差。

假设在一次编队调整过程中, 参考卫星的轨道半径为8000 km, 伴随卫星相对于参考卫星的初始状态为

积分步长为τ =500 s, 经过 4 小时的运行, 轨迹图如图4所示。的结果, 与其他曲线均有明显差距。由图 4 可知,经过 4 小时运行后, 线性化方程的计算误差为 102m 量级。从图 4 的局部放大图中可看到, IF4_4 的计算轨迹与RK3_6和RK8_6的计算轨迹基本上重合, 而RK4_4的误差为 10 m 量级。对比IF4_4与RK4_4 的计算结果可知, 在卫星编队飞行动力学方程的求解中, 精细指数积分方法能提供更好的计算精度。

图 4 中PIM即为精细积分法求解线性化方程

4.3卫星编队保持动力学仿真

当参考卫星的轨道半径为 8000 km 时, 伴随卫星相对于参考卫星的初始状态为即可使伴随卫星相对于参考卫星的运动轨道为周期轨道。分别用RK3_6, RK4_4和IF4_4方法求解动力学方程, 取积分步长为τ = 600 s, 积分时间为100小时, 得到伴随卫星相对于参考卫星的运动轨迹图(图5)。

由图 5 可知, 六阶的 RK3_6 方法和四阶的IF4_4方法能精确地模拟周期轨道, RK4_4方法求解的周期轨道则出现较大的误差。由此可知, 在弱非线性问题的求解上, 与同阶的 Runge-Kutta 方法对比, 指数积分方法能提供较高的精度, 因此在求解时可以采用较大的步长。

5　结论

本文将精细积分法用于指数积分方法中的指数矩阵的求解, 得到精细指数积分方法, 并将其用于弱非线性方程的求解。针对弱非线性问题, 例如卫星编队飞行的动力学方程, 精细积分法不再适用。如果忽略其弱非线性的性质, 直接用 Runge-Kutta方法求解, 则线性部分也会引入误差。

精细指数积分方法能精确求解方程的线性部分, 并用合适的算法求解非线性部分, 保证了算法的精度。与同阶的 Runge-Kutta 方法对比, 方程的非线性越弱, 精细指数积分方法在精度上越有优势。在卫星编队飞行的动力学仿真中, 精细指数积分方法取得了比同阶的Runge-Kutta方法更精确的结果。

[1] 李俊峰, 高云峰, 宝音贺西, 等. 卫星编队飞行动力学与控制研究. 力学与实践, 2002, 24(2): 1-6

[2] Clohessy W H, Wiltshire R S. Terminal guidance system for satellite Rendezvous. Journal of the Aerospace Sciences, 1960, 27(9): 653-658

[3] 周建平. 天宫一号/神舟八号交会对接任务总体评述. 载人航天, 2012, 18(1): 1-5

[4] 王博, 邓子辰, 李文成, 等. Hill’s方程的Magnus积分. 西北工业大学学报, 2011, 29(6): 988-991

[5] 林来兴. 交会对接动力学模型和动力学特性. 中国空间科学技术, 1994, 14(3): 47-53

[6] 钟万勰. 结构动力方程的精细时程积分法. 大连理工大学学报, 1994, 34(2): 131-136

[7] Zhong Wanxie, Williams F W. Precise time step integration method. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 1994, 208(6): 427-430

[8] Zhong Wanxie. On precise integration method. Journal of Computational and Applied Mathematics, 2004, 163(1): 59-78

[9] Ashi H A, Cummings L J, Matthews P C. Comparison of methods for evaluating functions of a matrix exponential. Applied Numerical Mathematics,2009, 59(3): 468-486

[10] Hersch J. Contribution à la méthode des équations aux différences. Zeitschrift für angewandte Mathematik und Physik ZAMP, 1958, 9(2): 129-180

[11] Certaine J. The solution of ordinary differential equations with large time constants. Mathematical Methods for Digital Computers, 1960, 1: 128-132

[12] Hochbruck M, Ostermann A. Exponential integrators. Acta Numerica, 2010, 19(1): 209-286

[13] Rahrovani S, Abrahamsson T, Modin K. An efficient exponential integrator for large nonlinear stiff systems Part 1: Theoretical investigation.Conference Proceedings of the Society for Experimental Mechanics Series, 2014, 2: 259-268

[14] Xu Yang, Zhao Jingjun, Sui Zhenan. Exponential Runge-Kutta methods for delay differential equations. Mathematics and Computers in Simulation, 2010, 80(12): 2350-2361

[15] Frénod E, Hirstoaga S A, Sonnendrücker E. An exponential integrator for a highly oscillatory Vlasov equation. Discrete and Continuous Dynamical Systems: Ser S, 2015, 8(1): 169-183

[16] Luan V T, Ostermann A. Explicit exponential Runge-Kutta methods of high order for parabolic problems. Journal of Computational and Applied Mathematics, 2014, 256: 168-179

[17] Celledoni E, Marthinsen A, Owren B. Commutatorfree Lie group methods. Future Generation Computer Systems, 2003, 19(3): 341-352

[18] Miyatake Y. An energy-preserving exponentiallyfitted continuous stage Runge–Kutta method for Hamiltonian systems. BIT Numerical Mathematics, 2014, 54(3): 777-799

[19] Maset S, Zennaro M. Unconditional stability of explicit exponential Runge-Kutta methods for semilinear ordinary differential equations. Mathematics of Computation, 2009, 78: 957-967

[20] Maset S, Zennaro M. Stability properties of explicit exponential Runge-Kutta methods. IMA Journal of Numerical Analysis, 2013, 33(1): 111-135

[21] Berland H, Owren B, Skaflestad B. B-Series and order conditions for exponential integrators. SIAM Journal on Numerical Analysis, 2005, 43(4): 1715-1727

[22] Lawson J D. Generalized Runge-Kutta processes for stable systems with large Lipschitz constants. SIAM J Numer Anal, 1967, 4: 372-380

[23] 马振华. 现代应用数学手册: 计算与数值分析卷.北京: 清华大学出版社, 2005

Precise Exponential Integrator and Its Application in Dynamics of Spacecraft Formation Flying

DENG Zichen†， LI Qingjun

Department of Engineering Mechanics, Northwestern Polytechnical University, Xi’an 710072; † E-mail: dweifan@nwpu.edu.cn

The dynamic equations of spacecraft formation flying are weakly nonlinear equations since the distance between spacecrafts is quite small compared with the orbital radius of the spacecrafts. To solve weakly nonlinear equations effectively, a precise exponential integrator (PEI) was proposed. Precise integration method (PIM) was applied to calculate exponential function in the formulas of exponential integrators (EI). Firstly, PEI was validated by solving a weakly nonlinear equation compared with Runge-Kutta method. Secondly, the dynamic equations of spacecraft formation flying were obtained through Lagrange equations, and then the equations were tansfered into semi-linear form. Ultimately, PEI and Runge-Kutta method were comparatively used to solve these equations. Through numerical analysis, PEI gave higher precision of the dynamic equations of spacecraft formation flying, indicating that PEI can be applied to other weakly nonlinear problems as well.

exponential integrator; precise integration method; spacecraft formation flying; Runge-Kutta method

O316

10.13209/j.0479-8023.2016.069

国家自然科学基金(11432010)资助

2015-10-07;

2016-02-03; 网络出版日期: 2016-07-14