基于傅里叶级数的小推力航天器快速轨迹设计

曾奎 耿云海 陈炳龙

基于傅里叶级数的小推力航天器快速轨迹设计

曾奎1耿云海1陈炳龙2

为了满足小推力航天器交会轨迹的快速性设计需求,基于形状逼近理论,设计了一种三维轨迹模型.将轨迹设计问题转换为求解傅里叶级数的系数问题,避免了轨迹运动方程非线性强、难以求解的难题,极大地提高了计算效率.考虑到推力加速度的限制,建立了加速度约束方程,并结合轨迹的运动方程,给出了傅里叶级数的求解过程.同时根据边界条件和最大推力加速度值,定性地分析了傅里叶系数的存在条件.仿真验证了该方法的正确性和可行性,并从计算效率上与高斯伪谱法进行了对比,结果表明本文的方法计算耗时仅为高斯伪谱法的0.67%.

快速性,小推力,轨迹设计,推力限制,傅里叶级数,航天器

航天器任务设计的初步阶段,需要快速地设计出一条合理的初始轨迹并评估燃料的消耗量,然而由于小推力轨迹模型的非线性较强,直接以解析解的形式求出运动轨迹非常困难,而且实际工程中推力具有最大限制,极大地增加了任务挑战性[1].传统的方法中一般将小推力轨迹设计问题转化为最优控制问题来求解,求解方法主要分为两类[2−4]:直接法和间接法.但是这两类方法对初始猜测值比较敏感,而且对于不同的初始值需要对全过程重新进行数值积分,不能满足任务初步阶段的快速性需求.因此,寻找新的快速性轨迹设计方法非常重要.

为了解决任务设计初步阶段的快速性需求, 1999年Petropoulos等[5]提出了一种形状逼近法,为小推力轨迹设计问题打开了一扇新的大门.形状逼近法是一种基于逆向设计的思想,它首先假设运动轨迹呈某一特定的曲线,然后用形状曲线拟合的方法设计出满足要求的转移轨迹.由于该方法简单、快速,迅速引起了学者的广泛关注.Petropoulos等[6]用正弦指数函数模型设计了小推力拦截轨迹.Izzo[7]研究了基于正弦指数函数模型的多圈Lambert问题.然而,对于交会问题正弦指数函数模型不能满足终端速度约束条件.随后,Wall等[8−9]提出了一种6阶逆多项式的方法,该方法虽然能满足边界条件约束,但是对于交会问题,只能通过增大转移时间或增加轨迹的圈数来减小最大推力加速度值.为了解决推力约束的问题,Taheri等[10−11]引入了傅里叶函数,设计了一种近似轨迹模型,但是该模型适合于转移圈数给定的情况,而且文中没有给出系数的存在条件,对于快速性轨迹设计需求存在一定的局限性.

针对空间小推力轨道交会问题,本文在前人的基础上引入傅里叶级数模型,设计了一种基于傅里叶级数的空间小推力轨迹设计方法.通过轨迹模型,将小推力轨迹设计问题转化成了求解傅里叶系数的问题.根据加速度约束方程和轨迹运动约束方程,结合序列二次规划算法给出了系数的求解过程,并定性分析了傅里叶系数解的存在条件.最后仿真验证了该方法的正确性和可行性.本文所设计的方法可为工程设计初步阶段提供一定的技术参考,或为更进一步的精确计算提供初值依据.

1 基于傅里叶级数的轨迹模型

1.1 问题的提出

在轨道交会过程中,航天器在空间的转移轨迹是一条满足约束条件的空间曲线.通常情况下,对于任意给定的一条空间曲线都可以近似用一个傅里叶级数的形式来表示,所以航天器的轨迹设计问题可以转化为傅里叶级数的设计问题,而傅里叶级数取决于多项式的系数,因此,只要确定了多项式的系数,就可以确定航天器在空间的运动轨迹.

为了便于分析和计算,本文以初始轨道平面为参考面,建立柱坐标系,将航天器的轨迹表示为极角θ、半径r和轨道面法向位置z的函数,图中x轴和y轴位于初始轨道面内,与z轴构成右手坐标系,如图1所示.

图1 柱坐标系Fig.1 Cylindrical coordinate system

根据牛顿引力定律,航天器的运动方程在柱坐标系下可以表示为[9,12]

式中,z表示oz轴方向的位置,Tain表示航天器在oxy平面内的推力加速度大小,Taz表示oz轴方向的推力加速度大小,α为oxy平面内的推力方向角,µ为地心引力常数,s表示地心距,其中

