一种无初值猜测的不规则小行星着陆轨迹优化方法与流程

本发明涉及航空航天技术领域，尤其涉及一种无初值猜测的不规则小行星着陆轨迹优化方法。

背景技术：

小行星种类数目众多，具有很高的科学价值。21世纪以来，丰富的小行星探测任务表明小行星探测已成为深空探测的主要方向之一。到目前为止，人类已开展了五次专门的小行星探测任务，其中多次任务涉及不规则小行星表面着陆或附着，如美国的舒梅克号探测器着陆了小行星433eros，日本的隼鸟号探测器附着到了小行星25143itokawa的表面。由于探测器着陆到小行星表面可以获得更高分辨率的小行星数据也可以采集小行星表面样本，在小行星探测任务中非常重要。而轨迹优化则可以用于设计低消耗的着陆轨迹。此外，小行星的形状也往往具有不规则性。因此，在小行星探测设计中，对不规则小行星着陆轨迹优化将是任务设计的重要基础。

不规则小行星着陆轨迹优化方法中，基于同伦法求解不规则小行星着陆燃料最优轨迹,需要随机猜测协态变量的初值。

无初值猜测轨迹优化方法中，基于双脉冲解推导了小推力转移轨迹燃料最优控制协态变量初值的近似解析解，避免了初值猜测问题。但是该方法推导过程是基于质点引力场开展的，对于不规则小行星引力场中着陆轨迹优化将不适用。

综上,现有技术中缺乏一种无需初值猜测的轨迹优化方法,能适用不规则的小行星的引力场。

技术实现要素：

本发明提供一种无初值猜测的不规则小行星着陆轨迹优化方法，能够适用不规则小行星着陆轨迹的优化。

为达到上述目的，本发明采用如下技术方案：

一种无初值猜测的不规则小行星着陆轨迹优化方法，包括：

s1、设定小行星着陆轨迹优化问题所涉及的参数；

s2、依据小行星着陆轨迹优化问题所涉及的参数建立探测器的动力学方程和最优控制方程；

s3、求解惯性坐标系下质点引力场兰伯特问题,得到初始速度；

s4、依据初始速度求解小行星本体坐标系下不规则引力场兰伯特问题,得到双脉冲解；

s5、利用双脉冲解估计着陆轨迹最优控制问题的协态变量初值与控制开关时间点；

s6、利用协态变量初值与控制开关时间点的估计值，通过打靶法得到最优控制解和最优着陆轨迹。

进一步的，在s1中，参数包括：

小行星的自转角速度ω、密度σ及多面体模型参数，根据密度σ和小行星的多面体模型，计算得到小行星的中心引力常数；

探测器初始质量m0、推力幅值tmax、比冲isp、探测器初始位置r0、出差速度v0、末端rf、速度vf、飞行时间；

轨迹优化数值求解过程中归一化计算的长度单位l，时间单位t和质量单位m。

进一步的，s2包括：

建立小行星本体坐标系，本体坐标系的原点位于小行星的质心，本体坐标系的x轴和z轴分别与最小惯量主轴和最大惯量主轴重合，y轴构成右手系，在小行星的本体坐标系下，探测器的动力学方程为

其中，r为探测器位置、v为探测器速度、m为探测器质量、u表征推力大小、α为推力单位方向向量、ω为小行星自转角速度、g0＝9.80665m/s为海平面处重力加速度常数、u为小行星的引力势函数，u对r处的导数即小行星的引力，小行星的引力通过多面体模型方法计算获得；

选取优化指标为其中，t表示时间，下标t0表示初始时刻和tf表示末端时刻；

根据庞特里亚金极小值原理，得到最优控制方向和大小为

其中，ρ为开关函数，协态变量方程组为

横截条件为

λm(tf)＝0(4)

其中，λr表示位置的协态变量、λv表示速度的协态变量、λm表示质量的协态变量；

双脉冲控制情形下，初末时刻之间的动力学方程组为

协态变量方程组为

速度协态变量与脉冲方向反向，即

其中，δv(t0)为初始时刻脉冲，δv(tf)为末端时刻脉冲，a和b为常数；

根据lawden理论，在脉冲施加前后瞬时速度协态变量保持不变，则

λv(t0)＝λv(t0^-)＝λv(t0⁺),λv(tf)＝λv(tf^-)＝λv(tf⁺)(8)

其中，上标-和上标+分别表示前瞬时和后瞬时，而开关函数为零，则

其中，m1和mf分别为施加第一次脉冲后和第二次脉冲后的质量；结合横截条件λm(tf⁺)＝0，并综合式(7)、(8)、(9)整理可得双脉冲情形下速度和质量协态变量的计算式，标记为最优控制方程

进一步的，s3包括：

建立惯性坐标系，设定初始时惯性坐标系的三轴分别与小行星本体坐标系三轴重合；

将小行星本体坐标系中初始时刻的位置和末端时刻的位置转换至惯性坐标系；

