CN107168056B - 一种自适应的月球软着陆轨道快速优化控制器 - Google Patents

Abstract

本发明公开了一种自适应的月球软着陆轨道快速优化控制器，该轨道快速优化控制器由仿真模块、优化模块、执行模块组成。在仿真模块中，月球软着陆轨道被分成若干个小段，利用分段处的各个节点的推力方向角及对应的时刻拟合出轨道的推力方向角变化曲线，进而由软着陆动力学模型仿真得到整条软着陆轨道。其中各个节点的最佳推力方向角及最佳软着陆终端时刻由优化模块寻优得到。优化模块采用改进的智能算法，改进的智能优化方法控制了子群的大小，添加了扰动因子，增加了搜索的多样性，同时根据进化状态自适应变化的惯性权重提高了算法的收敛性能。本控制器能够快速寻优得到一条使燃料消耗最少的月球软着陆轨道并最终实现轨道的最优控制。

Description

一种自适应的月球软着陆轨道快速优化控制器

技术领域

本发明涉及航空航天领域，具体地，涉及一种自适应的月球软着陆轨道快速优化控制器。

背景技术

在进行无人月面勘探或载人登月任务时，都需要使着陆器实现月球表面软着陆，以保证仪器设备以及航天员的安全。在典型的月球软着陆过程中，从环月轨道下降到月球表面主要可分为2个阶段：霍曼转移段和动力下降段。其中，霍曼转移段指着陆器从环月停泊轨道通过霍曼变轨进入环月椭圆轨道，并运行至近月点；动力下降段指着陆器从近月点处开始发动机制动抵消水平速度，并最终以较小速度着陆月面。由于霍曼变轨速度增量较小，着陆器的大部分燃料都消耗在动力下降段，故有必要对动力下降段的着陆轨迹进行优化设计以减少燃料消耗。

将常推力月球软着陆轨道离散化是求解月球软着陆轨道问题的一种方法，而离散点处的状态是求解的关键，可使用优化算法求解离散点的状态。但是不恰当的优化算法导致规划出非最优轨道，浪费大量的燃料。

发明内容

为了减少月球软着陆时消耗的燃料，本发明的目的在于提供一种自适应的月球软着陆轨道快速优化控制器。

本发明解决其技术问题所采用的技术方案是：一种自适应的月球软着陆轨道快速优化控制器，该轨道快速优化控制器由优化模块、仿真模块、执行模块组成；其中：

仿真模块把月球软着陆轨道分割成n个小段，并将n的值输入优化模块；每个节点的时刻由式(1)得到：

其中，t_k为第k+1个节点的时刻，k＝0,1,...,n，t₀为初始时间，记t₀＝0；

优化模块把包括初始节点和终节点在内的n+1个节点的推力方向角ψ与软着陆终端时刻t_f作为待优化的参数；初始化种群规模为N_s的粒子群，随机生成维度为n+2的粒子i的初始位置x_i＝(x_i1,x_i2,...,x_i(n+2))和初始速度v_i＝(v_i1,v_i2,...,v_i(n+2))，i＝1,2,...,N_s，并将粒子的位置信息传入仿真模块；定义维度变量d，d＝1,2,...,n+2；当d＝1,2,...,n+1时，x_id代表第d个节点的推力方向角，当d＝n+2时，x_id代表软着陆终端时刻，x_id∈[500,700]，v_id∈[-200,200]；种群规模N_s＝300～600；然后按以下方法进行迭代，初始时迭代计数T＝0：

(1)仿真模块中，月球软着陆过程中推力方向角表示成多项式(2)：

ψ(t)＝λ₀+λ₁t+λ₂t²+λ₃t³(2)

其中ψ(t)表示着陆轨道t时刻的推力方向角，λ₀,λ₁,λ₂,λ₃为系数；优化模块输入的粒子的位置代表n+1个节点的推力方向角及软着陆终端时刻，按照式(1)得到n+1个节点的对应时刻；采用函数逼近法，利用n+1个节点的推力方向角及其对应时刻，对式(2)进行拟合，可以求得多项式的系数λ₀,λ₁,λ₂,λ₃，进而得到整个着陆轨道各个时刻的推力方向角ψ(t)；

