沈紅新，李恒年

(西安衛(wèi)星測控中心宇航動力學國家重點實驗室，西安 710043)

靜止衛(wèi)星小推力多圈轉(zhuǎn)移軌道間接優(yōu)化

沈紅新，李恒年

(西安衛(wèi)星測控中心宇航動力學國家重點實驗室，西安 710043)

本文研究了電推進平臺地球靜止(GEO)衛(wèi)星轉(zhuǎn)移軌道優(yōu)化問題。在動力學模型中考慮了地球J2項攝動，基于龐特里壓金極大值原理(PMP)推導了時間最優(yōu)變軌控制律及其一階必要性邊界條件，采用推力幅值延拓技術(shù)和牛頓下山法對最優(yōu)控制兩點邊值問題求解。仿真算例表明本文方法的有效性，相對已有結(jié)果性能指標可以提升5%左右。分析了靜止衛(wèi)星小推力多圈時間最優(yōu)轉(zhuǎn)移軌道要素變化特性，可供工程實際參考。

電推進; 地球靜止衛(wèi)星；間接法；兩點邊值問題；軌道優(yōu)化

0 引 言

隨著電推進技術(shù)的成熟發(fā)展，電推進技術(shù)的應(yīng)用逐漸從位置保持向完成整個軌道轉(zhuǎn)移控制任務(wù)轉(zhuǎn)變。2012年3月，波音衛(wèi)星系統(tǒng)公司在一次商業(yè)通信衛(wèi)星競標中推出了全球首款全電推進平臺——BSS702SP平臺，采用該平臺，亞洲廣播衛(wèi)星-3A和歐洲通信衛(wèi)星-115B在2015年3月1日由獵鷹9號火箭以一箭雙星方式成功發(fā)射[1]。以往傳統(tǒng)的做法是火箭將靜止衛(wèi)星發(fā)射到地球同步轉(zhuǎn)移軌道(Geosynchronous transfer orbit, GTO)以后，通常采用化學推進發(fā)動機在遠地點進行3～5次變軌可以將地球同步轉(zhuǎn)移軌道變成地球靜止軌道(Geostationary orbit, GEO)。如果用電推進代替化學推進，由于電推進高比沖的特點，可以節(jié)省大量燃料。但電推進系統(tǒng)推力很小，一般只有幾百毫牛，遠遠小于化學推進的幾百牛，因而采用電推進從GTO變軌至GEO的時間將長達數(shù)月，且變軌策略與化學推進的策略有很大區(qū)別。

地球靜止衛(wèi)星軌道轉(zhuǎn)移是一類經(jīng)典問題，很多學者對此進行了大量研究，但以往的研究大多針對的情況是推力比較大，軌道轉(zhuǎn)移圈次只有幾十圈的情況，例如文獻[2-3]，或者利用月球近旁轉(zhuǎn)向技術(shù)以節(jié)省發(fā)射地球靜止衛(wèi)星的速度增量[4]。以目前電推進系統(tǒng)的水平，推力還不能達到幾牛的量級。由于初值估計偏差隨時間的傳播和累積，一般來說軌道轉(zhuǎn)移時間越長、圈次越多，優(yōu)化難度越大。因此很多學者將轉(zhuǎn)移軌道根數(shù)進行均化處理，主要的工作見文獻[5-7]，這樣得到的解對估計時間和燃料消耗固然有意義，但也會一定程度上降低軌道動力學模型的準確性。另外一種近似估計全電推進GEO衛(wèi)星轉(zhuǎn)移時間和燃料消耗的方法是解析法[8-9]，但是這類方法是基于小偏心率的假設(shè)，因此只適用于低地球軌道(LEO)向GEO轉(zhuǎn)移的情況。針對推力大小更接近實際的情況(1 N以下)，田百義等[10]基于李雅普諾夫優(yōu)化原理進行了GTO-GEO轉(zhuǎn)移軌道的優(yōu)化設(shè)計。Graham和Rao[11]采用偽譜法優(yōu)化多圈靜止衛(wèi)星轉(zhuǎn)移軌道，偽譜法本質(zhì)上是配點法的一種，解的最優(yōu)性只能事后驗證。

