一种智能的月球软着陆轨道控制器的制作方法

本发明涉及航空航天领域，具体地，涉及一种智能的月球软着陆轨道控制器。

背景技术：

我国的探月工程已于2004年启动，根据开展月球探测工程的基本原则，月球探测工程将分为“绕”、“落”、“回”三个阶段实施。在“落”和“回”这两个阶段里，在月球表面实施软着陆是一项关键技术，也是进行月球勘探的重要前提。由于月球表面没有大气，着陆器的速度必须完全由制动发动机抵消，才能安全实现软着陆，这一过程需要消耗大量燃料。所以有必要对这一过程进行优化设计。

求解月球软着陆轨道问题的方法，可分间接法和直接法两类。间接法比较少，主要有由pontryagin极大值原理得出最优着陆轨道。直接法比较多，有将常推力月球软着陆轨道离散化，利用离散点处状态连续作为约束条件，把常推力月球软着陆轨道优化问题归结为一个非线性规划问题；有使用伪光谱方法将软着陆轨道优化问题转为一个约束参数优化问题等。无论使用何种方法，月球软着陆轨道控制问题最后都将转化为一个优化问题，而优化方法的选择对最终的结果影响十分大。目前的优化方法并不能保证一定能够得到最优解，因此月球软着陆轨道控制问题还有很大的探索空间。

技术实现要素：

为了减少月球软着陆时消耗的燃料，本发明的目的在于提供一种智能的月球软着陆轨道控制器。

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

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

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

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

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

ψ(t)＝λ0+λ1t+λ2t²+λ3t³(2)

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

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

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

初始条件为：

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

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

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

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

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

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

fitness＝j-α[(r(tf)-r)²+v²(tf)+ω²(tf)](7)

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

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

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

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

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

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

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

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

|fitness_1-fitness_2|<tol(9)

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

x'qd＝xqd+η·rand(10)

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

η＝0.05(xmax-xmin)(11)

其中，xmax，xmin为粒子的搜索上下限，xmin＝0；

(7)按照式(12)(13)更新每个粒子的速度与位置：

vid(t)＝vid(t-1)+c1·rand·(pid-xid(t-1))+c2·rand·(pigd-xid(t-1))(12)

xid(t)＝xid(t-1)+vid(t)(13)

其中加速因子c1＝c2＝2，rand为0到1之间的随机数，pi＝(pi1，pi2,...,pi(n+2))为粒子xi的历史最优位置，pig＝(pig1，pig2,...,pig(n+2))为粒子i所在的子群的最优粒子的位置；更新后，当d＝1,2,...,n+1时，若xid<0，则令xid＝0，若则令当d＝n+2时，若xid<500，则令xid＝500，若xid>700，则令xid＝700；

(8)迭代计数累加，t＝t+1；

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

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

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

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

附图说明

图1是轨道离散化图；

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

图3是本发明的结构图；

图4是本发明的流程图。

具体实施方式

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

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

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

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

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

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

ψ(t)＝λ0+λ1t+λ2t²+λ3t³(2)

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

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

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

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

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

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

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

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

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

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

fitness＝j-α[(r(tf)-r)²+v²(tf)+ω²(tf)](7)

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

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

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

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

其中，xi＝(xi1,xi2,...,xi11)代表粒子i的位置，xj＝(xj1,xj2,...,xj11)代表粒子j的位置，i,j＝1,2,...,ns。若该粒子与某一个子群中心的欧几里得距离小于半径r，则将该粒子归为该子群中心所在的子群，并不再计算该粒子与剩下的子群中心的欧几里得距离；若该粒子与所有子群中心的距离都大于半径r，则将该粒子置为一个新的子群中心。根据搜索空间的大小，本实例中半径r＝20。

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

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

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

|fitness_1-fitness_2|<tol(9)

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

x'qd＝xqd+η·rand(10)

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

η＝0.05(xmax-xmin)(11)

其中，xmax，xmin为粒子的搜索上下限，xmin＝0。

(7)按照式(12)(13)更新每个粒子的速度与位置：

vid(t)＝vid(t-1)+c1·rand·(pid-xid(t-1))+c2·rand·(pigd-xid(t-1))(12)

xid(t)＝xid(t-1)+vid(t)(13)

其中加速因子c1＝c2＝2，rand为0到1之间的随机数，pi＝(pi1，pi2,...,pi11)为粒子xi的历史最优位置，pig＝(pig1，pig2,...,pig11)为粒子i所在的子群的最优粒子的位置。更新后，当d＝1,2,...,10时，若xid<0，则令xid＝0，若则令当d＝11时，若xi11<500，则令xi11＝500，若xi11>700，则令xi11＝700。

(8)迭代计数累加，t＝t+1。

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

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

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

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