(2)仿真模块存储了月球软着陆时着陆器的质心动力学方程，见式(3)：

式(3)中r为着陆器的月心距离，v为着陆器的径向速度，θ为着陆器极角，ω为着陆器极角角速度，μ为月球引力常数，m为着陆器质量，F为制动发动机推力，I_SP为制动发动机比冲；其中月球引力常数μ为常数，μ＝4902.75km³/s²，制动发动机推力F与制动发动机比冲I_SP与实际使用的发动机有关，也为常数；着陆器初始质量m₀根据实际确定；其他参数在着陆器着陆过程中发生变化；

初始条件为：

其中，r_p和r_a分别为霍曼转移段的近地点半径和远地点半径，r_p＝1753km，r_a＝1838km；

将步骤(1)中拟合得到的推力方向角ψ(t)、式(4)的初始条件以及着陆器初始质量m₀带入动力学方程(3)，所有数据单位统一，获得月球软着陆的轨道，并将获得的轨道信息输入给优化模块；

(3)优化模块中，优化目标为软着陆过程消耗燃料最少，即令式(5)中指标J最大：

同时，为实现软着陆，终端约束条件为：

其中，R为月球半径，R＝1738km；

在适应度函数中考虑约束条件，构造适应度函数fitness：

fitness＝J-α[(r(t_f)-R)²+v²(t_f)+ω²(t_f)](7)

其中r(t_f)、v(t_f)、ω(t_f)分别表示仿真模块输入的软着陆轨道终端时刻的月心距离、径向速度、极角角速度；α为罚因子，α＝10000；根据仿真模块输入的软着陆轨道，按照式(7)计算适应度函数值；适应度函数值最大的粒子为全局最优粒子，其位置为p_best＝(p_best1,p_best2,...,p_best(n+2))；

(4)在优化模块中对所有粒子进行分群操作，包括以下子步骤：

(4.1)将所有粒子按照适应度函数值大小从大到小排序，选取适应度函数值最大的粒子作为一个子群中心；

(4.2)在剩下的粒子中选取适应度函数值最大的粒子，依次计算该粒子与各个子群中心的欧几里得距离；粒子i与粒子j的欧几里得距离dist(i,j)定义为：

其中，xi＝(x_i1,x_i2,...,x_i(n+2))代表粒子i的位置，x_j＝(x_j1,x_j2,...,x_j(n+2))代表粒子j的位置，i,j＝1,2,...,N_s；若该粒子与某一个子群中心的欧几里得距离小于半径r，则将该粒子归为该子群中心所在的子群，并不再计算该粒子与剩下的子群中心的欧几里得距离；若该粒子与所有子群中心的距离都大于半径r，则将该粒子置为一个新的子群中心；半径r＝15～25；

(4.3)重复步骤4.2)，直到处理完所有粒子，则分群完成，且每个子群中心为该子群中适应度函数值最大的粒子；

(5)检查每个子群中的粒子数，若一个子群中的粒子数为S，且S>S_max，S_max为子群允许的最大粒子数，则将适应度最差的(S-S_max)个粒子的位置和速度进行重置；S_max＝8～10；

(6)在一个子群中，fitness_1为该子群中最大的适应度函数值，fitness_2为该子群中第二大的适应度函数值，tol为搜索精度，tol＝0.01，若满足式(9)：

fitness_1-fitness_2<tol(9)

则将第二大的适应度函数值的粒子q的位置按照式(10)处理：

x'_qd＝x_qd+η·rand(10)

其中，x_qd为粒子q原本的第d维位置，η为扰动因子，rand为0到1之间的随机数，x'_qd为粒子q扰动以后的第d维位置；扰动因子η的大小为：

η＝0.05(x_max-x_min)(11)

