高教社杯全国大学生数学建模竞赛A题

1、2014高教社杯全国大学生数学建模竞赛承诺书我们仔细阅读了全国大学生数学建模竞赛章程和全国大学生数学建模竞赛参赛规则（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则

2、的行为，我们将受到严肃处理。我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的报名参赛队号为（8位数字组成的编号）： 07033001 所属学校（请填写完整的全名）：吉林师范大学博达学院参赛队员 (打印并签名) ：1. 2. 3.指导教师或指导教师组负责人 (打印并签名)：（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。） 日期： 201

3、4 年 9 月 15 日赛区评阅编号（由赛区组委会评阅前进行编号）：2014高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：嫦娥三号软着陆轨道设计与控制策略摘要本文针对嫦娥三号软着陆轨道设计与控制策略问题，通过提取题目中的信息，利用拱点的概念、B样条函数逼近的统计定位方法、非线性规划问题及哈密尔顿函数为理论基础进行了完整的建模工作。首先，通过建立坐标系结合物理学运动公式求解出了近月点与远月点的位置及相应的速度；在此基础上，

4、利用B样条函数逼近的方法确定了嫦娥三号的着陆轨；最后通过分解着陆过程并利用非线性规划问题及哈密尔顿函数确定着陆阶段的最优控制策。针对问题一，利用拱点的概念及物理学中天体的运动的方程先求解出了嫦娥三号在近月点和远月点得速度大小分别为：1692.2m/s 、1613.9m/s及方向为运行轨道方向的切线方向，并利用复平面上建立坐标系的方法建立了极坐标，运用椭圆轨道方程进而得出了近月点与远月点的位置分别为：，针对问题二，采用B样条函数逼近的运动学统计定位方法确定了在着陆弧段上任意时刻的位置方程，从而刻画出了嫦娥三号的着陆轨道，并用matlab对轨迹进行了模拟。在6个阶段的最优控制策略上，先通过直角坐标

5、系得出质心的运动方程，再通过对6个阶段初始条件和终端状态的分解，利用非线性规划问题求解哈密尔顿函数，得出性能指标（耗燃量）的最小值为：382.6531kg，从而确定了最优控制策略。针对问题三，对于设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。通过对月球半径、陨石坑对嫦娥三号着陆所需燃料及轨道的影响、着陆过程中消耗燃料对嫦娥三号运行质量的影响三个方面进行了误差分析，合理解释了所求数据与实际数据之间的误差。又通过对问题二所采用的B样条函数逼近方法中的参数进行敏感性分析，得出模拟着陆轨道的最佳节点数。关键词：拱点复平面建立坐标系 B样条函数非线性规划哈密尔顿函数一、 问题重述1.1引言航空航

6、天技术对一个国家的综合国力发展具有重要的战略意义，一直是世界各国重点建设的内容。中国的航空航天技术创建于1956年，50年来，在中国人民的努力下取得了伟大成就，为国家的科技发展、经济建设和国家安全作出了巨大贡献。嫦娥三号是我国探月工程第二期的第二颗月球探测器，在此之前，我国已成功实施了探月的一期工程，发射了嫦娥一号，并发射了探月工程二期的嫦娥二号探测器，为嫦娥三号的成功发射打下了坚实的基础。嫦娥三号的软着陆，与嫦娥一号和嫦娥二号不同的是，嫦娥三号使用了先进大推力变推力发动机。嫦娥三号使用了专门研制的15007500N变推力发动机作为推进系统的主发动机，这台发动机不仅用于完成嫦娥三号直接奔月轨道

7、飞行期间的中途修正、近月制动和变轨任务，还用于嫦娥三号的软着陆。可见，嫦娥三号的软着陆过程非常复杂，变推力主发动机在其中起到了至关重要的作用，也是嫦娥三号研制的主要难点之一。1.2问题的提出要保证嫦娥三号能够准确的实现软着陆，其关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道从近月点至着陆点，其过程共经历6个阶段，要求满足每个阶段所在关键点的状态，并尽量减少着陆过程的燃料消耗。其着陆点大致位置：19.51W，44.12N，海拔-2641m。为了更好的找到嫦娥三号软着陆的轨道设计与控制策略，本文依次提出以下的几个问题：

