数学建模获奖论文嫦娥三号软着陆轨道设计与控制策略

1、2014高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了全国大学生数学建模竞赛章程和全国大学生数学建模竞赛参赛规则（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛

2、规则的行为，我们将受到严肃处理。我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的报名参赛队号为（8位数字组成的编号）： 所属学校（请填写完整的全名）： 华南理工大学 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： （论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。） 日期：2014 年 9 月

3、15 日赛区评阅编号（由赛区组委会评阅前进行编号）：2014高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：嫦娥三号软着陆轨道设计与控制策略摘 要嫦娥三号的核心任务是实施高可靠高安全的月面软着陆, 要求着陆器必需具备自主障碍识别与规避能力。本文将通过研究月球软着陆主减速段、快速调整段、避障段和缓速下降段的飞行动力学模型，在动力学方程的基础上对各阶段制导律进行优化设计后，结合工程实际确立嫦娥三号软着陆轨道的下降轨迹和

4、控制策略。针对第一问，本文基于开普勒定律求解着陆准备轨道的运动方程，得到近月点与远月点的速度。由于着陆点的位置是由近月点的位置和软着陆轨迹共同决定的。因此本文通过求解嫦娥三号在主减速段和粗避障段的经过的横向距离，利用着陆点的位置近似确定近月点的位置。实验结果表明，此近似求解的近月点位置误差在10%以内。同时，本文所建立的基于飞行器动力学模型的主减速段燃料最优制导律，通过蚁群算法的参数优化，也取得误差5%以内的良好结果。本文还将证明主减速段的燃料最优制导律等效于时间最优制导律。针对第二问，本文通过查阅相关资料，给出了软着陆各个阶段的动力学方程和制导律。主减速段是耗时最长，燃料消耗最多的阶段，因此

5、采用了燃料最优制导律。快速调整阶段是过渡阶段，拟采用推力线性变化制导律。粗避障段的主要目的是避开大的陨石坑，由于着陆器此时接近垂直下降，因而采用最优开关制导律。为实现粗避障，本文采用对光学图像进行分格，并且判断各格为障碍区（0）或非障碍区（1），最终通过在0-1图上寻找具有安全半径的连片非障碍区，找到合适的着陆区域，取得了良好的实验效果。精避障段则要构建三维高程图和平均坡面，计算平均坡度，估算出障碍高度，最终找出安全着陆区。缓速下降段主要考虑的是着陆的安全性，本文将利用变推力制导律和垂直软着陆模型进行分析。针对第三问，本文将在模型的分析与检验中给出相应的误差分析，实验结果表明，本文所提出的模型

6、与嫦娥三号实际结果相符合，具有实用性和准确性。关键词：制导律 蚁群算法 数字高程图 嫦娥三号 软着陆1.问题的重述1.1 背景实施在月球表面的软着陆是月球勘探计划的重要一步。所谓月球软着陆,是指着陆器在制动系统作用下以很小的速度准确降落到月面指定区域,以保证试验设备和宇航员的安全。由于月球表面没有大气，因此整个软着陆过程需要制导发动机的控制。正确地建立着陆器的飞行动力学模型，并在此基础对着陆器的下降轨迹和制导律进行优化设计，对于着陆器飞行程序的设计和燃料资源消耗等方面具有重要的指导意义。1.2 需要解决的问题根据题目信息和附件所示内容，分析嫦娥三号软着陆过程中的6个阶段，研究以下的问题：（1）

7、根据附件1与附件2研究嫦娥三号的着陆准备轨道，确定该椭圆轨道上近月点与远月点的位置，以及嫦娥三号相应的速度大小与方向。（2）建立合理的着陆器飞行动力学方程，根据软着陆过程中6个阶段不同的特点与要求，结合附件1与附件2的内容，确立嫦娥三号的下降轨迹和各个阶段的最优控制策略。（3）根据问题一和问题二建立的数学模型，对模型进行相应的误差分析和敏感性分析。2.模型假设（1）假设诸如月球引力非球项、日月引力摄动等影响因素均可忽略不计（2）月球表面附近没有大气，所以在飞行器的动力学模型中没有大气阻力项（3）嫦娥三号在软着陆全过程不发生意外碰撞和震荡3.符号说明符号符号说明嫦娥三号进行霍曼转移的初始时刻嫦娥