其中，x_max，x_min为粒子的搜索上下限，x_min＝0；

(7)确定种群的进化状态；首先，定义每个粒子与其所在子群的子群中心的距离的绝对值之和d_g：

其中，p_ig＝(p_ig1,p_ig2,...,p_ig(n+2))为粒子i所在子群的子群中心的位置；其次，定义每个粒子与其所在子群的子群中心的距离之和的绝对值D_g：

定义进化因子δ为：

由定义可知进化因子δ∈[0,1]；

(8)按照式(15)(16)更新每个粒子的速度与位置：

v_id(T)＝w·v_id(T-1)+c₁·rand·(p_id-x_id(T-1))+c₂·rand·(p_igd-x_id(T-1))(15)

x_id(T)＝x_id(T-1)+v_id(T)(16)

其中，w为惯性权重；c₁,c₂为加速因子，c₁＝c₂＝2；rand为0到1之间的随机数；p_i＝(p_i1，p_i2,...,p_i(n+2))为粒子x_i的历史最优位置；p_ig＝(p_ig1，p_ig2,...,p_ig(n+2))为粒子i所在的子群的最优粒子的位置；惯性权重按照式(17)变换：

更新后，当d＝1,2,...,n+1时，若x_id<0，则令x_id＝0，若则令当d＝n+2时，若x_id<500，则令x_id＝500，若x_id>700，则令x_id＝700；

(9)迭代计数累加，T＝T+1.

(10)重复步骤(1)～(9)，直到达到最大迭代次数T_max停止迭代，T_max＝100～2000；

种群全局最优粒子p_best＝(p_best1,p_best2,...,p_best(n+2))所在的位置即优化后的n+1个节点的推力方向角和终端时刻；将全局最优粒子的位置信息输入仿真模块，按照式(2)、(3)、(4)获得优化后的月球软着陆轨道，月球软着陆最优轨道规划完成；

仿真模块将最优软着陆轨道通过实时通讯传递给执行模块的执行元件执行，实现燃料消耗最少的月球软着陆最优控制。

本发明的有益效果主要表现在：用直接离散法将轨道控制问题转化为参数优化问题，使优化过程变得简单易行；改进的智能优化方法控制了子群的大小，并且添加了扰动因子，增加了搜索的多样性；根据进化状态自适应变化的惯性权重提高了算法的收敛性能。本控制器具有很强的搜索能力，能够寻找到软着陆轨道的最优参数并且实现使燃料消耗最少的月球软着陆最优控制。

附图说明

图1是轨道离散化图；

图2是月球软着陆极坐标系；

图3是本发明的结构图；

图4是本发明的流程图；

图5是本发明中惯性权重w随进化因子δ的变化图。

具体实施方式

下面根据附图以一个实例具体说明本发明。

一种自适应的月球软着陆轨道快速优化控制器，参照图3，该轨道快速优化控制器由优化模块、仿真模块、执行模块组成；其中：

参照图1，仿真模块把月球软着陆轨道分割成n个小段，本实例中令n取9，并将n的值输入优化模块。每个节点的时刻由式(1)得到：

其中，t_k为第k+1个节点的时刻，k＝0,1,...,9，t₀为初始时间，记t₀＝0。

优化模块把包括初始节点和终节点在内的10个节点的推力方向角ψ与软着陆终端时刻t_f作为待优化的参数。初始化种群规模为N_s的粒子群，随机生成维度为11的粒子i的初始位置x_i＝(x_i1,x_i2,...,x_i11)和初始速度v_i＝(v_i1,v_i2,...,v_i11)，i＝1,2,...,N_s，并将粒子的位置信息传入仿真模块。定义维度变量d，d＝1,2,...,11。当d＝1,2,...,10时，x_id代表第d个节点的推力方向角，由经验可知软着陆过程中推力方向角不会超过因此当d＝11时，x_i11代表软着陆终端时刻，根据实际情况软着陆所花时间为500秒至700秒，因此x_i11∈[500,700]，v_i11∈[-200,200]。种群规模N_s＝300～600。然后按以下方法进行迭代，初始时迭代计数T＝0：