8、（1） 确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小和方向。（2） 确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。（3） 对设计的着陆轨道和控制策略做出相应的误差分析和敏感性分析。二、 模型假设1.月球引力场均匀。2.月球自转对嫦娥三号的着陆无影响。3.嫦娥三号运行的轨道是以月球球心为一个焦点的椭圆。4、月球的形状扁率对嫦娥三号的着陆轨道无影响。三、 符号说明在近月点时的速度 在远月点时的速度 半主轴长 椭圆轨道离心率 月球标准重力参数值 远月点到月球球心的距离 近月点到月球球心的距离 万有引力常数月球质量 嫦娥三号运行质量四、 问题分析问题一：问题要求确定近月点与远月

9、点的位置，以及嫦娥三号相应速度大小和方向。首先，可以确定嫦娥三号近月的运动轨迹是一个椭圆，其中把月球球心作为椭圆长轴负半轴的焦点，则距离焦点较近的长轴交点为近月点，距离焦点较远的长轴交点为远月点。利用椭圆运动的天体方程，可求出嫦娥三号在近月点及远月点的相应速度。再以月球球心为坐标原点建立极坐标，通过建立常微分方程并通过matlab软件进行求解，确定了嫦娥三号在近月点及远月点的相应位置。问题二：由问题一解出的准备着陆的位置和着陆的6个阶段的运动特点，可以推算出这6个阶段关键点的运动状态及大致位置，再利用B样本函数逼近的统计定位方法得出着陆轨 迹上任意时刻探测器的位置，从而模拟出了嫦娥三号的着陆轨

10、道。在最优控制策略上，本文采用非线性规划求解月球软着陆最优控制中的两点边值问题，构造性能指标函数，由哈密尔顿函数取极小值的思想得出最优着陆曲线，进而对函数 进行积分求得最终消耗燃料，即最少燃料。问题三：通过对月球半径、陨石坑对嫦娥三号着陆所需燃料及轨道的影响、着陆过程中消耗燃料对嫦娥三号运行质量的影响三个方面进行了误差分析，合理解释了所求数据与实际数据之间的误差。同时，通过对问题二所采用的B样条函数逼近方法中的参数进行敏感性分析，得出模拟着陆轨道的最佳节点数。五、 模型的建立与求解5.1问题一的建模与解答 嫦娥三号在近月点及远月点相应速度的大小和方向 模型的建立通常，可以把月球作为椭圆的焦点，

11、把与焦点（月球）距离较近的交点叫做近月点，把与焦点（月球）距离较远的交点叫做远月点，如图1所示：1.拱点拱点在天文学上的意义是在椭圆轨道上运行的天体最接近或最远离它的引力中心，通常也就是系统的质量中心的点。最靠近引力中心的点称为近拱点（也就是本文中的近月点）而距离最远的点就称为远拱点（也就是本文中的远月点）连接近拱点和远拱点的直线称为拱点线，是椭圆的长轴，也是椭圆内最长的线段。近拱点（近月点）：在最短距离处有最快的速度即 （4）远拱点（远月点）：在最远距离上有最慢的速度即 （5）其中。2.轨道离心率 （6）在太空动力学上，一个天体的标准重力参数是万有引力常数和它质量之积即 （7） 其中：远月点

12、到月球中心的距离； 近月点到月球中心的距离。模型的求解联立（4）（6）（7）式可得近月点（近拱点）的速度： （8）联立（5）（6）（7）式可得远月点（远拱点）的速度： （9）通过matlab计算（程序详见附录）可得：于是，可得知嫦娥三号在近月点的速度为，方向是沿运动的切线方向；嫦娥三号在远月点的速度为 ，方向是沿运动切线方向。着陆准备轨道近月点与远月点的位置 模型的建立嫦娥三号经过轨道修正、减速制动、近月制动，从而进入环月轨道。环月轨道可以确定为是个椭圆轨道（如图2所示）。我们要确定近月点与远月点的位置，就必须研究嫦娥三号软着陆过程中的主减速过程。在此，我们要确定嫦娥三号着陆准备轨道近月点与远