求解惯性坐标系的质点引力场兰伯特问题，获得惯性坐标系的初始速度,将惯性坐标系的初始速度转换至小行星本体坐标系。

进一步的,s4包括:

建立小行星本体坐标系下不规则引力场兰伯特问题，即不规则引力场下无推力转移轨道两点边值问题；

转移轨道动力学方程为：

其中，r为探测器位置、v为探测器速度、ω为小行星自转角速度、u为小行星的引力势函数，u对r处的导数为小行星的引力，

设定转移轨道初始时刻位置为r(t0)＝r0，将s3中获得的初始速度标记为初始估计值，通过打靶法求解初始速度，使得积分转移轨道动力学方程至末端时刻后位置满足条件r(tf)＝rf；

根据转移轨道的初始速度与探测器初始速度之差、转移轨道的末端速度与探测器末端速度之差，计算得到双脉冲解，即转移轨道的初始速度与探测器初始速度之差δv(t0)和探测器末端速度与转移轨道的末端速度之差δv(tf)。

进一步的,s5包括：

根据双脉冲解δv(t0)和δv(tf)，计算双脉冲解中m1和mf，

根据最优控制方程

解析计算双脉冲情形下λv(t0)、λv(tf)和λm(t0^-)，λv(t0)和λm(t0^-)用作对应的最优控制协态变量初值的估计值；

以协态变量方程组(6)为常微分方程，初始时刻速度和末端时刻速度的协态变量为两点边值条件，选取位置协态变量初始估计值为零并采用打靶法获得位置协态变量初值解，用作最优控制情形的估计值，根据

估计控制开关时间点。

进一步的,s6包括：

根据最优控制需满足的末端位置、速度、质量协态变量条件和控制开关时间点开关函数为零的条件，建立最优控制打靶方程为

φ(λr0,λv0,λm0,ts1,ts2)＝[r(tf)-rf；v(tf)；λm(tf)；ρ(ts1)；ρ(ts2)]＝0(12)

在进行数值积分时，将轨迹按飞行时间分为[t0,ts1]、[ts1,ts2]、[ts2,tf]三段，并设置[t0,ts1]、[ts1,ts2]、[ts2,tf]各段推力依次为满推、关机和满推；

根据探测器的动力学方程

依次对[t0,ts1]、[ts1,ts2]、[ts2,tf]三段时间进行数值积分；

采用最优控制情形的估计值为初始值，求解最优控制打靶方程获得协态变量初值和控制开关时间点；

基于最优控制打靶方程解，通过数值积分获得最优控制解，并输出最优着陆轨迹。

本发明的有益效果是：

本发明提供的不规则小行星着陆轨迹优化方法，无需提供协态变量初值的猜测值，能适用不规则的小行星的引力场；

本发明为打靶方程中协态变量初值和控制开关点提供估计值，收敛性好；

本发明通过将轨迹分段依次积分，每段轨迹对应推力大小为常数，避免了被积动力学方程右端项不连续的问题。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1是本发明的流程示意图；

图2是本发明得到的用于不规则小行星着陆的最优控制解；

图3本发明得到的用于不规则小行星着陆的最优轨迹。

具体实施方式

为使本领域技术人员更好地理解本发明的技术方案，下面结合具体实施方式对本发明作进一步详细描述。

本发明实施例提供一种无初值猜测的不规则小行星着陆轨迹优化方法，以小行星433eros着陆轨迹优化为例，本发明的流程图如图1所示，具体步骤包括：

步骤一，设置小行星着陆轨迹优化问题参数

选择的目标小行星为433eros，设定其自转角速度ω＝[0,0,3.31458761013812×10^-4]rad/s、密度σ＝2.67g/cm³及多面体模型参数。小行星多面体模型参数选择near-a-5-collected-models-v1.0模型库中具有1708个面的多面体模型参数。根据小行星密度和小行星多面体模型，计算小行星中心引力常数，其值为4.4390168×10^-4km³/s²。设置探测器初始质量m0，推力幅值tmax及比冲isp的数值分别为1400kg、80n、225s。设置探测器初始位置r0、速度v0和末端rf、速度vf的数值分别为[22.167,-26.860,-0.029]km、[-6.200,-5.087,0.002]m/s、[9.302,-6.213,1.478]km、[0,0,0]m/s。设置飞行时间为2400s。轨迹优化数值求解过程中采用归一化以提升计算效率，设置归一化计算所需的长度单位l，时间单位t和质量单位m为30km、3016.966566s、1400kg。

步骤二，建立动力学方程和最优控制方程

建立小行星本体坐标系，原点位于小行星质心，x轴和z轴分别与最小和最大惯量主轴重合，y轴构成右手系。在该本体坐标下，探测器的轨道动力学方程组为

其中，r、v和m分别为探测器位置、速度和质量，u表征推力大小，α为推力单位方向向量，ω为小行星自转角速度，g0＝9.80665m/s为海平面处重力加速度常数，u为小行星的引力势函数。u对r处的导数即小行星的引力通过多面体模型方法计算获得。

