鄧子辰　李慶軍

西北工業(yè)大學工程力學系, 西安 710072; ? E-mail: dweifan@nwpu.edu.cn

?

精細指數(shù)積分法在衛(wèi)星編隊飛行動力學中的應用

鄧子辰?李慶軍

西北工業(yè)大學工程力學系, 西安 710072; ? E-mail: dweifan@nwpu.edu.cn

編隊飛行衛(wèi)星間的距離遠小于衛(wèi)星的軌道半徑, 其動力學方程表現(xiàn)為弱非線性。針對弱非線性方程的求解, 提出精細指數(shù)積分方法, 用精細積分法求解指數(shù)積分方法中的指數(shù)矩陣。用精細指數(shù)積分法和 Runge-Kutta方法, 在不同條件下求解弱非線性方程的算例, 驗證了精細指數(shù)積分法的有效性。通過Lagrange方程,建立衛(wèi)星編隊飛行動力學方程的半線性形式, 用精細指數(shù)積分方法與 Runge-Kutta 方法求解方程。數(shù)值計算結果表明, 與同階的 Runge-Kutta 求解弱非線性微分方程相比, 精細指數(shù)積分法具有更高的精度, 為衛(wèi)星編隊飛行動力學仿真提供了一種有效的數(shù)值算法。

指數(shù)積分方法; 精細積分法; 衛(wèi)星編隊飛行; Runge-Kutta方法

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

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

衛(wèi)星編隊飛行指若干個小衛(wèi)星按特定的編隊飛行, 各衛(wèi)星協(xié)同工作, 形成一個分布式的編隊衛(wèi)星群。與大衛(wèi)星相比, 編隊衛(wèi)星群成本低、可靠性強、靈活性高、編隊自由, 能完成大衛(wèi)星所不能完成的特殊任務。因此, 編隊衛(wèi)星技術在導航、通信、觀測、偵查等軍事和民用方面都扮演著重要角色, 具有廣闊的應用前景[1]。

編隊的控制是編隊衛(wèi)星群運行的重要前提, 包括編隊調整與編隊保持, 編隊控制的基礎則是衛(wèi)星間的相對運動動力學關系。1960 年, Clohessy 等[2]最早開展衛(wèi)星飛行的相對運動動力學關系研究, 并針對空間兩臨近飛行器交會問題, 提出 Clohessy-Wiltshire方程(C-W方程)。C-W方程是線性化、常系數(shù)的微分方程組, 形式簡單, 存在解析解[1], 在早期的衛(wèi)星編隊飛行和衛(wèi)星空間交會的動力學與控制中得到廣泛應用[3-4]。然而, 萬有引力與衛(wèi)星軌道半徑之間的關系是非線性的, 因此描述衛(wèi)星編隊飛行的動力學方程是非線性方程[5]。相對于衛(wèi)星的軌道半徑, 編隊衛(wèi)星群的衛(wèi)星間距很小, 若將萬有引力引起的非線性項做泰勒展開, 高階項部分數(shù)量級相對較小, 其線性部分占有主導地位, 因此, 衛(wèi)星編隊飛行的動力學方程呈現(xiàn)弱非線性的特點。在研究衛(wèi)星編隊飛行動力學時, 若采用 C-W 模型而忽略非線性項的影響, 勢必帶來較大的誤差, 其后果是消耗更多的燃料以維持衛(wèi)星編隊, 影響功能的充分發(fā)揮。對衛(wèi)星編隊飛行動力學方程的積分, 應當充分利用其弱非線性的特點, 首先將其做泰勒展開后寫成半線性的形式(即方程包含線性部分和非線性部分), 然后用精細積分法[6-7]求解線性部分,用合適的積分方法求解高階項的非線性部分。目的是精確計算數(shù)量級較大的線性部分, 而由于高階非線性部分的數(shù)量級較小, 誤差也比較小。

針對線性常系數(shù)微分方程, 經(jīng)典常微分方程理論的主要困難是指數(shù)矩陣的計算。1994年, 鐘萬勰等[6-7]針對指數(shù)矩陣的計算提出精細積分法(precise integration method, PIM), 用矩陣的加法代替矩陣乘法, 巧妙地避免了舍入誤差。精細積分法有多方面的優(yōu)勢, 如可以采用大步長, 計算精度高, 穩(wěn)定性好[8]。另外, Ashi 等[9]比較了 6 種常用的指數(shù)矩陣的求解方法, 其中矩陣分解(matrix decomposition)方法無論在大步長還是小步長都能精確計算指數(shù)矩陣, 但當矩陣獨立特征向量個數(shù)小于維數(shù)時, 此方法的應用比較困難。精細積分法則不受此限制, 無論步長大小都能精確計算指數(shù)矩陣, 是一種有效的計算方法。