13、月点的位置，就可以先确定嫦娥三号运动的轨迹方程。通过轨迹方程的求出，从而确定近月点与远月点的具体位置，进而得到嫦娥三号着陆准备轨道近月点与远月点的位置。设月球中心所在的位置为复平面的原点O，在时刻t，嫦娥三号位于所表示的点P。这里均为t的函数，分别表示的模和辐角。于是嫦娥三号的速度为：加速度为： （10）月球对嫦娥三号的引力依据物理学中的万有引力定律，大小为方向指向月球的中心O，故有，其中M为月球质量，m为嫦娥三号质量。 依据牛顿第二定律我们可得 （11）将公式（10）带入（11），然后比较其实部与虚部就得到：这是两个未知函数的二阶微分方程组。在确定某一卫星轨道时，需要加上定解条件。假设当t=

14、0时，卫星正处于远月点，而远月点位于正实轴上，距原点O为卫星的速度为那么就有初值条件：因此问题转化为求解带初值问题的微分方程组为：又将两边同乘以r即得： （12）其中。这样有向线段在时间内扫过的面积等于这正是开普勒第二定律，从月球指向嫦娥三号的线段在单位时间内扫过的面积相等。将（12）式代入得：于是我们可以得到嫦娥三号运行较为简单形式的数学模型：对上式进行化简得到：其中 A和是待定参数；。上式即为嫦娥三号的轨道方程，是一条平面二次曲线。由于嫦娥三号绕月球运行，故必有01.由于r在t=0时取到最大值（远月点），意味着此时函数取到最大值1，于是有： （13）从而轨迹方程为： （14） 对于嫦娥三号

15、而言，又知嫦娥三号的近月点到月球的距离为： 模型的求解 由（13）式 化简可得：利用matlab求解得（程序详见附录）：将所得数值代入（14）式中， 即可得到与P之间的关系：又因为，表示远月点到月球中心的距离，表示近月点到月球中心的距离，那么即可求得：当时，解得，则；当时，解得，则。则在远月点的位置是；近月点的位置是。根据所求位置，利用matlab做出图像（程序详见附录）：图3 嫦娥三号绕月运行的轨道曲线5.2 问题二的建模与解答嫦娥三号的软着陆过程是从近月点到月面的全过程。嫦娥三号的软着陆过程共分为6个阶段，前3个阶段主要是利用火箭的反推力制动减速，消除探测器较大的水平速度，通过姿态调整发动

16、机调整探测器姿态为竖直向下，并使探测器在到达预定高度时的速度接近0；第四阶段，悬停避障段主要是根据图像敏感器对着陆区的成像选择安全着陆点；第五阶段，缓速下降段是保证探测器在离月面约4米时速度为0，从而关闭发动机；最后一阶段探测器以自由落体方式撞机月球表面着陆。（具体阶段详见图4） 嫦娥三号的着陆轨道 模型的建立B样条函数：在曲线拟合设计中，B样条曲线主要用于要求局部交互式修改的自由曲线设计和实验数据平滑。设为给定空间的n个点，称下列参数曲线为k阶（k-1次）的B样条曲线。称折线为的控制多变形（如图5）。B样条曲线及其控制多变形图5本文采用B样条函数逼近的运动学统计定位方法对落月的轨迹进行仿真刻

