2014数学建模A题嫦娥三号软着陆轨道设计与控制策略

.2014高教社杯全国大学生数学建模竞赛 承 诺 书我们仔细阅读了全国大学生数学建模竞赛章程和全国大学生数学建模竞赛参赛规则（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的报名参赛队号为（8位数字组成的编号）： 13003024 所属学校（请填写完整的全名）： 厦门理工学院 参赛队员 (打印并签名) ：1. 刘煌 2. 江泽鹏 3. 章芳敏 指导教师或指导教师组负责人 (打印并签名)： 王琛晖 （论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。） 日期： 2014 年 9 月 14日赛区评阅编号（由赛区组委会评阅前进行编号）：.2014高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）： 全国评阅编号（由全国组委会评阅前进行编号）：. 嫦娥三号软着陆轨道设计与控制策略 摘要 嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。由于月球上没有大气，嫦娥三号无法依靠降落伞着陆，只能靠变推力发动机，才能完成中途修正、近月制动、动力下降、悬停段等软着陆任务。本文在着陆轨道设计的基本要求下，建立最优控制模型满足了每个阶段嫦娥三号在关键点所处的状态，以尽量减少软着陆过程燃料消耗的原则，完成了对三个问题的分析探究。 针对问题1，本文将嫦娥三号作为质点，根据椭圆公式求出半焦距的长度，采用适用于一切二体问题的开普勒第三定律模型，计算出近月点和远月点相应的速度大小，近月点速度为1.69204 KM/S，远月点速度为1.61390 KM/S，之后再通过给出的着陆点的方向反推出嫦娥三号相应的方向。 针对问题2，本文研究了一种应用参数化控制求解月球探测器精确定点软着陆最优控制问题的方法。首先用约束变换技术将小等式约束进行了近似处理，而后利用若十个分段的常数去逼近最优解，再根据强化技术通过时间轴上的变换，将每一段参数的持续时间转变为一组新的参数，于是最优控制问题被转化为一系列参数优化问题。最后应用经典的参数优化方法即可求得最优控制函数的一个近似解，通过增加参数个数，重复优化得到逼近连续最优解的参数化解。同时在优化过程中考虑了制动初始点的选取对结果的影响。运用matlab软件绘制着陆轨道的曲线，结果表明了所提设计方法足简单、有效的。得到最优初始点坐标为X0=83771 km，Y0=1423.9 km，Z0=586.26 km。 针对问题3，本文根据月球探测器向月飞行轨道动力学方程式得到了飞行轨道误差的迭代方程，采用协方差分析方法对轨道初始误差的误差源造成的轨道误差进行了分析，结合具体算例，给出了探测器初始轨道位置和速度误差引起的向月飞行轨道误差的时间历程和轨道终点误差。计算结果表明，在发射嫦娥三号卫星过程中，必须进行多次中途轨道修正。关键词：定点软着陆，最优控制，参数化控制，轨道误差，协方差分析方法 1 问题重述 嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t，其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m（见附件1）。 嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段（见附件2），要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。 根据上述的基本要求，请你们建立数学模型解决下面的问题： 问题一：确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。 问题二：确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。 问题三：对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。 2 模型假设1. 忽略太阳、地球对嫦娥三号卫星的引力2.将月球近视为一个质量均匀的标准球体3. 将嫦娥三号视为一个质点4. 主减速忽略动作调整所产生的燃料消耗5. 制动发动机推力可变6. 图片所给数据真实可靠 3 符号说明.近月点速度.远月点速度.惯性坐标系.月固坐标系.F发动机推力.m探测器质量.C制动火箭的比冲.月球半径.约束函数.指标函数.控制变量.调节参数.P与Ay1轴正向所成夹角.P在x1Az1平面上的投影与Ax1轴负向所成夹角.Ax1在xOz平面上的投影与Ox轴正向所成夹角.Ay1与Oy所成夹角.月球自转而产生的固坐标系相对惯性坐标系的转角 4 模型的建立与求解4.1 问题一确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。4.1.1 问题一的分析开普勒定律适用于宇宙中一切绕心的天体运动。在宏观低速天体运动领域具有普遍意义。对于高速的天体运动，开普勒定律提供了其回归低速状态的方程。根据附件1和附件2中给出的数据资料，近月点和远月点的相对距离都是已知的。本问题的模型解答可以根据椭圆公式求出半焦距的长度，然后利用适用于一切二体问题的开普勒第三定律计算出近月点和远月点相应的速度大小，之后再通过给出的着陆点的方向反推出嫦娥三号相应的方向。 4.1.2 问题一的解答 图1 开普勒第三定律模型 在图1中，A，B分别为嫦娥三号运动的近月点和远月点，以 和 分别表示嫦娥三号在该点的速度，由于速度沿轨道切线方向，可见 和 的方向均与此椭圆的长轴垂直，则嫦娥三号在此两点时对应的面积速度分别为 1 根据开普勒第二定律，应有 ，因此得 2嫦娥三号运动的总机械能E等于其动能与势能之和，则当它经过近月点和远月点时，其机械能应分别为 3 根据机械能守恒，应有 ，故得 4由24两式可解得 5 A,B的值题目已经给出 C2=A2+B2 6由56两式可解得： =1.69204 KM/S =1.61390 KM/S 图2 嫦娥三号软着陆曲线模型由图2可知，NS平面为着陆的那个平面，A点为着陆点，A点在NS的平面坐标为19.51W，44.12N，B点为嫦娥三号的近月点，根据题意近月点到月球表面的距离为15km，又有着陆点的海拔为-2641m，所以有OB的距离为17641m,因此可以得出近月点的速度方向为（19.51W，44.12N，17641），根据近月点和远月点处在同一轨道，并且共线，所以远月点的方向和近月点刚好相反，高度都为17641m，因此可得出远月点的速度方向为（19.51E，44.12S，17641）。近月点速度为1.69204 KM/S，远月点速度为1.61390 KM/S。4.2 问题二确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。4.2.1 问题二的分析 探测器经过环月轨道的软着陆方式因其具有较长的软着陆准备时间、对着陆位置的限制比较小以及减少着陆舱部分的燃料消耗等优点故而被广泛采用。该方式的关键环节就是从距离月面15 km的近月点至月面的动力下降过程。在本问题的分析中，模型考虑了月球自转，针对三维空间内精确定点软着陆问题，利用参数化控制解决了变推力软着陆最优控制问题，此外还针对仅知制动初始点到月心距离而具体位置未知的情况，对初始点的选取进行了研究。4.2.2 问题二的解答1. 动力学模型建立与控制律设计 探月飞行器首先进行霍曼变轨，从圆形环月轨道进入一条近月点高度为15 km的椭圆轨道；当到达近月点时，制动发动机点火，探测器进人动力下降段，最终以很小的相对速度(小于6 ms)降落到月面指定位置。 图3坐标示意图 如图3坐标示意图所示，定义惯性坐标系Oxyz，原点在月心，参考平面是月球赤道面，Ox轴指的是月球赤道相对于白道的升交点，Oy轴指向月球自转角速度方向，Oz轴按右手坐标系确定。再定义月固坐标系OxLyLzL, 以月球赤道面为参考面,OxL轴指向赤道面与起始子午面的交线方向，OyL指向月球自转角速度方向，OzL轴按右手坐标系确定。Ax1y1z1为原点在探测器质心的轨道坐标系，Ay1指向从月心到着陆器的延伸线方向，Ax1垂直Ay1指向运动方向，Az1按右手坐标系确定。制动发动机推力P的方向与探测器纵轴重合，为P与Ay1轴正向所成夹角，为P再x1Az1平面上的投影与Ax1轴负向所成夹角。为Ay1与Oy所成夹角，为Ax1在xOz平面上的投影与Ox轴正向所成夹角。为月球自转而产生的固坐标系相对惯性坐标系的转角，不妨假设初始时刻月固坐标系与惯性坐标系重合。 显然有轨道坐标系到惯性坐标系转换矩阵 惯性坐标系到月固坐标系的转换矩阵为 根据牛顿第二定律，结合科氏定律整理可以得到探测器在月固坐标系中的运动方程为 其中VxL,VyL和VzL为探测器速度矢量在月固坐标系各轴上的投影，F为发动机推力，m为探测器质量，gxL,gyL和gzL为该高度月球重力加速度在月固坐标系各轴上的投影，wL为月球自转角速度。 因此，在月固坐标系中探测器耳朵运动方程可表示如下（7）其中 GM为月球引力常数，C为制动火箭的比冲，是一个常值。取为系统状态变量，为控制变量，则式(7)可以简记为 X=f（x，u，t） 按照耗燃最优的要求，取性能指标为 在实际情况下，通常没必要令探测器着陆速度严格等于零，只要能保证探测器以很小的相对速度降落到月面就足可以接受的。因此，考虑到这一点，本文将软着陆的末速度要求以惩罚冈子的形式加入到指标中如下式所示，主要目的是降低最优控制问题求解的复杂度，该惩罚因子可以通过反复的数值仿真运算，按经验设定。 （8）此外，显然有约束条件（9）其中x，y，z为预定着陆点在月固坐标系中的坐标；为着陆点到月心距离，即月球半径。 对于含有形如 这类关于状态变量在连续时间上都要满足的不等式约束最优化问题，至今还是最优化领域的一个难 点。参考资料中2中给出一种约束变换技术，使得该类问题得到解决。 显然等价于 （10） 但上式显然在时不可微，因此用如下不等式去近似上式 （11） 其中 0，是调节参数。文献2证明了当足够小的时候，存在，使得对任何满足 的能够令(11)对(10)达到满足要求的近似。不妨记为用(11)式替换 后得到的新的约束函数。因此本文所讨论的软着陆耗燃最优问题转化为： 问题1在系统(7)满足约束函数的情况下，求取适当的控制变量 u 使指标函数(8)达到最小。2. 参数化控制求解耗燃最优问题 假定初始时刻为0，终端时刻tf为待定参数。 选取满足的序列和三组参数，构造形如 的参数化分段常数控制器。其中 用控制器替换系统（7）中的，则问题1转变为：问题2 寻找三组参数 ， 来最小化指标函数（8），并且满足约束函数。 显然，对于每个给定的P，这都是一个有限维的参数优化问题。文献4中第六章已经证明了当P一+时，问题2的最优解收敛于问题1的最优解。不过文献5已经证明了在数值计算中，求解问题2的参数梯度时难度很大甚至求不出真实解，因而本模型引入强化技术来解决这一问题。 从到构造如下变换 上式中0, ,序列为0,1区间上预先给定的分段点，并且满足。 将上式两边对s求导可得 其中 不妨令 则得到如下增广系统 即 （12） 其中，与， ， 分别为O,P,Q与， ，经过变换后的形式， ，；指标函数变为 （13）约束条件变为（14）其中 由于仅仅已知探测器在软着陆始点到月心的距离和探测器的初始速度以及初始质量，而软着陆起始点另两个空间位置信息角与角的初始值与未知，因而令与为系统待定参数，则系统初始状况可以表示为 （15） 那么问题2转化为如下问题： 问题3在系统(12)满足约束并且初始条件如式(15)的情况下，求取适当的控制变量u(s)使指标函数(13)达到最小。再由可知，问题3将最初的探月飞行器软着陆最优控制问题转化成了优化静态控制参数， ，以及系统参数， 的问题，利用经典的参数优化算法即可求出登月飞行器的软着陆最优控制的一组逼近解和软着陆最优初值点位置以及终端时刻。利用此算法，增加时间的分段点个数叮以重新优化，经过多次优化后即町得到满意精度的参数化解。 此外，假如令系统(7)中的推力F为已知的恒定推力，令控制变量，则本文问题变为恒定推力下软着陆最优控制问题，依然可以利用本文模型方法解决，而依据极大值原理结合传统的打靶法则只能解决恒定推力的情况，因I而相比之下本文方法适用性更广。3. 6个阶段的最优控制策略（1）问题分析由于月球表面没有大气，着陆器的速度必须完全由制动发动机抵消才能实现安全软着陆，所以减少燃料消耗是增加有效载荷的关键所在，将嫦娥三号软着陆问题，分为6个阶段依次为主减速、快速调整、粗避障、精避障、缓慢下降、自由下降。根据题意，将6个阶段分为4个阶段，依次为第一阶段（主减速和快速调整）、第二阶段（粗避障）、第三阶段(精避障)、第四阶段（缓慢下降和自由下降）。在第一阶段，利用最优控制过程推导，减少燃料消耗，这一阶段已在上一问得到解决。在第二阶段，粗避障阶段，嫦娥三号悬停在月球表面约2400米上方，对星下月表进行二维和三维成像，利用遗传算法的思想，从图像中先随机选取部分点，能直接从三维图像中得知该点的海拔高度，再分别扫描这些点附近的地貌，找出一些地势平坦的区域，我们用区域内所有点与中心点海拔的均方差作为地势判断依据之一，保留这些坐标，并进行重新组合，并改变某些坐标以便能获得其他新区域的坐标，再次搜索地势平坦的区域，重复进行多次搜索，直到没有出现崎岖地势的时候，我们将此时地势最平坦的地方作为全局最优降落地点。在第三阶段，精避障阶段，嫦娥三号悬停在月球表面约100米上方，由于在第二阶段已选定了平坦区域，也是用第二阶段采取遗传算法的思想，对采取数据进行高分辨率三维成像，对选取的区域进行更进一步的判断，选取地势更为平坦的区域作为最优降落地点。在第四阶段，卫星先是水平速度为0米/s，推力方向向下，进行缓慢降落，然后在离地面4m时，卫星将暂时处于悬停状态，关闭推进器，卫星呈自由落体降落，确保软着陆成功。（2）问题解决1. 对于第一阶段，在大约距离月球15公里到3公里，此时，反推发动机就要点火工作，着陆器的大部分燃料消耗在这个阶段。利用最优控制过程推导，减少燃料消耗。详细分析已在上问得到解决。2. 对于第二阶段，利用附件3高程图用matlab软件分析得到矩阵数据避开障碍物。该高程图的水平分辨率是1m/像素，其数值的单位是1m。例如数字高程图中第1行第1列的数值是102，则表示着陆区域最左上角的高程是102米。 图4 三维粗障碍物高程图3. 对于第三阶段，利用附件4高程图用matlab软件分析得到矩阵数据避开障碍物。该数字高程的水平分辨率为0.1m/像素，高度数值的单位是0.1m。例如数字高程图中第1行第1列的数值是100，则表示着陆区域最左上角的高程是10米。 图5三维精障碍物高程图4. 运用matlab计算结果 已知嫦娥三号在着陆准备轨道上的运行质量为2400kg，安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s由题一算出近月点的速度为：1.69204 KM/S，远月点的速度为：1.61390 KM/S。月球平均半径1737.013km、赤道平均半径1737.646km和极区半径1735.843km，月球的形状扁率为1/963.7256，月球质量是7.34771022kg。月球与地球距离最远（远地点）：406610km，最近（近地点）：356330km，平均距离为384400km。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m。利用最优控制软件miser32，通过计算机仿真运算令=0.01,=0.0025,np=30即可得到符合精度的最优解，最终利用本文的参数化控制得到软着陆末时刻tf=538.2s，末时刻探测器质量1334.88 kg，燃料消耗为1065.12 kg，最终燃料消耗为飞行器总质量的44. 38%。最后探测器以3.9811ms的对月速度精确降落到指定登月点。此外可得0=34.9860，0=35.6810,从而有最优初始点坐标为X0=83771 km，Y0=1423.9 km，Z0=586.26 km。 图6 软着陆最优控制律 图6为利用本文方法得出的最优控制律，由于最优控制律是分段常值函数因而为阶梯图；图7为软着陆速度曲线，可以看出探测器着陆时相对月面速度足够小，软着陆成功实现；图8为软着陆最优轨线，显示了角与角以及探测器距月心的距离随时间变化的曲线；图9是探测器质量变化曲线。 图7 软着陆速度曲线图8 软着陆最优轨线 图9 质量变化曲线4.3 问题三 对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。4.3.1 问题三的分析 卫星向着陆轨道动力学方程式得到了着陆轨道误差的迭代方程，采用协方差分析方法对轨道初始误差误差源造成的轨道误差进行分析，结合具体算例，给出了卫星初始轨道位置和速度误差引起的着陆轨道误差的时间历程和轨道终点误差。4.3.2 问题三的解答1.对卫星向着陆轨道的详细分析的基础上，采用协方差分析方法对其轨道初始误差源造成的轨道误差进行分析，并给出仿真算例，分析初始轨道位置误差和初始速度误差对轨道终点的位置和速度误差的影响，及从近地点轨道入轨点开始至着陆为止轨道位置的相应的轨道位置和速度总误差（3）的时间历程。2. 协方差分析方法的基本原理对于如下非线性函数关系 （16）可以使用一阶泰勒级数展开对其进行线性化，有 （17）其中，o(x1xn )为x1xn的高阶项。从而得到线性化方程 （18）或表示为 （19） 这里P 是偏导数矩阵： （20）若自变量 是随机变量，则线性化方程的函数 y的协方差矩阵为： （21） （22）式中 CX 是自变量的协方差矩阵；Cy是函数Y 的协方差矩阵。协方差矩阵中对角线元素是方差，非对角线元素为协方差。显然，只要求出传递矩阵P，便可确定源误差与欲求量误差之间的关系。若给定各种源误差，如发动机安装误差、敏感器测量误差或发动机推力和点火时间等误差时，便可以分析其对目标轨道误差的影响以及对控制系统精度的影响，进一步对各系统及元部件提出适当的精度要求。3. 计算向月飞行轨道误差的协方差迭代方程考虑到轨道参数的误差之相对于轨道参数的标称值是小量，因此可以将轨道运动方程进行线性化，从而得到能够反映轨道参数偏差量的传播关系的误差方程。在应用双二体模型且在地球影响球范围内时，对轨道运动产生摄动影响的各项，如月球引力摄动、太阳引力摄动、大气阻力摄动和太阳光压摄动等对误差方程的影响很小，因此在误差方程中将它们忽略掉。反映轨道位置和速度误差的线性化方程如下： （23）式中 ，其中 为地球引力常数。 （24）写成状态方程形式： （25）式中 （26）则式（24）变为X=F.X （27）下面推导矩阵F 的表达式： （28）式中 rx，ry 和rz 是探测器在地心惯性坐标系里的轨道位置坐标。则 （29） （30）将式（29）、（30）代入（25），得： （31） 积分式（26），得到： （32）式中 （33）取前6 阶截断，即： （34） 得到计算误差方程的迭代方程： （35）eFt 相当于式（19）中的 P 阵，由于误差方程是时变方程，因此每一步迭代都需要重新计算P 阵，计算P 阵需要利用标称轨道参数数据。进一步根据式（22），得到协方差矩阵的迭代方程： （36）4 向月飞行轨道误差的协方差分析引起轨道误差的误差源主要是导航误差，包括位置误差和速度误差。其中位置误差：r=（rx，ry，rz） ，rx，ry，rz 分别为在地心惯性坐标系中X 轴、Y 轴、Z 轴的分量。速度误差：v=(vx ，vy ，vz) ，vx ，vy ，vz 分别是在地心惯性坐标系X 轴、Y 轴、Z 轴的分量。向月飞行轨道的初始轨道位置和速度误差由运载火箭的发射入轨精度决定，若探测器在飞行途中进行轨道修正，则经过轨道修正以后的轨道位置误差将由导航误差决定，速度误差将由姿态误差和制导误差决定。上述误差决定了轨道误差协方差分析的计算初始条件，图10 给出了在不进行中途轨道修正情况下，在地心惯性坐标系里，初始轨道位置误差和初始速度误差对轨道终点的位置和速度误差的影响。项目计算一计算二计算三轨道初始误差位置误差/kmRX1.00.01.0RY1.00.01.0RZ1.00.01.0速度误差m/sVX0.01.01.0VY0.01.01.0VZ0.01.01.0轨道终点误差位置误差/kmRX7.379.5112.05RY6.258.3510.24RZ73.061.496.3速度误差m/sVX4.045.296.65VY4.326.157.71VZ0.660.720.83 图10 初始轨道位置和速度误差对轨道终点误差的影响 5 模型的优缺点1. 本文利用MATLAB处理了大量的公式计算和繁琐的坐标运算并且运用MATLAB 软件绘制大多数图形，而MATLAB软件具有较强的仿真性，功能也十分强大，精确度也很高，从这个角度可证明建立模型结果的可靠性与方法的合理性。2. 本文考虑了月球自转，针对三维空间内精确定点软着陆问题利用参数化控制解决了变推力软着陆最优控制问题，此外还针对仅知制动初始点到月心距离而具体位置未知的情况，对初始点的选取进行了研究。3. 本文在考虑月球自转及制动发动机推力可变前提下，针对月球探测器在三维空间内精确定点软着陆的实际问题，基于燃耗最优的原则，利用参数化控制将动态最优控制问题转化为静态参数最优化问题，通过数值计算，设计了软着陆最优控制律，进而得到了探测器软着陆的最优轨线，并给出了最终的燃料消耗值，此外还得到了最优制动初始点的坐标。相比较传统的打靶法而言，本文方法无需猜测毫无物理意义的共轭变量初值，既可以解决变推力软着陆问题，又适用于恒定推力的情况，同时优化速度更快。 4. 由于忽略了太阳、地球对嫦娥三号卫星的引力，导致计算会产生一定的误差。5. 由于最优控制模型一般很难准确估计次优解偏离最优解的程度，而且次优解的次优程度往往依赖于初始点的选取，导致结果稍微有点偏差。 参考文献 1 单永正，段广仁，张烽等，月球精确定点软着陆轨道设计及初始点选取，宇航学报， 30(6)，2009。2 Teo K L，Rehbock V，Jennings L S等 ，A new computational algorithm for functional inequality constrained optimization problemsJ.Auto-matica，789792，1993。3 王吉力，李俊峰，崔乃刚等；刘暾，登月琶行器软着陆轨道的遗传算法优化J，清华大学学报(自然科学版)，43(8)，10561059，2003。4 Tex,KI，Cola C J，WongKH等，A UnifiedComputationalApproachtoOptimal Control ProblemsM，England，Longman Scientific and Technical，1991。5 k H WJ，Teo K L，Jennings L S等，Control pammetrzation enhancing technique for time optimal control problemLJjDynamical Syst，6，243261，1997。6 孙宝，荣思远等，向月飞行轨道误差分析，导弹与航天运载技术，20067 苏叶果洛夫，关于向月球飞行的若干动力学问题，物理译刊: 有关人造地球卫星的运动与科学研究的若干问题，北京：科学出版社，1958:5892 8 周净扬，周获等，月球探测器软着陆精确建模及最优轨道没汁J，宇航学报， 28(6)，14621466，2007。 附录-MATLAB求解程序如下：%显示高程矩阵i=imread(附件3 距2400m处的数字高程图.tif);-MATLAB求解程序如下：%距2400m处二维和三维成像z=imread(附件3 距2400m处的数字高程图.tif);z=double(z);x=1:1:100;y=1:1:100;mesh(x,y,z);colorbar;gridxlabel(x),ylabel(y),zlabel(z)figure(2),contour(x,y,z,k);grid-MATLAB求解程序如下：%粗避障程序function Xe = roughAvoid( X, R, k )% 粗避障,对矩阵X的元素按边长为2R+1的区域进行避障搜索% 先在X中以某个点(m.,n)为中心取出(2R+1)*(2R+1)矩阵X(m-R:m+R,n-R:n+R)% 求斜率% X矩阵% R控制搜索区域边长% k最大斜率控制mX, nX=size(X);Xe=ones(mX, nX);for m=R+1:R:mX-R %以步长为R进行纵向扫描 for n=R+1:R:nX-R %以步长为R进行横向扫描 Xe(m-R:m+R,n-R:n+R) = isFlat( X(m-R:m+R,n-R:n+R), k ); endendendX=imread(附件3 距2400m处的数字高程图.tif);R=10; %扫描步长k=tan(pi/6); %斜率控制，倾斜角度不能超过pi/6 Xe = roughAvoid( X, R, k ); %趋避障函数，可得矩阵Xeimwrite(Xe,Xe.bmp,bmp); %将Xe矩阵存储成Xe.bmp文件function C = isFlat( Y, k )% 判断Y矩阵是否平坦，如果平坦返回1，如果不平坦，返回0m,n=size(Y);C=ones(size(Y);% Yl=Y(:,1);% Yr=flipud(Y(:,m);for i=1:m j=m-i+1; kc=(Y(j,m)-Y(i,1)/sqrt(j-i)2+m2); if abs(kc)k C=zeros(size(Y); return; endendfor i=1:n j=n-i+1; kc=(Y(n,j)-Y(1,i)/sqrt(j-i)2+n2); if abs(kc)k C=zeros(size(Y); return; endendend-MATLAB求解程序如下：%三维月球曲面图mstruct.origin=19.51,44.12,3;mstruct=defaultm(mstruc