基于hp自适应伪谱法的太阳帆星际转移轨迹优化方法

摘要

本发明提出了一种基于hp自适应伪谱法的太阳帆星际转移轨迹优化方法，适用于太阳帆从某行星轨道航行至另一行星轨道时，寻求一条最优转移轨迹，使其轨道转移的时间最小，用于解决现有技术中存在的收敛速度慢的技术问题，实现步骤包括：建立极坐标系下的太阳帆星际转移运动方程；构建太阳帆星际转移轨迹的最优控制问题；采用hp自适应伪谱法对太阳帆星际转移轨迹的最优控制问题进行优化，得到太阳帆星际转移的最优飞行轨迹。本发明根据计算精度的要求，自适应地调整相应区间的网格个数或多项式的阶次，从而具有更少的配点个数和更合理的配点分布，减少了计算量，有效地提高了收敛速度。

说明书

技术领域

本发明属于深空探测技术领域，涉及一种太阳帆转移轨迹优化方法，具体涉及一种基于hp自适应伪谱法的太阳帆星际转移轨迹优化方法，适用于太阳帆从某行星轨道航行至另一行星轨道时，寻求一条最优转移轨迹，使其轨道转移的时间最小。

背景技术

随着深空探测任务的复杂性及空间距离的增加，传统的推进方式已经很难满足其要求。太阳帆是一种利用太阳光子在大面积轻质薄膜表面反射产生推力的新型航天器，其动力来源是天然的太阳能，不需要消耗化学燃料。因此，太阳帆在太空中的航行寿命不受有限燃料的制约。此外，连续小推力、结构轻便、成本低廉、发射风险小等特点，使得太阳帆在深空探测、星际航行等长周期的航天任务方面具有很大的优势。

太阳帆涉及到的关键技术包括：结构设计、材料选取、太阳帆的折叠储存和展开、姿态控制、轨迹优化等。其中，航天器的轨迹优化能够实现减少航天器在轨运行的燃料消耗，延长在轨寿命，增大执行任务能力等最优性能指标。对于执行深空探测任务的太阳帆而言，由于其转移至目标轨道的时间跨度很大，因此十分有必要对其航行的轨迹进行优化，缩短其转移时长。

太阳帆星际转移最优轨迹的生成本质上是求解非线性、带有状态约束和控制约束的动态最优控制问题，其求解非常困难。目前，太阳帆星际转移轨迹优化方法主要分为间接法和直接法。间接法存在包括推导最优性条件困难、收敛半径小、以及协态变量初值猜测困难的缺点，因此，间接法的操作比较困难。直接法的基本思想是将连续的太阳帆星际转移最优控制问题离散化并转换为NLP问题，再通过优化算法对性能指标直接寻优。直接法具有初估值敏感度低、不需推导一阶最优性条件、易于程序化的优点，因此相对于间接法的应用更为广泛。例如：2013年，Naiming Qi在Journal of Astronautics的34卷第5期发表的《Electric Sail Trajectory Optimization and Performance Analysis》，公开了一种采用Gauss伪谱法对太阳帆从地球至火星转移轨迹进行优化的方法。该方法首先建立了极坐标系下的太阳帆运动方程，然后构建了太阳帆地球—火星转移轨迹的最优控制问题，最后采用Gauss伪谱法对该最优控制问题进行了优化求解。Gauss伪谱法也可以实现太阳帆星际转移轨迹的优化问题，但是由于其配点分布形式固定，往往需要较高维的插值多项式才能得到理想的近似解，计算量较大，导致收敛速度较慢。

发明内容

本发明的目的在于克服上述已有技术的不足，提出了一种基于hp自适应伪谱法的太阳帆星际转移轨迹优化方法，用于解决现有技术中存在的收敛速度慢的技术问题。

为实现上述目的，本发明采取的技术方案包括如下步骤：

(1)建立极坐标系下的太阳帆星际转移运动方程：

(1a)假设太阳帆为完全反射太阳光子的理想平面，计算太阳光子作用在太阳帆面上的太阳光压力F_s；

(1b)沿太阳帆星际转移轨道的径向和切向对太阳光子作用在太阳帆面上的太阳光压力F_s进行分解，得到太阳光压在轨道径向方向的加速度大小a_r和切向方向的加速度大小a_θ；

(1c)将太阳和太阳帆作为质点，通过步骤(1b)得到的径向光压加速度a_r和切向光压加速度a_θ，建立极坐标系下的太阳帆星际转移运动方程；

