2014全国大学生数学建模竞赛A题论文示范

仅供参考

仅供参考

嫦娥三号软着陆轨道设计与控制策略

摘要

本文针对嫦娥三号软着陆轨道设计与控制策略的实际问题，以理论力学（万有引力、开普勒定律、万能守恒定律等）和卫星力学知识为理论基础，结合微分方程和微元法，借助MATLAB软件解决了题目所要求解的问题。

针对问题（1），在合理的假设基础上，利用物理理论知识、解析几何知识和微元法，分析并求解出近月点和远月点的位置，即139.1097 。再运用能量守恒定律和相关数据，计算出速度v1（近月点的速度）=1750.78m/s,v2（远月点的速度）=1669.77m/s，，最后利用曲线的切线方程，代入点（近月点与远月点）的坐标求值，计算出方向余弦即为相应的速度方向。

针对问题（2）

关键词：模糊评判，聚类分析，流体交通量，排队论，多元非线性回归

一、问题重述

嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t，其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m（见附件1）。

嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段（见附件2），要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。

根据上述的基本要求，请你们建立数学模型解决下面的问题： （1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。

（2）确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。

仅供参考

（3）对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。

二、问题分析

2.1问题（1）的分析

首先根据问题的假设、题目中所提供的数据及图片分析，可以知道嫦

娥三号绕月球的轨道是由圆形轨道变为椭圆形轨道，借助开普勒定律、能量守恒定律求解出近月点的速度。

为了确定近月点和元月点的精确位置及相应的速度方向，我们建立以赤道（月球的赤道）平面为xoy平面、月心为原点、月心与零度经线和零度纬线交线的交点的连线为坐标轴的坐标系和赤道（月球的赤道）平面为xoy平面，为极轴（月球的极轴）为z轴建立空间直角坐标系，x轴与极坐标系的轴相重合。

首先根据着陆点的经度、纬度及月球的半径求解出着陆点和近月点（带参数 )的空间直角坐标。

其次利用两点间的距离公式，并借助MATLAB软件求解出近月点与着陆点最短距离。从而计算出 （近月点的经度）=。

最后利用卫星的轨迹是以月心为其中一个焦点，以近月点与远月点的距离为长轴的椭圆，从而求解出卫星的轨迹方程，再运用隐函数求导的应用的知识，求解出在近月点和远月点的方向导数，进而求解近月点和远月点方向余即为近月点和远月点的速度的方向。

2.2问题（2）的分析

首先在根据题意，将嫦娥三号软着陆问题，分为6个阶段依次为主减速、

快速调整、粗避障、精避障、缓慢下降、自由下降，我们先将6个阶段分为4个阶段，依次为第一阶段（主减速和快速调整）、第二阶段（粗避障）第三阶段(精避障),第四阶段（缓慢下降和自由下降）。 其次在第一阶段

粗避障阶段，嫦娥三号悬停在月球表面约2400米上方，对星下月表进行二维和三维成像，利用遗传算法的思想，从图像中先随机选取部分点，能直接从三维图像中得知该点的海拔高度，再分别扫描这些点附近的地貌，找出一些地势平坦的区域，我们用区域内所有点与中心点海拔的均方差作为地势判断依据之一，保留这些坐标， 并进行重新组合，并改变某些坐标以便能获得其他新区域的坐标，再次搜索地势平坦的区域，重复进行多次搜索，直到没有出现崎岖地势的时候， 我们将此时地势最平坦的地方作为全局最优降落地点

仅供参考

三、模型假设

1、不考虑空间飞行器上各点因燃料消耗而产生的位移；

2、在对卫星和空间飞行器进行轨道估计时，认为作用于其上的所有外力都通过其质心；

3、卫星和空间飞行器的运动是在真空中进行的；

4、卫星只受重力影响，空间飞行器除自身推力外只受重力影响； 5、卫星的观测图片及数据精准； 6、

四、变量与符号说明

C0

一条车道的基本通行能力 连续车流的车头间距 n 条车道的基本通行能力 排队长度 车流量

