«上一篇
文章快速检索    
下一篇»
  深空探测学报  2016, Vol. 3 Issue (1): 51-55  DOI: 10.15982/j.issn.2095-7777.2016.01.008
0

引用本文 [复制中英文]

袁旭, 朱圣英. 基于伪谱法的小天体最优下降轨迹优化方法[J]. 深空探测学报, 2016, 3(1): 51-55. DOI: 10.15982/j.issn.2095-7777.2016.01.008.
[复制中文]
Yuan X, Zhu S Y. Small Body Descent Trajectory Optimization Based on Pseudospectral Method[J]. Journal of Deep Space Exploration, 2016, 3(1): 51-55. DOI: 10.15982/j.issn.2095-7777.2016.01.008.
[复制英文]

基金项目

国家重点基础研究发展计划“973”项目(2012CB720000);国家自然科学基金资助项目(61374216,61304226,61304248);教育部博士点基金(20121101120006)

文章历史

收稿日期:2014-11-29
修回日期:2015-10-09
基于433 Eros的多面体引力模型精度与运行时间研究
袁旭1, 2, 3, 朱圣英1, 2, 3    
1. 北京理工大学 深空探测技术研究所, 北京 100081;
2. 深空自主导航与控制工信部重点实验室, 北京 100081;
3. 飞行器动力学与控制教育部重点实验室, 北京 100081
摘要: 针对小天体着陆探测下降阶段的多约束轨迹优化问题,基于Gauss伪谱法进行了燃耗最优下降轨迹优化设计,得出了燃耗最优下降轨迹。建立了小天体下降轨迹优化问题的最优控制问题模型,采用Gauss伪谱法进行离散化,转化为非线性规划问题进行了求解。数学仿真结果显示:优化结果符合各项约束条件,以零速度到达了目标着陆点,且符合燃耗最优的优化目标。利用Gauss伪谱法进行小天体最优下降轨迹优化,计算速度快,求解精度高。
关键词: 小天体    下降    着陆    轨迹优化    Gauss伪谱法    
Small Body Descent Trajectory Optimization Based on Pseudospectral Method
YUAN Xu1, 2, 3, ZHU Shengying1, 2, 3    
1. School of Aerospace Engineering, Beijing Institute of Technology, Beijing 100081, China;
2. Key Laboratory of Autonomous Navigation and Control for Deep Space Exploration, Ministry of Industry and Information Technology, Beijing 100081, China;
3. Key Laboratory of Dynamics and Control of Flight Vehicle, Ministry of Education, Beijing 100081, China
Abstract: Aimed at the multiple-constraint trajectory optimization problem in the descent phase of small body landing exploration, Gauss pseudospectral method is used for the optimization design and the fuel-optimal descent trajectory is derived. The optimal control model of small body descent trajectory optimization problem is established and discretized using Gauss pseudospectral method. The problem is then solved by transforming into a nonlinear programming problem. Mathematical simulation results show that all the constraints are satisfied, the optimization index of fuel consumption is minimized and the spacecraft reaches the target landing site with zero velocity. Using Gauss Pseudospectral method, the computation process of small body descent trajectory optimization is fast and has high solution accuracy.
Key words: small body    descent    landing    trajectory optimization    Gauss pseudospectral method    
0 引言

小天体探测是人们认识和研究太阳系的起源与演化的重要手段,是21世纪深空探测活动的重要内容。随着小天体探测活动的不断发展,探测手段从飞跃探测等简单形式逐渐向撞击、着陆、采样返回等更为复杂、科学成果更加丰富的形式转变[1]。迄今,人类2次成功完成了着陆任务,一次为2001年2月,NASA发射的NEAR探测器着陆于Eros小行星[2],另一次为2005年11月,日本JAXA发射的Hayabusa探测器着陆于Itokawa小行星并采样返回[3]。此外,欧洲空间局(ESA)发射的Rosetta探测器于2014年11月12日对Churyumov-Gerasimenko彗星进行了着陆探测。

下降与着陆阶段是小天体着陆或采样返回探测的关键阶段,对能否安全、准确到达预设的具有科学探测价值的目标区域起着决定性的作用。小天体探测器的下降过程可以转化为一个轨迹优化问题以及对标称轨迹的跟踪控制问题。设定的标称轨迹需要能够安全、准确地到达指定着陆点,满足始末状态约束、路径约束、控制约束等多重约束,同时使某项重要的性能指标(如燃耗、时间等)最优化。