(2)构建太阳帆星际转移轨迹的最优控制问题：

将步骤(1)所述太阳帆星际转移运动方程转化为状态微分方程约束，并将太阳帆发射位置的状态量x₀作为初始条件，在状态变量x满足状态微分方程约束及终端入轨约束时，优化控制变量u，使得性能指标J取极小值；

其中，状态变量x取太阳帆的径向速度v_r和切向速度v_θ，以及步骤(1)所述太阳帆星际转移运动方程中的太阳帆与太阳的距离r和太阳帆相对于太阳的角位置θ；控制变量u取步骤(1)所述太阳帆星际转移运动方程中的俯仰角α；性能指标J取太阳帆星际转移的时间；

(3)采用hp自适应伪谱法对步骤(2)中太阳帆星际转移轨迹的最优控制问题进行优化，得到太阳帆星际转移的最优飞行轨迹：

(3a)设置太阳帆星际转移轨迹的总区间个数为K，当前区间为k，{k|1≤k≤K}，每个区间初始的配点数为M，M≥1，误差门限为ε，相对曲率阈值为r_max；

(3b)在每个区间的配点上将最优控制问题离散化为NLP问题，并对其进行求解，得到一条飞行轨迹对应的状态变量值x(τ)和控制变量值u(τ)；

(3c)通过状态变量值x(τ)和控制变量值u(τ)，计算每个区间的最大残差并根据和ε的大小关系，判断步骤(3b)得到的飞行轨迹是否为最优轨迹，若所有区间的则步骤(3b)得到的飞行轨迹为最优轨迹，优化结束，若存在某些区间的计算这些区间的相对曲率并执行步骤(3d)；

(3d)更新的区间的配点个数：

判断这些区间和r_max的大小关系，若则重新确定对应区间的配点个数N_k，若则将对应区间划分为n_k个子区间，并将每个子区间内的配点个数设置为M，得到这些区间新的配点个数及分布，并执行步骤(3b)，进行下一次迭代。

本发明与现有技术相比，具有如下优点：

本发明由于在求解太阳帆星际转移轨迹的最优控制问题时，采用hp自适应伪谱法作为求解方法，该方法根据计算精度的要求，自适应地调整相应区间的网格个数或多项式的阶次，从而具有更少的配点个数和更合理的配点分布，减少了计算量，与现有技术相比，有效地提高了收敛速度。

附图说明

图1为本发明的实现流程图；

图2为太阳光压力沿太阳帆地球—火星转移轨道径向和切向分解的示意图。

具体实施方式

以下结合附图和太阳帆地球—火星转移轨迹优化的具体实施例，对本发明作进一步详细说明。

参照图1、一种基于hp自适应伪谱法的太阳帆地球—火星转移轨迹优化方法，包括如下步骤：

步骤1、建立极坐标系下的太阳帆地球—火星转移运动方程：

步骤1a)假设太阳帆为完全反射太阳光子的理想平面，计算太阳光子作用在太阳帆面上的太阳光压力F_s，其表达式为：

其中，m为太阳帆的质量；β为太阳帆的光照度；μ为太阳引力系数；r为太阳帆与太阳的距离；r为入射光线单位向量；n为太阳帆单位法向量。

步骤1b)沿太阳帆地球—火星转移轨道的径向和切向对太阳光子作用在太阳帆面上的太阳光压力F_s进行分解，如附图2所示，太阳光压力F_s的方向垂直于理想帆面，α为太阳入射光线单位向量r与太阳帆单位法向量n之间的夹角，得到太阳光压在轨道径向方向的加速度大小a_r和切向方向的加速度大小a_θ，其表达式分别为：

其中，β为太阳帆的光照度；μ为太阳引力系数；r为太阳帆与太阳的距离；α为太阳帆的俯仰角。

步骤1c)将太阳和太阳帆作为质点，通过步骤1b)得到的径向光压加速度a_r和切向光压加速度a_θ，建立极坐标系下的太阳帆地球—火星转移运动方程，其表达式为：

其中，r为太阳帆与太阳的距离；θ为太阳帆相对于太阳的角位置；μ为太阳引力系数；a_r和a_θ分别为太阳光压在地球—火星转移轨道径向和切向方向的加速度大小。

步骤2、构建太阳帆地球—火星转移轨迹的最优控制问题：