目前，在非線性最優(yōu)控制領(lǐng)域主要有兩大類數(shù)值方法：直接法和間接法[12-14]。直接法一般采用參數(shù)化方法將連續(xù)空間的非線性最優(yōu)控制問題轉(zhuǎn)換成離散空間的非線性規(guī)劃問題，雖然它具有容易實現(xiàn)和收斂半徑較大的優(yōu)點，但是直接法不能保證所獲的解的最優(yōu)性。相對地，間接法主要是構(gòu)造和求解由最優(yōu)控制的一階必要條件得到的Hamilton邊值問題，其優(yōu)點是：1)構(gòu)造邊值問題而不是參數(shù)優(yōu)化問題，理論上能夠準確快速獲得最優(yōu)解；2)由于采用最優(yōu)控制理論(這里主要指變分法和極大值原理)，間接法能夠提供優(yōu)化問題必要的理論信息。但間接法也有其缺點，獲得最優(yōu)性條件的過程復雜繁瑣，而且收斂半徑較小[14]。本文利用間接法研究靜止衛(wèi)星最短時間小推力多圈轉(zhuǎn)移軌道，推導了一階最優(yōu)性必要條件，采用牛頓下山法和推力幅值延拓方法求解兩點邊值問題，求解結(jié)果和文獻中相同算例對比結(jié)果表明，性能指標可以提升5%左右。

1 動力學模型

電推進系統(tǒng)可以提供的推力加速度在10-4m/s2的量級，與自然攝動加速度的量級相當，因此可以將電推進加速度和攝動加速度統(tǒng)一處理。本文采用高斯攝動方程作為軌道控制模型，形式如下[3]：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

式中自變量定義為：

(8)

ex=ecos(ω+Ω)

(9)

ey=esin(ω+Ω)

(10)

ix=tan(i/2)cosΩ

(11)

iy=tan(i/2)sinΩ

(12)

L=Ω+ω+θ

(13)

為了簡化公式形式，定義中間變量

(14)

X=ixsinL-iycosL

(15)

(16)

式中：a、e、i、Ω、ω和θ分別表示軌道半長軸、偏心率、傾角、升交點赤經(jīng)、近地點幅角和真近點角，μ為地球引力常數(shù)，m表示衛(wèi)星質(zhì)量，c表示推力器比沖Isp和地面重力加速度ge的乘積。fR、fT和fN分別表示徑向、切向和法向加速度，其中包含發(fā)動機推力加速度和地球J2項攝動加速度，即fR=fr+frJ2，fT=ft+ftJ2，fN=fn+fnJ2，其中發(fā)動機推力加速度為

(17)

(18)

(19)

式中：F為發(fā)動機推力大小，α和β表示推力方向角，定義推力方向角β為推力矢量與軌道面的夾角，推力方向角α為推力矢量在軌道面內(nèi)投影與地心矢徑方向的夾角。地球J2項攝動加速度在徑向、切向和法向的分量如下：

(20)

(21)

(22)

2 最優(yōu)控制律及邊界條件

假設(shè)初始軌道為GTO，初始軌道根數(shù)已知，目標軌道為GEO，要求目標軌道的偏心率和傾角為0，半長軸為42164 km。性能指標是轉(zhuǎn)移時間最短，寫為最大化性能指標為

φ=-tf

(23)

利用式(1)～(7)定義Hamilton函數(shù)

(24)

整理式(24)可得

H=λR(fr+frJ2)+λT(ft+ftJ2)+

(25)

式中：

(26)

(27)

λN= -λexnXey+λeynXex+

(28)

根據(jù)極大值原理[15]，最優(yōu)控制律應(yīng)使得H最大，因此最優(yōu)推力控制為

(29)

(30)

(31)

式中：Fmax表示推力最大值。據(jù)此Hamilton函數(shù)可以再次重寫為

(32)

由式(23)得

Hf=1

(33)

如果終端相位不作約束，則

λL f=0

(34)

此外，終端質(zhì)量不作約束，因此

λm f=0

(35)

3 兩點邊值問題求解算法

常用的邊值問題求解方法是打靶法。打靶法通過求解一系列初值問題來獲得邊值問題的解。如果給定了積分初值，積分常微分方程(Ordinary differential equation,ODE)后檢查邊界條件的滿足程度，再利用迭代算法逐步減小偏差。積分算法采用變步長Adams積分，具體算法見文獻[16]。Colasurdo和Pastrone[17]設(shè)計了利用解析梯度的Newton法。本文的邊值問題求解器綜合了Colasurdo和Pastrone的理論成果和Newton下山法[18]。

首先對自變量進行無量綱化，可以降低迭代算法的收斂難度。不失一般性，第j段的時刻t通過下式進行轉(zhuǎn)換

(36)