針對半線性微分方程, Hersch[10]于 1958 年首次提出指數(shù)積分的思想。Certaine[11]于 1960 年用常數(shù)變異公式, 首次給出指數(shù) Adams-Moulton 方法。Hochbruck 等[12]于 2010 年綜述了指數(shù)積分方法的發(fā)展, 并系統(tǒng)地總結了指數(shù)積分法的計算格式。由于指數(shù)積分方法的思想是精確地求解微分方程的線性部分, 而方程的剛性和高振蕩性通常表現(xiàn)在其線性部分, 所以在指數(shù)積分方法提出的初期, 主要用于解決剛性問題或高振蕩問題[13]。在指數(shù)積分方法的計算格式中, 指數(shù)矩陣往往難以精確計算, 尤其在維數(shù)較大的情況下。這限制了指數(shù)積分方法的進一步發(fā)展, 因此在指數(shù)積分方法提出后的很長時間都沒有得到重視。21世紀初, 由于新的針對指數(shù)矩陣的計算方法的提出[9], 指數(shù)積分方法得到快速發(fā)展。目前指數(shù)積分方法在各個領域得到廣泛應用, 如延時問題[14]、高振蕩問題[15]、拋物型方程的求解[16]、李群算法[17]以及 Hamilton 系統(tǒng)的保辛算法的構造[18]等。針對此方法的研究重點由計算格式的構造逐漸轉移到其計算性能的研究, 如穩(wěn)定性[19-20]、收斂性[21]等。

本文在精細積分法與指數(shù)積分方法的基礎上,提出精細指數(shù)積分方法, 并將其用于弱非線性方程的求解, 為弱非線性方程的求解提供一種有效的計算方法, 為衛(wèi)星編隊飛行動力學仿真提供一種精確的算法。

1　衛(wèi)星編隊飛行動力學方程

衛(wèi)星編隊飛行動力學主要研究伴隨衛(wèi)星(follower satellite)相對于參考衛(wèi)星(reference satellite)的相對運動與受力之間的關系。如圖 1 所示, 假設參考衛(wèi)星在圓形軌道上運行, 軌道半徑為 r = |r|, O2X2Y2Z2為軌道坐標系, 其中O2Y2軸的方向與地心到參考衛(wèi)星的方向一致, O2X2方向為參考衛(wèi)星運動的反方向, O2Z2方向由右手螺旋定則確定。在軌道坐標系中, 伴隨衛(wèi)星的位置矢量為r2= [x, y, z]T, 在不考慮攝動力和控制力的情況下, 伴隨衛(wèi)星相對于參考衛(wèi)星運動的Lagrange函數(shù)為

其中,μ為地球引力常數(shù), m為伴隨衛(wèi)星的質量, θ˙為參考衛(wèi)星的定常軌道角速度, 等號右邊第一項為伴隨衛(wèi)星的動能, 第二項為勢能。由Lagrange方程可導出伴隨衛(wèi)星相對于參考衛(wèi)星的動力學方程:

方程(2)的非線性項在于萬有引力引起的方程右端的分母。衛(wèi)星編隊飛行的距離通常比參考衛(wèi)星的軌道半徑小得多, 即, 故方程(2)右端的分母部分表現(xiàn)為弱非線性。為了簡化方程, 尋求解析解, 一般可將其做泰勒展開并保留線性項, 可得

由于方程的線性化會引入誤差, 當誤差積累到一定程度時需要對衛(wèi)星施加控制, 從而消耗更多的燃料以維持編隊, 甚至影響編隊衛(wèi)星群功能的實現(xiàn)。在求解衛(wèi)星編隊飛行的非線性方程時, 為了方便精細指數(shù)積分方法的使用, 將方程(2)寫成半線性的形式:

方程(4)右端的第一項為方程(2)的線性部分; 第二項和第三項合起來為方程(2)的高階非線性部分, 由于方程表現(xiàn)為弱非線性, 該部分的數(shù)值比線性部分小若干個數(shù)量級。

2　精細積分法

對于衛(wèi)星編隊飛行的線性化動力學方程(式(3)), 可以利用精細積分法, 計算得到伴隨衛(wèi)星相對于參考衛(wèi)星的近似運動規(guī)律。精細積分法的詳細推導過程[6-7]如下。

將方程(3)降階成一階微分方程組:

其中,

若能精確求解矩陣 T, 則可以精確求解常系數(shù)齊次線性微分方程(5)。精細積分法的核心正在于矩陣T的精確求解, 求解思路是利用指數(shù)函數(shù)的加法原理, 將時間步長τ分成2n等分, 令則有

由于tΔ是一個很小的時間段, 所以在 exp(LΔt)的泰勒展開式中直接取其前 5 項, 即可達到足夠的精度, 即

