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

嫦娥三号软着陆轨道设计与控制策略摘要嫦娥三号的核心任务是实施高可靠高安全的月面软着陆,要求着陆器必需具备自主障碍识别与规避能力。本文将通过研究月球软着陆主减速段、快速调整段、避障段和缓速下降段的飞行动力学模型，在动力学方程的基础上对各阶段制导律进行优化设计后，结合工程实际确立嫦娥三号软着陆轨道的下降轨迹和控制策略。针对第一问，本文基于开普勒定律求解着陆准备轨道的运动方程，得到近月点与远月点的速度。由于着陆点的位置是由近月点的位置和软着陆轨迹共同决定的。因此本文通过求解嫦娥三号在主减速段和粗避障段的经过的横向距离，利用着陆点的位置近似确定近月点的位置。实验结果表明，此近似求解的近月点位置误差在10%以内。同时，本文所建立的基于飞行器动力学模型的主减速段燃料最优制导律，通过蚁群算法的参数优化，也取得误差5%以内的良好结果。本文还将证明主减速段的燃料最优制导律等效于时间最优制导律。针对第二问，本文通过查阅相关资料，给出了软着陆各个阶段的动力学方程和制导律。主减速段是耗时最长，燃料消耗最多的阶段，因此采用了燃料最优制导律。快速调整阶段是过渡阶段，拟采用推力线性变化制导律。粗避障段的主要目的是避开大的陨石坑，由于着陆器此时接近垂直下降，因而采用最优开关制导律。为实现粗避障，本文采用对光学图像进行分格，并且判断各格为障碍区（0）或非障碍区（1），最终通过在0-1图上寻找具有安全半径的连片非障碍区，找到合适的着陆区域，取得了良好的实验效果。精避障段则要构建三维高程图和平均坡面，计算平均坡度，估算出障碍高度，最终找出安全着陆区。缓速下降段主要考虑的是着陆的安全性，本文将利用变推力制导律和垂直软着陆模型进行分析。针对第三问，本文将在模型的分析与检验中给出相应的误差分析，实验结果表明，本文所提出的模型与嫦娥三号实际结果相符合，具有实用性和准确性。关键词：制导律蚁群算法数字高程图嫦娥三号软着陆1.问题的重述1.1背景实施在月球表面的软着陆是月球勘探计划的重要一步。所谓月球软着陆,是指着陆器在制动系统作用下以很小的速度准确降落到月面指定区域,以保证试验设备和宇航员的安全。由于月球表面没有大气，因此整个软着陆过程需要制导发动机的控制。正确地建立着陆器的飞行动力学模型，并在此基础对着陆器的下降轨迹和制导律进行优化设计，对于着陆器飞行程序的设计和燃料资源消耗等方面具有重要的指导意义。1.2需要解决的问题根据题目信息和附件所示内容，分析嫦娥三号软着陆过程中的6个阶段，研究以下的问题：（1）根据附件1与附件2研究嫦娥三号的着陆准备轨道，确定该椭圆轨道上近月点与远月点的位置，以及嫦娥三号相应的速度大小与方向。（2）建立合理的着陆器飞行动力学方程，根据软着陆过程中6个阶段不同的特点与要求，结合附件1与附件2的内容，确立嫦娥三号的下降轨迹和各个阶段的最优控制策略。（3）根据问题一和问题二建立的数学模型，对模型进行相应的误差分析和敏感性分析。2.模型假设（1）假设诸如月球引力非球项、日月引力摄动等影响因素均可忽略不计（2）月球表面附近没有大气，所以在飞行器的动力学模型中没有大气阻力项（3）嫦娥三号在软着陆全过程不发生意外碰撞和震荡3.符号说明符号符号说明嫦娥三号进行霍曼转移的初始时刻嫦娥三号到达近月点的时刻嫦娥三号在主减速段的终止时刻着陆器的月心距月球半径远月点与月面的距离着陆器在近月点的月心距近月点与月面的月心距月球的引力常数在时刻的横向速度在时刻的径向速度在时刻的横向速度在时刻的径向速度横向速度径向速度着陆器的月心距着陆器的极角推力加速度发动机的推力大小着陆器的初始质量推力方向角的系数，待优化参量燃料消耗率燃料消耗质量比冲终止时刻合速度城市和的距离t时刻蚂蚁k由城市转移到城市的概率信息量表示蚂蚁k下一步允许走过的城市集合信息量对蚂蚁的影响力大小由城市i转移到城市j的期望程度终止时刻，待优化的时间参量的作用信息增量接近段的水平距离主减速段的水平距离近月点与着陆点的距离粗避障段的初始时刻粗避障段的终止时刻月面高度下降速度与垂直方向的夹角，称为下降角月球的重力加速度制动力的开关量等效推力加速度的斜率近月点与着陆点的水平距离相对误差近月点速度相对误差主减速仰幅角的相对误差主减速段时间相对误差单格内障碍点数单格内障碍点阈值过暗灰度临界值过亮灰度临界值注：未列出的符号及重复的符号以出现处为准。4.问题一的建模与求解4.1问题分析根据题目信息，着陆准备轨道为近月点15km，远月点100km的椭圆轨道，根据开普勒定律，可以建立椭圆轨道上的运动方程，求解出近月点与远月点的速度的大小与方向。着陆点的位置是由近月点在月心坐标系中的位置和软着陆形态共同决定的。因此可以通过求解软着陆的下降轨迹，结合着陆点的位置，确定近月点在月心坐标系中的切向距离。此外，利用月球的半径参数，可以求得近月点的径向距离，进而确定近月点的位置。相应地，远月点的位置可以结合近月点的位置，利用椭圆的轨迹方程进行确定。在本问题中，软着陆的下降轨迹的求解是一个难点，考虑到着陆器在近月点处离月面距离较高，因此可以在不引起较大的误差的情况下，建立简单的模型进行求解。根据相关参考资料[1]有：在接近段中，着陆器采用接近直线下降方式逐步接近着陆区，利用附件1与2提供的接近段的高度信息，可以求解出对应的水平距离。而主减速段的轨迹可以建立飞行动力学模型进行求解。考虑到主减速段是软着陆过程用时最长、推进剂消耗最多的任务段。该段的主要任务是消除较大的动力下降段初始水平速度(约1.7km/s),因此推进剂消耗优化是该段制导律的主要设计目标。因此，求解基于动力学模型的燃料消耗最低制导律方程，可以确定主减速段的水平距离。综合上述求解出的距离信息，就可以在一定的误差允许范围内给出近月点的位置，进而得到远月点的位置。总体而言，问题一的流程图如下所示：图4-1 问题一的流程图4.2模型的建立4.2.1着陆准备轨道的运动学方程假设嫦娥三号在时刻开始进行霍曼转移，在时刻到达近月点，则在惯性坐标系中，以月心为原点的极坐标形式的着陆准备轨道的运动学方程[2]为：在时刻： （4-1）其中的为着陆器的月心距，为月球半径，可取值为月球的平均半径1737.013km，为远月点与月面的距离，取值为100km，为着陆器在近月点的月心距，为近月点与月面的距离，取值为15km，为月球的引力常数，其值为：，为初始横向速度，即着陆器从远月点开始运动的速度，为初始径向速度，其值为0.在时刻： （4-2）其中的为初始横向速度，即着陆器从近月点开始运动的速度，为初始径向速度，其值为主减速段的动力学方程在惯性坐标系中，以月心为原点的极坐标形式受控飞行器动力学方程为[2]： （4-3）其中为月球引力常数，为径向速度，为横向速度，为着陆器的月心距，为着陆器的极角。其中为推力加速度，其值可以由下式计算： （4-4）其中为发动机推力的大小，查阅相关资料[1]可得：主减速段的主要任务是消除较大的动力下降段初始水平速度，所以发动机一直工作在最大推力模式即7500N。为着陆器的初始质量，其值为2400kg，为比冲（即单位质量的推进剂产生的推力），其值为2940m/s。为推力方向角，即推力方向与当地水平线的夹角，查阅相关资料[3]可得： （4-5）其中,,,为待优化的参量。4.2.3基于蚁群算法[3]的主减速段燃料最优制导律燃料最优制导律设主减速段的初始时刻为，终止时刻为，对于推力幅值恒定的飞行器，燃料消耗最省的性能指标为： （4-6）其中为燃料消耗率，为了使消耗的燃料质量最小，由于初始时刻，发动机推力和比冲均为定值，因此只能使主减速段的终止时刻最小。所以，对于主减速段而言，燃料最优制导律等同于时间最优制导律。查阅相关资料[1]可得，主减速末期的俯仰角为，因此可列出如下等式： （4-7）其中为终止时刻的合速度，为57m/s.因此，只要确定了的函数表达式或函数图像，就可以用解析法或作图法求得。而由式（4-3）所确定，式（4—3）的唯一参量是，又由参量,,,所确定，因此优化的问题就转换为优化参量,,,的问题。根据式(4-5)，参量,,,的优化可以用多个离散值和值来多阶拟合求得。其中： （4-8）用函数逼近法进行参数化的相关参数设置为：将轨迹离散化化为9段，一共10个离散点，相应时刻的值可以由来确定： （4-9）因此待优化参数共11个，即10个推力方向角和1个终端时刻。对于终端时刻，根据齐奥尔科夫斯基公式和软着陆初始条件，可由下式估计： （4-10）式中和分别表示着陆器的终端速度和初始速度，经计算确定搜索范围为(单位秒)： （4-11）我们将用蚁群算法来找到这11个参数。蚁群算法[3]设有n个城市,(i,j=1,2,…,n)表示城市i和j间的距离,表示在t时刻城市i和j之间的信息量。我们以此来模拟实际蚂蚁的分泌物。设共有m只蚂蚁,用表示在t时刻蚂蚁k由城市i转移到城市j的概率,则 (4-12)其中,表示蚂蚁k下一步允许走过的城市的集合,表示路径上的信息量对蚂蚁选择路径所起的作用大小,为由城市i转移到城市j的期望程度,例如,可以取＝。表示的作用.当时,算法就是传统的贪心算法;而当时,就成了纯粹的正反馈的启发式算法.经过n个时刻,蚂蚁可走完所有的城市,完成一次循环.每只蚂蚁所走过的路径就是一个解.此时,要根据下面公式对各路径上的信息量作更新: (4-13)其中,表示信息量随时间的推移而衰减的程度.信息增量可表示为 (4-14)表示蚂蚁k在本次循环中在城市i和j之间留下的信息量。下面给出蚁群算法应用过程的具体步骤。(1)设置初始参数，包括蚂蚁数N，循环次数K，挥发系数Rho，调节系数Lambda、Q，所有路径信息素量初值Tau。，蚂蚁初始位置等。(2)根据式(4-8)计算每只蚂蚁的转移概率，然后依据赌轮原则为每只蚂蚁选择下一个路径。重复上述操作直至所有蚂蚁均完成一次循环。(3)将每只蚂蚁的路径解码为优化参数值和和1个终端时刻，用这10个数据拟合求出，从而得到参量的表达式，根据式(4-6)计算目标函数值求得,对比当前最优值，判断找出当前最好的路径。然后根据式(4-9)、(4-10)更新相关路径上的信息素。(4)判断是否满足终止条件，不满足则重复(1)至(3)步过程；满足则结束计算输出结果。本节模型综述综上所述，4.2.3节的模型可以总结为下式所述的规划问题：目标函数：约束条件： （4-15）4.2.4远月点的位置嫦娥三号的运行轨道近似一个椭圆，月球则是该椭圆的焦点之一，若设椭圆半长轴为a，半短轴为b，焦点为c。则焦点（即月球圆心）距离椭圆轨道（即嫦娥三号的绕月轨道）长轴两端的距离分别为：。也就是远地点和近地点的高度（未减去月球半径）。根据椭圆方程： （4-16）可以求得近月点和远月点的位置关系。4.3模型的求解4.3.1近月点和远月点速度利用式（4-1）可求得远月点的速度为： 利用式（4-2）可求得近月点的速度为：4.3.2近月点的位置接近段的水平距离查阅相关资料[1]可得，处于接近段的着陆器采用接近直线下降方式逐步接近着陆区，因此接近段的水平距离可以近似等于接近段的垂直距离。主减速段的水平距离距离可以根据主减速段最终时刻的最优参量，以及的表达式积分求得，如下所示： （4-17）根据蚁群算法的优化的优化结果，求得: 从而得到的表达式为： （4-18）运行结果图如下：注：时间大于交点处时的线段不在主减速阶段内，故无效图4-2 的函数图像图4-3 的函数图像由图可得，近月点与着陆点的距离在一定的误差允许范围内，可以将近月点与着陆点的距离做如下近似计算：近月点的速度方向速度方向为，即主减速段的俯仰角为5.问题二的建模与求解5.1问题分析嫦娥三号软着陆过程大致可以划分为6个阶段，着陆准备阶段、主减速段、快速调整段、粗避障段、精避障段、缓速下降段。问题一中已经确立了着陆准备阶段和主减速段的运动方程和动力学模型，主减速段的主要目的是通过推力发动机抵消近月点处较大的初始速度，所以燃料消耗最优是这一阶段的主要优化目标。问题一中已给出了相应的数学模型。对于快速调整段，其介于主减速段和粗避障段之间，主要目的是缓冲作用，即使着陆器平缓地从主减速段过渡到接近段。查阅相关资料[1]，为了满足主发动机推力和着陆器姿态匀速过渡的要求,我们利用推力大小和方向线性变化的制导律,其中制导参数利用实际的主减速段末端状态和接近段初始状态约束确定。对于粗避障段，主要任务是要求避开大的陨石坑，实现在设计着陆点上方100m处悬停，并初步确定落月地点。查阅相关资料[4],为了满足制导目标的位置、速度、姿态以及初始高度和速度等多项约束，我们构建平面月球模型，并利用最优开关制导律模拟。在粗避障过程中，利用月球岩石坑的图像特征，采用基于光学图像的成像，通过对图像进行分析和处理，判断出过亮和过暗障碍，识别出当前下降区域上的大障碍，从而确定安全区域。对于精避障段，主要目的是在粗避障选取的较安全区域内进行精确的障碍检测。由于精避障时已极接近月球表面，务必识别并剔除危及安全的小尺度障碍，确保落点安全。为此需要进行高分辨率的三维成像，查阅相关资料，运用基于三维图像的精避障识别和安全着陆区选取算法处理三维图像数据。对于缓速下降段，主要考虑的是着陆的安全性，查阅相关资料[5],我们利用变推力制导律和垂直软着陆模型进行分析。问题二的流程图如下所示： 图5-1 问题二的流程图5.2模型的建立5.2.1快速调整阶段为了满足主发动机推力和着陆器姿态匀速过渡的要求，建立快速调整阶段的推力模型如下所示： （5-1）其中为推力的大小，为推力方向角，即推力方向与当地水平线的夹角。假设快速调整阶段的初始时刻为，终止时刻为，本阶段的动力学方程可以将式（5-1）、式（4-3）与式（4-4）联立，结合下述初值条件： （5-2）来构造。其中,,与之前的定义相同。5.2.2粗避障段最优开关制导律假设粗避障段的初始时刻为，终止时刻为，基于平面月球模型[6]的运动方程如下所示： （5-3） 其中为月面高度，为下降速度与垂直方向的夹角,称为下降角，为月球的重力加速度，为制动力的开关量，其定义如下所示： （5-4）式（5-4）表明粗避障段的着陆器将先做自由落体运动，然后依靠发动机制动下降到接近段终端。为了简化过程，在一定的误差允许范围内，引入以下两个简化条件：（1）简化条件1：由于着陆器在本阶段接近垂直下降，因此可以将视为小量，因而有： （5-5）（2）简化条件2：根据相关资料[1]，本阶段历时时间较短，故引入质量不变假设： （5-6）其中下标表示粗避障段初始状态，下标表示粗避障段终止状态。最终，我们得到了简化后的运动方程如下所示： （5-7）基于光学图像的粗障碍检测算法通过观察，月球表面的岩石坑图像特征如下：1）岩坑表面石头具有高亮度；2）岩坑存在明显的阴影区3）亮区与阴影之间存在强对比度。这些亮暗区都是我们剔除的目标。（1）光学图像的分格与障碍区的识别将光学图像分格为个等大方格，A题提供的TIF图为分辨率，将其灰度值分格为，每格像素为。然后对每格的灰度值进行统计。每格共有个灰度，而灰度过大则表示过亮，即岩坑表面石头；过小表示岩坑，为暗区。当一个格的灰度过大和过小的点过多，这个格将被认定为障碍区，不能降落。因此，需要对灰度过大和过小的点进行统计，统计方法如下：设障碍点数为初值为0，障碍点数阈值为，过暗灰度临界值为,过亮灰度临界值为，对于一格内个灰度点，若： 则表示该灰度点是障碍点，统计值加一。当一格统计完成时，若：即障碍点数超过阈值，此时判断该格为障碍区。为方便起见，标记障碍区所在格为0，非障碍区为1。（2）采用螺旋搜索算法寻找安全着陆区通过螺旋搜索算法（图5-2），从上步得到的的表格中间开始搜索非障碍区。当找到非障碍区后，为其设置一个安全半径。继续搜索，若在安全半径内无障碍区，则将该半径内区域设置为安全着陆区。图5-2 螺旋搜索算法示意图[7]5.2.3精避障段基于三维图像的精避障识别和安全着陆区选取[7]设计的精障碍识别和安全着陆区选取算法为：1）数据预处理，构建测量坐标系下的地形三维高程图；2）平均坡面构建，采用最小二乘法拟合一定单元区域的平均坡面；3）平均坡度计算，根据平均坡面计算该区域的平均坡度；4）障碍高度计算，根据平均坡面估算该区域内每个单元格的障碍高度；5）安全着陆区选取，采用从着陆区中心开始顺时针陆炫前进搜索的方法（图），直至找到符合着陆要求的区域位置。图5-3 中心螺旋式安全着陆区搜索方法[4]5.2.4缓速下降段垂直软着陆模型[6]接近段结束后,着陆器距离月面很近,且下降速度与当地水平面几乎垂直。采用分段控制方式,着陆器依次经过悬停、匀加速、匀减速和关机降落几个过程,各段加速度分别如下所示： （5-8）等效连续变推力制导律[6]假设等效推力的加速度是线性变化的，取向上为正方向，有： （5-9）其中,表示等效推力加速度的斜率,表示推力逐渐增大,表示推力逐渐减小；表示等效推力加速度相对于月心重力加速度；的初始偏移量,表示初始推力加速度大于,表示初始推力加速度小于。结合关机条件： （5-10）对式（5-9）等式两边连续积分两次，得： （5-11）其中，下标和分别表示着陆段初始和关机时刻状态,下标表示着陆段,下标表示关机时刻。5.3模型的求解5.3.1粗避障光学图像的分格与障碍区的识别通过MATLAB的imread函数，读取粗避障光学图像的灰度图。考虑到该灰度图有个像素，用MATLAB进行障碍区的判断会相当慢，于是将该灰度图录入C++下进行模型的求解。用C++进行分格和识别障碍区时，经过多次测试设定各值：阈值过暗灰度临界值过亮灰度临界值得到障碍区的0-1图如下，通过与原TIF图比较可验证其正确性:图5-4 障碍区0-1图与原tif图的比较显然，0区对应原图的岩石坑区，1区对应原图较平坦的区域，1区才是适合着陆的安全区。计算结果和实际TIF图表示的地形非常吻合。圆圈所包括的区域，是最后得出粗避障过程后适合安全降落的区域。5.3.2精避障三维高程图图5-5 精避障三维高程图6.模型的分析与检验6.1 着陆准备段与主减速段模型的分析与检验根据相关资料[5],着陆器动力下降初始状态为:经度19.0464°W,纬度28.9989°N,速度1695.7m/s,俯仰姿态角85°，主减速阶段经历时间为487s。根据经纬度距离计算公式，计算出实际近月点与着陆点的水平距离为： 解得：所以计算结果误差为： 其中，，，，依次为近月点与着陆点水平距离相对误差、近月点速度相对误差、主减速段仰幅角相对误差、主减速段时间相对误差。结论：由上述数据可知，与后三项的值相比，距离误差较大。从近月点与着陆点距离计算过程中可以知道，为了简化计算，我们只考虑了主减速段和粗避障段产生的水平距离，而忽略了其他阶段产生的水平距离，因而误差相对较大，但是仍然控制在10%以内，说明在一定误差允许范围内，此结果具有一定的参考性。而后三项的相对误差比较小，说明他们都与嫦娥三号实际飞行得到的结果比较符合。有力地论证了着陆准备阶段和主减速段模型的准确性与可靠性。6.2粗避障光学图像分格与障碍区识别的分析与检验粗避障的模型通过模型的求解检验了其正确性，在着陆器飞行时，若着陆器能及时计算出安全区域，该模型还是具有一定可行性的。该粗避障模型计算量偏大，对着陆器计算能力无疑是一大考验。该粗避障模型计算结果与实际地形误差较少，若要进一步减少误差，提高准确性，则需提高分格的格数，同时还要再度优化过暗灰度临界值、过亮灰度临界值和阈值。7.模型的创新点与不足之处7.1模型的创新点（1）将月球软着陆的轨迹优化问题转化为参数优化问题，将求燃料最优问题转化为求时间最优问题，使问题的求解变得简单。（2）使用蚁群算法求解复杂的角度变化参数优化问题，具有很好的试验效果。（3）利用光学图像的分格和设定阈值来寻找连片非障碍区作为安全着陆区域。7.2模型的不足之处蚁群算法比较适合处理离散的优化问题，从而限制了其在连续优化问题中的应用，本模型中将角度变化离散化，存在较大的误差。8.参考资料[1]张洪华,关轶峰,黄翔宇,等.嫦娥三号着陆器动力下降的制导导航与控制.中国科学:技术科学,2014,44:377–384[2]王劼，李俊峰，崔乃刚，刘暾，登月飞行器软着陆发动机持续工作初始点选择研究，工程力学，1000-4750，(2003)06-0145-04[3]段佳佳，徐世杰，朱建丰，基于蚁群算法的月球软着陆轨迹优化，宇航学报，第29卷第2期：1000-1328，(2008)02.0476.06[4]KlumppAR.Apollolunardescentguidance.Automatica,1974,10:133–146[5]张洪华,梁俊,黄翔宇,等.嫦娥三号自主避障软着陆控制技术.中国科学:技术科学,2014,44:559–568[6]王鹏基，张熇，曲广吉，月球软着陆飞行动力学和制导控制建模与仿真，中国科学E辑:技术科学2009年，第39卷第3期:521~527[7]张洪华,梁俊,黄翔宇,等.嫦娥三号自主避障软着陆控制技术.中国科学:技术科学,2014,44:559–5689.附录MATLAB源程序function[BESTX,BESTY,ALLX,ALLY]=ACOUCP(K,N,Rho,Q,Lambda,LB,UB)%%ACOUCP.m%在命令行中输入LB=[00000000000];UB=[11111111111];ACOUCP(50,30,0.95,1,0.5,LB,UB)%%输入参数列表%K迭代次数%N蚁群规模%Rho信息素蒸发系数，取值0～1之间，推荐取值0.7～0.95%Q信息素增加强度，大于0，推荐取值1左右%Lambda蚂蚁爬行速度，取值0～1之间，推荐取值0.1～0.5%LB决策变量的下界，M×1的向量%UB决策变量的上界，M×1的向量%%输出参数列表%BESTXK×1细胞结构，每一个元素是M×1向量，记录每一代的最优蚂蚁%BESTYK×1矩阵，记录每一代的最优蚂蚁的评价函数值%ALLXK×1细胞结构，每一个元素是M×N矩阵，记录每一代蚂蚁的位置%ALLYK×N矩阵，记录每一代蚂蚁的评价函数值%%参考设置%[BESTX,BESTY,ALLX,ALLY]=ACOUCP(50,30,0.95,1,0.5,LB,UB)%%第一步：初始化M=length(LB);%决策变量的个数%蚁群位置初始化X=zeros(M,N);fori=1:Mx=unifrnd(LB(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);%适应值初始化%%第二步：迭代过程whilek<=KYY=zeros(1,N);forn=1:Nx=X(:,n);YY(n)=FIT(x);endmaxYY=max(YY);temppos=find(YY==maxYY);POS=temppos(1);%蚂蚁随机探路forn=1:Nifn~=POSx=X(:,n);Fx=FIT(x);mx=GaussMutation(x);Fmx=FIT(mx);ifFmx<FxX(:,n)=mx;Y(n)=Fmx;elseifrand>1-(1/(sqrt(k)))X(:,n)=mx;Y(n)=Fmx;elseX(:,n)=x;Y(n)=Fx;endendendforn=1:Nifn~=POSx=X(:,n);Fx=FIT(x);mx=GaussMutation(x);Fmx=FIT(mx);ifFmx<FxY(n)=Fmx;elseifrand>1-(1/(sqrt(k)))X(:,n)=mx;Y(n)=Fmx;elseX(:,n)=x;Y(n)=Fx;endendend%朝信息素最大的地方移动forn=1:Nifn~=POSx=X(:,n);r=(K+k)/(K+K);p=randperm(N);t=ceil(r*N);pos=p(1:t);TempTau=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);ifFx3<FxX(:,n)=x3;Y(n)=Fx3;elseifrand>1-(1/(sqrt(k)))X(:,n)=x3;Y(n)=Fx3;elseX(:,n)=x;Y(n)=Fx;endendend%更新信息素并记录Tau=Tau*(1-Rho);maxY=max(Y);minY=min(Y);DeltaTau=(maxY-Y)/(maxY-minY);Tau=Tau+Q*DeltaTau;ALLX{k}=X;ALLY(k,:)=Y;minY=min(Y);pos4=find(Y==minY);BESTX{k}=X(:,pos4(1));BESTY(k)=minY;disp(k);k=k+1;end%从蚁群算法中得到：tf和十个角度值angle(单位：度)%计算各个时刻functiontf1=angle(M,tf)fori=1:10td(i)=0+(i-1)*(tf-0)/9;end%三次性项拟合globalA;A=polyfit(td,M,3)aL=1737.013*10^3;%月球半径%求解飞行器动力学方程[T,Y]=ode45('dip',[0750],[148.165,aL+15*10^3,1693.53,0]);%A=[-6.1321*10^(-8)8.0647*10^(-5)0.01635]%figure%plot(T,Y(:,3),'b-',0:750,51.66,'r-')%从数值解上选取坐标(tf,51.66)tf1=interp1(Y(:,3),T,51.66,'nearest','extrap');%[tf1,Vr]=ginput;functiondy=dip(t,y)u=0.4902802627*10^13;%月球引力常数globalA;%飞行器动力学方程%y1是径向速度，y2是月心距，y3是切向速度,y4是横向距离%[Vr,Vo,r,Q]=dsolve('DVr=-u/(r^2)+(Vo^2)/r+7500/(2400-7500*t/2940)*sin(a0+a1*t+a2*t^2+a3*t^3)','DVo=-Vr*Vo/r+7500/(2400-7500*t/2940)*cos(a0+a1*t+a2*t^2+a3*t^3)','Dr=Vr','DQ=Vo/r','(Vr(0))^2+(Vo(0))^2=2890000','r(0)=aL+15*10^3','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)*t^2+A(1)*t^3)*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)*t^2+A(1)*t^3)*pi/180);dy(4)=y(3);end%使优化参数反归一化，代入目标检测函数中functionk=FIT(x)M=[x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10)];M=90*M;M=sort(M);tf2=angle(M,300+x(11)*200);k=tf2;%高斯变异functionmx=GaussMutation(x)l=normrnd(0,1,11,1);k=1/max(l);mx=x.*(1+k*l);mx=mx-floor(mx);aL=1737.013*10^3;%月球半径globalA;%A=[-6.1321*10^(-8)8.0647*10^(-5)0.01635];A=[-6.4716*10^(-7)5.7154*10^-40.00582.2066][T,Y]=ode45('dip',[0750],[148.165,aL+15*10^3,1693.53,0]);figureplot(T,Y(:,3),'b-',0:1:750,51.66,'-');holdonfigureplot(T,Y(:,4),'b-');C++源程序：//粗避障分格、设置障碍区（0）和非障碍区（1）//先用MATLAB读取TIF图的灰度矩阵保存至TXT文件，再用C++通过读取TXT从而读取灰度矩阵#include<iostream>#include<fstream>usingnamespacestd;intmain(){ifstreamfile1("d:\\b.txt",ios::in);short*x;x=newshort[2300*2300];charc;for(inti=0;i<2300*2300;i++){file1>>x[i]>>c;}file1.close();inty[23*23];intthreshold=700,sum;shortDark=84;shortLight=105;intt0,t1;for(intxvec=0;xvec<23;xvec++){for(intyvec=0;yvec<23;yvec++){sum=0;t0=2300*100*xvec; t1=100*yvec; for(inti=0;i<100;i++){for(intj=0;j<100;j++) {if((x[t0+2300*i+t1+j]<Dark)||(x[t0+2300*i+t1+j]>Light)) sum++;}}if(sum<threshold)y[xvec*23+yvec]=1; elsey[xvec*23+yvec]=0;}}ofstreamfile2("d:\\testc.txt",ios::out);for(intxvec=0;xvec<23;xvec++){for(intyvec=0;yvec<23;yvec++){file2<<y[xvec*23+yvec]<<'';}file2<<'\n';}file2.close();}============================================================================基于C8051F单片机直流电动机反馈控制系统的设计与研究基于单片机的嵌入式Web服务器的研究MOTOROLA单片机MC68HC（8）05PV8/A内嵌EEPROM的工艺和制程方法及对良率的影响研究基于模糊控制的电阻钎焊单片机温度控制系统的研制基于MCS-51系列单片机的通用控制模块的研究基于单片机实现的供暖系统最佳启停自校正（STR）调节器单片机控制的二级倒立摆系统的研究基于增强型51系列单片机的TCP/IP协议栈的实现基于单片机的蓄电池自动监测系统基于32位嵌入式单片机系统的图像采集与处理技术的研究基于单片机的作物营养诊断专家系统的研究基于单片机的交流伺服电机运动控制系统研究与开发基于单片机的泵管内壁硬度测试仪的研制基于单片机的自动找平控制系统研究基于C8051F040单片机的嵌入式系统开发基于单片机的液压动力系统状态监测仪开发模糊Smith智能控制方法的研究及其单片机实现一种基于单片机的轴快流CO〈,2〉激光器的手持控制面板的研制基于双单片机冲床数控系统的研究基于CYGNAL单片机的在线间歇式浊度仪的研制基于单片机的喷油泵试验台控制器的研制基于单片机的软起动器的研究和设计基于单片机控制的高速快走丝电火花线切割机床短循环走丝方式研究基于单片机的机电产品控制系统开发基于PIC单片机的智能手机充电器基于单片机的实时内核设计及其应用研究基于单片机的远程抄表系统的设计与研究基于单片机的烟气二氧化硫浓度检测仪的研制基于微型光谱仪的单片机系统单片机系统软件构件开发的技术研究基于单片机的液体点滴速度自动检测仪的研制基于单片机系统的多功能温度测量仪的研制基于PIC单片机的电能采集终端的设计和应用基于单片机的光纤光栅解调仪的研制气压式线性摩擦焊机单片机控制系统的研制基于单片机的数字磁通门传感器基于单片机的旋转变压器-数字转换器的研究基于单片机的光纤Bragg光栅解调系统的研究单片机控制的便携式多功能乳腺治疗仪的研制基于C8051F020单片机的多生理信号检测仪基于单片机的电机运动控制系统设计Pico专用单片机核的可测性设计研究基于MCS-51单片机的热量计基于双单片机的智能遥测微型气象站MCS-51单片机构建机器人的实践研究基于单片机的轮轨力检测基于单片机的GPS定位仪的研究与实现基于单片机的电液伺服控制系统用于单片机系统的MMC卡文件系统研制基于单片机的时控和计数系统性能优化的研究基于单片机和CPLD的粗光栅位移测量系统研究单片机控制的后备式方波UPS提升高职学生单片机应用能力的探究基于单片机控制的自动低频减载装置研究基于单片机控制的水下焊接电源的研究基于单片机的多通道数据采集系统基于uPSD3234单片机的氚表面污染测量仪的研制基于单片机的红外测油仪的研究96系列单片机仿真器研究与设计基于单片机的单晶金刚石刀具刃磨设备的数控改造基于单片机的温度智能控制系统的设计与实现基于MSP430单片机的电梯门机控制器的研制基于单片机的气体测漏仪的研究基于三菱M16C/6N系列单片机的CAN/USB协议转换器基于单片机和DSP的变压器油色谱在线监测技术研究基于单片机的膛壁温度报警系统设计基于AVR单片机的低压无功补偿控制器的设计基于单片机船舶电力推进电机监测系统基于单片机网络的振动信号的采集系统基于单片机的大容量数据存储技术的应用研究基于单片机的叠图机研究与教学方法实践基于单片机嵌入式Web服务器技术的研究及实现基于AT89S52单片机的通用数据采集系统基于单片机的多道脉冲幅度分析仪研究机器人旋转电弧传感角焊缝跟踪单片机控制系统基于单片机的控制系统在PLC虚拟教学实验中的应用研究基于单片机系统的网络通信研究与应用基于PIC16F877单片机的莫尔斯码自动译码系统设计与研究基于单片机的模糊控制器在工业电阻炉上的应用研究基于双单片机冲床数控系统的研究与开发基于Cygnal单片机的μC/OS-Ⅱ的研究基于单片机的一体化智能差示扫描量热仪系统研究基于TCP/IP协议的单片机与Internet互联的研究与实现变频调速液压电梯单片机控制器的研究基