8、三号到达近月点的时刻嫦娥三号在主减速段的终止时刻着陆器的月心距月球半径远月点与月面的距离着陆器在近月点的月心距近月点与月面的月心距月球的引力常数在时刻的横向速度在时刻的径向速度在时刻的横向速度在时刻的径向速度横向速度径向速度着陆器的月心距着陆器的极角推力加速度发动机的推力大小着陆器的初始质量推力方向角的系数，待优化参量燃料消耗率燃料消耗质量比冲终止时刻合速度城市和的距离t时刻蚂蚁k由城市转移到城市的概率信息量表示蚂蚁k下一步允许走过的城市集合信息量对蚂蚁的影响力大小由城市i转移到城市j的期望程度终止时刻，待优化的时间参量的作用信息增量接近段的水平距离主减速段的水平距离近月点与着陆点的距离粗避障

9、段的初始时刻粗避障段的终止时刻月面高度下降速度与垂直方向的夹角，称为下降角月球的重力加速度制动力的开关量等效推力加速度的斜率近月点与着陆点的水平距离相对误差近月点速度相对误差主减速仰幅角的相对误差主减速段时间相对误差单格内障碍点数单格内障碍点阈值过暗灰度临界值过亮灰度临界值注：未列出的符号及重复的符号以出现处为准。4.问题一的建模与求解4.1 问题分析根据题目信息，着陆准备轨道为近月点15km，远月点100km的椭圆轨道，根据开普勒定律，可以建立椭圆轨道上的运动方程，求解出近月点与远月点的速度的大小与方向。着陆点的位置是由近月点在月心坐标系中的位置和软着陆形态共同决定的。因此可以通过求解软着陆

10、的下降轨迹，结合着陆点的位置，确定近月点在月心坐标系中的切向距离。此外，利用月球的半径参数，可以求得近月点的径向距离，进而确定近月点的位置。相应地，远月点的位置可以结合近月点的位置，利用椭圆的轨迹方程进行确定。在本问题中，软着陆的下降轨迹的求解是一个难点，考虑到着陆器在近月点处离月面距离较高，因此可以在不引起较大的误差的情况下，建立简单的模型进行求解。根据相关参考资料1有：在接近段中，着陆器采用接近直线下降方式逐步接近着陆区，利用附件1与2提供的接近段的高度信息，可以求解出对应的水平距离。而主减速段的轨迹可以建立飞行动力学模型进行求解。考虑到主减速段是软着陆过程用时最长、推进剂消耗最多的任务段

11、。该段的主要任务是消除较大的动力下降段初始水平速度(约 1.7 km/s), 因此推进剂消耗优化是该段制导律的主要设计目标。因此，求解基于动力学模型的燃料消耗最低制导律方程，可以确定主减速段的水平距离。综合上述求解出的距离信息，就可以在一定的误差允许范围内给出近月点的位置，进而得到远月点的位置。总体而言，问题一的流程图如下所示：图4-1问题一的流程图4.2 模型的建立4.2.1 着陆准备轨道的运动学方程假设嫦娥三号在时刻开始进行霍曼转移，在时刻到达近月点，则在惯性坐标系中，以月心为原点的极坐标形式的着陆准备轨道的运动学方程2为：在时刻：（4-1）其中的为着陆器的月心距，为月球半径，可取值为月球

12、的平均半径1737.013km，为远月点与月面的距离，取值为100km，为着陆器在近月点的月心距，为近月点与月面的距离，取值为15km，为月球的引力常数，其值为：，为初始横向速度，即着陆器从远月点开始运动的速度，为初始径向速度，其值为0.在时刻：（4-2）其中的为初始横向速度，即着陆器从近月点开始运动的速度，为初始径向速度，其值为0.4.2.2 主减速段的动力学方程在惯性坐标系中，以月心为原点的极坐标形式受控飞行器动力学方程为2：（4-3）其中为月球引力常数，为径向速度，为横向速度，为着陆器的月心距，为着陆器的极角。其中为推力加速度，其值可以由下式计算：（4-4）其中为发动机推力的大小，查阅相