横断面通行能力系数车流量 持续时间

L

C y

x1 x2 x3

五、模型建立与求解

5.1 问题(1)的分析、模型建立与求解 5.1.1建模准备

（1）开普勒定律

开普勒第一定律开普勒第一定律开普勒第一定律，也称椭圆定律：每一个行星都沿各自的椭圆轨道环绕太阳，而太阳则处在椭圆的一个焦点中。开普勒第二定律开普勒定律开普勒第二定律，也称面积定律：在相等时间内，太阳和运动着的行星的连线所扫过的面积都是相等的。 这一定律实际揭示了行星绕太阳公转的角动量守恒。用公式表示为开普勒定律开普勒第

仅供参考

三定律开普勒定律开普勒第三定律，也称调和定律：各个行星绕太阳公转周期的平方和它们的椭圆轨道的半长轴的立方成正比。由这一定律不难导出：行星与太阳之间的引力与半径的平方成反比。这是牛顿的万有引力定

a3

律的一个重要基础。用公式表示为2 K开普勒定律

T

这里，是行星公转轨道半长轴，是行星公转周期，是常数 。

（2）万有引力

万有引力：任意两个质点有通过连心线方向上的力相互吸引。该引力大小与它们质量的乘积成正比与它们距离的平方成反比，与两物体的化学组成和其间介质种类无关。即：

M1M2

， r2 11

其中M1，M2为两物体的质量，G 6.67 10Nm.2kg.2（牛顿每平方米二次方千

F G

克）

5.1.2 模型的建立

