1、2014 高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注嫦娥三号软着陆轨道设计与控制策略摘要在世界各国纷纷制定和实施“重返月球”的战略计划之际，我国的月球探测工程的展开，对于积极参与到月球资源的开发、维护我国对月球的权益，进而促进我国航天技术的创新与发展有着重要的意义。问题要求我们对嫦娥三号软着陆轨道与控制策略进行最优化设计。按照问题的要求，本文从以下三个方面进行了研究。针对问题一，在合理的假设基础上，利用物理理论知识、解析几何知识、微元法，且构建月心坐标系与舱体坐标系，利用坐标系之间的转换，结合最
2、优控制策略，确定了近月点和远月点的位置。且嫦娥三号绕月球的轨道是由圆形轨道变为椭圆形轨道，借助开普勒定律、能量守恒定律求解出近月点的速度。又因为我们建立了舱体坐标系，所以速度的方向也就求出来了。关键词： 坐标系转换 极大值原理 显式制导律 非线性规划一、问题重述1.1 引言嫦娥三号于 2013 年 12 月 2 日 1 时 30 分成功发射， 12 月 6 日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为 2.4t，其安装在下部的主减速发动机能够产生 1500N 到 7500N 的可调节推力，其比冲（即单位质量的推进剂产生的推力）为 2940m/s，可以满足调整速度的控制要求。在其四周安装
3、有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m。 嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆,关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点 15km，远月点 10 的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为 6 个阶段，要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。 1.2 问题的提出（1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。（2）确
4、定嫦娥三号的着陆轨道和在 6 个阶段的最优控制策略。（3）对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。二、问题分析问题一：问题要我们求着陆准备轨道近月点和远月点的位置，以及相应位置的速度大小与方向。首先我们要明确如何在月球上表示两个点的位置？速度的方向我们又应该如何表达？考虑到题目已给出我们一个着陆点 19.51W，44.12N，所以我们选择在月球上建立一个坐标系，再做一系列减小误差的措施，故上述问题就解决了。至于速度的计算，显然会与第二问的最优策略有关，我们通过对最优控制策略的计算，逆推就能算得近月点与远月点的速度。问题二：问题要我们确定嫦娥三号的着陆轨道和在六个阶段的最优
5、控制策略，嫦娥三号沿着陆准备轨道下降到距离月面一定高度时，嫦娥三号发动机点火工作，开始动力下降段。这个阶段的主要任务在于消除嫦娥三号速度的水平分量。由于着陆器在月面上软着陆只能依靠制动发动机的能量来实现，因此，在这一过程中，如何使能量最省成为所研究的关键问题，一般从以下两个方面考虑：一方面是制动火箭开始工作时机的选择；另一方面是在设计动力下降段的制导控制策略时，要在末端轨道参数满足下一阶段要求的前提下，根据燃料消耗最小的原则进行设计。对于动力下降段起始点的选择，从理论上讲，制动发动机开始工作时的高度越低，所消耗的推进剂量越少；制动发动机开始工作时的高度越高，制动行程越长，则重力损耗越大，所损耗
6、的推进剂就越多。这是因为制动发动机不仅要抵消己有的速度，还要防止下降速度在月球引力作用下进一步增长。但是，探测器下降到过于低的高度时才开始制动，则有可能制动过载太大，导致损坏仪器设备。本文对动力下降段起始点的选择不作具体研究，从兼顾节省能量和保证安全出发，将动力下降段的起始点选在自由下降段椭圆过渡轨道的近月点，其高度为 15km。问题三：问题要求我们对所设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。在这里考虑到测量误差和系统参数偏差是影响软着陆制导精度的主要原因，在存在测量误差和参数偏差的情况下实现软着陆，满足终端条件是对软着陆控制系统的一个基本要求，实际中测量误差可以通过滤波消除，对
7、闭环系统的影响不大。而系统参数如制动发动机的推力，比冲以及着陆舱的质量都是不可测的，它们是在发射以前在地面标定给出的，但在飞行过程中参数会由于一定的影响而产生一定的偏差考察这些偏差对制导过程的影响就显得十分重要。三、基本假设（1） 计算着陆准备轨道近月点和远月点的速度不考虑天体对嫦娥三号的万有引力。（2） 嫦娥三号在完成所有的着陆过程中本身固定机能不会发生改变。（ 3） 不 考 虑 空 间 飞 行 器 上 各 点 因 燃 料 消 耗 而 产 生 的 位 移 。（4） 在 对 卫 星 和 空 间 飞 行 器 进 行 轨 道 估 计 时 ， 认 为 作 用 于 其 上 的 所 有 外 力 都 通过
8、 其 质 心 。（5） 卫 星 和 空 间 飞 行 器 的 运 动 是 在 真 空 中 进 行 的 。（6） 卫 星 只 受 重 力 影 响 ， 空 间 飞 行 器 除 自 身 推 力 外 只 受 重 力 影 响， 不 考 虑 科 氏 加 速度（7） 卫 星 的 观 测 图 片 及 数 据 精 准 。四、符号及变量说明符号 符号说明V为纬度 处的卯酉圈曲率半径，5.02)sin1/(eav为坐标点的纬度为坐标点的经度E为椭球第一偏心率， 222/)(fabe推力控制角F主减速器推力eV是以米/秒为单位的比冲.m是单位时间燃料消耗的公斤数eg月球引力常量 着陆舱姿态角俯仰角偏航角 滚转角a 长半
9、轴b 短半轴G 万有引力常数五、模型的建立与求解5.1 问题一的建立与解答：5.1.1 坐标系的定义 （1）月心坐标系 zyxO,月心坐标系是和月心固连的坐标系。坐标系原点 在月心上， 垂直月球Oz赤道面指向月球北极。 轴在月球赤道面内。 轴和与其他两轴垂直并构成y x右手坐标系，如图 51 1 所示。(2) 月心极坐标系 ,rO这个坐标系是月心坐标系的极坐标形式， 为动力下降段开始时月心和rO着陆舱质心的连线，方向由月心指向着陆舱质心， 是 轴和 轴的夹角。rzO是 轴在赤道面投影和 轴的夹角。如图 51 1 所示。rOxO图 51 1 月心坐标系 与 zyxO,r（8） 参考坐标系 11,
10、rO 和 原点重合， = ，以 绕 轴旋转。如图 511,r, 1Z2。图 512 参考系 11,rO（4）月面坐标系 11,zyxO坐标系原点 在动力下降段开始时月心和着陆舱质心的连线与月面的交点，为动力下降段开始时月心和着陆舱质心的连线，方向由月心指向着陆舱质1zO心。 与 和月球交点处的经线相切， 轴和与其他两轴垂直并构成x1z 1yO右手坐标系。（5）舱体坐标系 22,zyxO坐标系的原点取在着陆舱的质心上， 轴与舱体纵轴重合，指向顶部为2y正， 轴位于舱体纵向对称面内指向上为正， 和其他两轴垂直并构成右2z xO手坐标系。舱体坐标系和舱体固连，是动坐标系。图 513 舱体坐标系（6）
11、 参考坐标系 33,zyxO坐标系原点 位于 与月面的交点处， 为月心和着陆舱质心的连线，1r3zO方向由月心指向着陆舱质心。 与 和月面交点处的经线相切， 轴与3x3 3yO其他两轴垂直并构成右手坐标系。如图 512 所示。5.1.2 坐标系之间的转换关系（1） 与极坐标系 之间的关系及其转换zyxO, ,rO和 分别为惯性坐标系 , 极坐标 , ，00,rzyxr两个坐标下的单位矢量：(511) 00 10cosincos0sin1i zyxr (512)0000 sinicoscosincoii zyxr (513) 00 cossin1i10sini rzyx 000cossinrz(
12、2) 月面坐标系 与舱体坐标系 之间的关系及其转换11,yxO22,zyxO定义着陆舱姿态角俯仰角 ，偏航角 滚转角 变换关系：先绕 以角速度1zo转 ，再绕 （过渡坐标系）以角速度 转 ，最后绕 以角速度 转 ，1xo1y图 514 12 0cosincosin01cossin0i zyxzy q(515) 12 coscsinosisncosinsicoci inii zyxzyx (516)图 514 月面坐标系与舱体坐标系之间的关系（3）月面坐标系 与动体系 之间的关系及其转11zyxO33zyxO变关系：先绕 转 ，再绕 （过渡体系）体系转 ，绕 转 。1yo0o11yo如图 515
13、。令 0图 515 月面坐标系与参考系之间的关系（3-6） 1003 cossin1ic10cosincos0sini zyxzyx 100003 coscsinsisnicoscosincoiniii zyxzyx （3-7）（4）动体系 与舱体体系 之间的关系及其转33zyxO22zyxO换定义 为推力方向角。变换关系如下：先绕 转 再绕 转, 3o),(3xo，最后绕 转 。)2(2yov图 516 参考系与舱体坐标系之间的关系（3- 32 10sincoisinco0icossin0ic zyxzyxvvvv 8） 33 sincosicosinzyxzyx 5.1.3 根据之前所建立
14、的坐标系可设月球长半轴为 a,短半轴为 b,扁率为1/f，那么：sin)1(co2hveZYX5.1.4 椭圆轨道初速度的数学表达式行星在椭圆轨道初速度的平方，等于引力常数与中心天体（为了简便本文以下均称为主星）质量及焦弦的积，除以主星中心与行星中心连线距离的平方。数学表达式为：2rGMDu式中 u 为椭圆轨道的初速度，G 为引力常数，M 为主星天体质量，d 为椭圆半焦弦，r 为主星中心与客星中心在拱点边线（即主星中心到行星近日点或远日点）距离。对椭圆轨道拱点速度数学表达式的理论推导：设：行星的质量为 m，行星椭圆轨道周期为 T，轨道半长轴为 a，轨道半短轴为b,面积速度为 ,角速度为 w，r
15、 为轨道半径。U 为轨道速度。u理论推导如下：根据面积速度公式面积速度 令：h 为 2 倍的面积速度，所以有rumu21h,其中 M 为太阳的质量；m 为行星的质量。但是一般行星质量)(423GTa远小于太阳质量（mM）故可以近似为 常数。也许当时牛顿是着眼于研究23Ta引力理论的需要才这样处理的，但是用于拱点速度求值时，出于对精度的要求就不能取近似意义。所以，必须理解开普勒第三定律的准确意义，那么上式就应该为：)(423mMGTa所以，以式上来推导，于是轨道（或拱点）速度公式为： 22)(rdu通过对坐标的描述以及转化，计算得出嫦娥三号动力下降初始状态（近月点）位置:经度 19.9989 N
16、，纬度 28.9989 W，高度 14.8km，速度 1695.7m/s，俯仰00姿态角 85 。05.2 问题二的建立与解答：5.2.1 动力学模型的化简为了得到燃耗最优制导律我们有必要对模型在一定的假设上进行化简：(1) 从 100km 左右的停泊轨道高度下降到月球表面的时间比较短，一般在几百秒的范围内，自转角速度也只有地球角速度的二十九分之一左右，因此，忽略由月球自转引起的科氏力的影响。即令 =0。(2) 由于燃耗最优，制动推力矢量、初始速度矢量以及由月球引力所产生的速度矢量这三个矢量应该在一个平面内。所以令侧滑角 ，以及滚转角 均为 0。在动力下降段开始时，经过姿态调整着陆舱主轴也处于
17、制动推力矢量、初始速度矢量以及由月球引力所产生的速度矢量所构成的平面之内。(3) 在上面两个假设的条件下，着陆舱的运动成为一个平面运动，不妨令=/ 2， = / 2。使得着陆舱的轨道在月球的赤道平面内。(4) 由于制动发动机的燃耗远大于姿控发动机，因此忽略姿态控制系统工作时对着陆舱质量所产生的影响。由以上 4 点假设，我们把一个三维坐标下的运动学方程化为平面中的运动学方程。5.2.2 性能指标的选取 选取优化的性能指标为月球探测器软着陆过程中燃料消耗最小，该性能指标等价于软着陆束 时探测器的质量最大，因此本文的性能指标可以表示为使下式最小。在建立动力学方程组中还需要补充描述着陆舱质量变化的方程
18、。 dtmvFe为发动机比冲； 为月球表面的重力加速度； 为着陆舱质量变化律。evegt( 1)约束条件 终端状态约束：为了保证探测器安全到达月球表面 ，对终端月心距有如下约束 ffrt对终端速度有如下约束： ffvt（2）控制变量约束。 在设计月球软着陆最优轨道时应限制控制量的变化范围。对推力攻角 a 有如下约束 maxin对发动机推力 F 有如下约束 maxinF FIgvmdtrdtdtrdtrrtdrtdrtvrmvreqr 1sinco00cos2sincosinis)( ssiin 11111122 121211121 111 有坐标系之间的几何关系我们容易得到， 和 之间的关系如
19、下：,r1,tr11,将上式化简，且取 为一组状态变量，可得在惯性坐Trmv标系下轨道动力学模型的状态空间描述：（521）Fgvmrvrvrvmvr err 1sincos00sin)(sin2由于忽略姿态控制系统工作时对着陆舱质量所产生的影响。在运用极大值原理求解燃耗最优轨迹时系统方程可以改写为如下形式：（522）errrgvFmvuvsinco2化为状态空间形式（523）)(1sinco002tuFgvmrvurvmvrerr 寻找一个控制律，实现燃耗最优着陆定义如下的性能指标（524）)()(00tmtdtJff 初始条件（525）TrT vrtx 00000)( 末端条件（526）),
20、(txGf 0)()(frfftvt引入拉格朗日乘子向量：（527）54321T哈密顿函数为：（528） )(/)()(),( tuvFtutxBtfuxHeTT 运用极大值原理有： x0uH（529）utxGtfTf)()(0tf对所有的容许控制 u(t)（5210）),(),(*tuxHtux同 可得控制角方程0H（5211）)(34arctg控制力矩为（5212）)(2xJM由 可得推力开关函数,(*tuHtu（5213）)(/()(), tuvFxBtxfx eTT 要 H（x, ,u）取得最大值（5214）0)(1)(tsttu不 定 )/()(eTgvtxts（5215）evFmF
21、ts )1(inco)( 543 满足共轭方程543,（5216）xH一般来说方程(1-15)很难得到解析解。上述最优着陆轨迹初值和终值（着陆点的值）已知，所以最优控制量-推力方向角以及推力开关函数可以通过求解两点边值问题得到。所以上述就是嫦娥三号的着陆轨道的最优控制策略以及相应的最优轨道。5.3 问题三的建立与解答：5.3.1 误差分析与敏感性分析测量误差和系统参数偏差是影响软着陆制导精度的主要原因，在存在测量误差和参数偏差的情况下实现软着陆，满足终端条件是对软着陆控制系统的一个基本要求，实际中测量误差可以通过滤波消除，对闭环系统的影响不大。而系统参数如制动发动机的推力，比冲以及着陆舱的质量
22、都是不可测的，它们是在发射以前在地面标定给出的，但在飞行过程中参数会由于一定的影响而产生一定的偏差考察这些偏差对制导过程的影响就显得十分重要。（a）(b)(c)图 531 推力偏差对软着陆过程的影响（a）（b）（c） 图 532 比冲偏差对软着陆过程的影响 （a）（b）（c）图 533 初始质量偏差对软着陆过程的影响显式制导律只和当前的状态信息，终端约束有关。从仿真曲线图可以看出，这些系统偏差对显式制导律的着陆效果并没有实质上的影响，都满足终端约束条件。从图 532 和图 533 可以看出：优化轨道变化平缓 ，能很好地 收敛 到终端约束值，且精度较高，说明探月器可以精确稳定地到达月球表面。整个
23、软着陆过程中的控制量变化曲线如图 531 所示，推力攻角和发动机推力变化较平滑，易于控制操作，且能很好满足控制量约束，即整个飞行 轨道是可控制的。软着陆初始阶段推力攻角一般取较大值 ，目的是为 了更快地降低探月器飞行高度。另外，仿真初始参数选取相对自由，因此 ，结果也说明了显示制导律对于初始状态量和控制量的取值不敏感，具有很好的鲁棒性。说明该方法具有一定的实时性。六、模型的科学性分析6.1 模型的优点及推广性通过仿真分析和误差分析的结果，我们可以看出这个显式制导律模型可以实现燃耗次优软着陆，而且无需迭代计算，对系统的误差不敏感，具有较好的鲁棒性。且所有的假设七、参考文献1 卢波。月球探测的势态
24、及发展意义。 国际太空, 1998, (4): 142 王景全。宋智。21 世纪国际与行星探测的发展趋势。国际太空, 1998, (4): 493 Stewart nozette 。The Clementine Mission: Past, Present, and Future。ActaAstronautica , 1995, (35): 1611694 郗晓宁, 曾国强, 任营, 赵汉元。月球探测器轨道设计。国防工业出版社,2001: 1121225 J。Wubbo Ockels, The EuroMoon 2000 Core Team. EuroMoon 2000-a plan for
25、aEuropean lunar south pole expedition 。Acta Astronautica, 1997: 146 崔乃刚, 刘墩。通过 “建立月球垂线”实现月球软着陆方法研究。导弹与航天运载技术, 2000: 45477 郗晓宁，曾国强，等。月球探测器轨道设计（M） 。长沙：国防科学技术大学出版社，2001 。 8 曾国强， 赵汉元，等月球探测器月面软着陆制动轨道研究（J） 国防科学技术大学学报， 1996 ，1 8( 4 ) ： 4O 一 43。9 徐敏 ， 李俊峰。月球探测器软着陆的最优控制（J） 清华大学学报（自然科学版） ，2001 ， 41( 8 ) ： 8 7
26、 8 9 。 10 王劫，李俊峰， 崔乃刚，等。登月飞行器软着陆轨道的遗传算法优化（ J ） 。清华大学学报（自然科学版)，2003 ，43（8）：10561059。 11 孙军伟，乔栋，等基于 S QP 方法的常推力月球软着陆轨道优化方法 E J - 1 。宇航学报。2006 ， 27(1) ：99102 。 八、附录附图 1 100m 高程图附图 1 2400m 高程图MATLAB 程序：子函数 1a=imread(附件 4 距月面 100m 处的数字高程图.tif);A=im2double(a);xa=1:1000;xb=1:1000;xx,yy=meshgrid(xa,xb);mesh(xx,yy,A);figure(gcf); 子函数 2a=imread(C:Documents and Settingsadmin桌面cumcm2014problemsA附件 3 距 2400m 处的数字高程图.tif);A=im2double(a);xa=1:2300;xb=1:2300;xx,yy=meshgrid(xa,xb);mesh(xx,yy,A);figure(gcf);