13、关资料1可得：主减速段的主要任务是消除较大的动力下降段初始水平速度，所以发动机一直工作在最大推力模式即7500N。为着陆器的初始质量，其值为2400kg，为比冲（即单位质量的推进剂产生的推力），其值为2940m/s。为推力方向角，即推力方向与当地水平线的夹角，查阅相关资料3可得：（4-5）其中,为待优化的参量。4.2.3 基于蚁群算法3的主减速段燃料最优制导律4.2.3.1 燃料最优制导律设主减速段的初始时刻为，终止时刻为，对于推力幅值恒定的飞行器，燃料消耗最省的性能指标为：（4-6）其中为燃料消耗率，为了使消耗的燃料质量最小，由于初始时刻，发动机推力和比冲均为定值，因此只能使主减速段的终止时

14、刻最小。所以，对于主减速段而言，燃料最优制导律等同于时间最优制导律。查阅相关资料1可得，主减速末期的俯仰角为，因此可列出如下等式：（4-7）其中为终止时刻的合速度，为57m/s.因此，只要确定了的函数表达式或函数图像，就可以用解析法或作图法求得。而由式（4-3）所确定，式（43）的唯一参量是，又由参量,所确定，因此优化的问题就转换为优化参量,的问题。根据式(4-5)，参量,的优化可以用多个离散值和值来多阶拟合求得。其中： （4-8）用函数逼近法进行参数化的相关参数设置为：将轨迹离散化化为9段，一共10个离散点，相应时刻的值可以由来确定： （4-9）因此待优化参数共11个，即10个推力方向角和1

15、个终端时刻。对于终端时刻，根据齐奥尔科夫斯基公式和软着陆初始条件，可由下式估计：（4-10）式中和分别表示着陆器的终端速度和初始速度，经计算确定搜索范围为(单位秒)：（4-11）我们将用蚁群算法来找到这11个参数。4.2.3.2蚁群算法3设有n个城市,(i,j=1,2,n)表示城市i和j间的距离, 表示在t时刻城市i和j之间的信息量。我们以此来模拟实际蚂蚁的分泌物。设共有m只蚂蚁,用表示在t时刻蚂蚁k由城市i转移到城市j的概率,则(4-12)其中,表示蚂蚁k下一步允许走过的城市的集合,表示路径上的信息量对蚂蚁选择路径所起的作用大小,为由城市i转移到城市j的期望程度,例如,可以取。表示的作用.当

16、时,算法就是传统的贪心算法;而当时,就成了纯粹的正反馈的启发式算法.经过n个时刻,蚂蚁可走完所有的城市,完成一次循环.每只蚂蚁所走过的路径就是一个解.此时,要根据下面公式对各路径上的信息量作更新:(4-13)其中,表示信息量随时间的推移而衰减的程度.信息增量可表示为(4-14)表示蚂蚁k在本次循环中在城市i和j之间留下的信息量。下面给出蚁群算法应用过程的具体步骤。(1)设置初始参数，包括蚂蚁数N，循环次数K，挥发系数Rho，调节系数Lambda、Q，所有路径信息素量初值Tau。，蚂蚁初始位置等。(2)根据式(4-8)计算每只蚂蚁的转移概率，然后依据赌轮原则为每只蚂蚁选择下一个路径。重复上述操作

17、直至所有蚂蚁均完成一次循环。(3)将每只蚂蚁的路径解码为优化参数值和和1个终端时刻，用这10个数据拟合求出，从而得到参量的表达式，根据式(4-6)计算目标函数值求得,对比当前最优值，判断找出当前最好的路径。然后根据式(4-9)、(4-10)更新相关路径上的信息素。(4)判断是否满足终止条件，不满足则重复(1)至(3)步过程；满足则结束计算输出结果。4.2.3.2 本节模型综述综上所述，4.2.3节的模型可以总结为下式所述的规划问题：目标函数： 约束条件：（4-15）4.2.4 远月点的位置嫦娥三号的运行轨道近似一个椭圆，月球则是该椭圆的焦点之一，若设椭圆半长轴为a，半短轴为b，焦点为c。则焦点

18、（即月球圆心）距离椭圆轨道（即嫦娥三号的绕月轨道）长轴两端的距离分别为：。也就是远地点和近地点的高度（未减去月球半径）。根据椭圆方程：（4-16）可以求得近月点和远月点的位置关系。4.3 模型的求解4.3.1 近月点和远月点速度利用式（4-1）可求得远月点的速度为：利用式（4-2）可求得近月点的速度为：4.3.2 近月点的位置4.3.2.1 接近段的水平距离查阅相关资料1可得，处于接近段的着陆器采用接近直线下降方式逐步接近着陆区，因此接近段的水平距离可以近似等于接近段的垂直距离。4.3.2.2 主减速段的水平距离距离可以根据主减速段最终时刻的最优参量，以及的表达式积分求得，如下所示：（4-17