在柱坐标系下,只要确定了航天器每个时刻的状态(r,θ,z),航天器的空间运动轨迹便可以唯一确定.下面采用傅里叶级数的设计思想,分别对oxy平面内的运动和z向的运动进行建模分析.

1.2 模型的建立

1.2.1 oxy平面内的运动

在平面内航天器的位置是轨道半径r和极角θ的函数,r和θ可以表示成两种函数形式[6,8],一种是以θ为自变量,用θ来表示函数r;另一种是将r和θ都表示成关于时间t的函数.这里选取第二种形式,航天器在oxy平面内的运动轨迹模型如下:

式中,T表示从离轨点到入轨点之间的飞行时间, a0,an,bn和c0,cn,dn为待定系数,nr,nθ∈N且nr≥2,nθ≥2.

假设平面内推力方向角等于飞行路径角,即

式中,γ表示飞行路径角,当推力方向与速度方向相同时b等于1,相反时b等于0.

对式(1)中的第二项进行整理可得

将上式代入式(1)中的第一项,整理可得

将式(3)代入式(5)可得

根据切向速度与径向速度之间的关系,可知推力方向角的表达式为

将式(7)代入式(6),可得oxy平面内轨迹约束方程为

1.2.2 平面外的运动

平面外的运动即oz轴方向的运动,考虑到转移轨迹具有一定的周期性,为了尽量减少优化参数,同时提高轨迹模型的逼近能力,这里选取极角θ作为自变量,将平面外的运动表示成关于余弦函数和高阶多项式的复合函数的形式

式中,az,bz,cz,dz为多项式的系数,q为大于等于3的正整数.由于平面外的运动边界条件已知,即初始和末端时刻航天器z向运动的位置和速度参数是确定的,所以可以直接求出式(9)中的多项式系数.

将式(9)代入式(1)中的第三项得z轴方向推力加速度值为

由式(4)和式(10)可得总的推力加速度大小为

由于实际工程中,发动机的推力是有限制的,所以在轨迹设计中还需要考虑如下约束方程

式中,Ta,max表示轨道机动过程中允许的最大推力加速度值.

2 傅里叶级数系数的求解

由于傅里叶级数各项系数的好坏直接影响着轨迹模型的逼近能力,因此求解满足边界条件和约束方程的最优系数是本文设计方法的关键.

将式(2)和(9)代入式(8)和(12),可以得到一组仅与未知系数和时间变量相关的非线性代数方程

式中,j=0,···,m,k=0,···,n.式(13)和(14)即为各个时刻航天器的运动轨迹需要满足的所有约束函数.求解未知系数至少需要2(nr+nθ+1)个方程.为了得出约束方程,需要先对轨道转移时间进行离散化,然后针对每一个离散点建立满足约束函数的方程.

通常情况下,nr和nθ取10以内的值即可达到较高的精度.对于带推力约束的轨迹设计问题,nr和nθ的取值大于等于3才能保证推力在约束范围以内,为了保证边界点处的位置和速度精度,未知系数分两部分求取.

2.1 由边界条件直接确定的系数

令初始时刻t0=0,将初始时刻t0和末端时刻tf的状态参数分别代入式(13),可得式(2)中函数r的可确定系数为

式中,λ1≥3,λ1为奇数,λ2≥4,λ2为偶数.r0和 rf分别为初始和末端时刻轨道半径,和分别为初始和末端时刻径向速度.

同理,将初始时刻t0和末端时刻tf的状态参数分别代入式(13),可得式(2)中函数θ的可确定系数为

式中,τ1≥3,τ1为奇数,τ2≥4,τ2为偶数.θ0和 θf分别为初始和末端时刻的极角,0和分别为初始和末端时刻极角的变化率.

同理,将t0和tf时刻的状态参数分别代入式(2)和式(9),整理得

z0和zf分别为初始时刻oz轴方向的位置,和分别为末端时刻oz轴方向的速度值.

图2 极角的定义Fig.2 Definitions of angles

2.2 待优化系数的求解

2.2.1 优化算法和求解步骤

为了使计算过程尽量简单,本文选用序列二次规划(Sequential quadratic programming,SQP)[10]的方法对待求系数进行优化求解.在计算的过程中为了避免有多个极值点和陷入局部极小值的情况,采用了Lagrange乘子法,对Lagrange函数取二次逼近,将带有非线性约束的优化问题转换为二次规划问题,通过求解指标函数的最小值来获得待求系数的最优值.本文以燃料消耗作为优化指标函数

对于给定的轨道转移时间T和极角θf、傅里叶级数nr和nθ,以及离散点DP的个数,待优化多项式系数的计算步骤如下:

2)结合离散点DP的个数,将轨道交会时间区间0～tf,离散成(DP−1)个区间,利用步骤1)中给出参考轨迹模型,计算出每个节点i处的状态参数将这些状态参数代入式(13)中,通过等式约束方程初步估算出多项式系数的猜测值;

3)在猜测值的基础上,根据每个节点处的约束方程,采用序列二次规划的方法,优化求解出满足所有约束方程(13)和(14)的一组最优解,所得的最优解即为所求的待优化多项式系数.

2.2.2 待优化系数的初值猜测

在用SQP方法求解待优化系数时,首先要给出一组初始猜测值.为了避免自主猜测的盲目性,本文采用初始值自主给定的方式,首先假设存在一条满足等式(13)约束的参考轨迹,然后根据参考轨迹上的离散点,拟合出初始猜测值.文献[8,10]给出了逆多项式法和指数正弦函数两种参考轨迹模型,但是文献中的模型都是以极角θ作为自变量,而本文的模型是以时间为自变量,为了能够简单、快速地得出初始参考轨迹,这里选用了立方多项式的方法.

r和θ的参考轨迹模型如下

式中,a,b,c,d和e,f,g,h均为常值,其值可由初始和末端边界状态参数直接得出.

2.3 待优化系数的存在条件

待优化系数的存在性,决定着是否存在满足设计要求的轨迹曲线,在有推力约束轨迹的设计任务中,轨迹设计的难点是寻找满足推力加速限制的曲线.对于本文的方法,影响轨迹曲线性能的因素,除了轨迹模型的逼近能力外,主要的影响因素就是转移轨迹的圈数Nrev.本节将定性给出Nrev的确定方法.

对于给定的轨道交会任务,交会时间是确定的,轨道交会时间T满足以下约束

式中Pt0和Ptf分别表示初始和目标轨道的周期(这里假设Pt0＜Ptf).

为了保证算法的收敛、避免轨迹发生奇异现象,由式(20)可得转移轨迹的圈数Nrev取值范围为

此外,为了保证在推力约束范围以内存在满足边界条件的最优解,除了θf≥以外,还应满足以下条件

式中v0和vf表示离轨点和入轨点的速度大小,a0和af表示初始轨道和目标轨道的半长轴.

对于不同的轨道交会任务,根据给定的设计参数Nrev可能存在多个可行值,由于Nrev的值越大,所需离散点的个数就越多,直接会导致计算量增大,为了满足快速性设计需求,在计算过程中选取Nrev的最小值作为最优值.

3 数值分析

为了论证本文所设计的方法的可行性和有效性,选取如下算例进行仿真验证.仿真电脑为Windows 7系统,内存(RAM)4GB,仿真平台为Matlab2014a.

假设航天器初始轨道参数和目标轨道参数如表1所示,表中的a、e、i、ω、Ω和f分别为轨道半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经和真近点角.输入参数nr=4,nθ=5, q取9,离散点数DP 取22,整个轨道转移时间为17449s,初始质量为4000kg,发动机的比冲为 3000m/s.根据式(21)和(23)计算可得1.9306＜Nrev＜2.8830,由于Nrev为正整数,所以Nrev取2.为了便于计算,仿真中将参数进行无量纲化,选取参考长度为DU,1DU=6378.1km,参考时间为发动机的推力限制设为Ta,max=0.014DU/TU2.

表1 输入状态参数Table 1 Input boundary parameters

由前文可知,本文方法中设计参数有q、离散点个数DP、nr和nθ,为了分析各设计参数对计算效率和精度的影响,这里结合仿真算例,在给定的设计参数基础上,分别以某一设计参数作为变量进行仿真.从仿真结果来看,初始和末端的速度、位置计算精度均在10−16量级,由于篇幅有限,表2～4给出了部分仿真结果,其中∆v表示燃料消耗,Tmax表示最大推力加速度值,∆t表示仿真时间.

表2 q取不同值时的仿真结果Table 2 Results with different q

表3 DP取不同值时的仿真结果Table 3 Results with different DP

表4 nr,nθ取不同值时的仿真结果Table 4 Results with different nrand nθ

