1、2014 高教社杯全国大学生数学建模竞赛 承 诺 书 我们仔细阅读了全国大学生数学建模竞赛章程和全国大学生数学建模竞赛参赛规则 （ 以下简称为 “ 竞赛章程和参赛规则 ” ，可从全国大学生数学建模 竞赛网站下载 ） 。 我们完全明白，在竞赛开始后参赛队员不能以任何方式 （ 包括电话、电子邮件、网上咨询等 ） 与队外的任何人 （包括指导教师 ） 研究、讨论与赛题有关的问题。 我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料 （包括网上查到的资料 ） ，必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。 我们郑重承诺，严格遵守竞赛章程和参赛规则
2、，以保证竞赛的公正、公平性。 如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。 我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示 （ 包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等 ） 。 我们参赛选择的题号是（从 A/B/C/D 中选择 一 项填写 ） ： A 我们的报名参赛队号为（ 8 位数字组成的编号 ） ： 10009072 所属学校（请填写完整的全名 ） ： 东南大学 参赛队员 (打印并签名 ) ： 1. 吉张鹤轩 2. 杨升 3. 陈同广 指导教师或指导教师组负责人 (打印并签名 )： （论文纸质版与电子版中的以上信息必须一致，只是电
3、子版中无需签名。以 上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。 ） 日期： 2014 年 09 月 15 日 赛区评阅编号（由赛区组委会评阅前进行编号）： 2014 高教社杯全国大学生数学建模竞赛 编 号 专 用 页 赛区评阅编号（由赛区组委会评阅前进行编号）： 赛区评阅记录（可供赛区评阅时使用）： 评 阅人 评 分 备注 全国统一编号（由赛区组委会送交全国前编号）： 全国评阅编号（由全国组委会评阅前进行编号）： 1 嫦娥三号软着陆轨道设计与控制策略 摘要 本题要求我们以嫦娥三号登月为背景，分析登月轨道参数，重点探讨了登月过程最具难 度的着陆轨道设计优化
4、，并对所使用的优化方案进一步作出了误差分析与灵敏度分析。 对于第一问，由于正面求解条件有限，难以从已有的条件中得到近月点和远月点的位置 以及准备轨道参数，因此巧妙的使用了逆推思路，通过已知条件求解主减速阶段运动过程， 通过水平位移量反推近月点位置。首先，我们由已知的端点约束条件，结合所学物理知识， 并查阅相关资料，通过线性正切制导率等工程经验判断减速过程中推力控制方案应为推力大 小始终为最大值，推力与速度反方向夹角也为恒量，由此建立微分方程模型。但在求解的过程中我们发现，正面求解难度十分大，于是对微分方程离散化，转化为差分方程组，继而通 过计算机模拟行为，拟合出最接近解的轨迹，求得水平位移量为
5、 385.21m，由此得到近月点 位置为 19.51W ，31.50N ，距月面 15km，速度为 1692.46m/s，平行月面指向预期落点方位； 远月点位置为 160.49E， 31.50S ，距月面 100km，速度为 1612.15m/s，方向与近日点反相平 行。 对于第二问，依据题意将其分为五个阶段。对于前两个阶段，建立天体力学分析中常用的二体模型及力学方程组，并对约束条件进行归一化处理。为了后续优化，列出全部归一化 约束条件并建立目标函数，从而获得非线性规划模型。之后采用了序列化遗传算法对其进行 筛选，从而逼近该阶段全局最优解，最终此阶段燃耗为 1055.39kg。第二阶段同样属于
6、复杂 多变量优化问题，同阶段一建立优化模型并通过遗传算法求解该段全局最优解，最终燃耗为 26.71kg。第三第四阶段主要在于图 像处理与统计，第三阶段，通过对高程图进行 K 均值聚 类分析，将图中像素点分为安全点与危险点，在对地图栅格化，并对方格内点类型统计取整， 对方格二元化为安全格与危险格，在通过扩大寻找最大安全半径，综合考虑水平偏移量，建 立合理的落点评价体系，最终找出最优点坐标 （ 1275,1000） ，燃耗为 86.97kg。第四阶段同样 做聚类分析，并根据嫦娥三号实际体积选取合适的栅格大小，并对栅格通过最小二乘法对空 间进行线性统计回归，求出平均坡面与平均坡度，结合最大安全半径建
7、立最优落点评价体系， 最终获得最优落点坐标为（ 88,56） ，燃耗为 20.68kg。第五阶段最优燃耗为 8.09kg。最后还 讨论了简单运动的局部最优模型，简化了后几个阶段的运动学分析与计算。最终综合各段最 优解，获得最优着陆轨道与控制策略。 对于第三问，首先总结了优化模型中引入的一些误差因素，并针对主要因素做了数值上 的相对误差分析，证明了误差对于优化方案并未产生很大影响。其次从初始变量和约束条件 入手，分析了这些变量的波动对于结果产生的影响，最终发现角度控制向量的灵敏性较高， 而其他因素的灵敏性普遍处在较低水平，侧面说明了优化方案的对于初值的不敏感性与方案 对于全局 最优解的逼近程度较
8、高。 关键词 ： 非 线 性 规 划 模型 序 列 化 遗 传 算 法 K 均 值 聚类 空 间 线 性 回归 二 体模型 2 一问题的提出 1.1 背景介绍 根据计划，嫦娥三号将在北京时间 12 月 14 号在月球表面实施软着陆。嫦娥三号如何实 现软着陆以及能否成功成为外界关注焦点。目前，全球仅有美国、前苏联成功实施了 13 次无人月球表面软着陆。 北京时间 12 月 10 日晚，嫦娥三号已经成功降轨进入预定的月面着陆准备轨道，这是嫦 娥三号 “落月 ”前最后一次轨道调整。在实施软着陆之前，嫦娥三号还将在椭圆轨道上继续飞 行，做最后准备。 嫦娥三号着陆地点选在较为平坦的虹湾区。但由于月球地形
9、的不确定性，最终 “落月 ”地点的选择仍存在一定难度。在整个 “落月 ”过程中， “动力下降 ”被业内形容为最惊心动魄的环节。在这个阶段，嫦娥三号要完全依靠自主导航控制，完成降低高 度、确定着陆点、实施软 着陆等一系列关键动作，人工干预的可能性几乎为零。在距月面 100 米处时，嫦娥三号要进 行短暂的悬停，扫描月面地形，避开障碍物，寻找着陆点。 1.2 问题重述 嫦娥三号在着陆准备轨道上的运行质量为 2.4t，其安装在下部的主减速发动机能够产生的推力可调节，变化范围为 1500N 到 7500N，其比冲为 2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，能够自动通过多个
10、发动机的脉冲组合实现各种姿态的 调整控制。嫦娥三号的预定着陆点为 19.51W， 44.12N，海拔为 -2641m。 嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题 是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点 15km， 远月点 100km 的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为 6 个阶段，要求满足每个阶段在关键点所处的状态，尽量减少软着陆过程的燃料消耗。 根据上述的基本要求，建立数学模型解决下面的问题： （ 1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。 （ 2）确定嫦娥三号
11、的着陆轨道和在 6 个阶段的最优控制策略。 （ 3）对设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。 二问题的分析 2.1 问题一 由于题目已经明确给出准备轨道的形状参数，可以通过已有的物理知识与几何关系明确 计算出近月点与远月点处速度大小和相对于月面的速度方向。为了预期落点与着陆轨道在同 一个平面内，且准备轨道过月心，可以大致确定轨道所在平面有无数多个，无法确定近月点和远月点的位置，因此需要其他条件来推测。本题现有的条件下，只有通过着陆轨道逆推近 月点，并结合地月轨道制动的实际情况综合考虑，才能得到。对于着陆过程一，由相关报道 以 及 NASA 在 1976 年提出的线性正切制导率
12、1，得知主减速阶段通常都是恒推力作用在轨 道切线上，且嫦娥三号主减速阶段实际也是保持着最大推力依照这一定律进行制导。我们由 此出发，通过二体模型，结合已知条件，建立微分方程组，通过计算机模拟降落轨迹即可求 出降落弧线距离，从而反推近月点，对称得出远月点。 2.2 问题二 由问题一已经得出近月点，即开始降落点位置，也知道每一阶段的状态，因此，降落轨 3 道大致范围基本确定，但六个过程的精确路径是要通过策略优化来控制的。由于燃料消耗表 现在推力在时间上的积累量，即减小推力作用的冲量，即可优化燃料消耗。 对第一个过程，由于推力很大且历时较长，因此燃料消耗主要体现在这一阶段，对应的， 优化策略也应重点
13、体现，由于有二体模型，建立微分方程模型，并由初值条件以及阶段限定 条件，可以写出非线性约束条件，本问题及转化为轨道优化中的非线性规划问题，一般可通 过成熟的 SQP 算法可以得到全局最优解，但本题采用了更为常见也相对传统的遗传算法， 逼近全局最优解，得出最优方案。 对于第二阶段，仅仅是为了是水平速度将为 0，且推力迅速减小，由上一阶段的优化结 果，得出末速度水平分量，由于冲量可分解，则此阶段分为水平方向和竖直方向分别优化， 即分为了两个变速直线运动模型，简化了优化模型，可以由这两个局部最优解加和得到该阶 段的全局最优解。 第三个阶段水平速度初始为 0，经过对月面成像分析，制定平坦度评价体系，选
14、择距离中心点最近且满足平坦度要求的区域中心为粗调整目标降落点。由于这一段终点悬停，速度 减为 0，因此可以对该段推力进行优化，从而局部燃料最优。 第四阶段悬停，精细成像并分析，同样制定平坦度评价体系并选择距离中心点最近且满 足平坦度要求的区域中心作为最终目标降落点 ，修正轨道。结束时水平速度依然为 0，因此 同样存在优化过程。 第五阶段与第六阶段是减速至 0 然后自由落体的过程，针对减速阶段也可考虑优化，可经过简单讨论得到结果。 最后根据各个阶段的最优方案，模拟出嫦娥三号着陆轨道即可。 2.3 问题三 为了分析设计轨道和控制策略的误差与敏感性，有必要制定误差指标并考虑各部分误差 对于结果的最大
15、影响，敏感性也同样需要这一思路，各个阶段的细微变化会对结果产生影响 的衡量。另一方面，或许还有必要寻找参考物以显示该方案的好坏。 三模型假设 1. 由于月球自转速度为 27.3d，十分缓慢，而着陆过程仅有十几分钟，因此在本题中月球不考虑自转； 2. 由于月球扁率很小，可认为月球为球体，半径以平均半径为准，并且引力场分布均匀； 3. 由于侧面姿态调整喷射装置对燃料影响很小，为简化模型，认为飞行器变换姿态的过程不消耗燃料； 4. 认为飞行器变换姿态是瞬间完成的； 5. 由于着陆时间较短，所以诸如月球引力非球项、日月引力摄动等影响因素均可忽略不计； 6. 推力大小可瞬间改变； 7. 燃料除供给推力外
16、，无任何其他耗散方式； 8. 月球空气稀薄，不考虑任何摩擦力； 9. 由于预定着陆点海拔 -2641m，因此在着陆过程中所使用的高度均不应是 海拔高度， 而是相对于着陆点海拔的高度。 4 四符号说明 符号 符号意义 G 引力常量 M 月球质量 m 嫦娥三号质量 比冲 gm 月球表面重力加速度 F 推力 J 燃耗 P 平坦度评价指标 H 综合评价指标 平均坡度评价指标 五模型建立与求解 5.1 问题一 5.1.1 模型建立 为了确定近月点和远月点，正常的思路是求出椭圆轨道所在平面以及椭圆长轴在空间中 的位置，但本题仅给出了轨道两点的高度信息以及平均半径，还有根据常识推断出的预定着 陆点在椭圆平面
17、内，除此之外并无其他信息，因此从现有条件准确判断两点的位置是不可能 的。因此，本题应采用逆推思路，由着陆轨道的第一个阶段反推近月点。 先对轨道参数进行分析，已知近月点高度 Hc=15km，远月点高度 Hf=100km，月球平均半径为 R=1737.013km，轨道为椭圆。如图 5.1.1.1. 图 5.1.1.1 近月轨道示意图（为了方便示意，本图不符合比例） 由开普勒定律，任何椭圆天体轨道的中心天体一定在椭圆的一个焦点上。则由椭圆几何 性质，可以近似得出方程： a+c=Hf+R a-c=Hc+R 联系万有引力定律与牛顿第二定律，可以列出嫦娥三号在近月点和远月点的运动学方程： Mm v2 G
18、0 m f (a c)2 0 f Mm v2 G 0 m c (a c)2 0 c 5 0 其中 M 7.34771022 kg 为月球质量， G=6.67210-11 为引力常量， m 2.4t 为在近月轨道 上飞行器的质量， f 与 c 为远月点和近月点的曲率半径。由椭圆的几何性质，在长轴两端点处的曲率半径分别为： c b = a c1 e a a2 f 联立以上各式，带入参数，得： = a c1 eb 近月点速度大小 vc 1692.46m / s ； 远月点速度大小 vf 1614.15m / s 。 由假设 2，月球视为球体，则近月点处速度方向方向应平行于月面且为了着陆方便，方向指向
19、预定着陆点所在方向，远月点处同样平行月面但速度方向与之相反。 下面开始研究主减速过程，该过程由近月点开始，切向速度 v 0 vc ，径向速度 vr 等于 0，由切线正切制导率及相关报道可知，该阶段推力始终保持最大推力，且始终调整使推力 与速度方向相反，即 F 7500N ， F v ，最终应大致达到预定着陆点目标上空，且竖直方向 速度为 57m/s，高度共下降 12000m。由于有端点限制条件，可以通过物理知识建立微分方程组，即采用微分方程模型进行分析求解。 由于飞行器绕月球表面飞行，且切向速度较大，在月心极坐标系下必须考虑向心力和科 里奥利力，且涉及转动变量，较为复杂。这里采用参考系转化，将
20、极坐标系转化为非惯性参 考系，并对该系中所有质点提供一个向上的离心力，大小由水平速度与质点到月心距离决定： v2 F m xd r 这样，问题就转化到了平面直角坐标系下，并且可以对速度和受力进行正交分解，如图 5.1.1.2。 图 5.1.1.2 非惯性系下力学分析图示 由此可得微分方程模型： 2 6 m0 mt r r dvy GM v2 x a sin dt r 2 r dvx a cos vx vy dt r dr v dt y a F vy tan vx 其中 r 为飞行器到月心的距离， 为推力与水平方向的夹角， a 为推力产生的加速度， 整个过程中质量是匀速减少的，减少系数由比冲与推
21、力定义，即 F m ，本题中 =2940N / kg 。 至此，本问题已经完全抽象为数学语言，并且有初始参数： vy 0 ， vx 1692.46m / s ， r 1749.37km ， =0 理论上可求得当高度下降 12000m 时速度方向大小以及飞行距离，但这个方程组要想解出一个描述运动轨迹的函数是十分困难的，因此，需要另寻它路，逼近它的解。 5.1.2 模型求解 本题由方程组反推运动方程是十分困难的，但这类问题就像是解决 NP 难题一样，无法由问题得到结论，但可以 “猜测 ”结论从而验证问题的正确性，继而得到想要的数据。 回顾本题，之所以无法解决是因为微分方程的连续性，这种连续性计算机
22、很难求解，但 若是让该组方程做某种近似，使之转化为差分方程，从而利用差分方程的离散型导出差分方 程组。输入初始条件，使用计算机模拟其运动过程，最后打点画出折线图。由微元思想可知， 当时间被微分成足够小的时候，模拟的运动轨迹可认为无限逼近真实运动轨迹。 微元化后的差分方程依然可以有运动学规律得出： v v (a cos vxn vyn ) t xn1 xn n n n GM v2vyn1 vyn ( 2 n an sin n xn ) t rn an m m F mn m t n1 n rn1 rn yn 1 2 xn1 xn vxn t 2 an cos n t yn1 yn vyn t 1
23、a 2 n sin n t 2 其 中 为时 间 微 元 小 量 ， 由于 减速 段 大 约 几 百 秒 ， 为 了精 确描 绘 运 动 曲 线 一般取 0.1s。带入初始条件： F 7500N ， m0 2400kg ， 0 =0 ， x0 0 ， y0 0 ， vx0 1692.46m / s ， vy0 0 即可通过迭代法，求出在 y=3000m 时的 x 的值，该值即为飞行器转过的角度所对应的 7 月面弧长。但是仿真结果表明，当高度降至 3000m 时，速度并不为 57m/s，而且远大于这一数值，显然，是我们采用的控制方案不合题意。 重新考虑这一过程，应该是竖直方向分力不够，导致竖直速
24、度分量增加过快。之后在了 解了嫦娥三号实际登月过程后，我们发现这一过程中推力并非与速度方向相反，而是始终与 速度反方向保持一个角度 ，并向曲线内侧倾斜。 于是我们通过改变 的大小，试图输出大量运动轨迹簇，这一过程由计算机模拟生成。通过计算机仿真模拟，得到非惯性参考系下的运动轨迹簇。 其中可以发现，在 取时，高度 3000m 处竖直方向速度为，在误差允许的范围内可以接受，由此得出符合题意的主减速阶段飞行轨迹图，如图 5.1.1.4. 图 5.1.1.4 主减速阶段飞行轨迹图（非惯性参考系） 图 5.1.1.5 水平与竖直速度变 化 示意图 图 5.1.1.6 夹角变化示意图 得到弧长为，对应转过
25、角度为，所用时间为。由于推力大小不变，由时间即可得出此阶 段燃耗。该阶段内水平方向与竖直方向 速度变化如图 5.1.1.5. 可以看出，水平速度不断下降，而竖直方向速度先增大后减小，合速度 57.12m/s，符合题意。水平总位移为 385.21km。 虽然该阶段末并不是准确到达预定着陆点上空，需要经过快速调整阶段才能准确到达， 但题目表述中明确指出该处以基本到达目标上空，再加之段末水平速度分量相比初始速度已 经十分小，快速调整阶段又很迅速，即快速调整阶段对到目标点上空距离的影响十分小，相 比主减速段弧长可以忽略，即认为主减速末段已经到达目标上空。 由此，我们结合弧长与转角可以计算近月点所在范围
26、。由于该角度可在以预定落点为圆 心的球面圆上任意选取，考虑到当天月球月面范围，降落过程应基本暴露在有光月面一侧， 以便拍摄或信息收集，并且一般选择轨道时为了方便计算控制，轨道平面会与落点所在经线 8 平面基本重合，即近月点经度也为 19.51W，纬度应更加靠近赤道，从而保证降落过程全程在光侧面，最终确定纬度位置为 31.50N 。 由月球的球对称性，远月点纬度南北对调，经度东西对调，角度互补。 最终结论，近月点位置为 19.51W ， 31.50N ，距月面 15km，速度为 1692.46m/s；远月点位置为 160.49E ， 31.50S ，距月面 100km，速度为 1612.15m/
27、s。 5.2 问题二 本题为典型的多因素优化问题，优化目标是减少燃料损耗，并且通过图像分析选择最为 平坦的地区着陆。 由于燃料损耗直接与推力有关，即 F m ，其中 =2940m / s 为燃料比冲，也就有 Fdt mdt 该式左边为推力对时间的积累量，即推力所做的冲量；右边为燃料损耗对时间的积累， 即一定时间内的燃料损耗总量。由于比冲为常数，因此推力冲量与燃耗成正比例关系，从而 使优化目标与力学量直接联系，便于分析。 另外，由动量定理 mdv Fdt ，在近似处理下，也可将燃耗与速度变化联系在一起。 5.2.1 过程一模型建立与求解 由 5.1 计算过程可知，主减速段是着陆过程用时最长，燃耗
28、最多的阶段。该阶段的主要任务是消除较大的初始水平速度 , 因此推进剂消耗优化是该阶段优化的主要设计目标，也是整个着陆轨道燃耗优化的重中之重。 嫦娥三号的主减速阶段是从近月点开始下降到距离预定地点 3 km 处，这一过程的时间比较短，仅有几百秒的时间，所以可以不考虑月球引力摄动。月球自转速度比较小 , 也可忽略 。所以，可以仅仅对月球和嫦娥三号进行分析，将问题转化为二体模型，如图 5.2.1.1。以月心为原点，建立平面直角坐标系，设嫦娥三号的月心距为 r ，极角为 ，角速度为 ，质量为 m ；又设 v 为嫦娥三号沿 r 方向上的速度， F 为主减速发动机的推力（ 在 本题中为固定值 ） ，发动机推力与当地水平线的夹角即为推力的方向 ，比冲为 ISP ；设月 球的引力常数为 。 图 5.2.1.1 二体模型力学分析示意图 首先，从运动学角度进行分析。由动力学基本方程可得 dr tdv ，所以 dr v dt 同理 ， 可 以 得 到 角 度 变 化 的关 系 式为 d dt