19、）根据蚁群算法的优化的优化结果，求得:从而得到的表达式为：（4-18）运行结果图如下：注：时间大于交点处时的线段不在主减速阶段内，故无效图4-2 的函数图像图4-3 的函数图像由图可得，4.3.2.3 近月点与着陆点的距离在一定的误差允许范围内，可以将近月点与着陆点的距离做如下近似计算：4.3.2.4 近月点的速度方向速度方向为，即主减速段的俯仰角为5.问题二的建模与求解5.1 问题分析嫦娥三号软着陆过程大致可以划分为6个阶段，着陆准备阶段、主减速段、快速调整段、粗避障段、精避障段、缓速下降段。问题一中已经确立了着陆准备阶段和主减速段的运动方程和动力学模型，主减速段的主要目的是通过推力发动机抵

20、消近月点处较大的初始速度，所以燃料消耗最优是这一阶段的主要优化目标。问题一中已给出了相应的数学模型。对于快速调整段，其介于主减速段和粗避障段之间，主要目的是缓冲作用，即使着陆器平缓地从主减速段过渡到接近段。查阅相关资料1，为了满足主发动机推力和着陆器姿态匀速过渡的要求,我们利用推力大小和方向线性变化的制导律, 其中制导参数利用实际的主减速段末端状态和接近段初始状态约束确定。对于粗避障段，主要任务是要求避开大的陨石坑，实现在设计着陆点上方100m处悬停，并初步确定落月地点。查阅相关资料4,为了满足制导目标的位置、速度、姿态以及初始高度和速度等多项约束，我们构建平面月球模型，并利用最优开关制导律模

21、拟。在粗避障过程中，利用月球岩石坑的图像特征，采用基于光学图像的成像，通过对图像进行分析和处理，判断出过亮和过暗障碍，识别出当前下降区域上的大障碍，从而确定安全区域。对于精避障段，主要目的是在粗避障选取的较安全区域内进行精确的障碍检测。由于精避障时已极接近月球表面，务必识别并剔除危及安全的小尺度障碍，确保落点安全。为此需要进行高分辨率的三维成像，查阅相关资料，运用基于三维图像的精避障识别和安全着陆区选取算法处理三维图像数据。对于缓速下降段，主要考虑的是着陆的安全性，查阅相关资料5,我们利用变推力制导律和垂直软着陆模型进行分析。问题二的流程图如下所示：图5-1问题二的流程图5.2 模型的建立5.

22、2.1 快速调整阶段为了满足主发动机推力和着陆器姿态匀速过渡的要求，建立快速调整阶段的推力模型如下所示：（5-1）其中为推力的大小，为推力方向角，即推力方向与当地水平线的夹角。假设快速调整阶段的初始时刻为，终止时刻为，本阶段的动力学方程可以将式（5-1）、式（4-3）与式（4-4）联立，结合下述初值条件：（5-2）来构造。其中,与之前的定义相同。5.2.2 粗避障段5.2.2.1最优开关制导律假设粗避障段的初始时刻为，终止时刻为，基于平面月球模型6的运动方程如下所示：（5-3）其中为月面高度，为下降速度与垂直方向的夹角,称为下降角，为月球的重力加速度，为制动力的开关量，其定义如下所示：（5-4

23、）式（5-4）表明粗避障段的着陆器将先做自由落体运动，然后依靠发动机制动下降到接近段终端。为了简化过程，在一定的误差允许范围内，引入以下两个简化条件：（1）简化条件1：由于着陆器在本阶段接近垂直下降，因此可以将视为小量，因而有：（5-5）（2）简化条件2：根据相关资料1，本阶段历时时间较短，故引入质量不变假设：（5-6）其中下标表示粗避障段初始状态，下标表示粗避障段终止状态。最终，我们得到了简化后的运动方程如下所示：（5-7）5.2.2.2 基于光学图像的粗障碍检测算法通过观察，月球表面的岩石坑图像特征如下：1）岩坑表面石头具有高亮度；2）岩坑存在明显的阴影区3）亮区与阴影之间存在强对比度。这