将步骤1)所述太阳帆地球—火星转移运动方程转化为状态微分方程约束，并将太阳帆在地球轨道的状态量x₀作为初始条件。在状态变量x＝[r,θ,v_r,v_θ]^T满足状态微分方程约束及终端入轨约束时，优化太阳帆的俯仰角α，使得性能指标J取极小值。其中，状态微分方程约束的表达式为：

初始条件x₀为：

r₀＝1,v_r0＝0,v_θ0＝1,θ₀＝0.1

终端入轨约束为太阳帆进入火星轨道的状态量x_f：

性能指标J为太阳帆地球—火星转移的时间，其表达式为：

其中，r为太阳帆与太阳的距离；θ为太阳帆相对于太阳的角位置；v_r和v_θ分别为太阳帆的径向速度和切向速度；a_r和a_θ分别为太阳光压在地球—火星转移轨道径向和切向方向的加速度大小；μ为太阳引力系数；t₀为太阳帆航行的初始时刻；t_f为太阳帆进入火星轨道的时刻。

步骤3、采用hp自适应伪谱法对步骤2)中太阳帆地球—火星转移轨迹的最优控制问题进行优化，得到太阳帆地球—火星转移的最优飞行轨迹：

步骤3a)设置太阳帆地球—火星转移轨迹的总区间个数为K，当前区间为k，{k|1≤k≤K}，每个区间初始的配点数为M，M≥1，误差门限为ε，相对曲率阈值为r_max；

步骤3b)在每个区间的配点上，采用Radau伪谱法将最优控制问题离散化为NLP问题，并对其进行求解，得到一条飞行轨迹对应的状态变量值x(τ)和控制变量值u(τ)；

步骤3c)通过状态变量值x(τ)和控制变量值u(τ)，计算每个区间的最大残差其表达式为：

根据和ε的大小关系，判断步骤3b)得到的飞行轨迹是否为最优轨迹，若所有区间的则步骤3b)得到的飞行轨迹为最优轨迹，优化结束；若存在某些区间的计算这些区间的相对曲率并执行步骤3d)，其表达式为：

其中，为第s个采样点处的第i个状态变量的一阶微分残差；为第s个采样点处的第j个路径约束残差；和分别为所有采样点中的最大曲率和平均曲率。第s个采样点处的状态变量一阶微分残差路径约束残差以及曲率κ^(k)的表达式分别为：

其中，f为状态微分方程约束；C为路径约束；为第s个采样点处的状态变量值；为第s个采样点处的控制变量值。

步骤3d)更新的区间的配点个数：

判断这些区间和r_max的大小关系，若说明该段轨迹较平滑，则重新确定对应区间的配点个数N_k，其表达式为：

若说明该段轨迹过于震荡，则将对应区间划分为n_k个子区间，其表达式为：

并将每个子区间内的配点个数设置为M，得到这些区间新的配点个数及分布，并执行步骤3b)，进行下一次迭代。其中，M为每个区间初始的配点数；ceil(·)表示向上取整；ε为误差门限；为每个区间的最大残差；A为大于0的整数，为一个可调因子；B为大于0的整数，为一个可调因子。

以下结合仿真实验对本发明的技术效果进一步说明：

1、仿真条件：

本发明和现有的采用Gauss伪谱法对太阳帆地球—火星转移轨迹的优化仿真都是在主频2.6GHZ的Intel(R)Core(TM)i5-3230M处理器、内存4GB的硬件环境和MATLAB R2014a的软件环境下进行的。在本发明的仿真中，太阳帆的光照度β＝0.17；太阳帆姿态调整是连续的，其俯仰角α满足：配点精度要求设置为10^-6。

2、仿真结果分析：

本发明和现有的采用Gauss伪谱法对太阳帆地球—火星转移轨迹的优化仿真，所使用的配点个数和计算耗时如表1所示：

表1

由表1中的对比结果可以看出，本发明使用了52个配点，耗时7.611s，而Gauss伪谱法使用了60个配点，耗时18.402s。本发明使用的配点个数和计算耗时更少，优化的收敛速度也更快，这是由于本发明将太阳帆地球—火星转移轨迹的震荡区域划分为多个较为平滑的区间，在每个区间内以低阶插值多项式来逼近最优轨迹，相对于Gauss伪谱法使用高阶全局多项式逼近，在配点个数和计算量上有所减少。