17、画，对着陆轨迹确定的策略进行了分析。1. 计算原理：三阶B样条具有良好的二阶光滑度，拟合的灵活性比较强，适合拟合曲率变化大、拐曲严重的任意形状函数，而且具有较好的稳定性和收敛性。本文选择三阶标准B样条，函数形式如下： 将样本数据处理区间为，P为节点数 则：其中：。以上可以得出：该弧段上任意时刻探测器的位置便可用B样条函数表示。2. 采用B样条函数逼近方法的本质：从统计定位的思想出发，把需要定位的弧段分成若干个小区间，在区间上用B样条函数逼近，然后利用最小二乘法求解各节点的参数。可调参数P越大，样本数据区间就会分的越小，求解出的参数也会越多。 模型的求解 用matlab程序作图得到嫦娥三号的着陆

18、轨道如图6（matlab程序详见附录）。 图6 定点直线拟合与模拟仿真模型5.2.2 6个阶段的最优控制策略 模型的建立建立着陆坐标系如图7所示：取月球球心为坐标原点；oy指向着陆准备轨道的近月点；r为探测器到球心的距离；为oy与or的夹角；(t)为推力方向与垂线的夹角；F为制动火箭的推力大小，F取值为1500-7500或0。根据题意及运动轨迹得出嫦娥三号的质心运动方程： （1）其中，v为探测器在r方向上的速度；为探测器方位角的角速度；m为探测器的质量；为月球的引力常数；C为制动火箭的排气速度，是一个常值。 嫦娥三号的软着陆分为6个阶段，根据近月点的位置和着陆的轨迹，可以推算出每个阶段的初始状

19、态和终端状态：初始条件终端条件 主减速057快速调整5757粗粗避障5700精精避障000 缓速下降003.61440 自由落体3.614403.61440约束条件变量阶段表1：着陆各阶段状态在上表中，r为月球半径； （，其中）探为测器到达月面是的速度。通过matlab编程得到着陆过程中月心距与时间的关系（程序详见附录）图8 月心距变化曲线在计算过程中，由于状态变量的量级相差很大，在积分过程中会导致有效位数的丢失，所以本文采取归一化处理来提高计算的精准度，这样也可以使变量保持在同一量级。因此，做以下处理：令，及：则质心运动方程（1）可以改写为 （2）最优耗燃问题：运动方程（1）的状态方程为，其

20、中状态变量，控制变量。终端性能指标 （3）若求耗然最少，则需构造哈密尔顿函数 （4）其中，满足 （5）即： （6）终端约束为：由（6）式知，沿最优轨迹有，得出：由于不受约束，因此由得最优制导律 （7）将（7）式代入（2）、（6）式积分可得最优轨道。 模型的求解根据非线性规划解两点边值问题的方法，可以把性能指标J看作与的函数，即可把问题转化为：优化目标函数：；优化参数：，；约束：即:其中为以米/秒为单位的比冲。利用matlab编程解得:J= 382.6531kg(程序详见附录)5.3问题三5.3.1 误差分析本文在计算15km处速度时采用的月球半径为月球的平均半径，而此位置的半径应该1735.8

21、43kmR1737.646km，所以导致计算速度与实际速度有微小的偏差。本文在燃料最优计算时所用的时间为750s，而实际的时间应该比750s还要小，如：嫦娥三号在离地面4m时就关闭了发动机，之后就做自由落体运动。根据：有：（）代入数据得：实际用时要比750s小。并且随着嫦娥三号的燃料消耗，会导致其质量随时间的增加而减小，从而影响发动机的推力大小，导致计算所需的燃料与实际所用燃料的偏差。嫦娥三号进行实施软着陆中，时间作为主要影响因素。除此之外，还受到地形因素的影响。所用燃料的实际值还与嫦娥三号软着陆位置的地形有关；因为软着陆的地形会影响嫦娥三号姿态调整发动机所需燃料的多少；当遇到的陨石坑面积大时

