基于傅立葉平均法下的連續(xù)小推力動力學分析①

馮維明，李 源，苗 楠

(山東大學 工程力學系，濟南 250062)

基于傅立葉平均法下的連續(xù)小推力動力學分析①

馮維明，李 源，苗 楠

(山東大學 工程力學系，濟南 250062)

通過將小推力展開為偏近點角的傅立葉級數(shù)，并對高斯攝動方程在一個軌道周期上的平均，將原方程的推力轉化為僅由14個傅立葉系數(shù)表示的控制變量。仿真計算表明，平均化后的高斯方程使計算量與牛頓積分相比顯著減少，且對小推力而言有足夠的精度。對利用平均化后的高斯方程計算軌道根數(shù)時產(chǎn)生誤差的原因進行了研究，并進一步分析小推力的范圍和小推力近似表達式對上述誤差的影響，為今后小推力下非開普勒軌道動力學分析提供了理論依據(jù)和參數(shù)。

傅立葉級數(shù);連續(xù)小推力;平均法;高斯方程;動力學分析;非開普勒軌道

0 引言

小推力推進系統(tǒng)為許多星際航行和地球軌道任務提供了一種高效的新選擇，但是這種系統(tǒng)卻對最優(yōu)化控制提出了新的挑戰(zhàn)［1］對于某些特殊的小推力作用的情況下，軌道轉移的最優(yōu)控制問題已經(jīng)出現(xiàn)了解析法和數(shù)值近似解的方法［2-3］，但是對于一般的小推力問題，需要對每一個初始條件和推力值進行完整的數(shù)值積分，軌道參數(shù)的確定往往對這些變量很敏感。因此，普通的最優(yōu)控制法則對于確定螺旋軌道的數(shù)十甚至上百圈軌道是十分困難的［4］。

在一些特殊情形的軌道轉移中，數(shù)值解法已經(jīng)得到了挖掘。目前較多的解法應用了變分法或直接優(yōu)化法［5-7］，以確定特定約束條件下的小推力的最優(yōu)控制律。此外，文獻［8］運用李亞普諾夫反饋控制的方法，用以解決開環(huán)軌道最優(yōu)時間控制和最優(yōu)軌道問題。然而實踐證明，把平均法和其他方法結合的思路在解決對初始軌道和推力變量預測敏感的問題上是很有效的［9-10］，但所有這些方法都局限于特定的推力和軌道參數(shù)范圍。

本文討論了一種有效解決航天器在小推力下軌道動力學問題的新方法。推力加速度的各分量以偏近點角展開為傅立葉級數(shù)，而高斯變分方程則在一個軌道周期上進行平均并且通過正交條件進行簡化，從而定義出了一組久期方程。這組久期方程是一個含有14個推力傅立葉系數(shù)的函數(shù)(無論初始傅立葉級數(shù)的階數(shù)是多少)。因此，該方法將普遍形式的推力分布情況簡化為僅有14個參數(shù)的形式，這使計算量與牛頓積分相比大大減少。并以此分析了影響傅立葉系數(shù)表示的平均久期方程求解精度的主要因素。

1 數(shù)學建模

考慮一個在軌航天器，其質量對于其環(huán)繞的中心體可忽略不計。將航天器認為是一個質點。它受一個大小和方向可能隨時間變化的連續(xù)推力加速度的作用。其軌道可用牛頓運動方程來表述:

力加速度矢量F可以沿著徑向FR、法向FW和周向FS分解，即

式中r和w分別為徑向和法向的單位矢量。

牛頓方程可分解為拉格朗日行星方程，高斯形式的拉格朗日行星方程可表示為

式中 a為半長軸;e為偏心率;i為軌道傾角;Ω為升交點赤經(jīng);ω為近地點幅角;f為真近點角;E為偏近點角，而ε1+∫ndt=l為平黃經(jīng)。

平近點角是平黃經(jīng)和近地點幅角的差值:

根據(jù)傅立葉理論，對于任何一個在(0，L)上只存在有限個跳躍間斷點的分段光滑的函數(shù)f(θ)，都可以表示為一系列周期性延伸的傅立葉級數(shù)，這個級數(shù)無限逼近于函數(shù)本身，當跳躍間斷點存在時，傅立葉級數(shù)將逼近于左右極限的平均值。因此，這種表示方法可以應用于幾乎任何一種普通的小推力航天器的控制律中。對于給定的任意加速度矢量F，其各個分量都可以展開為任意時間間隔上的傅立葉級數(shù)。實際上，傅立葉級數(shù)可按時間展開，也可按隨時間變化的軌道參數(shù)展開，如真近點角、偏近點角或平近點角?？紤]一個軌道周期L=2π，以θ來表示這一任意的軌道參數(shù):

