基于多項式設計法的小推力軌跡設計*

曾 奎，耿云海，謝成清

(哈爾濱工業大學 衛星技術研究所， 黑龍江 哈爾濱 150001)

基于多項式設計法的小推力軌跡設計*

曾 奎，耿云海，謝成清

(哈爾濱工業大學 衛星技術研究所， 黑龍江 哈爾濱 150001)

形狀逼近法是小推力軌跡設計中的一種有效方法，然而現有的方法大都假定運動軌跡為某一特定的形狀，而且沒有考慮推力加速度的約束限制。針對小推力軌道交會問題，提出一種基于多項式的軌跡設計方法。結合極坐標系，建立基于多項式的三自由度軌跡運動模型，將軌跡設計問題轉化為求解多項式的系數問題；根據運動模型推導軌跡的動力學特性，建立約束方程，并以消耗燃料最少作為性能指標，采用序列二次規劃的方法對多項式的系數進行尋優計算。該方法不僅能求解多個形狀設計參數不確定性問題，而且還能滿足推力加速度的約束限制。仿真驗證了該方法的正確性和可用性，該方法可為任務設計初步階段的軌跡設計和燃料消耗預估提供一定的技術參考。

多項式設計；軌跡設計；序列二次規劃；小推力；推力限制

隨著空間活動的增加和任務復雜性的提升，小推力發動機具有體積小、比沖高、推力大小可控性強等優點，受到了工程技術人員越來越廣泛的重視[1]。由于小推力極大地增強了運動方程的非線性，求解非常困難?，F有的小推力軌跡設計方法中最常見的方法是最優控制[2-3]，按照求解方式的不同一般可以分為間接法[4]和直接法[5]。然而直接法和間接法對初值均具有一定的敏感性，間接法雖然計算精度較高，但是計算精度對初始值有較強的依賴性，因此需要反復迭代來猜測初始協態變量，計算量大而且耗時。因此，尋找新的小推力軌跡設計方法非常有意義。

形狀逼近法是近些年新提出的一種新方法，由于該方法在初始軌跡設計方面效果良好，受到人們越來越廣泛的關注。但是目前的方法大都假設運動軌跡成某一特定的曲線，在一定程度上逼近函數模型的選取極大地限制了逼近能力，而且對于推力約束問題常常捉襟見肘。Petropoulos[6-7]采用指數函數的方法建立了軌跡模型，但是該方法不能滿足速度邊界條件。Wall[8]設計了六階逆多項式(Inverse Polynomial，IP)的軌道交會算法，將半徑表示為關于極角的逆多項式模型，但該方法不能滿足發動機推力約束的限制。隨后，尚海濱[9]等對6階逆多項式算法進行了進一步的深入研究，在前人的基礎上推導了關鍵系數的取值范圍。CHENG[10]在Wall的基礎上進行了改進，建立了三維軌跡運動模型，但是文中仍然沒有考慮推力約束的限制。Abdelkhalik[11]引入了傅里葉級數的逼近方法，但是文中并未給出傅里葉函數各項系數的存在條件。本文主要針對軌道交會任務中的小推力軌跡設計問題，提出了一種基于多項式的軌跡設計方法。

1 多項式設計法

一般情況下，空間任何一條曲線都可以用多項式函數的形式來表示。因此，如果根據軌道運動學方程和約束條件求出多項式的系數，便可以確定航天器的轉移軌跡。

對于共面情況下的軌跡設計，在極坐標系下，航天器的位置僅與軌道半徑r和極角θ有關，確定r和θ有兩種方式，一種是將r表示為θ的函數，通過極角θ的變化來確定r；另一種是將r和θ表示為時間t的函數。本文采用第二種方式，將軌跡模型表示為如下形式：

(1)

式中，ar0，ar1，…，arm和aθ0，aθ1，…，aθn為待確定系數，其中m和n為大于3的正整數，多項式設計法的關鍵就是根據任務參數和約束條件確定多項式的系數。

極坐標系下根據牛頓引力定理，航天器的運動方程為[12]：

(2)其中，r為航天器的位置矢量r的模，即軌道半徑，α為發動機的推力方向角，Ta為發動機的推力加速度大小，μ為引力常數，具體幾何關系如圖 1所示，圖中v表示航天器的速度，γ為飛行路徑角。

圖1 航天器軌跡Fig.1 Spacecraft trajectory

考慮到本文研究的是小推力軌跡設計問題，切向推力可以最快地改變航天器的動能，并且便于發動機姿態機動，在工程上可實現性好，這里假設推力沿著速度線方向，即：

α=γ+nπ

(3)

其中，當推力方向與速度方向相同時n取1，相反時n取0。

根據式(2)中的第二項可得：

(4)

將式(4)代入式(2)中的第一項可得：

(5)

根據切向速度與徑向速度的關系，可計算出推力方向角為：

(6)

將式(6)代入式(5)可得軌跡的等式約束方程為：

(7)

對于發動機推力大小有限制的情況，還需要考慮以下約束