从表2～4中的数据分析可知,单独改变q和DP的值对最大推力加速度值Tmax的影响不是很明显.当q较小或过大时,轨迹模型的逼近能力较差,计算耗时长,而且精度下降.q直接影响的是轨迹模型的逼近能力,根据多项式的特性可知,q取7附近的值时,式(9)的逼近能力最强,所以q一般取7附近的值.对于转移轨迹,当每圈离散点的个数在10以上时,可以达到较高的计算精度,对于本文的算例当每圈离散点的个数超过60时,计算精度基本保持不变,但是由于计算量增加,计算效率明显下降.nr和nθ的取值对Tmax影响比较明显,随着nr、nθ的增大,Tmax和∆v逐渐下降,并趋于稳定,当nr、nθ增大到一定值时,计算精度趋于饱和,由于模型中多项式的系数增加,计算时间快速上升.同时,由表2和表4对比分析可知,在轨迹设计过程中选定q,将nr和nθ作为优化变量,更容易快速得出精确的计算结果.

令柱坐标系下初始轨道面内的x轴、y轴与地心惯性坐标系在赤道面内的两个坐标轴重合,算例的仿真结果如图3～5所示.考虑到逆多项式设计法(Inverse polynomial,IP method)是形状逼近法中最经典的方法,为了直观地分析本文方法(Proposed method)的合理性,仿真中同时给出了逆多项式设计法的仿真结果.

图3为轨道交会过程中航天器的轨迹变化,从图中可以看出采用本文的设计方法,在离轨第一圈时间内,轨迹变化较平缓,基本沿着初始轨道附近运动,随着时间的推进,z方向的运动逐渐变快,随后平滑地切入到目标轨道交会点,进入目标轨道.仿真得到的离轨点状态值为(1.1254, 1.2034E–18,2.0311E–18,–3.0251E–19,0.8376, 4.3610E–18),入轨点的状态值为(1.4842,3.1415,–0.0518,–1.5261E–18,0.5501, 1.3892E–19),与初始给定的状态值误差均小于10−17量级,表明所得轨迹曲线能很好地满足边界约束条件.

为了更直观地分析轨迹变化情况,图4给出了轨道交会过程中航天器地心距的变化曲线,从曲线的变化趋势可以看出在轨迹末端实线和虚线完全重合,两种方式都能保证航天器精确地到达交会点,但是相比于逆多项式设计法,本文方法在离轨和入轨区间段轨迹变化更为平缓.

图4 地心距的变化Fig.4 Geocentric distance profiles

图5给出了推力加速的变化曲线,由于本文方法在算法中考虑了推力加速度的限制,推力加速度的最大值保持在给定的范围以内.而逆多项式方法假定运动轨迹为某一特定的曲线,限制了模型的逼近能力,最大推力加速度值超出了允许的范围.

图5 推力加速度变化曲线Fig.5 Thrust acceleration profiles

考虑到高斯伪谱法是最优控制方法中精度非常高的一种算法,最后采用高斯伪谱法对算例进行了仿真.由表5中的数据可知本文的方法燃料消耗非常接近高斯伪谱法,说明本文方法燃料预估精度较高,而且明显高于逆多项式方法.同时本文方法计算得出的Tmax值小于0.014,进一步表明本文的方法具有满足推力加速度约束的能力.从仿真时间来看,由于逆多项式是以解析解的形式给出的结果,计算时间最短,但是相比于逆多项式算法仅多了0.2789s;相比于高斯伪谱法,本文方法避免了反复迭代寻找初始协态变量的过程,计算时间减少了61.488s,仅为高斯伪谱法的0.67%,极大提高了计算效率.

表5 三种方法计算结果比较Table 5 Comparison of the results using three methods

4 结论

本文针对小推力航天器空间交会轨迹设计问题,提出了一种基于傅里叶级数的快速轨迹设计方法,并通过仿真验证了该方法的正确性和可行性.本文的方法主要有以下3个优点:1)模型简单、计算快速:通过形状曲线拟合运动轨迹,克服了非线性运动方程难以求解的困难,直接将轨迹设计问题简化成了求解傅里叶级数系数的问题,提高了计算效率;2)求解精度较高:避免了建模时假设运动轨迹成某一特定形状曲线的限制,增强了形状曲线的逼近能力,改善了求解精度;3)可用性好:设计过程中考虑了发动机的推力限制,所得结果不仅能满足边界条件,而且还能满足推力加速度约束条件.本文方法可为小推力航天器任务设计初步阶段的交会轨迹设计和燃料预估提供一定的参考,或为更进一步的精确计算提供可靠的初始值.

本文方法旨在为任务设计初步阶段小推力轨迹设计问题提供一种新思路,对于如何定量地选择设计参数(q、nr、nθ和DP)、考虑摄动干扰和常值小推力等轨迹设计问题,将在后续的工作中做进一步的深入研究.

1 Laipert F E,Longuski J M.Low-thrust trajectories for human missions to Ceres.Acta Astronautica,2014,95: 124−132

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