24、些亮暗区都是我们剔除的目标。（1）光学图像的分格与障碍区的识别将光学图像分格为个等大方格，A题提供的TIF图为分辨率，将其灰度值分格为，每格像素为。然后对每格的灰度值进行统计。每格共有个灰度，而灰度过大则表示过亮，即岩坑表面石头；过小表示岩坑，为暗区。当一个格的灰度过大和过小的点过多，这个格将被认定为障碍区，不能降落。因此，需要对灰度过大和过小的点进行统计，统计方法如下：设障碍点数为初值为0，障碍点数阈值为，过暗灰度临界值为,过亮灰度临界值为，对于一格内个灰度点，若：则表示该灰度点是障碍点，统计值加一。当一格统计完成时，若：即障碍点数超过阈值，此时判断该格为障碍区。为方便起见，标记障碍区所在格

25、为0，非障碍区为1。（2）采用螺旋搜索算法寻找安全着陆区通过螺旋搜索算法（图5-2），从上步得到的的表格中间开始搜索非障碍区。当找到非障碍区后，为其设置一个安全半径。继续搜索，若在安全半径内无障碍区，则将该半径内区域设置为安全着陆区。图5-2螺旋搜索算法示意图75.2.3 精避障段5.2.3.1基于三维图像的精避障识别和安全着陆区选取7设计的精障碍识别和安全着陆区选取算法为：1）数据预处理，构建测量坐标系下的地形三维高程图；2）平均坡面构建，采用最小二乘法拟合一定单元区域的平均坡面；3）平均坡度计算，根据平均坡面计算该区域的平均坡度；4）障碍高度计算，根据平均坡面估算该区域内每个单元格的障碍高

26、度；5）安全着陆区选取，采用从着陆区中心开始顺时针陆炫前进搜索的方法（图），直至找到符合着陆要求的区域位置。图5-3中心螺旋式安全着陆区搜索方法45.2.4 缓速下降段5.2.4.1 垂直软着陆模型6接近段结束后,着陆器距离月面很近,且下降速度与当地水平面几乎垂直。采用分段控制方式, 着陆器依次经过悬停、匀加速、匀减速和关机降落几个过程,各段加速度分别如下所示：（5-8）5.2.4.2 等效连续变推力制导律6假设等效推力的加速度是线性变化的，取向上为正方向，有： （5-9）其中,表示等效推力加速度的斜率,表示推力逐渐增大,表示推力逐渐减小；表示等效推力加速度相对于月心重力加速度；的初始偏移量,

27、 表示初始推力加速度大于,表示初始推力加速度小于。结合关机条件：（5-10）对式（5-9）等式两边连续积分两次，得：（5-11）其中，下标和分别表示着陆段初始和关机时刻状态,下标表示着陆段,下标表示关机时刻。5.3 模型的求解5.3.1 粗避障光学图像的分格与障碍区的识别通过MATLAB的imread函数，读取粗避障光学图像的灰度图。考虑到该灰度图有个像素，用MATLAB进行障碍区的判断会相当慢，于是将该灰度图录入C+下进行模型的求解。用C+进行分格和识别障碍区时，经过多次测试设定各值：阈值 过暗灰度临界值 过亮灰度临界值 得到障碍区的0-1图如下，通过与原TIF图比较可验证其正确性:图5-4

28、障碍区0-1图与原tif图的比较 显然，0区对应原图的岩石坑区，1区对应原图较平坦的区域，1区才是适合着陆的安全区。计算结果和实际TIF图表示的地形非常吻合。圆圈所包括的区域，是最后得出粗避障过程后适合安全降落的区域。5.3.2 精避障三维高程图图5-5精避障三维高程图6.模型的分析与检验6.1着陆准备段与主减速段模型的分析与检验根据相关资料5, 着陆器动力下降初始状态为: 经度 19.0464°W, 纬度 28.9989°N,速度 1695.7 m/s, 俯仰姿态角 85°，主减速阶段经历时间为487s。根据经纬度距离计算公式，计算出实际近月点与着陆点的水平距离