狀態(tài)和協(xié)變量的ODE也隨自變量作如下變換

(37)

式中：τj=tj-tj-1。為了簡潔，把所有的因變量寫在單獨的向量中，同時為了給出一般表達式，把常因變量也包含在內(nèi)

(38)

為了方法的通用性引入常因變量y，對應(yīng)ODE有dy/dε=0。

邊值問題的微分方程可以表示為

(39)

用矢量p表示未知初值，設(shè)p的初始解為p0，令z0=p0。給定初始解后，第r次迭代參數(shù)變化為

pr+1=pr+Δp

(40)

式中參數(shù)變化的大小Δp可以通過下式求得

(41)

式中：Ψ表示狀態(tài)或者協(xié)變量的等式邊界條件，注意到在迭代過程中Ψ可能不等于0，通過搜尋初值p使Ψ等于0。對[?Ψ/?p]可以用數(shù)值或解析的方式求逆?？紤]到可以分解矩陣[?Ψ/?p]為

(42)

式中：矩陣[?Ψ/?s]可通過對邊界條件求導獲得。注意到s對應(yīng)邊界點的z，[?s/?p]即是邊界點的[?z/?p]，而[?z/?p]滿足

(43)

積分式(43)的齊次微分方程可以獲得[?z/?p]，進而得到[?s/?p]。

經(jīng)典的Newton迭代法具有二次收斂速度，但是十分依賴于初始值的選取。參考Newton下山法的思想，采取兩種措施來降低迭代發(fā)散的幾率：1)當Ψr+1>2Ψr時，令

pr+1=pr+Δp/2

(44)

2)在迭代格式中引進一個小于等于1的系數(shù)，即

(45)

式中:K稱為下山因子，選取0.01≤K≤1。隨著迭代過程逐步接近最優(yōu)解時，可逐步提高K的值。Newton下山法的實質(zhì)是以降低收斂速度為代價來提高算法魯棒性，盡管如此，其收斂速度也比較快。

為了進一步降低多圈轉(zhuǎn)移兩點邊值問題求解難度，先從推力大的軌道開始求解，這樣軌道轉(zhuǎn)移圈次少，轉(zhuǎn)移時間較短，采用牛頓下山迭代法邊值問題更容易收斂。以大推力軌道解作為初值，將推力幅值不斷減小逐步求解(即推力幅值延拓)，每次延拓的邊值問題都利用牛頓下山迭代法求解，通常經(jīng)過十步左右可以求解得到實際推力的電推進轉(zhuǎn)移軌道。采用i5處理器主頻2.5 GHz CPU個人電腦編寫的Fortran程序，每次推力延拓迭代時間大約幾分鐘。

4 仿真校驗

衛(wèi)星采用全電推進方式由GTO向GEO轉(zhuǎn)移，優(yōu)化的性能指標為時間最短。假設(shè)衛(wèi)星推力大小為800 mN或者350 mN，比沖3500 s。初始軌道GTO和目標軌道GEO參數(shù)如表1所示。此算例來自文獻[10]。

表1 GTO-GEO初始軌道和目標軌道參數(shù)Table 1 Initial and final orbits for GTO-GEO transfer

當最大推力為800 mN時，初始加速度2×10-4m/s2，剩余質(zhì)量3721 kg，飛行時間138.4天，比文獻[10]中148.8天的結(jié)果節(jié)省10天，性能指標提升約6.7%。軌道參數(shù)變化如圖1所示。近地點高度全程增加，遠地點高度先增加再降低；半長軸增加的速率先快后慢，傾角減小速率先快后慢，偏心率減小速率則先慢后快。特別值得指出的是傾角的變化速率曲線是向下凹的形式，這有些類似最速下降線路徑的變化趨勢，這可能是由于隨著軌道能量越來越大，改變軌道面的效率也逐漸下降，仿真中發(fā)現(xiàn)一些次優(yōu)的結(jié)果往往呈現(xiàn)傾角速率曲線向上凸的變化。

由于未對動力學方程作均化處理，軌道要素變化圖放大觀察可以見到其短周期變化，見圖2。最大推力800 mN變軌的三維軌道位置變化如圖3所示，根據(jù)L的變化可知大約經(jīng)過197圈的變軌，遠地點高度仍然有一個先增加再降低的過程。該圖是將優(yōu)化后的軌道位置速度數(shù)據(jù)輸入STK軟件后得到的分辨率較高的三維仿真演示畫面截圖。

