A10009072_吉张鹤轩_杨升_陈同广(1)【公众号：数模加油站】.pdf

2014 高教社杯全国大学生数学建模竞赛高教社杯全国大学生数学建模竞赛 承承 诺诺 书书 我们仔细阅读了全国大学生数学建模竞赛章程和全国大学生数学建模竞赛参赛规则（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。我们参赛选择的题号是（从 A/B/C/D 中选择一项填写）：A 我们的报名参赛队号为（8 位数字组成的编号）：10009072 所属学校（请填写完整的全名）：东南大学 参赛队员(打印并签名)：1.吉张鹤轩 2.杨升 3.陈同广 指导教师或指导教师组负责人 (打印并签名)：（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。）日期：2014 年 09 月 15 日 赛区评阅编号（由赛区组委会评阅前进行编号）：获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576572014 高教社杯全国大学生数学建模竞赛高教社杯全国大学生数学建模竞赛 编编 号号 专专 用用 页页 赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评 阅 人 评 分 备 注 全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576571 嫦娥三号软着陆轨道设计与控制策略 摘要 本题要求我们以嫦娥三号登月为背景，分析登月轨道参数，重点探讨了登月过程最具难度的着陆轨道设计优化，并对所使用的优化方案进一步作出了误差分析与灵敏度分析。对于第一问，由于正面求解条件有限，难以从已有的条件中得到近月点和远月点的位置以及准备轨道参数，因此巧妙的使用了逆推思路，通过已知条件求解主减速阶段运动过程，通过水平位移量反推近月点位置。首先，我们由已知的端点约束条件，结合所学物理知识，并查阅相关资料，通过线性正切制导率等工程经验判断减速过程中推力控制方案应为推力大小始终为最大值，推力与速度反方向夹角也为恒量，由此建立微分方程模型。但在求解的过程中我们发现，正面求解难度十分大，于是对微分方程离散化，转化为差分方程组，继而通过计算机模拟行为，拟合出最接近解的轨迹，求得水平位移量为 385.21m，由此得到近月点位置为 19.51 W，31.50 N，距月面 15km，速度为 1692.46m/s，平行月面指向预期落点方位；远月点位置为 160.49 E，31.50 S，距月面 100km，速度为 1612.15m/s，方向与近日点反相平行。对于第二问，依据题意将其分为五个阶段。对于前两个阶段，建立天体力学分析中常用的二体模型及力学方程组，并对约束条件进行归一化处理。为了后续优化，列出全部归一化约束条件并建立目标函数，从而获得非线性规划模型。之后采用了序列化遗传算法对其进行筛选，从而逼近该阶段全局最优解，最终此阶段燃耗为 1055.39kg。第二阶段同样属于复杂多变量优化问题，同阶段一建立优化模型并通过遗传算法求解该段全局最优解，最终燃耗为26.71kg。第三第四阶段主要在于图像处理与统计，第三阶段，通过对高程图进行 K 均值聚类分析，将图中像素点分为安全点与危险点，在对地图栅格化，并对方格内点类型统计取整，对方格二元化为安全格与危险格，在通过扩大寻找最大安全半径，综合考虑水平偏移量，建立合理的落点评价体系，最终找出最优点坐标（1275,1000），燃耗为 86.97kg。第四阶段同样做聚类分析，并根据嫦娥三号实际体积选取合适的栅格大小，并对栅格通过最小二乘法对空间进行线性统计回归，求出平均坡面与平均坡度，结合最大安全半径建立最优落点评价体系，最终获得最优落点坐标为（88,56），燃耗为 20.68kg。第五阶段最优燃耗为 8.09kg。最后还讨论了简单运动的局部最优模型，简化了后几个阶段的运动学分析与计算。最终综合各段最优解，获得最优着陆轨道与控制策略。对于第三问，首先总结了优化模型中引入的一些误差因素，并针对主要因素做了数值上的相对误差分析，证明了误差对于优化方案并未产生很大影响。其次从初始变量和约束条件入手，分析了这些变量的波动对于结果产生的影响，最终发现角度控制向量的灵敏性较高，而其他因素的灵敏性普遍处在较低水平，侧面说明了优化方案的对于初值的不敏感性与方案对于全局最优解的逼近程度较高。关键词关键词：非线性规划模型 序列化遗传算法 K 均值聚类 空间线性回归 二体模型 获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576572 一问题的提出 1.1 背景介绍 根据计划，嫦娥三号将在北京时间 12 月 14 号在月球表面实施软着陆。嫦娥三号如何实现软着陆以及能否成功成为外界关注焦点。目前，全球仅有美国、前苏联成功实施了 13 次无人月球表面软着陆。北京时间 12 月 10 日晚，嫦娥三号已经成功降轨进入预定的月面着陆准备轨道，这是嫦娥三号“落月”前最后一次轨道调整。在实施软着陆之前，嫦娥三号还将在椭圆轨道上继续飞行，做最后准备。嫦娥三号着陆地点选在较为平坦的虹湾区。但由于月球地形的不确定性，最终“落月”地点的选择仍存在一定难度。在整个“落月”过程中，“动力下降”被业内形容为最惊心动魄的环节。在这个阶段，嫦娥三号要完全依靠自主导航控制，完成降低高度、确定着陆点、实施软着陆等一系列关键动作，人工干预的可能性几乎为零。在距月面 100 米处时，嫦娥三号要进行短暂的悬停，扫描月面地形，避开障碍物，寻找着陆点。1.2 问题重述 嫦娥三号在着陆准备轨道上的运行质量为 2.4t，其安装在下部的主减速发动机能够产生的推力可调节，变化范围为 1500N 到 7500N，其比冲为 2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为 19.51W，44.12N，海拔为-2641m。嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点 15km，远月点 100km 的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为 6 个阶段，要求满足每个阶段在关键点所处的状态，尽量减少软着陆过程的燃料消耗。根据上述的基本要求，建立数学模型解决下面的问题：（1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。（2）确定嫦娥三号的着陆轨道和在 6 个阶段的最优控制策略。（3）对设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。二问题的分析 2.1 问题一 由于题目已经明确给出准备轨道的形状参数，可以通过已有的物理知识与几何关系明确计算出近月点与远月点处速度大小和相对于月面的速度方向。为了预期落点与着陆轨道在同一个平面内，且准备轨道过月心，可以大致确定轨道所在平面有无数多个，无法确定近月点和远月点的位置，因此需要其他条件来推测。本题现有的条件下，只有通过着陆轨道逆推近月点，并结合地月轨道制动的实际情况综合考虑，才能得到。对于着陆过程一，由相关报道以及 NASA 在 1976 年提出的线性正切制导率1，得知主减速阶段通常都是恒推力作用在轨道切线上，且嫦娥三号主减速阶段实际也是保持着最大推力依照这一定律进行制导。我们由此出发，通过二体模型，结合已知条件，建立微分方程组，通过计算机模拟降落轨迹即可求出降落弧线距离，从而反推近月点，对称得出远月点。2.2 问题二 由问题一已经得出近月点，即开始降落点位置，也知道每一阶段的状态，因此，降落轨获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576573 道大致范围基本确定，但六个过程的精确路径是要通过策略优化来控制的。由于燃料消耗表现在推力在时间上的积累量，即减小推力作用的冲量，即可优化燃料消耗。对第一个过程，由于推力很大且历时较长，因此燃料消耗主要体现在这一阶段，对应的，优化策略也应重点体现，由于有二体模型，建立微分方程模型，并由初值条件以及阶段限定条件，可以写出非线性约束条件，本问题及转化为轨道优化中的非线性规划问题，一般可通过成熟的 SQP 算法可以得到全局最优解，但本题采用了更为常见也相对传统的遗传算法，逼近全局最优解，得出最优方案。对于第二阶段，仅仅是为了是水平速度将为 0，且推力迅速减小，由上一阶段的优化结果，得出末速度水平分量，由于冲量可分解，则此阶段分为水平方向和竖直方向分别优化，即分为了两个变速直线运动模型，简化了优化模型，可以由这两个局部最优解加和得到该阶段的全局最优解。第三个阶段水平速度初始为 0，经过对月面成像分析，制定平坦度评价体系，选择距离中心点最近且满足平坦度要求的区域中心为粗调整目标降落点。由于这一段终点悬停，速度减为 0，因此可以对该段推力进行优化，从而局部燃料最优。第四阶段悬停，精细成像并分析，同样制定平坦度评价体系并选择距离中心点最近且满足平坦度要求的区域中心作为最终目标降落点，修正轨道。结束时水平速度依然为 0，因此同样存在优化过程。第五阶段与第六阶段是减速至 0 然后自由落体的过程，针对减速阶段也可考虑优化，可经过简单讨论得到结果。最后根据各个阶段的最优方案，模拟出嫦娥三号着陆轨道即可。2.3 问题三 为了分析设计轨道和控制策略的误差与敏感性，有必要制定误差指标并考虑各部分误差对于结果的最大影响，敏感性也同样需要这一思路，各个阶段的细微变化会对结果产生影响的衡量。另一方面，或许还有必要寻找参考物以显示该方案的好坏。三模型假设 1.由于月球自转速度为 27.3d，十分缓慢，而着陆过程仅有十几分钟，因此在本题中月球不考虑自转；2.由于月球扁率很小，可认为月球为球体，半径以平均半径为准，并且引力场分布均匀；3.由于侧面姿态调整喷射装置对燃料影响很小，为简化模型，认为飞行器变换姿态的过程不消耗燃料；4.认为飞行器变换姿态是瞬间完成的；5.由于着陆时间较短，所以诸如月球引力非球项、日月引力摄动等影响因素均可忽略不计；6.推力大小可瞬间改变；7.燃料除供给推力外，无任何其他耗散方式；8.月球空气稀薄，不考虑任何摩擦力；9.由于预定着陆点海拔-2641m，因此在着陆过程中所使用的高度均不应是海拔高度，而是相对于着陆点海拔的高度。获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576574 四符号说明 符号 符号意义 G 引力常量 M 月球质量 m 嫦娥三号质量 比冲 mg 月球表面重力加速度 F 推力 J 燃耗 P 平坦度评价指标 H 综合评价指标 平均坡度评价指标 五模型建立与求解 5.1 问题一 5.1.1 模型建立 为了确定近月点和远月点，正常的思路是求出椭圆轨道所在平面以及椭圆长轴在空间中的位置，但本题仅给出了轨道两点的高度信息以及平均半径，还有根据常识推断出的预定着陆点在椭圆平面内，除此之外并无其他信息，因此从现有条件准确判断两点的位置是不可能的。因此，本题应采用逆推思路，由着陆轨道的第一个阶段反推近月点。先对轨道参数进行分析，已知近月点高度 Hc=15km，远月点高度 Hf=100km，月球平均半径为 R=1737.013km，轨道为椭圆。如图 5.1.1.1.图 5.1.1.1 近月轨道示意图（为了方便示意，本图不符合比例）由开普勒定律，任何椭圆天体轨道的中心天体一定在椭圆的一个焦点上。则由椭圆几何性质，可以近似得出方程：a+c=Hf+R a-c=Hc+R 联系万有引力定律与牛顿第二定律，可以列出嫦娥三号在近月点和远月点的运动学方程：20022002(a c)(a c)ffccvMmGmMmvGm 获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576575 其中227.3477 10Mkg为月球质量，-11=6.672 10G为引力常量，02.4mt为在近月轨道上飞行器的质量，f与c为远月点和近月点的曲率半径。由椭圆的几何性质，在长轴两端点处的曲率半径分别为：22=1=1cfbaceaaaceb 联立以上各式，带入参数，得：近月点速度大小1692.46/cvm s；远月点速度大小1614.15/fvm s。由假设 2，月球视为球体，则近月点处速度方向方向应平行于月面且为了着陆方便，方向指向预定着陆点所在方向，远月点处同样平行月面但速度方向与之相反。下面开始研究主减速过程，该过程由近月点开始，切向速度 0cvv，径向速度rv等于0，由切线正切制导率及相关报道可知，该阶段推力始终保持最大推力，且始终调整使推力与速度方向相反，即7500FN，Fv，最终应大致达到预定着陆点目标上空，且竖直方向速度为 57m/s，高度共下降 12000m。由于有端点限制条件，可以通过物理知识建立微分方程组，即采用微分方程模型进行分析求解。由于飞行器绕月球表面飞行，且切向速度较大，在月心极坐标系下必须考虑向心力和科里奥利力，且涉及转动变量，较为复杂。这里采用参考系转化，将极坐标系转化为非惯性参考系，并对该系中所有质点提供一个向上的离心力，大小由水平速度与质点到月心距离决定：2xdvFmr 这样，问题就转化到了平面直角坐标系下，并且可以对速度和受力进行正交分解，如图5.1.1.2。图 5.1.1.2 非惯性系下力学分析图示 由此可得微分方程模型：获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576576 220sincostanyxxyxyyxdvvGMadtrrv vdvadtrdrvdtFammtvv 其中 r 为飞行器到月心的距离，为推力与水平方向的夹角，a 为推力产生的加速度，整个过程中质量是匀速减少的，减少系数由比冲与推力定义，即Fm，本题中=2940/N kg。至此，本问题已经完全抽象为数学语言，并且有初始参数：0yv，1692.46/xvm s，1749.37rkm，=0 理论上可求得当高度下降 12000m 时速度方向大小以及飞行距离，但这个方程组要想解出一个描述运动轨迹的函数是十分困难的，因此，需要另寻它路，逼近它的解。5.1.2 模型求解 本题由方程组反推运动方程是十分困难的，但这类问题就像是解决 NP 难题一样，无法由问题得到结论，但可以“猜测”结论从而验证问题的正确性，继而得到想要的数据。回顾本题，之所以无法解决是因为微分方程的连续性，这种连续性计算机很难求解，但若是让该组方程做某种近似，使之转化为差分方程，从而利用差分方程的离散型导出差分方程组。输入初始条件，使用计算机模拟其运动过程，最后打点画出折线图。由微元思想可知，当时间被微分成足够小的时候，模拟的运动轨迹可认为无限逼近真实运动轨迹。微元化后的差分方程依然可以有运动学规律得出：1212112121(cos)(sin)1cos21sin2xnynxnxnnnnxnynynnnnnnnnnnnnnnxnnnnnynnnv vvvatrvGMvvatrrFammmm trryxxvtatyyvtat 其中t为时间微元小量，由于减速段大约几百秒，为了精确描绘运动曲线，t一般取 0.1s。带入初始条件：7500FN，02400mkg，0=0，00 x，00y，01692.46/xvm s，00yv 即可通过迭代法，求出在 y=3000m 时的 x 的值，该值即为飞行器转过的角度所对应的获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576577 月面弧长。但是仿真结果表明，当高度降至 3000m 时，速度并不为 57m/s，而且远大于这一数值，显然，是我们采用的控制方案不合题意。重新考虑这一过程，应该是竖直方向分力不够，导致竖直速度分量增加过快。之后在了解了嫦娥三号实际登月过程后，我们发现这一过程中推力并非与速度方向相反，而是始终与速度反方向保持一个角度，并向曲线内侧倾斜。于是我们通过改变 的大小，试图输出大量运动轨迹簇，这一过程由计算机模拟生成。通过计算机仿真模拟，得到非惯性参考系下的运动轨迹簇。其中可以发现，在 取时，高度 3000m 处竖直方向速度为，在误差允许的范围内可以接受，由此得出符合题意的主减速阶段飞行轨迹图，如图 5.1.1.4.图 5.1.1.4 主减速阶段飞行轨迹图（非惯性参考系）图 5.1.1.5 水平与竖直速度变化示意图 图 5.1.1.6 夹角变化示意图 得到弧长为，对应转过角度为，所用时间为。由于推力大小不变，由时间即可得出此阶段燃耗。该阶段内水平方向与竖直方向速度变化如图 5.1.1.5.可以看出，水平速度不断下降，而竖直方向速度先增大后减小，合速度 57.12m/s，符合题意。水平总位移为 385.21km。虽然该阶段末并不是准确到达预定着陆点上空，需要经过快速调整阶段才能准确到达，但题目表述中明确指出该处以基本到达目标上空，再加之段末水平速度分量相比初始速度已经十分小，快速调整阶段又很迅速，即快速调整阶段对到目标点上空距离的影响十分小，相比主减速段弧长可以忽略，即认为主减速末段已经到达目标上空。由此，我们结合弧长与转角可以计算近月点所在范围。由于该角度可在以预定落点为圆心的球面圆上任意选取，考虑到当天月球月面范围，降落过程应基本暴露在有光月面一侧，以便拍摄或信息收集，并且一般选择轨道时为了方便计算控制，轨道平面会与落点所在经线获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576578 平面基本重合，即近月点经度也为 19.51 W，纬度应更加靠近赤道，从而保证降落过程全程在光侧面，最终确定纬度位置为 31.50 N。由月球的球对称性，远月点纬度南北对调，经度东西对调，角度互补。最终结论，近月点位置为 19.51 W，31.50 N，距月面 15km，速度为 1692.46m/s；远月点位置为 160.49 E，31.50 S，距月面 100km，速度为 1612.15m/s。5.2 问题二 本题为典型的多因素优化问题，优化目标是减少燃料损耗，并且通过图像分析选择最为平坦的地区着陆。由于燃料损耗直接与推力有关，即Fm，其中=2940/m s为燃料比冲，也就有 Fdtmdt 该式左边为推力对时间的积累量，即推力所做的冲量；右边为燃料损耗对时间的积累，即一定时间内的燃料损耗总量。由于比冲为常数，因此推力冲量与燃耗成正比例关系，从而使优化目标与力学量直接联系，便于分析。另外，由动量定理mdvFdt，在近似处理下，也可将燃耗与速度变化联系在一起。5.2.1 过程一模型建立与求解 由 5.1 计算过程可知，主减速段是着陆过程用时最长，燃耗最多的阶段。该阶段的主要任务是消除较大的初始水平速度,因此推进剂消耗优化是该阶段优化的主要设计目标，也是整个着陆轨道燃耗优化的重中之重。嫦娥三号的主减速阶段是从近月点开始下降到距离预定地点 3km处，这一过程的时间比较短，仅有几百秒的时间，所以可以不考虑月球引力摄动。月球自转速度比较小,也可忽略。所以，可以仅仅对月球和嫦娥三号进行分析，将问题转化为二体模型，如图5.2.1.1。以月心为原点，建立平面直角坐标系，设嫦娥三号的月心距为r，极角为，角速度为，质量为m；又设v为嫦娥三号沿r方向上的速度，F为主减速发动机的推力（在本题中为固定值），发动机推力与当地水平线的夹角即为推力的方向，比冲为SPI；设月球的引力常数为。图 5.2.1.1 二体模型力学分析示意图 首先，从运动学角度进行分析。由动力学基本方程可得drtdv，所以drvdt 同理，可以得到角度变化的关系式为ddt 获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：5444576579 图 5.2.1.2 受力分析示意图 再从动力学角度进行分析。如图，对嫦娥三号进行受力分析可知。在径向上，受到万有引力22GM=mmFrr月引以及主减速发动机在径向上的分力Fsin。由此，根据牛顿第二定律可以得到：22mmvrr Fsindvmdt 联立，可得22Fsinmmvdvmrdtr，利用公式vr化简该式可得：22sindvFrdtmr 再对切向上的受力进行分析，此处需要注意到我们选择参考系月球自身的自转，是一个非惯性参考系。因此，嫦娥三号除了受到主减速发动机推动力的分量Fcos的作用，还需要计入大小为2v的科里奥利力，根据右手定则可以确定其方向为切向。由牛顿运动学第二定律可以得到Fcos2ydvvmdt。同样，利用yvr化简可以得到：*cos2*Fvdmdtr .由于飞船在运行过程利用燃料反冲制动，因此飞船得而质量是不断变化的。考虑比冲的定义为“火箭发动机单位质量推进剂产生的冲量”可得：SPdmFdtI 。综上所述，联立式，即可得到完整的描述嫦娥三号与月球这一二体运动模型的方程组：获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765710 drvdt 22sindvFrdtmr ddt cos2Fvdmdtr SPdmFdtI 对于嫦娥三号燃料使用指标的衡量，我们可以反映到其动量上去，即将动量作为衡量燃料消耗的性能指标。由质量变化的推导，我们可以得到0ftSPtFJdtI。在轨道优化过程中,由于各状态变量的量级相差较大,寻优过程中可能会导致有效位数的丢失2。归一化处理可以克服这一缺点，提高计算精度。另外，由于对轨道的优化也要求优化变量尽可能地保持在相同的量级，故作以下处理：令：0refrr,0refmm.则：refrrr,refvvv,refrefvr,refSPSPrII,refFFF,2refrefrefrefm vFr,refmmm,3refr,refttt,refrefrefrtv,.那么，嫦娥三号的动力学方程可改写为：drvdt 221sindvFrdtmr ddt cos2Fvdmdtr SPdmFdtI 又由第一问可以得到飞行器的初始条件和终端条件分别为:00001692.46/1749.372rvvm srkm 057/rffvvm s 注：此处的终端条件中，径向速度和切向速度为理想情况下的结果，实际只能得到最优获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765711 化的值。使用归一化条件可以得到初始条件和终端条件分别变为：00001692.46/1749.372rrefrefvvm svrkmr 057/rffrefvvm sv 在考虑燃料最优的情况，并要保障水平速度几乎减到0/m s 的前提下，可以选择推力较大的发动机完成软着陆任务，即7500FN。由此可以得到优化的目标函数为：0=ftfSPSPtFFJdttII 其中，推力F和比冲SPI 都是定值，所以即要求时间最短的方案。最优控制问题的数值方法主要包括直接法和间接法，对上述最优控制问题，使用一种直接方法参数化方法进行求解3。即作如下的假设：推力方向角 可以表示成一个多项式的形式：60iiia t 由此，对于所描述的问题可由一个有约束条件的非线性规划问题描述。由于确定了嫦娥三号的推力，那么在确定变量为0a到6a七个变量后，根据所建立的二体模型方程即可确定运行的轨迹，所以此处需要优化的参数为0a到6a七个变量，问题的约束条件为归一化后的动力学方程以及初值条件和终端条件。由此需要求出燃料最优化方案，即性能指标=fSPFJtI的最小值。对以上转化出来的问题，我们利用浮点数编码的遗传算法（FGA）进行求解。其步骤如下：1.将 7 个取值范围给定的优化参量按一定的浮点数编码排列在一起成为一个个体，随机产生 2000 个这样的个体作为初始种群；2.通过编写的适应函数计算每一个个体的性能指标；3.使用轮盘赌法作为选择算子并对这 N 个个体进行排序；4.选择出若干个性能指标取值较优的个体保留，并将其遗传到下一代；5.将个体随机两两配对，按照既定的交叉概率进行交叉操作；6.对每一个个体中的每一个参数，按照既定的变异概率进行编译操作；7.若满足收敛条件则输出最优解，否则继续进行编码、评价、选择、交叉和变异等操作。通过对约束条件的分析，可以选择交叉概率为 0.6，变异概率为 0.05。通过计算机仿真模拟，由遗传算法规划出的角度随时间控制函数参数向量的取值为：2568101215067.58 107.01 104.30 108.12 106.03 101.60 101.40 10aa 运动轨迹见图 5.2.1.2.获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765712 图 5.2.1.3 运动轨迹图 图 5.2.1.4 速度对比图 图 5.2.1.5 角度变化图 当下降 12000m 时速度 v=57.12m/s，水平方向与垂直方向分别为48.63/xvm s，29.96/yvm s，符合题意。整 个 过 程 耗 时 413.2s，由 于 推 力 F 恒 为 7500N，因 此 该 阶 段 燃 料 消 耗 为11055.39FtJkg，剩余质量 1344.61kg，水平位移为 385.21km。5.2.2 过程二模型建立与求解 快速调整阶段主要任务有两个：一是在主减速阶段完成的基础上，将嫦娥三号水平方向的速度减到0/m s；二是快速调整嫦娥三号的姿态，使推力方向和所收到的引力方向基本一致。在该阶段，考虑让F的方向逐渐调整到竖直方向，并且在调整的过程中（即下落的600m距离内），将水平方向的速度减少到0/m s。图 5.2.1.6 受力分析示意图 获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765713 由此，首先进行运动学分析。在水平方向上，速度最终减小为有0/m s，所以有 xxdvadt 在竖直方向上，由加速度的定义，有这样的关系：yydvadt 然后，在进行动力学分析。在水平方向上，根据牛顿第二定律可以得到：cosxFma 在竖直方向上，不妨设月球的重力加速度为mg，其值由黄金代换公式2mGMg R可以确定。那么同理可以得到：sinmymgFma 至此，受力分析完毕。而另一方面，考虑到减速过程中燃料的消耗对质量产生影响，有：SPdmFdtI 综合，式可以得到描述该阶段运行的方程：cossinxymSPdvFmdtdvmgFmdtdmFdtI 由主减速阶段的终值以及快速调整阶段的要求可以得到该方程初始条件和终端条件分别为：00048.63/29.96/1344.61xyvm svm smkg 0/9090 xfofofvm s 而力F的取值范围为15007500NFN。最优控制问题的数值方法主要包括直接法和间接法，对上述最优控制问题，使用一种直接方法参数化方法进行求解3。即作如下的假设：推力F 可以表示成一个多项式的形式：60iiiFa t。另外，在快速调整阶段，燃料的性能指标同样为0ftSPtFJdtI。综上所述，即可将求解快速调整阶段的问题转换成这样的一个非线性规划问题：需要优化的参数为0a到6a以及、九个变量，问题的约束条件为方程组以及初值条件和终端条件。由此需要求出燃料最优化方案，即性能指标J的最小值。同样使用浮点数编码的遗传算法来解决该非线性规划问题，由于具体步骤与主减速阶段基本一致，故此处不再赘述。在这样的算法下，求得的最优解为：获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765714 20/29.96/90903928.11520.02726.71xyoovm svm sFNtsJkg 5.2.3 过程三模型建立与求解 粗避障段的范围是距离月面 2.4km 到 100m 区间，其主要是要求避开大的陨石坑，实现在设计着陆点上方 100m 处悬停，并初步确定落月地点。嫦娥三号在距离月面 2.4km 处对正下方月面 2300 2300m 的范围进行拍照，获得数字高程，嫦娥三号在月面的垂直投影位于预定着陆区域的中心位置。该高程图的水平分辨率是 1m/像素，其数值的单位是 1m。由附件 3 所提供的灰度图，获取各像素点高程，统计各高程数据出现的个数，绘制高程数量分布曲线图，如图 5.2.3.1。图 5.2.3.1 高程数量分布曲线图 由高程图大致可以发现，在高程约为 95m 时有最高点，120m 以上或 80m 以下数量明显下降，据观测可以认为平坦的区域的高程大致在 80m 到 120m 之间，于是应大致将高程划分为三个组，低于 80m 的区域较暗，高于 120m 的区域较亮，输出的图像上大致可以分辨出较大的陨石坑，但是小陨石坑仍然无法与平地区分，因此需要从统计学角度合理划分。通过对全部高程数据进行 K 均值聚类分析，并设置聚落数，经过尝试发现当聚落数为 5时，可以较为清楚的区分出陨石坑与平地，如图 5.2.3.2。获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：544457657获取更多数学建模相关资料关注【公众号：数模加油站】国赛交流分享群：54445765715 图 5.2.3.2 聚落色阶图 聚类中心值分别为 92.4595，113.8753，138.4552，168.0559，44.6447，平坦区域聚类编号为 1。图中可以发现，对于平坦度造成影响的障碍地形大致分为三类：过高障碍区域，过低障碍区域，条纹障碍区域。通过排除 1 以外的聚类编号，可以准确排除前两种障碍区域，但是第三种条纹障碍区域无法去除，地图上会留下数个环形地带。至此为止，图像分析基本完成，但由于此处降落点的调整同时受制于落点平坦度与落点偏移量，即落点到图像中心的直线距离，因为这一距离会导致飞行器水平冲量的积累，燃耗增加，偏离最优解。因此，有必要制定数字评价手段，通过地图中现有参数，对平坦度以及偏移量做出数字化的评价，综合权重分配即可得出落点的综合评价指标。下面先制定平坦度评价指标P。由图像聚类结果，可以粗略将图中像素点分为安全点与危险点。但考虑到 2300 2300 点数量实在太大，为了减少数据处理量，对地图进行栅格化，每 5m 5m，即 25 像素点作为一个栅格单元，根据格内安全点与危险点数量的多少定义该格为安全格或危险格，这样，地图便被分为了 460 460 格的二元化栅格网络。以任意一个安全格为中心，搜索其外围 8 个方格，确认其中无危险格后，继续扩大一周考察范围，如此循环，直到出现危险格停止，记录无危险格的最大安全半径SR（以格数衡量，若第一圈就出现危险格，则认为半径为 0 格）。SR越大，则认为该格所在区域越安全，P越高。由于平坦度仅与该因素有关，因此，若要将SR转化为实际距离长度，正比例系数应取 5m，即5SPR。条纹障碍由于P非常小，可以有效排除。接下来考虑落点对中心的偏移量的评估指标。由于最终优化目标是燃耗最小，因此