轨迹优化方法主要包括基于参数化方法的直接法和基于庞特里亚金极小值原理的间接法。Gauss伪谱法(Gauss pseudospectral method,GPM)是直接法的一种,是近年来应用广泛的轨迹优化方法之一。Gauss伪谱法的原理是利用全局插值多项式的有限基来近似状态量和控制量,且在一系列的配点上满足动力学方程约束,将微分方程约束转化为代数约束,从而将原连续最优控制问题转化为离散的非线性规划问题。Gauss伪谱法对初始猜测值不敏感,配点后的状态和控制量采用多项式拟和,无需数值积分,大大节省了计算时间,且能满足一阶必要性条件,因此成为轨迹优化领域的研究热点[4]。本文采用Gauss伪谱法对小天体最优下降轨迹问题进行优化求解。

本文首先建立小天体最优下降轨迹优化问题的数学模型,给出小天体下降过程的动力学方程,进而建立原始最优控制问题的各项约束和目标性能函数;然后利用Gauss伪谱法对连续最优控制问题进行离散化,转化为静态的非线性规划问题;最后以Eros 433小行星为目标小天体,对上述过程进行仿真求解,对仿真结果进行了分析和总结。

1 小天体最优下降轨迹优化问题 1.1 问题描述

小天体最优下降轨迹优化问题可描述为探测器从初始的探测轨道上某处,经一段时间后以零速度到达目标着陆点,实现软着陆。轨迹优化的目的是在满足始末状态、动力学和控制约束的前提下,使在此过程中的某项指标达到最优,如时间最优或燃耗最优。由于小天体探测的目标星体距地球非常遥远,任务持续时间可达数年,尤其采样返回任务还需保留足够的燃料以将小天体表面样本带回地球,因此燃耗是任务设计中的关键环节。这里将燃耗最优作为下降轨迹的优化目标。

1.2 动力学模型

在进行下降轨迹优化前,首先建立原始最优控制问题的数学模型,包括探测器在小天体下降过程中的动力学模型、受到的各项约束和需要优化的指标函数。

探测器的下降过程采用小天体固连坐标系∑a进行建模:坐标原点oa位于小天体的质量中心,za轴沿小天体最大惯量轴即小天体自旋轴方向,xa轴沿小天体最小惯量轴方向,ya轴的定义使xayaza三轴满足右手法则。

假设小天体密度均匀并绕其最大惯量主轴匀速转动,忽略其他干扰力,在小天体固连坐标系下,着陆器的动力学方程可表示为