选取优化指标为其中t表示时间，下标0和f分别表示初始和末端时刻。

据庞特里亚金极小值原理，可以得到最优控制方向和大小为

其中，ρ为开关函数。协态变量方程组为

横截条件为

λm(tf)＝0(4)

其中，λr、λv和λm分别表示位置、速度和质量的协态变量。

双脉冲控制情形下，初末时刻之间的动力学方程组为

协态变量方程组为

速度协态变量与脉冲方向反向，即

其中，δv(t0)为初始时刻脉冲，δv(tf)为末端时刻脉冲，a和b为常数。

根据lawden理论，在脉冲施加前后瞬时速度协态变量保持不变，则

λv(t0)＝λv(t0^-)＝λv(t0⁺),λv(tf)＝λv(tf^-)＝λv(tf⁺)(8)

其中，上标-和+分别表示前后瞬时。而开关函数为零，则

其中，m1和mf分别为施加第一次和第二次脉冲后的质量。结合横截条件λm(tf⁺)＝0，并综合(7)、(8)、(9)整理可得双脉冲情形下速度和质量协态变量的计算式为

步骤三，求解惯性坐标系下质点引力场兰伯特问题

建立惯性坐标系，假设初始时刻该坐标系三轴分别与小行星本体坐标系三轴重合。将小行星本体坐标系下初末位置转换至惯性坐标系。求解质点引力场兰伯特问题即经典兰伯特问题，获得初始速度并转换至小行星本体坐标系。获得的归一化初始速度值为[-1.285598,0.429747,6.723504×10^-2]。

步骤四，求解小行星本体坐标系下不规则引力场兰伯特问题

建立小行星本体坐标系下不规则引力场兰伯特问题，即不规则引力场下无推力转移轨道两点边值问题。选择转移轨道动力学方程为式(5)。设定转移轨道初始时刻位置r(t0)＝r0，选取步骤三中获得的初始速度为初始估计值并通过打靶法求解初始速度，使得积分转移轨道动力学方程至末端时刻后位置满足条件r(tf)＝rf。获得的归一化初始速度值为[-1.299562,0.413655,6.800592×10^-2]。根据所得转移轨道初末速度与探测器初末速度之差计算得到所需双脉冲即δv(t0)和δv(tf)。获得的双脉冲的归一化值分别为[-0.676056,0.925232,6.780480×10^-2]、[-0.316602,-1.102242,-3.397051×10^-2]。

步骤五，基于双脉冲解，估计最优控制协态变量初值和控制开关时间点

根据步骤四所得δv(t0)和δv(tf)，计算双脉冲解中m1和mf，其归一化值分别为0.994840，0.989710。根据式(10)解析计算双脉冲情形下λv(t0)、λv(tf)和λm(t0-)，其值分别为[2.626812×10^-3,-3.594987×10^-3,-2.634553×10^-4]、[1.230801^-3,4.284994^-3,1.320613^-4]、1.029032^-2，λv(t0)和λm(t0^-)用作对应的最优控制协态变量初值的估计值。以动力学方程组(6)为常微分方程，初末时刻速度协态变量为两点边值条件，选取位置协态变量初始估计值为零并采用打靶法获得位置协态变量初值解，用作最优控制情形的估计值。该值为[2.632955×10^-3,-1.209140×10^-2,-4.339755×10^-4]。根据

估计控制开关时间点，其归一化数值分别为6.603973×10^-2和0.729836。

步骤六，计算最优控制解，并输出燃料最优轨迹

根据最优控制需满足的末端位置、速度、质量协态变量条件和控制开关时间点开关函数为零的条件，建立最优控制打靶方程为

φ(λr0,λv0,λm0,ts1,ts2)＝[r(tf)-rf；v(tf)；λm(tf)；ρ(ts1)；ρ(ts2)]＝0(12)

在进行数值积分时，将轨迹按飞行时间分为[t0,ts1]、[ts1,ts2]、[ts2,tf]三段并设置各段推力依次为满推、关机和满推。然后，根据动力学方程组(1)依次对每段时间进行数值积分。采用步骤五中获得的估计值为初始值并求解最优控制打靶方程获得协态变量初值和控制开关时间点，所得值分别为[3.096729×10^-3,-1.457667×10^-2,-4.915701×10^-4]、[3.166020×10^-3,-4.423136×10^-3,-3.041063×10^-4]、1.258645×10^-2、7.266694×10^-2、0.721790。对应的最优控制解如图2所示。基于最优控制打靶方程解，通过数值积分获得最优控制解，并输出最优轨迹，如图3所示。

本发明的有益效果是：

本发明提供的不规则小行星着陆轨迹优化方法，无需提供协态变量初值的猜测值，能适用不规则的小行星的引力场；

本发明为打靶方程中协态变量初值和控制开关点提供估计值，收敛性好；

本发明通过将轨迹分段依次积分，每段轨迹对应推力大小为常数，避免了被积动力学方程右端项不连续的问题。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。