(1)仿真模块中，月球软着陆过程中推力方向角表示成多项式(2)：

ψ(t)＝λ₀+λ₁t+λ₂t²+λ₃t³(2)

其中ψ(t)表示着陆轨道t时刻的推力方向角，λ₀,λ₁,λ₂,λ₃为系数，无实际意义。优化模块输入的粒子的位置代表10个节点的推力方向角及软着陆终端时刻，按照式(1)得到10个节点的对应时刻。采用函数逼近法，利用10个节点的推力方向角及其对应时刻，对式(2)进行拟合，可以求得多项式的系数λ₀,λ₁,λ₂,λ₃，进而得到整个着陆轨道各个时刻的推力方向角ψ(t)。

(2)仿真模块存储了月球软着陆时着陆器的质心动力学方程，见式(3)：

参照图2，式(3)中r为着陆器的月心距离，v为着陆器的径向速度，θ为着陆器极角，ω为着陆器极角角速度，μ为月球引力常数，m为着陆器质量，F为制动发动机推力，I_SP为制动发动机比冲。其中月球引力常数μ为常数，μ＝4902.75km³/s²，制动发动机推力F与制动发动机比冲I_SP与实际使用的发动机有关，也为常数。着陆器初始质量m₀根据实际确定。本实例中F＝1350N，I_SP＝2940m/s，m₀＝600kg。其他参数在着陆器着陆过程中发生变化。

由于起点处于霍曼转移轨道的近地点，故初始条件为：

其中，r_p和r_a分别为霍曼转移段的近地点半径和远地点半径，r_p＝1753km，r_a＝1838km。

将步骤(1)中拟合得到的推力方向角ψ(t)、式(4)的初始条件以及着陆器初始质量m₀带入动力学方程(3)，所有数据单位统一，获得月球软着陆的轨道，并将获得的轨道信息输入给优化模块。

(3)优化模块中，优化目标为软着陆过程消耗燃料最少，即令式(5)中指标J最大：

同时，为实现软着陆，终端约束条件为：

其中，R为月球半径，R＝1738km。

在适应度函数中考虑约束条件，构造适应度函数fitness：

fitness＝J-α[(r(t_f)-R)²+v²(t_f)+ω²(t_f)](7)

其中r(t_f)、v(t_f)、ω(t_f)分别表示仿真模块输入的软着陆轨道终端时刻的月心距离、径向速度、极角角速度。α为罚因子，罚因子应足够大，令α＝10000。适应度函数值越大，说明消耗的燃料越少，规划的轨道越优。若规划的轨道不满足终端约束，则适应度函数值将会变得很小。根据仿真模块输入的软着陆轨道，按照式(7)计算适应度函数值。适应度函数值最大的粒子为全局最优粒子，其位置为p_best＝(p_best1,p_best2,...,p_best11)。

(4)在优化模块中对所有粒子进行分群操作，包括以下子步骤：

(4.1)将所有粒子按照适应度函数值大小从大到小排序，选取适应度函数值最大的粒子作为一个子群中心；

(4.2)在剩下的粒子中选取适应度函数值最大的粒子，依次计算该粒子与各个子群中心的欧几里得距离。粒子i与粒子j的欧几里得距离dist(i,j)定义为：

其中，x_i＝(x_i1,x_i2,...,x_i11)代表粒子i的位置，x_j＝(x_j1,x_j2,...,x_j11)代表粒子j的位置，i,j＝1,2,...,N_s。若该粒子与某一个子群中心的欧几里得距离小于半径r，则将该粒子归为该子群中心所在的子群，并不再计算该粒子与剩下的子群中心的欧几里得距离；若该粒子与所有子群中心的距离都大于半径r，则将该粒子置为一个新的子群中心。根据搜索空间的大小，本实例中半径r＝20。