現(xiàn)在開始一階平均化分析，先假設一個加速度矢量，它可在一個軌道周期(L=2π)中表示，這個加速度矢量的數(shù)值足夠小，以使軌道的大小和形狀在一圈中的改變不是很明顯。因此，可將高斯方程相對于平近點角在一個軌道周期上平均來得到平均軌道根數(shù)。

式中 γ代表任意的軌道根數(shù)。

在對推力加速度矢量的各分量進行傅立葉展開時，所使用的軌道參數(shù)的選擇是很重要的，有的參數(shù)會使由此得出的久期方程將變得冗長而復雜(如按真近點角展開成傅立葉級數(shù)的話)，而有的參數(shù)會使久期方程會被簡化。經(jīng)比較，加速度矢量的每一個分量都是按偏近點角進行傅立葉級數(shù)展開是最佳選擇，并且各軌道參數(shù)的平均化也是按偏近點角進行，由dM=(1-ecosE)dE，式(13)可寫為

將式(11)代入式(4)～式(9)中，并由式(14)可得到關于偏近點角的平均化高斯方程。由傅立葉級數(shù)的正交性，可消掉推力加速度傅立葉展開級數(shù)中第二階以上及個別第二階傅立葉系數(shù)。因此，無論初始推力加速度傅立葉級數(shù)的階數(shù)如何，軌道根數(shù)a、e、i、Ω、ω和ε1的平均變化率只與這14個傅立葉系數(shù)有關，即(k=0，1，2)(k=0，1，2)(k=0，1，2)，(k=1，2)(1，2)則高斯平均方程為

2 計算分析

為在計算時數(shù)據(jù)的簡便起見，將方程做歸一化處理。在此，將重力參數(shù)μ歸一化常數(shù)，μ=3.986×105km3/s2=3.986×1014m3/s2，因此相當于各變量的長度(以m為單位)除以。初始的軌道根數(shù)由表1給出。

表1 初始軌道根數(shù)Table 1 initial orbit elements

設控制率為周向加速度FS=1×10-8×(1-sinE)、徑向加速度FR=1×10-8×sinE、法向加速度FW=1×10-8×(1+cosE)，上述推力的加速度分量為歸一化后的值，如航天器為1 000 kg，推力的最大值幅值為2.208 N，可算作小推力范圍。用龍格-庫塔法進行積分來估算軌道根數(shù)的變化情況。該情況下由牛頓方程和平均變分方程求解出的軌道根數(shù)分別用實線和虛線在圖1中給出。

圖1 航天器受到3個方向推力加速度作用時軌道根數(shù)吻合情況Fig.1 Osculating orbital elements of spacecraft subject to acceleration in three directions

由圖1可看到，所有的軌道根數(shù)均給出了非常吻合的結果。即用傅立葉系數(shù)表示的平均久期方程的求解結果和直接對拉格朗日方程進行積分得到的結果具有緊密的一致性。

再考慮一個簡單的周向分步加速度的問題，F(xiàn)R=FW=0，F(xiàn)S按如圖2中所示變化。每個周期中包含2次點火和2個滑行弧，初始軌道同表1所示。切向分步推力加速度作用時軌道根數(shù)吻合情況見圖3。

圖2 周向分步加速度Fig.2 Step tangential acceleration

圖3 切向分步推力加速度作用時軌道根數(shù)吻合情況Fig.3 Osculating orbital elements of spacecraft subject to step tangential acceleration

由圖3知，平均久期方程和牛頓方程給出的結果具有完全相同的趨勢，數(shù)值上也相差不大。

另外，基于傅立葉平均化后的高斯方程最顯著優(yōu)勢之一就是計算量大規(guī)模減少，尤其在進行非開普勒軌道優(yōu)化計算時，計算精度和計算耗時是困擾著眾多學者的問題，考慮在表2所列初始條件下，歸一化后的周向推力為FS=5×10-9(實際推力加速度為7.36×10-5m/s2)，分別用牛頓方程和平均變分方程求解出的軌道根數(shù)及誤差列入表2。

推力表達式為如下形式:

表2 連續(xù)小推力作用下牛頓方程和平均方程10圈后的計算結果比較Table 2 Comparison of the calculated results determined by the two methods on low-thrust continuous controls after running 10 laps

航天器運行時間為65 973.4 s，約運行10圈。由表2計算結果看出，由平均化后的高斯方程計算得到的10圈后的軌道的6個根與牛頓方程精確積分所得到的結果相比，除近地點幅角ω相對誤差略大些，其他相對誤差非常之小。利用牛頓法進行精確積分計算耗時(CPU耗時)為13.556 407 s，而利用基于傅立葉平均化后的高斯方程計算耗時為0.450 808 s，后者比前者快30倍。

3 影響計算精度的主要因素

平均久期方程能夠精確有效解決小推力螺旋軌道問題，且與通常的牛頓問題的積分方法相比能大大減少計算量。但每一種近似計算方法都有其局限性，對于平均久期方程也是如此，這也是很多文獻所未曾涉及到的問題，因此很有必要對利用該方法進行軌道計算可能導致的精度進行分析。

3.1 推力幅值大小的影響

由該方法的限制條件可知，推力幅值過大，會導致周期內軌道變形太大，從而使平均法的結果產(chǎn)生累計誤差，于是，關心在保證精度條件下，推力的極限值是多少。為討論推力增大對計算結果帶來的影響，研究推力加速度大小在切向分步加速度形式(如圖2示)的控制率中對計算結果的影響。令周向加速度FS的幅值分別增大為 1 ×10-7(歸一化后)的 1.0、1.5、2.0和2.5倍，限于篇幅僅討論對軌道形狀影響最大的根數(shù)解半長軸a和偏心率e，圖4給出其吻合情況(幅值為1×10-7見圖3)。

圖4 不同切向分步推力加速度作用時軌道根數(shù)吻合情況Fig.4 Osculating orbital elements of spacecraft subject to different step tangential acceleration

為更清晰地表現(xiàn)牛頓方程與平均方程的計算結果誤差，現(xiàn)將不同幅值的推力作用下航天器運行20圈2種方法計算得到的半長軸a、偏心率e及相對誤差列入表3。

由圖4和表3可看到，隨著加速度幅值的增大，a和e的誤差都在增大，a的誤差增幅較小，而e的誤差顯著增大。在FS達到2.5×10-7時，偏心率的2種求解結果已有顯著偏離，此時對應質量為1 000 kg的航天器，推力大小約為18.4 N，約20圈后，偏心率的結果已經(jīng)完全不可靠了。另外，由圖4也可看到，即使推力較大時，在前12圈之內，平均久期方程和牛頓方程給出的結果仍然保持著良好的一致性。由表3可見，半長軸變化很大，因此，若在一定圈數(shù)內完成變軌任務，平均方程仍有足夠的可靠度，即小推力的上限應根據(jù)變軌任務確定。

表3 推力加速度大小對a、e計算誤差的影響Table 3 Influence of acceleration on the accuracy of a and e

3.2 近似推力表達式的影響

由式(15)～式(20)可知，平均化后的高斯方程僅與14個傅立葉系數(shù)有關，推力的傅立葉級數(shù)高階項會因平均過程而被消去，平均化后非零的傅立葉系數(shù)是否能精確的描述原推力形式將成為關鍵。在本文第1個算例中，由傅立葉系數(shù)表達式(12)不難得到各推力加速度矢量分量的傅立葉表達式中的系數(shù)為=1 ×10-8=1 ×10-8=-1 ×10-8=1 ×10-8和=1×10-8，而其他系數(shù)均為零。將上述非零系數(shù)帶入式(11)可得

高斯方程經(jīng)平均化后得到的推力加速度分量傅立葉級數(shù)近似表達式與原初始推力表達式完全一致，因此傅立葉級數(shù)近似表達式對計算精度沒有影響。下面再來考慮3個推力分量皆為分步加速度的形式，即

在圖5中可表示加速度分量近似表達式與原表達式的差異。圖5中用點表示的線是高斯方程平均后3個加速度分量變化曲線，注意到此時徑向推力FR已變成常數(shù)，為分段推力的平均值;而周向推力FS和法向推力FW近似為正弦函數(shù)形式，從圖形上看比較接近原推力形式。這也就是在前例中(僅有周向推力FS)計算精度較高的原因之一。平均方程與解析解得到的計算結果如圖6所示。

圖5 分段切向加速度的傅立葉級數(shù)表示Fig.5 Fourier series for step circumferential acceleration