$\left\{ {\begin{array}{*{20}{l}} {\ddot x = {\omega ^2}x + 2\omega \dot y + {V_x} + \frac{{{T_x}}}{m}} \\ {\ddot y = {\omega ^2}y - 2\omega \dot x + {V_y} + \frac{{{T_y}}}{m}} \\ {\ddot z = {V_z} + \frac{{{T_z}}}{m}} \\ {\dot m = - \frac{T}{{{I_{sp}}{ \cdot _0}}}} \end{array}} \right.$ (1)
其中:x,y,z为探测器的三轴位置坐标;$\dot x,\dot y,\dot z$为探测器三轴速度大小;m为着陆器质量,共同组成探测器状态向量X = {\text{ }}{\left[{\begin{array}{*{20}{c}} x&y&z&{\dot x}&{\dot y}&{\dot z}&m \end{array}} \right]^{\text{T}}}$;ω为小天体的自旋角速度;Tx,Tv,Tz分别为三轴的控制力大小;Isp为发动机比冲;g0为地球海平面引力加速度;Vx,Vv,Vz为小天体作用在探测器上的三轴引力加速度大小。

小天体引力势函数采用球谐函数形式,即[5]

$\left\{ \begin{array}{l} V = \frac{{GM}}{a}\mathop \sum \limits_{n = 0}^\infty \mathop \sum \limits_{m = 0}^n {\left( {\frac{a}{r}} \right)^{n + 1}}({C_{nm}}{V_{nm}} + {S_{nm}}{W_{nm}})\\ {V_{nm}} = {P_{nm}}(\sin \varphi )\cos m\theta \\ {W_{nm}} = {P_{nm}}(\sin \varphi )\sin m\theta \end{array} \right.$ (2)
其中:a为小天体名义半径; θ,φ,r分别为探测器所处的经度、纬度和半径;Pnm为缔结勒让德多项式;CnmSnm为小天体对应的引力场球谐函数各阶系数。

1.3 约束与性能指标函数

探测器在下降过程中受到的约束包括动力学约束,初始、末端状态约束和路径约束。其中,动力学约束即为下降过程的动力学方程式(1)。

初始状态约束是指探测器在研究的初始时刻,位置、速度等状态变量满足一定的约束关系,如处于某一绕飞轨道,或处于某一悬停状态等,可表示为

$\psi (x,y,z,\dot x,\dot y,\dot z,m) = 0$ (3)

当初始状态为已知时即为

$X(0) = {\left[ {{x_0},{y_0},{z_0},{{\dot x}_0},{{\dot y}_0},{{\dot z}_0},{m_0}} \right]^{\text{T}}}$ (4)

末端状态约束包括到达预定的着陆位置、着陆速度为零,以及末端质量不小于探测器的干重

$\left\{ \begin{gathered} x({t_f}) = {x_f},y({t_f}) = {y_f},z({t_f}) = {z_f},\hfill \\ \dot x({t_f}) = 0,\dot y({t_f}) = 0,\dot z({t_f}) = 0,\hfill \\ m({t_f}) \geqslant {m_{dry}} \hfill \\ \end{gathered} \right.$ (5)
其中:tf为到达目标着陆点的时刻;${[{x_f},{y_f},{z_f}]^{\rm{T}}}$为目标着陆点;mdrv为探测器干重,即全部燃料耗尽后的质量。

此外,这里考虑探测器推力的限制为路径约束

$T \leqslant {T_{\max }}$ (6)

探测器下降轨迹优化的目标是燃料消耗最小,即末端质量最大。性能指标函数可表示为

$J = - m({t_f})$ (7)

动力学约束式(1)、边界约束式(4)、式(5)、路径约束式(6)和性能指标函数式(7)即组成了小天体最优下降轨迹优化问题的数学模型。这样,最优轨迹优化问题转换为最优控制问题,下面利用Gauss伪谱法对此连续最优控制问题进行求解。

2 基于伪谱法的最优下降轨迹优化

Gauss伪谱法首先求解Legendre多项式一阶导数在[-1,1]区间内的根作为时间节点值,将节点上的状态和控制变量作为参数,用Lagrange插值多项式拟合轨迹上各时刻状态和控制变量。对于动力学方程的微分方程约束,对Lagrange插值多项式在各节点上求导,令其等于动力学方程在各节点上的值,这样将原始的动力学微分方程约束转化为代数方程约束,再结合边界约束、路径约束以及性能指标函数共同构成有静态的限维非线性规划问题进行求解[6]

利用Gauss伪谱法解最优控制问题,首先需要将$t \in \left[{{t_0},{t_f}} \right]$映射到$\tau \in \left[{ - 1,1} \right]$上

$\tau = - 1 + \frac{{2(t - {t_0})}}{{{t_f} - {t_0}}}$ (8)

这里令t0 = 0 ,则式(8)化简为

$\tau = \frac{{2t}}{{{t_f}}} - 1$ (9)

Gauss伪谱法利用Lagrange插值多项式对状态变量和控制变量进行近似。在区间(–1,1)内设置N个Legendre-Gauss(LG)点

${\tau _i} = Roots({\dot P_{N - 1}}(\tau )) ({\rm i} = 1,\cdots,N) $ (10)
其中:${\dot P_{N - 1}}(\tau )$表示N–1阶Legendre多项式的导数,LG点为该多项式在区间(–1,1)内的根。则[–1,1]上状态变量和控制变量的近似值$\hat X(\tau )$和$\hat U(\tau )$可写作
$\hat X\left( \tau \right) \approx \hat X(\tau ) = \sum\limits_{i = 0}^N {\hat X({\tau _i})} {L_i}(\tau )$ (11)
$U\left( \tau \right) \approx \hat U(\tau ) = \sum\limits_{i = 1}^N {\hat U({\tau _i})} L_i^*(\tau )$ (12)
其中:$\hat X({\tau _i})$为在第i个LG点(i=0时为初始点)上的状态变量;$\hat U({\tau _i})$为在第i个LG点上的控制变量;Lagrange插值多项式${L_i}(\tau )$和$L_i^*(\tau )$的定义为
${L_i}(\tau ) = \prod\limits_{j = 0,j \ne i}^N {\frac{{\tau - {\tau _j}}}{{{\tau _i} - {\tau _j}}},\,i = 0,\cdots ,N}$ (13)
$L_i^ * (\tau ) = \prod\limits_{j = 1,j \ne i}^N {\frac{{\tau - {\tau _j}}}{{{\tau _i} - {\tau _j}}},\,i = 1,\cdots ,N} $ (14)
末端状态的近似值则需要通过动力学方程$\dot X\left( \tau \right) = f({\text{ }}X\left( \tau \right),u(\tau ),\tau )$求得。将
$X\left( {{\tau _f}} \right) = X\left( {{\tau _0}} \right) + \int_{ - 1}^1 {f({\text{ }}X\left( \tau \right),u(\tau ),\tau )} {\text{d}}\tau $
离散化可得
$X\left( {{\tau _f}} \right) = {\text{ }}X\left( {{\tau _0}} \right) - \frac{{{t_f}}}{2}\sum\limits_{i = 1}^N {{\omega _i}} f(\hat X\left( {{\tau _i}} \right),\hat U({\tau _i}),{\tau _i};{t_0},{t_f})$ (15)
其中:${\omega _i} = \int_{ - 1}^1 {{L_i}(\tau )} {\rm d}\tau $为Gauss求积权重;τi表示各LG点。这样原最优控制问题的状态和控制约束可表示为
${l_h}{\text{ }} \leqslant h(\hat X(\tau ),\hat U(\tau )){\text{ }} \leqslant {u_h}$ (16)
其中:uhlh分别表示这些约束的上、下界。

为使状态变量始终满足动力学方程约束,动力学方程的左端通过在LG点处对状态的近似值求导进行近似

$\dot X\left( \tau \right) \approx \frac{2}{{{t_f}}}\sum\limits_{i = 0}^N {\hat X({\tau _i})} {\dot L_i}(\tau )$ (17)
其中${\dot L_i}(\tau )$为
${\dot L_i}(\tau ) = \sum\limits_{l = 0}^N {\frac{{\prod\limits_{j = 0,j \ne i,l}^N {(\tau - {\tau _j})} }}{{\prod\limits_{j = 0,j \ne i}^N {({\tau _i} - {\tau _j})} }}} $ (18)

于是,动力学微分方程约束就转换为代数形式

$\sum\limits_{i = 0}^N {{{\dot L}_i}} (\tau )\hat X({\tau _i}) - \frac{{{t_f}}}{2}f(\hat X\left( {{\tau _i}} \right),\hat U({\tau _i}),{\tau _i};{t_0},{t_f}) = 0$ (19)

这样,经过在LG点进行离散化,并利用Lagrange插值多项式对状态变量和控制变量进行全局近似后,原最优控制问题的边值约束和路径约束转化为式(16)的形式,微分方程形式的动力学约束由式(1)的动态形式转化为式(19)的静态代数约束形式。原动态连续最优控制问题转化为静态的非线性规划问题,可以使用非线性规划问题的求解工具方便地进行解算。

3 数学仿真与分析

以小行星Eros 433为目标小天体,利用Gauss伪谱法进行下降和着陆阶段燃耗最优轨迹优化仿真。仿真中用到的Eros 433小行星的物理参数如表 1所示,引力场采用4阶模型。参照文献[7],探测器的发动机比冲设为300 s,最大推力设为22 N,仿真中采用的边界条件如表 2所示[8]

表 1 433 Eros小行星物理参数 Table 1 Physical parameters of 433 Eros

表 2 仿真采用的边界条件 Table 2 Initial and final states in the simulation

仿真得到的优化结果如图 1图 5所示。其中,图 1为探测器三轴位置的时间历程,图 2为三轴速度大小的时间历程,图 3为三轴推力大小的时间历程,图 4为总推力大小的时间历程,图 5为得到的燃耗最优下降轨迹。优化得到的末端时刻探测器的质量为793.7 kg,即消耗燃料6.3 kg;tf = 1 237 s,即下降和着陆过程持续1 237 s.

图 1 探测器三轴位置曲线 Fig. 1 Spacecraft position curves

图 2 探测器三轴速度曲线 Fig. 2 Spacecraft velocity curves

图 3 探测器三轴推力大小曲线 Fig. 3 Spacecraft thrust curves in each axes

图 4 探测器总推力大小曲线 Fig. 4 Spacecraft thrust magnitude curve

图 5 探测器的燃耗最优下降轨迹 Fig. 5 Fuel-optimal descent trajectory of the spacecraft

从以上数学仿真结果可以看出:

1)优化结果满足各项约束条件。经过1 237 s,探测器以零速度软着陆于预设着陆点,着陆时刻探测器质量大于探测器干重,同时优化得到的控制力始终在设定的发动机推力上限之内。

2)从图 3图 4可以看到,发动机在下降和着陆过程的初始和后半阶段处于点火状态,在中间阶段有近400 s的时间内不产生推力,以降低燃耗。在发动机工作的时段,虽然发动机的三轴推力大小不断变化,但总推力的大小基本保持在22 N,即发动机在点火工作时处于最大推力状态,只是推力方向不断变化。

3)从图 4显示的最优轨迹解算结果看,发动机在下降的初始阶段点火初步降低探测器速度,但探测器基本沿原轨迹上向小行星表面运行,在着陆阶段再次点火对探测器进行持续制动,并使其改变方向落向目标着陆点,最终实现软着陆。整个下降与着陆过程符合燃耗最优目标的预期。