29、为：解得：所以计算结果误差为：其中，依次为近月点与着陆点水平距离相对误差、近月点速度相对误差、主减速段仰幅角相对误差、主减速段时间相对误差。结论：由上述数据可知，与后三项的值相比，距离误差较大。从近月点与着陆点距离计算过程中可以知道，为了简化计算，我们只考虑了主减速段和粗避障段产生的水平距离，而忽略了其他阶段产生的水平距离，因而误差相对较大，但是仍然控制在10%以内，说明在一定误差允许范围内，此结果具有一定的参考性。而后三项的相对误差比较小，说明他们都与嫦娥三号实际飞行得到的结果比较符合。有力地论证了着陆准备阶段和主减速段模型的准确性与可靠性。6.2粗避障光学图像分格与障碍区识别的分析与检验粗

30、避障的模型通过模型的求解检验了其正确性，在着陆器飞行时，若着陆器能及时计算出安全区域，该模型还是具有一定可行性的。该粗避障模型计算量偏大，对着陆器计算能力无疑是一大考验。该粗避障模型计算结果与实际地形误差较少，若要进一步减少误差，提高准确性，则需提高分格的格数，同时还要再度优化过暗灰度临界值、过亮灰度临界值和阈值。7.模型的创新点与不足之处7.1 模型的创新点（1）将月球软着陆的轨迹优化问题转化为参数优化问题，将求燃料最优问题转化为求时间最优问题，使问题的求解变得简单。（2）使用蚁群算法求解复杂的角度变化参数优化问题，具有很好的试验效果。（3）利用光学图像的分格和设定阈值来寻找连片非障碍区作为

31、安全着陆区域。7.2 模型的不足之处蚁群算法比较适合处理离散的优化问题，从而限制了其在连续优化问题中的应用，本模型中将角度变化离散化，存在较大的误差。8.参考资料1张洪华, 关轶峰, 黄翔宇, 等. 嫦娥三号着陆器动力下降的制导导航与控制. 中国科学: 技术科学, 2014, 44: 3773842王劼，李俊峰，崔乃刚，刘暾，登月飞行器软着陆发动机持续工作初始点选择研究，工程力学，1000-4750，(2003)06-0145-043段佳佳，徐世杰，朱建丰，基于蚁群算法的月球软着陆轨迹优化，宇航学报，第29卷第2期：1000-1328，(2008)02.0476.064Klumpp A R.

32、Apollo lunar descent guidance. Automatica, 1974, 10: 1331465张洪华, 梁俊, 黄翔宇, 等. 嫦娥三号自主避障软着陆控制技术. 中国科学: 技术科学, 2014, 44: 5595686 王鹏基，张熇，曲广吉，月球软着陆飞行动力学和制导控制建模与仿真，中国科学 E 辑: 技术科学 2009 年，第39 卷 第3 期: 521 5277 张洪华, 梁俊, 黄翔宇, 等. 嫦娥三号自主避障软着陆控制技术. 中国科学: 技术科学, 2014, 44: 5595689.附录MATLAB源程序function BESTX,BESTY,ALLX,

33、ALLY=ACOUCP(K,N,Rho,Q,Lambda,LB,UB)% ACOUCP.m%在命令行中输入LB = 0 0 0 0 0 0 0 0 0 0 0 ;UB = 1 1 1 1 1 1 1 1 1 1 1;ACOUCP(50,30,0.95,1,0.5,LB,UB)% 输入参数列表% K 迭代次数% N 蚁群规模% Rho 信息素蒸发系数，取值01之间，推荐取值0.70.95% Q 信息素增加强度，大于0，推荐取值1左右% Lambda 蚂蚁爬行速度，取值01之间，推荐取值0.10.5% LB 决策变量的下界，M×1的向量% UB 决策变量的上界，M×1的向量%

34、输出参数列表% BESTX K×1细胞结构，每一个元素是M×1向量，记录每一代的最优蚂蚁% BESTY K×1矩阵，记录每一代的最优蚂蚁的评价函数值% ALLX K×1细胞结构，每一个元素是M×N矩阵，记录每一代蚂蚁的位置% ALLY K×N矩阵，记录每一代蚂蚁的评价函数值% 参考设置% BESTX,BESTY,ALLX,ALLY=ACOUCP(50,30,0.95,1,0.5,LB,UB)% 第一步：初始化M=length(LB);%决策变量的个数%蚁群位置初始化X=zeros(M,N);for i=1:M x=unifrnd(LB

35、(i),UB(i),1,N); %unidrnd()离散分布 X(i,:)=x;end%输出变量初始化ALLX=cell(K,1);%细胞结构，每一个元素是M×N矩阵，记录每一代的个体ALLY=zeros(K,N);%K×N矩阵，记录每一代评价函数值BESTX=cell(K,1);%细胞结构，每一个元素是M×1向量，记录每一代的最优个体BESTY=zeros(K,1);%K×1矩阵，记录每一代的最优个体的评价函数值k=1;%迭代计数器初始化Tau=ones(1,N);%信息素初始化Y=zeros(1,N);%适应值初始化% 第二步：迭代过程while k