3 Sun C,Yuan J P,Fang Q.Continuous low thrust trajectory optimization for preliminary design.Proceedings of the Institution of Mechanical Engineers.Part G:Journal of Aerospace Engineering,2016,230(5):921−933

4 Topputo F,Zhang C.Survey of direct transcription for lowthrust space trajectory optimization with applications.Abstract and Applied Analysis,2014,2014:851720,DOI: 10.1155/2014/851720

5 Petropoulos A E,Longuski J M,Vinh N X.Shape-based analytic representations of low-thrust trajectories for gravityassist applications.In:Proceedings of the 2000 Advances in Astronautical Sciences.Girdwood,USA:AAS,2000. 563−581

6 Petropoulos A E,Longuski J M.Shape-based algorithm for the automated design of low-thrust,gravity assist trajectories.Journal of Spacecraft and Rockets,2004,41(5): 787−796

7 Izzo D.Lambert's problem for exponential sinusoids.Journal of Guidance,Control,and Dynamics,2006,29(5): 1242−1245

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

9 Wall B J.Shape-based approximation method for lowthrust trajectory optimization.In:Proceedings of the 2008 AIAA/AAS Astrodynamics Specialist Conference and Exhibit.Honolulu,Hawaii:AIAA,2008.

10 Taheri E,Abdelkhalik O.Shape based approximation of constrained low-thrust space trajectories using Fourier series.Journal of Spacecraft and Rockets,2012,49(3): 535−546

11 Taheri E,Abdelkhalik O.Initial three-dimensional lowthrust trajectory design.Advances in Space Research,2016, 57(3):889−903

12 Gondelach D J,Noomen R.Hodographic-shaping method for low-thrust interplanetary trajectory design.Journal of Spacecraft and Rockets,2015,52(3):728−738

曾 奎 哈尔滨工业大学卫星技术研究所博士研究生.主要研究方向为小推力航天器轨迹设计.

E-mail:zenghit@126.com

(ZENG Kui Ph.D.candidate at the Research Center of Satellite Technology,Harbin Institute of Technology. His research interest covers low-thrust spacecraft trajectory design.)

耿云海 哈尔滨工业大学航天学院教授.主要研究方向为飞行器动力学与控制.本文通信作者.

E-mail:gengyh@hit.edu.cn

(GENG Yun-HaiProfessor at the Astronautics School,Harbin Institute of Technology. His research interest covers astrodynamics and control.Corresponding author of this paper.)

陈炳龙 博士,上海微小卫星工程中心工程师.主要研究方向为航天器动力学与控制.

E-mail:chenbinglonghit@163.com

(CHEN Bing-Long Ph.D.,engineer at Shanghai Engineering Center for Micro-Satellite.His research interest covers spacecraft dynamic and control.)

Rapid Design of Low-thrust Rendezvous Trajectory with Fourier Series

ZENG Kui1GENG Yun-Hai1CHEN Bing-Long2

Trajectory design with low-thrust propulsion needs a method for quickly approximating the spacecraft′s trajectory during the preliminary stage phase.To address this requirement,a 3D rendezvous trajectory design model is proposed based on the shape-based theory,by which the trajectory design problem is simplified into solving the coefficients of Fourier series.In the process of solving coefficients,thrust constraint is considered,and the existence conditions of the coefficients are qualitatively analyzed.Finally,the correctness and feasibility of this approach are verified through a numerical example.Results show that the proposed method can not only satisfy the boundary conditions and thrust constraint,but also has a good computational efficiency.Its computation time is only 0.67%of that of the Gauss pseudospectral method.

Quickly approximating,low-thrust,trajectory design,thrust constraint,Fourier serious,spacecraft

曾奎,耿云海,陈炳龙.基于傅里叶级数的小推力航天器快速轨迹设计.自动化学报,2016,42(11):1641−1647

Zeng Kui,Geng Yun-Hai,Chen Bing-Long.Rapid design of low-thrust rendezvous trajectory with Fourier series.Acta Automatica Sinica,2016,42(11):1641−1647

2015-12-21 录用日期2016-04-28

Manuscript received December 21,2015;accepted April 28, 2016

国家自然科学基金(61473096)资助

Supported by National Natural Science Foundation of China (61473096)

本文责任编委崔平远

Recommended by Associate Editor CUI Ping-Yuan

1.哈尔滨工业大学卫星技术研究所 哈尔滨150001 2.上海微小卫星工程中心上海201210

1.Research Centre of Satellite Technology,Harbin Institute of Technology,Harbin 150001 2.Shanghai Engineering Center for Micro-Satellite,Shanghai 201210

DOI 10.16383/j.aas.2016.c150859