(8)

其中，Ta,max表示發動機允許的最大推力加速度值。

2 多項式系數的確定

前面通過多項式法，將小推力軌跡設計問題轉換成了求多項式的系數問題，將式(1)分別代入式(7)和式(8)，可得僅與多項式系數相關的約束方程為：

f(ar0,ar1，…，arm,aθ0,aθ1，…，aθn,t)=0

(9)

f(ar0,ar1，…，arm,aθ0,aθ1，…，aθn,t)≤1

(10)

為了準確地求解未知系數，至少需要(m+n+2)個約束方程，這就需要對轉移軌跡進行離散化，建立每一個離散點的約束方程。為了保證運動軌跡滿足約束方程(9)和方程(10)，多項式的系數可以分兩部分求取。

2.1 可確定的多項式系數

將初始時刻t0的狀態參數代入式(9)，可得：

(11)

將末端時刻tf的狀態參數代入式(9)，可解出：

(12)

(13)

圖2 極角示意圖Fig.2 Definitions of angles

2.2 待優化系數的求解

2.2.1 優化算法

考慮到本文算法中含有等式和不等式約束方程，采用序列二次規劃(Sequence Quadratic Program, SQP)方法對參數進行優化求解，SQP優化算法是一種非常適合于求解帶有等式和不等式約束問題的有效方法，具體算法流程可參考文獻[13]。為了避免最優解的數值不穩定情況，采用基于Lagrange乘子法，對Lagrange函數取二次近似，將非線性約束優化問題轉化為二次規劃問題來獲得最優解。本文選取的優化目標函數為消耗燃料最小，即：

(14)

2.2.2 求解步驟

對給定的軌道轉移時間tf，多項式的數目m和n，多項式系數的確定步驟如下：

步驟2：將軌道轉移時間0～tf離散成N個區間，利用步驟1中的參考軌跡計算每個節點處的狀態參數(ri,θi)，并通過這些狀態參數解算出多項式系數的初始值；

步驟3：采用序列二次規劃優化算法，根據步驟2中的初始值，求解滿足式(9)和式(10)的一組多項式系數，得出優化后的轉移軌跡。

2.3 多項式系數的初始值猜測

文獻[11]提供了逆多項式方法和指數正弦函數法，但是所建立的函數模型都是將極角θ作為自變量，而本文中的多項式系數法采用的是將時間作為自變量，為了能快速地計算出初始猜測值，這里有兩種可供選用的方法：雙曲正切線函數法和立方多項式法。

雙曲正切函數法的r和θ近似模型為：

(15)

立方多項式法的r和θ的近似模型為：

(16)

其中，a，b，c，d和e，f，g，h為常數，可以通過邊界條件直接確定。

2.4 多項式系數的存在條件

為了避免奇異和算法不收斂的現象產生，任務設計中軌道交會時間需要滿足以下約束：

(17)

式中H0和Hf分別為初始和目標軌道的周期，這里假設H0

(18)

其中，v0和vf分別為離軌點和入軌點的速度大小，a0和af為初始軌道和目標軌道的半長軸。θf越大所需的最大推力加速度值越小。

3 數值仿真

表1 輸入參數和邊界條件Tab.1 Input parameters and boundary conditions

圖3 航天器轉移軌跡Fig.3 Spacecraft trajectories

地心距變化曲線如圖4所示，圖中兩條曲線的兩端點處完全重合，說明兩種方法均能滿足端點處的邊界約束條件，但是從變化趨勢來看，本文方法的半徑變化在中間段要比IP方法更為平緩一些，這也與圖5中的推力加速度曲線變化相對應。由于IP方法在算法中沒有考慮推力加速度值的約束，所以所需的最大推力加速度值超出了允許的約束范圍，而本文方法一直保持在約束范圍內，很好地滿足了設計要求。

圖4 地心距的變化曲線Fig.4 Geocentric distance profile

圖5 加速度變化曲線Fig.5 Thrust acceleration profile

為了更充分地分析本文方法的正確性，采用最優控制(高斯偽譜法)對算例進行了仿真，表2給出了三種不同算法的部分仿真結果，從表中數據可以看出，本文算法不僅保證了推力加速度最大值在約束范圍內，而且燃料更接近于最優控制算法，由于IP算法只是從尋求解析解的角度給出的仿真結果，假設轉移軌跡為6階逆多項式的形式，在一定程度上限制了曲線的逼近能力，并沒對轉移軌跡進行擇優求解，因此燃料消耗較多，相對于IP算法，本文方法的可用性更好。

表2 不同的算法仿真結果對比Tab.2 Comparisons between different types of simulation results

4 結論