圖4和圖5給出了最優(yōu)軌道的平經(jīng)度協(xié)變量和質(zhì)量協(xié)變量隨時間變化圖，根據(jù)式(34)和式(35)所描述的最優(yōu)性條件，在終端時刻λL f=0，λm f=0，由圖4～5可知，最優(yōu)軌道滿足上述條件。

當推力為350 mN時，初始加速度0.875×10-4m/s2，剩余質(zhì)量3715 kg，飛行時間322.7天，比文獻[10]中338.9天結(jié)果節(jié)省15天，性能指標提升約4.5%。根據(jù)L的變化可知整個過程大約經(jīng)過了452圈的變軌，大約是800 mN轉(zhuǎn)移軌道圈次的兩倍多。半長軸、傾角和偏心率等參數(shù)變化規(guī)律與圖1中800 mN情況類似，限于篇幅這里不再給出軌道參數(shù)變化曲線。推力方向角變化如圖6所示，推力方向角變化存在周期性規(guī)律，這點從局部放大圖中更容易觀察，從圖6的局部放大圖可以看出,推力方向角的變化周期接近半天。

需要指出的是，目前的動力學模型中只考慮了地球J2項攝動，雖然在衛(wèi)星多圈轉(zhuǎn)移過程中日月引力攝動、太陽光壓攝動對衛(wèi)星半長軸、傾角和偏心率都有不可忽視的影響，然而在長時間轉(zhuǎn)移過程中考慮全攝動模型從工程角度來說并無必要，這固然是由于采用全攝動模型設(shè)計長時間多圈轉(zhuǎn)移軌道十分困難，同時也因為控制執(zhí)行過程中推力大小和方向的標定等也不可避免地存在誤差累積，因此目前的多圈軌道設(shè)計的任務(wù)是將衛(wèi)星送到近同步軌道，對衛(wèi)星重新進行精確軌道確定之后，在之后的定點捕獲階段將衛(wèi)星準確控制到指定位置。定點捕獲階段轉(zhuǎn)移時間一般只有幾天，采用解析攝動模型分析即可，具體模型可參見文獻[19]。

5 結(jié) 論

對于利用電推進的小推力多圈地球靜止軌道衛(wèi)星轉(zhuǎn)移軌道優(yōu)化問題，本文利用極大值原理給出了最優(yōu)控制律所應(yīng)滿足的微分方程和邊界條件，按照這種方法設(shè)計，可以在收斂的前提下從理論上保證得到的解是最優(yōu)解。文中給出的仿真實例說明了該方法的有效性，并且指標優(yōu)于已有結(jié)果。

GTO-GEO多圈轉(zhuǎn)移一般先增加遠地點高度，然后再降低，這樣可以在軌道速度更低的地方調(diào)整軌道面，降低調(diào)整軌道面的燃料消耗。最短時間的多圈轉(zhuǎn)移軌道全程開機，主要軌道根數(shù)全程調(diào)整，描述軌道形狀和大小的軌道根數(shù)不存在持續(xù)不變的情況。

[1] 段傳輝，陳荔瑩. GEO衛(wèi)星全電推進技術(shù)研究及啟示[J]. 航天器工程, 2013, 22(3): 99-104. [Duan Chuan-hui, Chen Li-ying. Research and inspiration of all- electric propulsion technology for GEO satellite[J]. Spacecraft Engineering, 2013, 22(3): 99-104.]

[2] Yue X, Yang Y, Geng Z. Indirect optimization for finite-thrust time-optimal orbital maneuver[J]. Journal of Guidance, Control, and Dynamics, 2010, 33(2): 628-634.

[3] Gao Y. Advances in low-thrust trajectory optimization and flight mechanics[D]. Columbia: University of Missouri-Columbia, 2004.

[4] 曾國強，郗曉寧，任萱. 月球近旁轉(zhuǎn)向技術(shù)研究[J]. 宇航學報, 2000, 21(4): 107-110. [Zeng Guo-qiang, Xi Xiao-ning, Ren Xuan. A study on lunar swing-by technique[J]. Journal of Astronautics, 2000, 21(4): 107-110.]

[5] Kluever C A. Geostationary orbit transfers using solar electric propulsion with specific impulse modulation[J]. Journal of Spacecraft and Rockets, 2004, 41(3): 461-466.

[6] Kluever C A, Oleson S R. Direct approach for computing near-optimal low-thrust Earth-orbit transfers[J]. Journal of Spacecraft and Rockets, 1998, 35(4): 509-515.

