曾 奎，耿云海，謝成清

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

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

曾 奎，耿云海，謝成清

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

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

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

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

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

1 多項式設計法

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

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

(1)

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

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

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

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

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

α=γ+nπ

(3)

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

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

(4)

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

(5)

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

(6)

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

(7)

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

(8)

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

2 多項式系數(shù)的確定

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

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)

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

2.1 可確定的多項式系數(shù)

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

(11)

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

(12)

(13)

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

2.2 待優(yōu)化系數(shù)的求解

2.2.1 優(yōu)化算法

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

(14)

2.2.2 求解步驟

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

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

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

2.3 多項式系數(shù)的初始值猜測

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

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

(15)

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

(16)

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

2.4 多項式系數(shù)的存在條件

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

(17)

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

(18)

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

3 數(shù)值仿真

表1 輸入?yún)?shù)和邊界條件Tab.1 Input parameters and boundary conditions

圖3 航天器轉(zhuǎn)移軌跡Fig.3 Spacecraft trajectories

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

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

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

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

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

4 結(jié)論

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

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] 雍恩米， 陳磊， 唐國金. 飛行器軌跡優(yōu)化數(shù)值方法綜述[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] 尚海濱, 崔平遠, 喬棟. 基于形狀的行星際小推力轉(zhuǎn)移軌道初始設計方法[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] 尚海兵， 崔平遠， 欒恩杰. 地球-火星的燃料最省小推力轉(zhuǎn)移軌道的設計與優(yōu)化[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