其中, I為單位矩陣。與I相比, T1是一個很小的矩陣。為防止計算過程中“大數(shù)吃小數(shù)”, 要避免T1與I直接相加。由方程(7)和(8)可知:

若有矩陣序列:

則根據(jù)式(9)有如下遞推公式:

通過求此矩陣序列得到 Tn+1, 最后與單位矩陣I相加, 便可計算矩陣 T 的值。在此過程中, 巧妙地運用指數(shù)函數(shù)的性質, 避免了嚴重的舍入誤差, 所以這種算法有很高的精度, 一般取n=20可保證算法良好的性能。由于計算過程中在每個時間步長τ 內插入220個點, 所以即使步長比較大, 也不影響算法的精確性。精細積分法求解一般的常系數(shù)齊次線性微分方程, 能夠達到計算機意義上的精度。

3　精細指數(shù)積分方法

精細積分法只能處理衛(wèi)星編隊飛行的線性化動力學方程。為了更加精確地求解伴隨衛(wèi)星相對于參考衛(wèi)星的動力學規(guī)律, 還需考慮其非線性因素。將半線性微分方程(4)降階, 且寫成如下形式:

其中, L與方程(5)中的L相同, N(u, t)為方程的非線性部分, 當方程為弱非線性時有

對時間離散后, 方程(11)的精確解可寫為如下形式:

其中, 右端的積分項稱為Duhamel積分。用不同的方法求解, 可得到不同的指數(shù)積分方法, 目前應用較廣的是指數(shù) Runge-Kutta 方法。s 級的指數(shù) Runge-Kutta方法格式如下:

其中, aij(Lτ)和bi(Lτ)均為待定矩陣, 且為 L 和τ 的函數(shù), ci為待定常數(shù)。當N(u, t)=0時, 指數(shù)Runge-Kutta 方法退化為常系數(shù)齊次線性微分方程的精確解, 當 L=0 時, 此方法退化為經(jīng)典的 Runge-Kutta方法。

根據(jù) aij(Lτ)和 bi(Lτ)給出的方法不同, 指數(shù)Runge-Kutta方法可以分成不同的類型, 其中包括兩種最常用的類型: 一種是積分因子法(integration factor, IF), 也稱為 Lawson 方法[22], 這種方法的待定矩陣完全根據(jù)已有的 Runge-Kutta 方法確定, 構造和應用都比較簡單; 另一種是指數(shù)時間微分方法(exponential time differencing, ETD), 此方法通過常數(shù)變異公式, 將非線性部分用代數(shù)插值多項式代替, 需要定義φ函數(shù)[12,19]。本文只對積分因子法做簡單介紹。

積分因子法中的 aij(Lτ)和 bi(Lτ)按如下規(guī)律[19]給出:

員工是企業(yè)發(fā)展的關鍵，企業(yè)發(fā)展的動力是廣大員工，企業(yè)員工的管理工作十分重要。當然，在知識經(jīng)濟的大環(huán)境下，對員工管理所存在的各種問題，人們逐漸有了全新的認識。當前社會，很多企業(yè)都十分關心企業(yè)客戶的滿意度。隨著整個市場競爭活動的不斷加劇，越來越多的企業(yè)認識到完善的經(jīng)營理念并非只有外部因素。因此，企業(yè)經(jīng)營者不僅需要關注外部環(huán)境的要求，還需要關注客戶與企業(yè)內部的員工管理活動，探索適合當前企業(yè)經(jīng)營的新道路。

其中, aij, bi, ci與經(jīng)典Runge-Kutta方法的系數(shù)相同,即給定一種經(jīng)典的Runge-Kutta方法就有一種積分因子法與之對應, 且當 N(u, t)=0 時積分因子法退化為相應的經(jīng)典的Runge-Kutta方法。例如, 經(jīng)典四級四階的Runge-Kutta方法所對應的積分因子法如下:

4　衛(wèi)星編隊飛行數(shù)值仿真

在本節(jié)中, 先通過一個弱非線性的算例, 驗證本文提出的積分因子法的有效性, 然后將此方法用于衛(wèi)星編隊飛行動力學方程的求解。

4.1積分因子法有效性的驗證

考慮以下弱非線性微分方程:

方程有解析解:

當 a 取值較小時, 方程的非線性部分數(shù)量級比較小, 此處取 a=0.1。分別采用經(jīng)典的四級四階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 方法具有相同的階數(shù)。在相同的步長時, IF4_4比RK4_4、IF3_6比RK3_6的計算誤差都低4個數(shù)量級左右, 說明在弱非線性情況下, 指數(shù) Runge-Kutta 方法比 Runge-Kutta 方法在精度上有比較明顯的優(yōu)勢。將 IF4_4與 RK3_6 進行比較, 在較大步長時, IF4_4 方法比RK3_6 更精確; 在中等步長時, 兩種方法精度相當;而在較小步長時, RK3_6更有優(yōu)勢。