22、，水平移动的距离就大，姿态调整发动机所需的燃料就多；当遇到的陨石坑面积小或地势平坦时，水平移动的距离就小，姿态调整发动机所需的燃料就少（见图9,、图10）。图9 嫦娥三号着陆过程中距月面2400m处月面情况变化趋势图图10 嫦娥三号着陆过程中距月面100m处月面情况变化趋势图1.像素像素是指基本原色素及其灰度的基本编码。像素是构成数码影像的基本单元，通常以像素每英寸PPI为单位来表示影像分辨率的大小。通常，一个像素被视为图像的最小完整采样。根据附录3、附录4所给图片。我们发现图片阴影部分较深的地方为陨石坑，所以对应着色较深。图片中成像较浅的部分，可认为月球表面且较平坦。把两张图片所对像素点作为

23、图片的基本单位，我们发现像素所对应单位可以做成一个的矩阵记为。然后再依据每个像素点作为单位，以每一行的像素点（地面平坦程度）变化作为周期，做出整张图片的像素变化趋势。5.3.2 敏感性分析在问题二中，确定嫦娥三号着陆轨道时用到了B样条函数，而本文只利用题中涉及的7个节点确定了着陆轨迹，事实上，当节点个数增多时，会解得更多参数，从而更利于确定轨迹。下表示为取不同节点数时得到的位置和速度的数据：表2 不同节点数所对应的位置和速度不同节点数所对应的位置和速度节点数位置RMS/m速度RMS/mRTN总计RTN总计10节点78.743.729.394.60.8490.6080.8521.34720节点6

24、8.444.129.586.61.7340.6241.6012.43730节点63.440.129.680.72.2510.8261.7772.98440节点70.948.829.991.32.5881.0562.4983.78450节点75.737.433.390.81.8381.1681.7022.763表中，RMS是均方根，也就是有效值，是一组统计数据的平方的平均值的平方根。统计学计算公式为：从这些数据可以看出，可调参数（节点）从10增加到30是，位置误差的RMS再逐渐减小，P取到30是，位置误差的RMS为80m左右，但是，继续增加节点数，P增大到50时，位置误差的RMS并没有继续减小，

25、而是增大到了90m左右。这是因为如果P增加，需要求解的参数也会怎多，而求解参数N与P的关系为,所以，P值太大，反而会影响计算效率，因此，取为宜。六、 模型的评价6.1模型的优点1. 模型一.将嫦娥三号的运行轨道方程化，可以直观的求出运行轨道上的每一点的具体坐标。将近拱点（近月点）与远拱点（远月点）公式化。试用于运动轨道为椭圆的天体之间的运动。2模型二利用B样条差值法，处理嫦娥三号软着陆过程的拟合。使试验结果更加逼真，准确。3.通过使用matlab，excel，Word作图，来处理相关数据，使抽象的具体数值变成更加详尽的图像，更加直观的观察出数据变化趋势与规律。并且，利用matlab来进行大量的

26、数据处理。6.2模型的缺点1. 模型具有局限性，不适用于运动轨道为圆形的天体间的运动。2本文所建模型，只适用于软着陆过程参考文献1 钟玉泉复变函数论第四版北京：高等教育出版社.2 昌胜骐，黄勇，宋叶志，嫦娥三号动力落月段轨迹确定策略飞行器测控学报 2014,33（3）.3 单永正，段广仁，吕世良月球探测器软着陆的最优控制光学精密工程第17卷第9期 2009年九月.4 姜启源，谢金星，叶俊. 数学模型. 第三版. 北京：高等教育出版社，20035 东北师范大学微分方程教研室常微分方程第二版北京：高等教育出版社，2012年12月 P157-1696 刘卫国 MATLAB程序设计教程第二版北京：中国

27、水利水电出版社，2005附录：1.求clearG=6.67e-011;M=7.3477e+022;ra=1837.013;rp=1752.013;a=1794.513e+003vp=sqrt(ra*G*M)/(rp*a)ansvp =1.6922e+032.求clearG=6.67e-011;M=7.3477e+022;ra=1837.013;rp=1752.013;a=1794.513e+003va=sqrt(rp*G*M)/(ra*a)ansva = 1.6139e+033.求所对应的角度clearra=1837013;rp=1752013;e=(ra-rp)/(ra+rp);p=(1-e)*ra;ra=1837013;e=(r