(4.3)重复步骤4.2)，直到处理完所有粒子，则分群完成，且每个子群中心为该子群中适应度函数值最大的粒子。

(5)检查每个子群中的粒子数，若一个子群中的粒子数为S，且S>S_max，S_max为子群允许的最大粒子数，则将适应度最差的(S-S_max)个粒子的位置和速度进行重置。这样做的目的是防止一个子群中包含过多的粒子而导致搜索多样性下降。增加搜索过程的多样性能增大算法搜索到全局最优的概率从而规划出更优的轨道。在一个子群范围内，粒子数不需要太多，本实例中令S_max＝10。

(6)在一个子群中，fitness_1为该子群中最大的适应度函数值，fitness_2为该子群中第二大的适应度函数值，tol为搜索精度，tol＝0.01，若满足式(9)：

|fitness_1-fitness_2|<tol(9)

则说明这两个粒子的适应度函数值过分接近，将第二大的适应度函数值的粒子q的位置按照式(10)处理：

x'_qd＝x_qd+η·rand(10)

其中，x_qd为粒子q原本的第d维位置，η为扰动因子，rand为0到1之间的随机数，x'_qd为粒子q扰动以后的第d维位置。添加扰动因子可以增加搜索的多样性，从而规划出更优的轨道。但是扰动因子过大又会影响粒子群的正常更新，因此设扰动因子η的大小为：

η＝0.05(x_max-x_min)(11)

其中，x_max，x_min为粒子的搜索上下限，x_min＝0。

(7)确定种群的进化状态。随着粒子的更新，种群共经历四种进化状态，即探索期、开拓期、聚合期以及跳出期。下面利用进化因子来表示进化状态。首先，定义每个粒子与其所在子群的子群中心的距离的绝对值之和d_g：

其中，p_ig＝(p_ig1,p_ig2,...,p_ig11)为粒子i所在子群的子群中心的位置。其次，定义每个粒子与其所在子群的子群中心的距离之和的绝对值D_g：

在进化初始阶段，D_g取值略小于d_g；在进化收敛阶段，D_g取值远小于d_g；在跳出阶段，D_g取值接近于d_g。因此，定义进化因子δ为：

由定义可知进化因子δ∈[0,1]。

(8)按照式(15)(16)更新每个粒子的速度与位置：

v_id(T)＝w·v_id(T-1)+c₁·rand·(p_id-x_id(T-1))+c₂·rand·(p_igd-x_id(T-1))(15)

x_id(T)＝x_id(T-1)+v_id(T)(16)

其中，w为惯性权重；c₁,c₂为加速因子，c₁＝c₂＝2；rand为0到1之间的随机数；p_i＝(p_i1，p_i2,...,p_i11)为粒子x_i的历史最优位置；p_ig＝(p_ig1，p_ig2,...,p_ig11)为粒子i所在的子群的最优粒子的位置。

惯性权重w越大，算法的搜索能力越强，反之亦然。在探索期，希望惯性权重大一些，在聚合期，希望惯性权重小一些。由于进化因子可以反映进化状态，参照图5，惯性权重按照式(17)变换：

其中，δ为进化因子。由于进化因子δ∈[0,1]，因此惯性权重w∈[0.4,0.9]。进化因子大，表示为初始阶段，大的惯性权重能扩大搜索范围；进化因子小，表示为收敛阶段，小的惯性权重能精确搜索。跟随进化状态而变化的惯性权重能够根据实际情况随时调整，提高了算法的搜索能力与收敛速度。

更新后，当d＝1,2,...,10时，若x_id<0，则令x_id＝0，若则令当d＝11时，若x_i11<500，则令x_i11＝500，若x_i11>700，则令x_i11＝700。

(9)迭代计数累加，T＝T+1.

(10)重复步骤(1)～(9)，直到达到最大迭代次数T_max停止迭代，本实例中T_max＝500。

种群全局最优粒子p_best＝(p_best1,p_best2,...,p_best11)所在的位置即优化后的10个节点的推力方向角和终端时刻。将全局最优粒子的位置信息输入仿真模块，按照式(2)、(3)、(4)获得优化后的月球软着陆轨道，月球软着陆最优轨道规划完成。