為了分析方程非線性部分的數(shù)值大小與各種方法求解相對誤差之間的關系, 取步長為0.1τ=, a在0.01~1 之間變化, 同樣用上述5種方法求解方程(16), 結果如圖3所示。

由圖 3 可知, 對于普通的 Runge-Kutta 方法,方程非線性的強弱對求解結果的相對誤差影響不大, 然而, 精細指數(shù)積分方法求解結果的相對誤差隨非線性部分的數(shù)值增大而增加。方程的非線性越弱, 精細指數(shù)積分方法求解結果越準確。與同階的Runge-Kutta 方法相比, 精細指數(shù)積分方法在精度上有優(yōu)勢。

4.2衛(wèi)星編隊調整動力學仿真

在衛(wèi)星編隊飛行的非線性動力學方程中, 對于Runge-Kutta 方法, 方程(2)與方程(4)沒有本質上的差別, 為了便于計算, 可直接積分方程(2); 對于指數(shù)Runge-Kutta方法, 先將方程寫成半線性形式(即方程(4)), 再進行計算。由于非線性動力學方程沒有解析解, 所以用RK3_6和RK8_6兩種高階算法的計算結果作為參考, 主要比較RK4_4與IF4_4的計算結果。另外, 用精細積分法(PIM)求解線性化的動力學方程(方程(3)), 研究方程的線性化帶來的誤差。

假設在一次編隊調整過程中, 參考衛(wèi)星的軌道半徑為8000 km, 伴隨衛(wèi)星相對于參考衛(wèi)星的初始狀態(tài)為

積分步長為τ =500 s, 經(jīng)過 4 小時的運行, 軌跡圖如圖4所示。的結果, 與其他曲線均有明顯差距。由圖 4 可知,經(jīng)過 4 小時運行后, 線性化方程的計算誤差為 102m 量級。從圖 4 的局部放大圖中可看到, IF4_4 的計算軌跡與RK3_6和RK8_6的計算軌跡基本上重合, 而RK4_4的誤差為 10 m 量級。對比IF4_4與RK4_4 的計算結果可知, 在衛(wèi)星編隊飛行動力學方程的求解中, 精細指數(shù)積分方法能提供更好的計算精度。

圖 4 中PIM即為精細積分法求解線性化方程

4.3衛(wèi)星編隊保持動力學仿真

當參考衛(wèi)星的軌道半徑為 8000 km 時, 伴隨衛(wèi)星相對于參考衛(wèi)星的初始狀態(tài)為即可使伴隨衛(wèi)星相對于參考衛(wèi)星的運動軌道為周期軌道。分別用RK3_6, RK4_4和IF4_4方法求解動力學方程, 取積分步長為τ = 600 s, 積分時間為100小時, 得到伴隨衛(wèi)星相對于參考衛(wèi)星的運動軌跡圖(圖5)。

由圖 5 可知, 六階的 RK3_6 方法和四階的IF4_4方法能精確地模擬周期軌道, RK4_4方法求解的周期軌道則出現(xiàn)較大的誤差。由此可知, 在弱非線性問題的求解上, 與同階的 Runge-Kutta 方法對比, 指數(shù)積分方法能提供較高的精度, 因此在求解時可以采用較大的步長。

5　結論

本文將精細積分法用于指數(shù)積分方法中的指數(shù)矩陣的求解, 得到精細指數(shù)積分方法, 并將其用于弱非線性方程的求解。針對弱非線性問題, 例如衛(wèi)星編隊飛行的動力學方程, 精細積分法不再適用。如果忽略其弱非線性的性質, 直接用 Runge-Kutta方法求解, 則線性部分也會引入誤差。

精細指數(shù)積分方法能精確求解方程的線性部分, 并用合適的算法求解非線性部分, 保證了算法的精度。與同階的 Runge-Kutta 方法對比, 方程的非線性越弱, 精細指數(shù)積分方法在精度上越有優(yōu)勢。在衛(wèi)星編隊飛行的動力學仿真中, 精細指數(shù)積分方法取得了比同階的Runge-Kutta方法更精確的結果。

[1] 李俊峰, 高云峰, 寶音賀西, 等. 衛(wèi)星編隊飛行動力學與控制研究. 力學與實踐, 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積分. 西北工業(yè)大學學報, 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] 馬振華. 現(xiàn)代應用數(shù)學手冊: 計算與數(shù)值分析卷.北京: 清華大學出版社, 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; 網(wǎng)絡出版日期: 2016-07-14