圖6 三向分步推力加速度作用時軌道根數(shù)吻合情況Fig.6 Osculating orbital elements of spacecraft subject to step circumferential acceleration in three directions

圖6表示的是航天器運行了10圈的根數(shù)時程圖。長半軸a和偏心率e平均方程的解與牛頓積分法的解吻合的較好，軌道傾角i，升交點赤經(jīng)Ω，近地點幅角ω在10圈后誤差開始增大，較明顯的為Ω的值，平均結果開始偏離高斯方程積分結果。

4 結論

(1)將推力加速度分量展開為偏近點角表示的傅立葉級數(shù)，而高斯變分方程則在一個軌道周期上進行平均并且通過正交條件進行簡化。該方法將整個連續(xù)控制問題的參數(shù)減少為14個(無論初始傅立葉級數(shù)的階數(shù)是多少)，這使計算量與牛頓積分相比大大減少。

(2)通過平面變軌和空間變軌的計算結果分析，用該方法計算連續(xù)和非連續(xù)小推力下的軌道轉移是正確的和準確的。

(3)對造成計算誤差的主要因素進行了定性和定量分析，推力超出“小推力”的范圍和傅立葉級數(shù)下推力的近似表達式(高斯方程平均后)是影響精度的主要因素，即便如此，在一定范圍內平均法計算結果仍與牛頓方程的計算結果吻合良好。

［1］Gao Y.Advances in low-thrust trajectory optimization and flight mechanics ，dissertation thesis［R］.University of Missouri-Columbia，2003，64-12(B):6178.

［2］Craig A Kluever.Direct approach for computing near-optimal low-thrust earth-orbit transfers［J］.Journal of Spacecraft and Rockets，1998，29(1):45-61.

［3］Akella M R，Broucke R A.Anatomy of the constant radial thrust problem［J］.Journal of Guidance，Control，and Dy namics，2002，25(3):563-570.

［4］李俊峰，龔勝平.非開普勒軌道動力學與控制［J］.宇航學報，2009，30(1):47-53.

［5］Kluever C A.Optimal low-thrust interplanetary trajectories by direct method techniques［J］.Journal of the Astronautical Sciences，1997，45(3):162-247.

［6］Petropoulos A E.Some analytic integrals of the averaged variational equations for a thrusting spacecraft［R］.Interplanetary Network Progress Rept.42-150，Jet Propulsion Lab.，California Inst.of Technology，Pasadena，CA，Aug.2002:1-29.

［7］John T Betts.Very low-thrust trajectory optimization using a direct SQP method［J］.Journal of Computational and Applied Mathematics，2000，120:27-40.

［8］Gurfill P.Nonlinear feedback control of low-thrust orbital transfer in a central graviational field［J］.Acta Astronautica，2007，60(8-9):631-648.

［9］Jennifer S Hudson，Daniel J Scheeres.Reduction of lowthrust continuous controls for trajectory dynamics［J］.Journal of Guidance，Control，and Dynamics，2009，32(3):780-787.

［10］尚海濱，崔平遠，欒恩杰.基于平均法的小推力轉移軌道優(yōu)化研究［C］//25屆中國控制論文集.2006.

Dynamic analysis of continuous low-thrust based on fourier average method

FENG Wei-ming，LI Yuan，MIAO Nan
(Department of Engineering Mechanics，Shandong University，Jinan 250062，China)

Each component of the thrust vector was expanded as Fourier series in eccentric anomaly and Gauss variational equations were averaged over one orbit period，then the thrust vector was translated to a variable controlled by fourteen Fourier＇s parameters.Simulation results show that these secular equations are sufficient to accurately determine a low-thrust spiral trajectory with significantly reduced computation as compared with integration of the full Newtonian problem.In addition，error causes of orbit elements witch were calculated by the averaged Gauss equations were studied，and influence of low-thrust range and approximate expressions on the error was further analyzed，providing theoretical basis and parameters for dynamic analysis of the Non-Keplerian orbits of lowthrust.

Fourier series;continuous low-thrust;average method;Gauss equations;dynamic analysis;non-Keplerian orbits

V412

A

1006-2793(2012)03-0285-05

2011-08-18;;

2011-10-10。

國家863項目。

馮維明(1957—)，男，教授，主要研究方向為非線性動力學和軌道動力學。E-mail:fwm@sdu.edu.cn

(編輯:呂耀輝)