4 结论

本文针对小天体着陆或采样返回探测中的下降和着陆阶段,基于伪谱法进行了燃耗最优下降轨迹优化,得到了符合初始和末端状态约束、动力学约束和控制约束的燃耗最优下降轨迹。首先对最优下降轨迹问题进行数学描述,给出了小天体探测下降阶段的动力学模型,建立了包括边值约束和路径约束在内的各项约束条件和燃耗最优指标函数,得出轨迹优化的原始连续最优控制问题。然后利用Gauss伪谱法对此最优控制问题进行离散化,将其转化为离散的非线性规划问题。Gauss伪谱法避免了间接法的推导复杂,以及协态变量无物理意义不易猜测的问题,对初始值不敏感,节约了计算时间,是一种快速、高精度的轨迹优化方法。数学仿真结果显示,采用Gauss伪谱法得到的最优下降轨迹以零速度到达预定着陆点,符合各项约束条件,并通过发动机分阶段工作的方式达到优化燃耗的目标。

参考文献
[1] 崔平远, 乔栋. 小天体附近轨道动力学与控制研究现状与展望[J]. 力学进展, 2013, 43(5): 526-539.Cui P Y, Qiao D. State-of-the-art and prospects for orbital dynamics and control near small celestial bodies[J]. Advances in Mechanics, 2013, 43 (5): 526-539.(1)
[2] Dunham David W, Farquhar Robert W, McAdams James V, et al. Implementation of the first asteroid landing[J]. Icarus, 2002, 159(2): 433-438.(1)
[3] Uo M, Shirakawa K, Hashimoto T, et al. Hayabusa's touching-down to Itokawa-autonomous guidance and navigation[J]. The Journal of Space Technology and Science, 2006, 22(1): 41.(1)
[4] 雍恩米, 陈磊, 唐国金. 飞行器轨迹优化数值方法综述[J]. 宇航学报, 2008, 29(2): 397-406.Yong E M, Chen L, Tang G J. A survey of numerical methods for trajectory optimization of spacecraft[J]. Journal of Astronautics, 2008, 29(2): 397-406.(1)
[5] Stephen Broschart, Daniel Scheeres, Spacecraft descent and translation in the small-body fixed frame[M]. [S.L]: American Institute of Aeronautics and Astronautics, 2004.(1)
[6] 张万里. 轨道转移飞行器的轨迹优化与制导算法研究[D]. 哈尔滨: 哈尔滨工业大学, 2011.Zhang W L. Research on trajectory optimization and guidance of orbital transfer vehicle[D]. Harbin: Harbin Institute of Technology,2001.(1)
[7] 崔平远, 朱圣英, 崔祜涛. 小天体自主软着陆脉冲机动控制方法研究[J]. 宇航学报, 2008, 29(2): 511-516.Cui P Y, Zhu S Y, Cui H T. Autonomous impulse maneuver control method for soft landing on small bodies[J]. Journal of Astronautics, 2008, 29(2): 511-516.(1)