根据以上的分析，建立以月球赤道平面为xOy平面,月心为原点O、Ox为月心与零度经线和零度纬线交线的交点的连线，Oz为极轴（月球的极轴），Oy与Ox和Oz满足右手标架，建立空间直角坐标系(如图5-1所示）。

图5-1 卫星绕月轨迹及软着陆轨迹

由于着陆点在球面上且近月点与远月点是由月球的经度、纬度及高度唯一确定，在此为了便于计算 将极坐标转化为空间直角坐标，并代数题中相关数据，反解出经度 。 极坐标转化为空间直角坐标

x rsin cos

即： y rsin sin

z rcos

（5.1.1）

仅供参考

x' rsin(90- ）cos(- ) '

y rsin(90- ）cos(- ) （5.1.2） z' rcos(90- ）

距离公式：

d （5.1.3）

其中： 为纬度； 为经度；r为嫦娥三号距月心的距离；d为嫦娥三号距着陆点的距离；根据能量守恒、开普勒第二定律（面积定律），建立以下模型 即：

r1v1 r2v2

（5.1.4） 1122

mv1 mgh mv2 mgH 22

则近月点的速度，近月点的速度：

v1

（5.1.5）

v 2

其中：m为卫星的质量，h1为海拔高度，h近月点距月球表面的距离；

r1 h r0 h1，r2 H r0 h1，r0月球半径，H远月点距月球表面的距离， g月球重力加速度，v1 近月点的速度，v2 近月点的速度。

5.1.3模型的求解

5.1.3.1 近月点与远月点的位置

根据题目所给数据以上分析，可知：

0,h 15000m,r0 1737013m,h1 2641m

将以上数据代入（5.1.1）式可得，着陆点及近月点的空间直角坐标分别为：

x0 r0sin(90 )cos r0sin(90 19.51)cos 44.12

y0 r0sin(90 )sin r0sin(90 19.51)sin 44.12 （5.1.6）

z rcos(90 ) r0cos(90 19.51) 00

x' rsin(90- ）cos(- )=(r0 h)cos '

y rsin(90- ）sin(- )=-(r0 h)sin z' rcos(90- ）=0

（5.1.7）

仅供参考

再将（5.1.6）式和（5.1.7）式代入（5.1.3）式可得关于 与d（近月点和着陆点距离）的函数，？利用Mathematica 5.0编程求解可得： -139.107

5.1.3.2近月点与远月点的速度大小及方向

近月点与远月点的速度方向，即为相应速度在x轴与y轴方向上的投影（如图5-2所示）

图5-2 近月点与远月点的速度方向示意图

由图易知：

5.2 模型二的建立 5.2.1模型准备 5.2.1.1系统模型

1、着陆器的动力下降段一般从15km左右的轨道高度开始，下降到月球表面的时间比较短，在几百秒范围内，所以可以不考虑月球引力摄动。月球自转速度比较小，也可忽略。因此，可以利用二体模型描述系统的运动。建立图5-2所示的着陆坐标系，并假设着陆轨道在纵向平面内，令月心为坐标原点，Oy指向动力下降段的开始制动点，Ox 指向着陆器的开始运动方向。则着陆器的质心动力学方程可描述如下:

r v

v (F/m)sin /r2 r 2

[(F/m)cos 2v ]/r ⑴

m F/ISP

式中:r, , 和m分别为着陆器的月心距、极角、角速度和质量；v为着陆器

沿r 方向上的速度；F为制动发动机的推力(固定的常值或0）；ISP为其比

为月球引力常数； 为发动机推力与当地水平线的夹角即推力方向角。冲；

仅供参考

图5-3 月球软着陆坐标系

动力下降的初始条件由霍曼变轨后的椭圆轨道近月点确定，终端条件为

着陆器在月面实现软着陆。令初始时刻t0 0，终端时刻tf不定，则相应的初始条件为

r0

终端约束为

rf rL,vf 0, f 0 ⑶

rL h0,v0 0, 0 o ⑵

式中:rL为月球半径；h0为初始轨道高度； o为轨道角速度。

月球软着陆的最优轨道设计就是要在满足上述初始条件和终端约束的前提下，调整推力大小和方向9使得着陆器实现燃料最优软着陆，即要求以下性能指标达最大。

J mdt

tf

5.2.1.2模型归一化

在轨道优化过程中，由于各状态变量的量级相差较大，寻优过程中可能会导致有效位数的丢失。通过归一化处理可以克服这一缺点[9]，提高。计算精度。令rref r0,mtef

m0，则 r/rref, v/vref,vref ISp I

仅供参考

2 F/Fref,Fref mrefvref/rref, m/mref, t/tref

, rref/vref, 。那么，着陆器的动力学方程可改为：

v

22 (F/m)sin /

[(F/)cos 2]/ F/ISP

相应的初始条件和终端约束变为：

1, 0, 000/

f r1/r0,vf 0, f 0

性能指标改写为： 0

5.2.4模型评判

由以上计算可知，多车道道路通行能力从中心至边缘车道依次递减.视频一中撞车位置在距道路中心一、二条车道上，因而可行车道为第三条车道；视频二中撞车位置在距道路中心二、三条车道上，因而可行车道为第一条车道.而从计算中可得，可上述结论，即视频二的事故所处断面实际通行能力要比视频一要强.这与实际情况比较吻合。 5.3 问题（3）的模型建立与求解

当上游交通需求量大于事发路段现有的通行量，到达车流在事故地点陆续减慢速度甚至停车而集结成密度较高的队列，事故接触后，由于数段通过能力的恢复，排队车辆又陆续加速而疏散成一列具有适当密度的车队，排队消散完毕后，车流就会恢复顺畅的交通状态.

f

仅供参考

5.3.1模型推导

由车流波动理论可知，波速公式为：

WX,Y (QX QY)/(KX KY).

式中： Wx，y 为集散波的波速，Km/h；

（1）

QX、QY为前后两种车流状态的流量，辆/h； KX、KY为前后两种车流状态的密度，辆/Km.

根据交通流模型可知，交通量Q、行车速度v、车流密度K三者的关系为：

Q v K. （2） 速度-密度线性关系模型：

v vf(1 K/Kj）.

（3）

式中：vf为畅行速度，即车流密度为零时,车辆的最大速度；

Kj为阻塞密度，即车流密集到所有车辆无法移动时的密度；

由以上（1）、（2）、（3）式可以推导出波速与密度的关系：

K/K

WX,Y vf(1 x

y）.

Kj

（4）

5.3.2模型的建立与求解

事故发生后排队长及消散时间的计算

图为事故发生后累计车辆-时间图，实线表示交通需求流量，点划线表示通过能力.为叙述简便，对所有符号说明如下：事故发生时堵塞了部分车道，该路段通行能力下降S1；相应密度上升Ks1；交通事故处理所需时间为Ｔ事故解除后到车队消散前通行能０；

力回升为S2；车流密度相应地下降为Ks2.其中路段的通行能力由图２中点划线的斜率来表示.路段上游交通需求流量为Q1、Q2…….由图２中实线斜率表示；持续时间为

T1、T2；相应车流密度为K1、K2…….

在图中，由车流波动理论可知， 波速公式为：

WX,Y (QX QY)/(KX KY).

首先假设两波相遇之前该路段需求量始终未变，OA与CB相交处表示排队向上游的延伸达到的最远处，设两波相遇时的时间为T，集结波波速为W12，消散波波速为W2,3，

仅供参考

则根据两波相遇时波传动的距离相等这一关系可知：

W12 T W23 (T T0). (5） 其中W12=|

K Ks1Q1 S1

| |Vf(1 1)|;

K1 Ks1Kf

K Ks2S S

W23=|12| |Vf(1 s1)|.

Ks1 Ks2Kf

KS1 KS2kx

) (T1 (s2 T0).KfS2 Q

则: T=x Vf(1

若T> T1，则说明在车队消散之前该路段上游要求流量发生了变化，需求流量变为Q2，

相应的密度为K2.所以（5）式改写为

W T W (T T) W (T T). (6)

12

1

34

1

23

其中W34 则T

K KS1Q2 S1

Vf(1 2

K2 KS1Kf

(K1 KSI KS2） T0 (K2 K1) T1

.

K1 KS2

KS1 KS2

) (T T0))KS2/S2可解出本次事故引起的排队Kf

根据公式；T1 T (Vf(1 长.

由资料可知车队消散时间为:

T1 T x/V3.

其中V3为路段通行能力为S2时的行车速度，

vm S2/KS2;

WX,Y (QX QY)/(KX KY);

Q v K.

所以，根据以上推导可以得到排队长度（x）关于横断面实际通行能力（S2）、拥堵持续时间（消散时间T1）、路段上游车流量(Q)的关系式如下：

x Vf(1

KS1 KS2kx

) (T1 (s2 T0). Kfs2 Q

5.4 问题（4）的模型建立与求解

5.4.1数值分析确立各个参数值

在此问中我们采用非线性回归法：非线性回归是指因变量y对于 1... m（不是自变量）是非线性. MATLAB统计工具箱的命令nlinfit等不仅可以给出拟合的回归系数及其置信区间，而且可以给出预测值及置信区间等.

多元线性回归分析模型

:

仅供参考

其中都是与无关的未知参数，，其中

为y的观测值，

称为线性回归系数.分别为

现在得到n个独立观测数据的观测值，

由上式得

:

在进行非线性回归的时候，道路通行能力以道路通行系数表示. 我们选取了视频1情况下的六次排队情况（由于每次持续时间较短，因此可以假设排队长度为120，）道

第一步：以回归系数和自变量为输入变量，将要拟合的模型写成匿名函数:y= ( 1 ( 1 2x1)/x2) ( 3 x3)；

第二步：用nlinfit计算回归系数，用nlparci 计算回归系数的置信区间，用nlpredci计算预测值及其置信区间；

第三步：确立各参数，得出 1、 2、 3的参考值为[0.0001,1.0243,-0.3253] 所以y= (00001. (0.0001 1.0243x1)/x2) (-0.3253 x3).

5.4.2模型建立与求解

基于第三问的模型建立求解过程： 由第三问我们可以得出：y= (00001. (0.0001 1.0243x1)/x2) (-0.3253 x3). 在该方程中，有三个变量，分别为道路上游车流量、通行系数以及堵车时间，假设从开始排队到排队位置到达上游路口的时间内上游车流量保持不变，即为1500pcu/h，假设该车流为连续流，于是设定x2 为 0.5 . 解该等式：可得出，堵车时间间隔大约为8分钟，即从事故发生开始，经过大约8分钟后排队长度将达到上游

六、模型的分析与评价

在模型一中，综合考虑了驾驶时间、车身长度、车道折算系数等多种因素，使得计算出来的畅通率更加的正确. 但是我们为了简化模型，我们没有考虑到模型的岔道对车流的影响.

在问题二的模型中，用模糊评判法对附件1和附件2两种状态下的通行能力进行定性的研究；用聚类分析模型对附件1和附件2两种状态下的通行能力进行定量上的分析，理论严谨.

在问题三中，通过建立基于车流波动理论的交通流模型，推导过程思路清晰. 但是

仅供参考

我们是把车流运输当做连续的流体处理，这个因素使得结果有些误差.

在问题四中的模型中，通过问题三所得关系式，进行参数分析，通过多元非线性回归，得出了数值表达式，进而求出最短时间，思路明朗. 但在模拟的多元非线性自身也存在一些误差.

七、模型的改进

由于交通问题是一个复杂的系统，可以考虑进行基于元胞自动机的交通流计算机模拟研究. 在问题一中，可以进一步把统计的区间段进行细分，进行差值拟合，是图像更加接近于真实情况；在问题二中，同理可以进行细分；在问题三中，引入相位变换等，可以考虑添加修正系数，获取相关数据，并用BP神经网络进行训练.

八、模型的应用与推广

为了更好的反映车道被占用对城市道路通行能力的影响，问题三的模型给出了车辆排队长度与车流量、间隔时间、通行能力的关系.然而在实际问题中，还有许多非确定性因素. 且车辆排队长度与车流量、间隔时间、通行能力以及其他非确定性因素都存在相互影响的关系，经过相关数据提取技术，在获取相关数据后，可以结合片最小二乘回归分析进行处理：

第一步：提取两变量组的第一对成分，使之相关性达到最大； 第二步：建立回归方程； 第三步：残差替换；

第四步：进行偏最小二陈回归分析法； 第五步：交叉有效性检验.

参考文献：

[1] 姜启源．数学模型（M）．第三版．北京：高等教育出版社，2003．8. [2] 张德丰．MATLAB数值计算方法（M）．北京：机械工业出版社，2010．1. [3] 汪晓银，周保平．数学建模与数学实验（M）．北京：科学出版社，2010．2.

[4] 孔惠惠，秦超，李新波，李引珍，交通事故引起的排队长度及消散时间的估算，第二十七卷，

第5期，2004-12-23

[5] 石磊等，《数学建模优秀论文》，北京：清华大学出版社，2011年9月第1版 [6] 姜启源，《数学模型（第二版）》，北京：高等教育出版社，1993年8月第2版 [7] 王炜，过秀成等编著．交通工程学．南京：东南大学出版社，2000 [8] 杨肇夏．计算机模拟及其应用．北京：中国铁道出版社，1999

[9] 谌红．模糊数学在国民经济中的应用．武汉：华中理工大学出版社，1994

仅供参考

附录

附录一： 程序

1．问题二中模型一：模糊评判模型 程序一（源代码：p1.m） clc, clear

a=[0.1 0.3 0.4 0.1 0.1 0.2 0.4 0.4 0 0 0.1 0.3 0.5 0.1 0 0 0.1 0.1 0.4 0.4 0 0 0.1 0.4 0.5 0.2 0.4 0.3 0.1 0 ];

w=[0.6 0.4];

w1=[0.5 0.25 0.25]; w2=[0.5 0.25 0.25]; b(1,:)=w1*a([1:3],:); b(2,:)=w2*a([4:6],:); c=w*b ;

程序二（源代码：p2.m） clc, clear

a=[0.1 0.3 0.4 0.2 0 0 0.4 0.2 0.2 0.2 0.2 0.3 0.4 0.1 0 0.1 0.1 0.1 0.4 0.3 0.1 0.1 0.1 0.4 0.3 0.2 0.4 0.2 0.1 0.1 ];

w=[0.6 0.4];

w1=[0.25 0.25 0.5]; w2=[0.25 0.25 0.5]; b(1,:)=w1*a([1:3],:); b(2,:)=w2*a([4:6],:); c=w*b

2．问题二中模型二：聚类分析

视频1程序代码（源代码：moxinger.m） clc,clear load hua.txt

gj=zscore(hua); %数据标准化

y=pdist(hua); %求对象间的欧氏距离,每行是一个对象

仅供参考

z=linkage(y,'average'); %按类平均法聚类 h=dendrogram(z); %画聚类图

set(h,'Color','k','LineWidth',1.3) %把聚类图线的颜色改成黑色，线宽加粗 for k=3:5

fprintf('划分成%d类的结果如下：\n',k)

T=cluster(z,'maxclust',k); %把样本点划分成k类 for i=1:k

tm=find(T==i); %求第i类的对象

tm=reshape(tm,1,length(tm)); %变成行向量

fprintf('第%d类的有%s\n',i,int2str(tm)); %显示分类结果 end if k==5 break end

fprintf('**********************************\n'); end

视频二程序代码（源代码：p4.m）： clc,clear load hu.txt

gj=zscore(hu); %数据标准化

y=pdist(hu); %求对象间的欧氏距离,每行是一个对象 z=linkage(y,'average'); %按类平均法聚类 h=dendrogram(z); %画聚类图

set(h,'Color','k','LineWidth',1.3) %把聚类图线的颜色改成黑色，线宽加粗 for k=3:5

fprintf('划分成%d类的结果如下：\n',k)

T=cluster(z,'maxclust',k); %把样本点划分成k类 for i=1:k

tm=find(T==i); %求第i类的对象

tm=reshape(tm,1,length(tm)); %变成行向量

fprintf('第%d类的有%s\n',i,int2str(tm)); %显示分类结果

仅供参考

end if k==5 break end

fprintf('**********************************\n'); end

3．问题四中模型：非线性回归法（源代码：p5.m）

xy0=[120 1200 0.6 10 120 930 0.5 300 120 840 0.4 172 120 760 0.5 65 120 1100 0.3 60 120 900 0.6 90 120 890 0.7 560 ];

x=xy0(:,[2:4]); y=xy0(:,1);

huaxue=@(beta,x)

(beta(1)-x(:,1).*beta(1).*beta(2)./x(:,2)).*(beta(3).*x(:,1)-x(:,3)); beta0=[0.05,0.005,0.02]';

[beta,r,j]=nlinfit(x,y,huaxue,beta0) betaci=nlparci(beta,r,'jacobian',j)

[yhat,delta]=nlpredci(huaxue,x,beta,r,'jacobian',j)

运行结果： beta =

0.0001 1.0243 -0.3253 r =

28.4765 -8.0898 13.1246 65.7630 -55.1525 54.3847 -3.4178 j =

1.0e+006 *

0.8056 0.0001 -0.0003 1.1274 0.0001 -0.0002 0.9407 0.0001 -0.0002 0.4774 0.0001 -0.0001 1.5416 0.0002 -0.0005

仅供参考

0.5775 0.0001 -0.0002 1.0863 0.0001 -0.0001

betaci =

1.0e+003 *

-0.0005 0.0005 -4.2144 4.2165 -0.0015 0.0008

yhat =

91.5235 128.0898 106.8754 54.2370 175.1525 65.6153 123.4178

delta =

1.0e+004 *

1.3405 1.8908 1.5770 0.7937 2.5870 0.9604

1.8150

本文来源：https://www.bwwdw.com/article/6hlj.html