[7] Edelbaum T N. Propulsion requirements for controllable satellites[J]. ARS Journal, 1961, 31: 1079-1089.

[8] Kechichian J A. Reformulation of Edelbaum’s low-thrust transfer problem using optimal control theory[J]. Journal of Guidance, Control, and Dynamics, 1997, 20(5): 988-994.

[9] Casalino L, Colasurdo G. Improved Edelbaum’s approach to optimize low Earth geostationary orbits low-thrust transfers[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(5): 1504-1510.

[10] 田百義，雪丹，黃美麗. 全電推進GEO衛(wèi)星的變軌策略研究[J]. 航天器工程, 2015, 24(2):7-13. [Tian Bai-yi, Xue Dan, Huang Mei-li. Orbit transfer strategies for GEO satellites using all-electric propulsion[J]. Spacecraft Engineering, 2015, 24(2): 7-13.]

[11] Graham K F, Rao A V. Minimum-time trajectory optimization of multiple revolution low-thrust Earth-orbit transfers[J]. Journal of Spacecraft and Rockets, 2015, 52(3): 711-727.

[12] 雍恩米, 陳磊, 唐國金. 飛行器軌跡優(yōu)化數(shù)值方法綜述[J]. 宇航學報, 2008, 29(2): 397-406. [Yong En-mi, Chen Lei, Tang Guo-jin. A survey of numerical methods for trajectory optimization of spacecraft[J]. Journal of Astronautics, 2008, 29(2): 397-406.]

[13] 李俊峰, 蔣方華. 連續(xù)小推力航天器探測軌道優(yōu)化方法綜述[J]. 力學與實踐, 2011, 33(3): 1-6. [Li Jun-feng, Jiang Fang-hua. Survey of low-thrust trajectory optimization methods for deep space exploration[J]. Mechanics in Engineering, 2011, 33(3): 1-6.]

[14] Betts J T. Survey of numerical methods for trajectory optimization[J]. Journal of Guidance, Control, and Dynamics, 1998, 21(2): 193-207.

[15] Bryson A E, Ho Y C. Applied optimal control[M]. Washington, DC: Hemisphere, 1975.

[16] Shampine L F, Gordon M K. Computer solution of ordinary differential equations: the initial value problems[M]. San Francisco: Freeman, 1975.

[17] Colasurdo G, Pastrone D. Indirect optimization method for impulsive transfer[R]. AIAA Paper 1994-3762, 1994.

[18] 沈紅新. 基于解析同倫的月地應(yīng)急返回軌跡優(yōu)化方法[D]. 長沙：國防科技大學. 2014. [Shen Hong-Xin. Optimization method for the Moon-Earth abort return trajectories based on analytic homotopic technique[D]. Changsha: National University of Defense Technology, 2014.]

[19] 李恒年. 地球靜止衛(wèi)星軌道與共位控制技術(shù)[M]. 北京：國防工業(yè)出版社, 2010:173-190.

IndirectOptimizationofLow-ThrustMulti-RevolutionOrbitTransfersforGeostationary-OrbitSatellites

SHEN Hong-xin, LI Heng-nian

(State Key Laboratory of Astronautic Dynamics, Xi’an Satellite Control Center, Xi’an 710043, China)

A method for the optimization of geostationary orbit (GEO) satellite orbit transfer with all-electric propulsion is proposed. The Earth’s oblateness (J2) perturbation is included in the dynamic model. The time-optimal control law of the GEO satellites as well as the first order optimality necessary conditions, based on the Pontryagin’s maximum principle (PMP), are deduced. The resulting two-point boundary value problem is solved based on thrust continuation approach and by means of the Newton’s downhill method. Numerical simulations have demonstrated the effectiveness of the approach proposed, with performance increasing about 5% compared to the existing results. From the numerical results, some characteristics of the time-optimal transfer orbit element changes are drawn, which could provide reference for aerospace engineering.

Electric propulsion; Geostationary-orbit satellite; Indirect method; Two-point boundary value problem; Trajectory optimization

V412.4

A

1000-1328(2017)10- 1041- 07

10.3873/j.issn.1000-1328.2017.10.004

2017- 05- 11;

2017- 08- 10

國家自然科學基金(11702330)

沈紅新(1986-)，男，博士，工程師，主要從事航天器軌道姿態(tài)最優(yōu)控制研究。

通信地址：西安505信箱28分箱(710043)

電話：(029)84762520

E-mail:shxnudt@163.com