針對軌道交會任務中的小推力軌跡設計問題，提出了一種基于多項式的軌跡設計方法。該方法主要具有以下三個優點：①模型簡單，通過多項式的方法有效地將軌跡設計問題簡化成了求解多項式的系數問題；②準確性高，本文方法不僅保證了軌跡滿足邊界條件約束，還考慮了推力加速度的限制；③可用性好，該算法避免了建模時假定運動軌跡成某一特定形狀曲線的限制，增強了形狀曲線的逼近能力，所得結果與最優控制算法非常接近。本文方法可用于任務設計初步階段的軌跡設計和燃料消耗預估，該方法為小推力航天器的任務設計提供了一種新方法。

References)

[1] 湯建國， 張洪波. 小推力軌道機動動力學與控制[M]. 北京:科學出版社, 2013: 1-20.

TANG Jianguo, ZHANG Hongbo. Dynamics and control of low-thrust orbital maneuver[M]. Beijing: Science Press, 2013: 1-20. (in Chinese)

[2] 雍恩米， 陳磊， 唐國金. 飛行器軌跡優化數值方法綜述[J]. 宇航學報, 2008, 29(2): 397-406. YONG Enmi, CHEN Lei, TANG Guojin.A survey of numerical methods for trajectory optimization of spacecraft[J]. Journal of Astronautics, 2008, 29(2): 397-406. (in Chinese)

[3] Hull D G. Conversion of optimal control problems into parameter optimization problems[J]. Journal of Guidance Control and Dynamics, 1997, 20(1): 57-60.

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

[5] Betts J T. Very low-thrust trajectory optimization using a direct SQP method[J]. Journal of Computational & Applied Mathematics, 2000, 120(1/2): 27-40.

[6] Petropoulos A E, Longuski J M. Shape-based algorithm for the automated design of low-thrust, gravity assist trajectories[J]. Journal of Spacecraft & Rockets, 2004, 41: 787-796.

[7] Petropoulos A E. Shape-based approach to automated, low-thrust, gravity-assist trajectory design[D]. USA:Purdue University, 2001.

[8] Wall B J, Conway B A. Shape-based approach to low-thrust rendezvous trajectory design[J]. Journal of Guidance Control & Dynamics, 2009, 32(1): 95-101.

[9] 尚海濱, 崔平遠, 喬棟. 基于形狀的行星際小推力轉移軌道初始設計方法[J]. 宇航學報,2010, 31(6): 1569-1574. SHANG Haibin, CUI Pingyuan, QIAO Dong. A shape-based design approach to interplanetary low-thrust transfer trajectory[J]. Journal of Astronautics, 2010, 31(6): 1569-1574. (in Chinese)

[10] Cheng X, Shi X N, Cui N G. Modified shape-based method for three-dimensional trajectory design[J]. Journal of Harbin Institute of Technology, 2012, 19(2): 1-4.

[11] Abdelkhalik O, Taheri E. Shape-based approximation of constrained low-thrust space trajectories using fourier series[J]. Journal of Spacecraft & Rockets, 2012, 49(3): 535-545.

[12] Miele A, Wang T, Williams P N. Computation of optimal Mars trajectories via combined chemical/electrical propulsion, part 1: baseline solutions for deep interplanetary space[J]. Acta Astronautica, 2004, 55(2): 95-107.

[13] 尚海兵， 崔平遠， 欒恩杰. 地球-火星的燃料最省小推力轉移軌道的設計與優化[J]. 宇航學報, 2006, 27(6): 1168-1171. SHANG Haibin, CUI Pingyuan, LUAN Enjie. Design and optimization of Earth-Mars optimal-fuel low-thrust trajectory[J]. Journal of Astronautics, 2006, 27(6): 1168-1171. (in Chinese)

Low-thrust trajectory design of rendezvous based on polynomial function

ZENG Kui, GENG Yunhai, XIE Chengqing

(Research Center of Satellite Technology, Harbin Institute of Technology, Harbin 150001, China)

Shape-based approximation is an effective method for the low-thrust trajectory design. However, the vast majority of methods assume that the motion trajectory is a specific shape, and there is no constraint on the thrust acceleration. In this representation, according to the issue of low-thrust spacecraft rendezvous, a new method for transfering trajectory design was proposed. Based on the polar coordinates, the trajectory design was successfully turned into solving polynomial coefficients problem through three degrees of freedom motion model built by introduced polynomial function. Meanwhile, the dynamic characteristics of the trajectory as well as the constraint equations were deduced. Subsequently, the optimal polynomial coefficients were solved by the approach of sequential quadratic programming. This method can not only solve problems with a greater number of free parameters, but also meet the thrust acceleration constraints. The simulation has confirmed that this method is accuracy and availability. It can provide a certain technical reference for the trajectory design and fuel consumption estimation during the preliminary stage.

polynomial design; trajectory design; sequential quadratic programming; low-thrust; thrust constrained

10.11887/j.cn.201606013

2015-06-24

國家自然科學基金資助項目(61473096)

曾奎(1987—)，男，湖北襄陽人，博士研究生，E-mail：zenghit@126.com； 耿云海(通信作者)，男，教授，博士，博士生導師，E-mail：gengyh@hit.edu.cn

V412.4+1

A

1001-2486(2016)06-077-05

http://journal.nudt.edu.cn