36、<=K YY=zeros(1,N); for n=1:N x=X(:,n); YY(n)=FIT(x); end maxYY=max(YY); temppos=find(YY=maxYY); POS=temppos(1); %蚂蚁随机探路 for n=1:N if n=POS x=X(:,n); Fx=FIT(x); mx=GaussMutation(x); Fmx=FIT(mx); if Fmx<Fx X(:,n)=mx; Y(n)=Fmx; elseif rand>1-(1/(sqrt(k) X(:,n)=mx; Y(n)=Fmx; else X(:,n)=x; Y(n)

37、=Fx; end end end for n=1:N if n=POS x=X(:,n); Fx=FIT(x); mx=GaussMutation(x); Fmx=FIT(mx); if Fmx<Fx Y(n)=Fmx; elseif rand>1-(1/(sqrt(k) X(:,n)=mx; Y(n)=Fmx; else X(:,n)=x; Y(n)=Fx; end end end %朝信息素最大的地方移动 for n=1:N if n=POS x=X(:,n); r=(K+k)/(K+K); p=randperm(N); t=ceil(r*N); pos=p(1:t); Tem

38、pTau=Tau(pos); maxTempTau=max(TempTau); pos2=find(TempTau=maxTempTau); pos3=pos(pos2(1); x2=X(:,pos3(1); x3=(1-Lambda)*x+Lambda*x2; Fx=FIT(x); Fx3=FIT(mx); if Fx3<Fx X(:,n)=x3; Y(n)=Fx3; elseif rand>1-(1/(sqrt(k) X(:,n)=x3; Y(n)=Fx3; else X(:,n)=x; Y(n)=Fx; end end end %更新信息素并记录 Tau=Tau*(1-Rho

39、); maxY=max(Y); minY=min(Y); DeltaTau=(maxY-Y)/(maxY-minY); Tau=Tau+Q*DeltaTau; ALLXk=X; ALLY(k,:)=Y; minY=min(Y); pos4=find(Y=minY); BESTXk=X(:,pos4(1); BESTY(k)=minY; disp(k); k=k+1;end%从蚁群算法中得到：tf和十个角度值angle(单位：度)%计算各个时刻function tf1 = angle(M,tf)for i=1:10 td(i) = 0+(i-1)*(tf-0)/9; end%三次性项拟合glob

40、al A;A = polyfit(td,M,3)aL = 1737.013*103; %月球半径 %求解飞行器动力学方程T,Y = ode45('dip',0 750,148.165,aL+15*103,1693.53,0);%A = -6.1321*10(-8) 8.0647*10(-5) 0.0163 5%figure%plot(T,Y(:,3),'b-',0:750,51.66,'r-')%从数值解上选取坐标(tf,51.66) tf1 = interp1(Y(:,3),T,51.66,'nearest','extr

41、ap');%tf1,Vr = ginput;function dy = dip(t,y)u = 0.4902802627*1013; %月球引力常数 global A;%飞行器动力学方程%y1是径向速度，y2是月心距，y3是切向速度,y4是横向距离%Vr,Vo,r,Q=dsolve('DVr = -u/(r2)+(Vo2)/r+7500/(2400-7500*t/2940)*sin(a0+a1*t+a2*t2+a3*t3)','DVo = -Vr*Vo/r+7500/(2400-7500*t/2940)*cos(a0+a1*t+a2*t2+a3*t3)'

42、,'Dr = Vr','DQ = Vo/r','(Vr(0)2+(Vo(0)2 = 2890000','r(0)=aL+15*103','t');dy = zeros(4,1);dy(1) = -u/(y(2)2+(y(3)2/y(2)-7500/(2400-7500*t/2940)*sin(A(4)+A(3)*t+A(2)*t2+A(1)*t3)*pi/180);dy(2) = y(1);dy(3) = -y(1)*y(3)/y(2)-7500/(2400-7500*t/2940)*cos(A(4)+A(3)*t+A(2)*t2+A(1)*t3)*pi/180);dy(