仿真模块将最优软着陆轨道通过实时通讯传递给执行模块的执行元件执行，实现燃料消耗最少的月球软着陆最优控制。

上述实施例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。

Claims (1)

1.一种自适应的月球软着陆轨道快速优化控制器，其特征在于：该轨道快速优化控制器由优化模块、仿真模块、执行模块组成；其中：

仿真模块把月球软着陆轨道分割成n个小段，并将n的值输入优化模块；每个节点的时刻由式(1)得到：

其中，t_k为第k+1个节点的时刻，k＝0,1,...,n，t₀为初始时间，记t₀＝0；

优化模块把包括初始节点和终节点在内的n+1个节点的推力方向角ψ与软着陆终端时刻t_f作为待优化的参数；初始化种群规模为N_s的粒子群，随机生成维度为n+2的粒子i的初始位置x_i＝(x_i1,x_i2,...,x_i(n+2))和初始速度v_i＝(v_i1,v_i2,...,v_i(n+2))，i＝1,2,...,N_s，并将粒子的位置信息传入仿真模块；定义维度变量d，d＝1,2,...,n+2；当d＝1,2,...,n+1时，x_id代表第d个节点的推力方向角，当d＝n+2时，x_id代表软着陆终端时刻，x_id∈[500,700]，v_id∈[-200,200]；种群规模N_s＝300～600；然后按以下方法进行迭代，初始时迭代计数T＝0：

(1)仿真模块中，月球软着陆过程中推力方向角表示成多项式(2)：

ψ(t)＝λ₀+λ₁t+λ₂t²+λ₃t³ (2)

其中ψ(t)表示着陆轨道t时刻的推力方向角，λ₀,λ₁,λ₂,λ₃为系数；优化模块输入的粒子的位置代表n+1个节点的推力方向角及软着陆终端时刻，按照式(1)得到n+1个节点的对应时刻；采用函数逼近法，利用n+1个节点的推力方向角及其对应时刻，对式(2)进行拟合，可以求得多项式的系数λ₀,λ₁,λ₂,λ₃，进而得到整个着陆轨道各个时刻的推力方向角ψ(t)；

(2)仿真模块存储了月球软着陆时着陆器的质心动力学方程，见式(3)：

式(3)中r为着陆器的月心距离，v为着陆器的径向速度，θ为着陆器极角，ω为着陆器极角角速度，μ为月球引力常数，m为着陆器质量，F为制动发动机推力，I_SP为制动发动机比冲；其中月球引力常数μ为常数，μ＝4902.75km³/s²，制动发动机推力F与制动发动机比冲I_SP与实际使用的发动机有关，也为常数；着陆器初始质量m₀根据实际确定；其他参数在着陆器着陆过程中发生变化；

初始条件为：

其中，r_p和r_a分别为霍曼转移段的近地点半径和远地点半径，r_p＝1753km，r_a＝1838km；

将步骤(1)中拟合得到的推力方向角ψ(t)、式(4)的初始条件以及着陆器初始质量m₀带入动力学方程(3)，所有数据单位统一，获得月球软着陆的轨道，并将获得的轨道信息输入给优化模块；

(3)优化模块中，优化目标为软着陆过程消耗燃料最少，即令式(5)中指标J最大：

同时，为实现软着陆，终端约束条件为：

其中，R为月球半径，R＝1738km；

在适应度函数中考虑约束条件，构造适应度函数fitness：

fitness＝J-α[(r(t_f)-R)²+v²(t_f)+ω²(t_f)] (7)

其中r(t_f)、v(t_f)、ω(t_f)分别表示仿真模块输入的软着陆轨道终端时刻的月心距离、径向速度、极角角速度；α为罚因子，α＝10000；根据仿真模块输入的软着陆轨道，按照式(7)计算适应度函数值；适应度函数值最大的粒子为全局最优粒子，其位置为p_best＝(p_best1,p_best2,...,p_best(n+2))；

