A数学建模优秀论文.docx
《A数学建模优秀论文.docx》由会员分享，可在线阅读，更多相关《A数学建模优秀论文.docx（63页珍藏版）》请在冰豆网上搜索。
A数学建模优秀论文
2014A数学建模优秀论文
2014高教社杯全国大学生数学建模竞赛
承诺书
我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。
如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。
我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）：
A我们的报名参赛队号为（8位数字组成的编号）：
10009072所属学校（请填写完整的全名）：
东南大学参赛队员（打印并签名）：
1.吉张鹤轩
2.杨升
3.陈同广指导教师或指导教师组负责人（打印并签名）：
（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。
以
上内容请仔细核对，提交后将不再允许做任何修改。
如填写错误，论文可能被取消评奖资格。
）
日期：
2014年09月15日
赛区评阅编号（由赛区组委会评阅前进行编号）：
2014高教社杯全国大学生数学建模竞赛
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：
赛区评阅记录（可供赛区评阅时使用）：
评
阅人
评
分
备注
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
对于第二问，依据题意将其分为五个阶段。
对于前两个阶段，建立天体力学分析中常用的二体模型及力学方程组，并对约束条件进行归一化处理。
为了后续优化，列出全部归一化约束条件并建立目标函数，从而获得非线性规划模型。
之后采用了序列化遗传算法对其进行筛选，从而逼近该阶段全局最优解，最终此阶段燃耗为1055.39kg。
第二阶段同样属于复杂多变量优化问题，同阶段一建立优化模型并通过遗传算法求解该段全局最优解，最终燃耗为
26.71kg。
第三第四阶段主要在于图像处理与统计，第三阶段，通过对高程图进行K均值聚类分析，将图中像素点分为安全点与危险点，在对地图栅格化，并对方格内点类型统计取整，对方格二元化为安全格与危险格，在通过扩大寻找最大安全半径，综合考虑水平偏移量，建立合理的落点评价体系，最终找出最优点坐标（1275,1000），燃耗为86.97kg。
第四阶段同样做聚类分析，并根据嫦娥三号实际体积选取合适的栅格大小，并对栅格通过最小二乘法对空间进行线性统计回归，求出平均坡面与平均坡度，结合最大安全半径建立最优落点评价体系，最终获得最优落点坐标为（88,56），燃耗为20.68kg。
第五阶段最优燃耗为8.09kg。
最后还讨论了简单运动的局部最优模型，简化了后几个阶段的运动学分析与计算。
最终综合各段最优解，获得最优着陆轨道与控制策略。
对于第三问，首先总结了优化模型中引入的一些误差因素，并针对主要因素做了数值上的相对误差分析，证明了误差对于优化方案并未产生很大影响。
其次从初始变量和约束条件入手，分析了这些变量的波动对于结果产生的影响，最终发现角度控制向量的灵敏性较高，而其他因素的灵敏性普遍处在较低水平，侧面说明了优化方案的对于初值的不敏感性与方案对于全局最优解的逼近程度较高。
关键词：
非线性规划模型序列化遗传算法K均值聚类空间线性回归二体模型
一．问题的提出
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问题二
由问题一已经得出近月点，即开始降落点位置，也知道每一阶段的状态，因此，降落轨
道大致范围基本确定，但六个过程的精确路径是要通过策略优化来控制的。
由于燃料消耗表现在推力在时间上的积累量，即减小推力作用的冲量，即可优化燃料消耗。
对第一个过程，由于推力很大且历时较长，因此燃料消耗主要体现在这一阶段，对应的，优化策略也应重点体现，由于有二体模型，建立微分方程模型，并由初值条件以及阶段限定条件，可以写出非线性约束条件，本问题及转化为轨道优化中的非线性规划问题，一般可通过成熟的SQP算法可以得到全局最优解，但本题采用了更为常见也相对传统的遗传算法，逼近全局最优解，得出最优方案。
对于第二阶段，仅仅是为了是水平速度将为0，且推力迅速减小，由上一阶段的优化结果，得出末速度水平分量，由于冲量可分解，则此阶段分为水平方向和竖直方向分别优化，即分为了两个变速直线运动模型，简化了优化模型，可以由这两个局部最优解加和得到该阶段的全局最优解。
第三个阶段水平速度初始为0，经过对月面成像分析，制定平坦度评价体系，选择距离中心点最近且满足平坦度要求的区域中心为粗调整目标降落点。
由于这一段终点悬停，速度减为0，因此可以对该段推力进行优化，从而局部燃料最优。
第四阶段悬停，精细成像并分析，同样制定平坦度评价体系并选择距离中心点最近且满足平坦度要求的区域中心作为最终目标降落点，修正轨道。
结束时水平速度依然为0，因此同样存在优化过程。
第五阶段与第六阶段是减速至0然后自由落体的过程，针对减速阶段也可考虑优化，可经过简单讨论得到结果。
最后根据各个阶段的最优方案，模拟出嫦娥三号着陆轨道即可。
2.3问题三
为了分析设计轨道和控制策略的误差与敏感性，有必要制定误差指标并考虑各部分误差对于结果的最大影响，敏感性也同样需要这一思路，各个阶段的细微变化会对结果产生影响的衡量。
另一方面，或许还有必要寻找参考物以显示该方案的好坏。
三．模型假设
1.由于月球自转速度为27.3d，十分缓慢，而着陆过程仅有十几分钟，因此在本题中月球不考虑自转；
2.由于月球扁率很小，可认为月球为球体，半径以平均半径为准，并且引力场分布均匀；
3.由于侧面姿态调整喷射装置对燃料影响很小，为简化模型，认为飞行器变换姿态的过程不消耗燃料；
4.认为飞行器变换姿态是瞬间完成的；
5.由于着陆时间较短，所以诸如月球引力非球项、日月引力摄动等影响因素均可忽略不计；
6.推力大小可瞬间改变；
7.燃料除供给推力外，无任何其他耗散方式；
8.月球空气稀薄，不考虑任何摩擦力；
9.由于预定着陆点海拔-2641m，因此在着陆过程中所使用的高度均不应是海拔高度，而是相对于着陆点海拔的高度。
四．符号说明
符号
符号意义
G
引力常量
M
月球质量
m
嫦娥三号质量
ν
比冲
gm
月球表面重力加速度
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+Ra-c=Hc+R
联系万有引力定律与牛顿第二定律，可以列出嫦娥三号在近月点和远月点的运动学方程：
⎧Mmv2
⎪G0=mf
⎪（a+c）20
⎨f
⎪Mmv2
G0=mc
⎪（a-c）20
0
其中M=7.3477⨯1022kg为月球质量，G=6.672⨯10-11为引力常量，m=2.4t为在近月轨道
上飞行器的质量，ρf与ρc为远月点和近月点的曲率半径。
由椭圆的几何性质，在长轴两端点处的曲率半径分别为：
⎧
⎪⎪c
⎨
=b=（a-c）（1+e）a
a2
⎪ρ=
⎪⎩f
联立以上各式，带入参数，得：
=（a+c）（1-e）
b
近月点速度大小vc=1692.46m/s；远月点速度大小vf=1614.15m/s。
由假设2，月球视为球体，则近月点处速度方向方向应平行于月面且为了着陆方便，方向指向预定着陆点所在方向，远月点处同样平行月面但速度方向与之相反。
下面开始研究主减速过程，该过程由近月点开始，切向速度vθ（0）=vc，径向速度vr等于
0，由切线正切制导率及相关报道可知，该阶段推力始终保持最大推力，且始终调整使推力
与速度方向相反，即F≡7500N，Fv，最终应大致达到预定着陆点目标上空，且竖直方向
速度为57m/s，高度共下降12000m。
由于有端点限制条件，可以通过物理知识建立微分方程组，即采用微分方程模型进行分析求解。
由于飞行器绕月球表面飞行，且切向速度较大，在月心极坐标系下必须考虑向心力和科里奥利力，且涉及转动变量，较为复杂。
这里采用参考系转化，将极坐标系转化为非惯性参考系，并对该系中所有质点提供一个向上的离心力，大小由水平速度与质点到月心距离决定：
v2
F=mxdr
这样，问题就转化到了平面直角坐标系下，并且可以对速度和受力进行正交分解，如图
5.1.1.2。
图5.1.1.2非惯性系下力学分析图示
由此可得微分方程模型：
⎧dvyGMv2
⎪=-+x+a⋅sinβ
⎪dtr2r
⎪dvx=a⋅cosβ-vxvy
⎪dtr
⎪dr=v
⎨dty
⎪
⎪a=F
⎪
⎪vy=tanβ
⎩vx
其中r为飞行器到月心的距离，β为推力与水平方向的夹角，a为推力产生的加速度，整个过程中质量是匀速减少的，减少系数由比冲与推力定义，即F=νm，本题中ν=2940N/kg。
至此，本问题已经完全抽象为数学语言，并且有初始参数：
vy=0，vx=1692.46m/s，r=1749.37km，β=0
理论上可求得当高度下降12000m时速度方向大小以及飞行距离，但这个方程组要想解出一个描述运动轨迹的函数是十分困难的，因此，需要另寻它路，逼近它的解。
5.1.2模型求解
本题由方程组反推运动方程是十分困难的，但这类问题就像是解决NP难题一样，无法由问题得到结论，但可以“猜测”结论从而验证问题的正确性，继而得到想要的数据。
回顾本题，之所以无法解决是因为微分方程的连续性，这种连续性计算机很难求解，但若是让该组方程做某种近似，使之转化为差分方程，从而利用差分方程的离散型导出差分方程组。
输入初始条件，使用计算机模拟其运动过程，最后打点画出折线图。
由微元思想可知，当时间被微分成足够小的时候，模拟的运动轨迹可认为无限逼近真实运动轨迹。
微元化后的差分方程依然可以有运动学规律得出：
⎧v=v
-（a
⋅cosβ
+
vxnvyn）⋅t
⎪xn+1
⎪
⎪
xnnn
n
r
GMv2
⎪vyn+1
⎪
=vyn+（2
n
-
an⋅sinβn
-
xn）⋅trn
⎪
⎪an=
⎪⎪
⎨m=m
F
mn
-m⋅t
⎪n+1n
⎪rn+1=rn-yn
⎪12
⎪xn+1=xn+vxn⋅t-
⎪
2
ancosβn⋅t
⎪
⎪yn+1
⎪
⎪⎩
=yn+vyn
⋅t-1a
2
nsinβn
⋅t2
其中
为时间微元小量，由于减速段大约几百秒，为了精确描绘运动曲线
一般取0.1s。
带入初始条件：
F=7500N，m0=2400kg，β0=0，x0=0，y0=0，vx0=1692.46m/s，vy0=0
即可通过迭代法，求出在y=3000m时的x的值，该值即为飞行器转过的角度所对应的
月面弧长。
但是仿真结果表明，当高度降至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。
虽然该阶段末并不是准确到达预定着陆点上空，需要经过快速调整阶段才能准确到达，但题目表述中明确指出该处以基本到达目标上空，再加之段末水平速度分量相比初始速度已经十分小，快速调整阶段又很迅速，即快速调整阶段对到目标点上空距离的影响十分小，相比主减速段弧长可以忽略，即认为主减速末段已经到达目标上空。
由此，我们结合弧长与转角可以计算近月点所在范围。
由于该角度可在以预定落点为圆心的球面圆上任意选取，考虑到当天月球月面范围，降落过程应基本暴露在有光月面一侧，以便拍摄或信息收集，并且一般选择轨道时为了方便计算控制，轨道平面会与落点所在经线
平面基本重合，即近月点经度也为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问题二
本题为典型的多因素优化问题，优化目标是减少燃料损耗，并且通过图像分析选择最为平坦的地区着陆。
由于燃料损耗直接与推力有关，即F=νm，其中ν=2940m/s为燃料比冲，也就有
⎰Fdt=ν⎰mdt
该式左边为推力对时间的积累量，即推力所做的冲量；右边为燃料损耗对时间的积累，即一定时间内的燃料损耗总量。
由于比冲为常数，因此推力冲量与燃耗成正比例关系，从而使优化目标与力学量直接联系，便于分析。
另外，由动量定理mdv=Fdt，在近似处理下，也可将燃耗与速度变化联系在一起。
5.2.1过程一模型建立与求解
由5.1计算过程可知，主减速段是着陆过程用时最长，燃耗最多的阶段。
该阶段的主要任务是消除较大的初始水平速度,因此推进剂消耗优化是该阶段优化的主要设计目标，也是整个着陆轨道燃耗优化的重中之重。
嫦娥三号的主减速阶段是从近月点开始下降到距离预定地点3km处，这一过程的时间比较短，仅有几百秒的时间，所以可以不考虑月球引力摄动。
月球自转速度比较小,也可忽略。
所以，可以仅仅对月球和嫦娥三号进行分析，将问题转化为二体模型，如图
5.2.1.1。
以月心为原点，建立平面直角坐标系，设嫦娥三号的月心距为r，极角为θ，角速度为ω，质量为m；又设v为嫦娥三号沿r方向上的速度，F为主减速发动机的推力（在本题中为固定值），发动机推力与当地水平线的夹角即为推力的方向ϕ，比冲为ISP；设月球的引力常数为μ。
图5.2.1.1二体模型力学分析示意图
首先，从运动学角度进行分析。
由动力学基本方程可得dr=tdv，所以dr=v①
dt
同理，可以得到角度变化的关系式为dθ=ω②
dt
图5.2.1.2受力分析示意图
再从动力学角度进行分析。
如图，对嫦娥三号进行受力分析可知。
在径向上，受到万有引力F
GMm
μm
=月=以及主减速发动机在径向上的分力Fsinϕ。
引r2r2
由此，根据牛顿第二定律可以得到：
μmmv2
③
r2r
Fsinϕ=mdv④
dt
2
m，利用公式vr化简该式可得：
r2rdt
dv=Fsinϕ-μ+rω2⑤
dtmr2
再对切向上的受力进行分析，此处需要注意到我们选择参考系月球自身的自转，是一个非惯性参考系。
因此，嫦娥三号除了受到主减速发动机推动力的分量Fcosϕ的作用，还需要计入大小为2vω的科里奥利力，根据右手定则可以确定其方向为切向。
由牛顿运动学
第二定律可以得到Fcosϕ+2vω=-mdvy。
同样，利用v
dty
=rω化简可以得到：
F*cosϕ+2*v*ω
dω=-m⑥.
dtr
由于飞船在运行过程利用燃料反冲制动，因此飞船得而质量是不断变化的。
考虑比冲的定义为“火箭发动机单位质量推进剂产生的冲量”可得：
dm=-F⑦。
dtISP
综上所述，联立式①，②，⑤，⑥，⑦即可得到完整的描述嫦娥三号与月球这一二体运动模型的方程组：
dr=vdt
dv=Fsinϕ-μ+rω2
dtmr2
dθ=ω
dt
Fcosϕ+2vω
dω=-m
dtr
dm=-FdtISP
对于嫦娥三号燃料使用指标的衡量，我们可以反映到其动量上去，即将动量作为衡量燃
t
fF
料消耗的性能指标。
由质量变化的推导，我们可以得到J=⎰
t0SP
dt。
在轨道优化过程中,由于各状态变量的量级相差较大,寻优过程中可能会导致有效位数的丢失[2]。
归一化处理可以克服这一缺点，提高计算精度。
另外，由于对轨道的优化也要求优化变量尽可能地保持在相同的量级，故作以下处理：
令：
rref
=r0
mref=m0.
rvF
mv2m
则：
r=,v=,v
=,I=I
F=,F
=refref
m=,
rref
vref
ref
SPSP
Fref
ref
rref
mref
r
ω=ω
t=
ttref
tref
=ref
vref
θ=θ.
那么，嫦娥三号的动力学方程可改写为：
dr=vdt
dv=Fsinϕ-1+rω2
dtmr2
dθ=ω
dt
Fcosϕ+2vω
dω=-m
dtr
dm=-FdtISP
又由第一问可以得到飞行器的初
升级会员 