(4)在优化模块中对所有粒子进行分群操作，包括以下子步骤：

(4.1)将所有粒子按照适应度函数值大小从大到小排序，选取适应度函数值最大的粒子作为一个子群中心；

(4.2)在剩下的粒子中选取适应度函数值最大的粒子，依次计算该粒子与各个子群中心的欧几里得距离；粒子i与粒子j的欧几里得距离dist(i,j)定义为：

其中，x_i＝(x_i1,x_i2,...,x_i(n+2))代表粒子i的位置，x_j＝(x_j1,x_j2,...,x_j(n+2))代表粒子j的位置，i,j＝1,2,...,N_s；若该粒子与某一个子群中心的欧几里得距离小于半径a，则将该粒子归为该子群中心所在的子群，并不再计算该粒子与剩下的子群中心的欧几里得距离；若该粒子与所有子群中心的距离都大于半径a，则将该粒子置为一个新的子群中心；半径a＝15～25；

(4.3)重复步骤4.2)，直到处理完所有粒子，则分群完成，且每个子群中心为该子群中适应度函数值最大的粒子；

(5)检查每个子群中的粒子数，若一个子群中的粒子数为S，且S＞S_max，S_max为子群允许的最大粒子数，则将适应度最差的(S-S_max)个粒子的位置和速度进行重置；S_max＝8～10；

(6)在一个子群中，fitness_1为该子群中最大的适应度函数值，fitness_2为该子群中第二大的适应度函数值，tol为搜索精度，tol＝0.01，若满足式(9)：

|fitness_1-fitness_2|＜tol (9)

则将第二大的适应度函数值的粒子q的位置按照式(10)处理：

x'_qd＝x_qd+η·rand (10)

其中，x_qd为粒子q原本的第d维位置，η为扰动因子，rand为0到1之间的随机数，x'_qd为粒子q扰动以后的第d维位置；扰动因子η的大小为：

η＝0.05(x_max-x_min) (11)

其中，x_max，x_min为粒子的搜索上下限，x_min＝0；

(7)确定种群的进化状态；首先，定义每个粒子与其所在子群的子群中心的距离的绝对值之和d_g：

其中，p_ig＝(p_ig1,p_ig2,...,p_ig(n+2))为粒子i所在子群的子群中心的位置；其次，定义每个粒子与其所在子群的子群中心的距离之和的绝对值D_g：

定义进化因子δ为：

由定义可知进化因子δ∈[0,1]；

(8)按照式(15)(16)更新每个粒子的速度与位置：

v_id(T)＝w·v_id(T-1)+c₁·rand·(p_id-x_id(T-1))+c₂·rand·(p_igd-x_id(T-1)) (15)

x_id(T)＝x_id(T-1)+v_id(T) (16)

其中，w为惯性权重；c₁,c₂为加速因子，c₁＝c₂＝2；rand为0到1之间的随机数；p_i＝(p_i1，p_i2,...,p_i(n+2))为粒子x_i的历史最优位置；惯性权重按照式(17)变换：

更新后，当d＝1,2,...,n+1时，若x_id＜0，则令x_id＝0，若则令当d＝n+2时，若x_id＜500，则令x_id＝500，若x_id＞700，则令x_id＝700；

(9)迭代计数累加，T＝T+1；

(10)重复步骤(1)～(9)，直到达到最大迭代次数T_max停止迭代，T_max＝100～2000；

种群全局最优粒子p_best＝(p_best1,p_best2,...,p_best(n+2))所在的位置即优化后的n+1个节点的推力方向角和终端时刻；将全局最优粒子的位置信息输入仿真模块，按照式(2)、(3)、(4)获得优化后的月球软着陆轨道，月球软着陆最优轨道规划完成；

仿真模块将最优软着陆轨道通过实时通讯传递给执行模块的执行元件执行，实现燃料消耗最少的月球软着陆最优控制。