2014高教社杯全国大学生数学建模竞赛
编 号 专 用 页
评 阅 人 评 分 备 注 赛区评阅编号（由赛区组委会评阅前进行编号）：
赛区评阅记录（可供赛区评阅时使用）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
嫦娥三号软着陆轨道设计与控制策略
摘要
1 问题重述
1.1 背景
嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t，其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m
嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段（见附件2），要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。
1.2 问题
（1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。 （2）确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。
（3） 对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。
2 基本假设
3 问题分析、模型的建立与求解
3.1.1 问题一分析
对于问题一，需确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。文章中给定了着陆点坐标，该段中,将主要减速阶段（即从15km到3km）作为制动段。制动终点到达着陆点。 采用逆推的方法，通过着陆点坐标和建立的动力学方程，逆推出近月点的坐标以及速度的大小方向 。然后通过
角动量守恒定理，算出远月点的速度。
着陆器距离月面相对较高, 且着陆器走过的月面距离比较长, 将月球视为平面建立模型会带来较大的偏差. 因此, 制动段有必要将月球视为球体来建立均匀球体下的三维软着陆模型. 制动段推进系统采用常值推力方式, 通过姿态控制来完成制动力方向的改变. 3.1.2均匀球体三维动力学模型
首先 定义几个坐标 系: 1)参考惯性坐标系OXrYrZr 原点O位于月球中心, Zr轴由月心指向初始软着陆点, Xr轴位于环月轨道平面内且指向前进方向, Yr轴与 构成直角坐标系. 该坐标系仅用于软着陆下降轨迹和制导律设计中; 2)下降轨道参考坐标系ox0y0z0. 原点 o 位于着陆器质心, zo轴由月心指向着陆器质心为正, xo轴位于当地水平面内且指向着陆器前进方向,yo轴与xo和zo轴构成直角坐标系; 3)着陆器oxbybzb体坐标系 . 原点 o 位于着陆器质心, xb轴在制动推力矢量延长线上, 沿推力方向为正,yb,zb轴分别根据着陆器上仪器设备的安装而定, 并与xb轴构成直角坐标系. 坐标系示意图及着陆器位置与推力矢量关系如图 2 所示： .
图1 软着陆坐标系定义与推力矢量空间关系
图2(a)给出了上面各坐标系的示意和着陆器在坐标系中的位置, 图2(b)给出了F在下降轨道参考坐标系中的位置. 其中, α为在 XrYr平面内的横向月心角; β为下降轨道平面内的纵向月心角; 推力F与坐标系 ox0y0z0之间的2个推力方向角分别为推力方位角ψ和推力仰角θ, 他们定义为: 推力方位角绕正zo轴旋转为正, 推力仰角绕负 yo轴旋转为正.
分别用U, V, W表示着陆器下降速度在坐标系ox0y0z0三轴上的分量,于是有
若不考虑摄动影响且忽略月球自转, 同时引入质量方程, 可利用球坐标系与直角坐标系的关系最终得到下降轨道参考坐标系下的软着陆动力学模型
=w, =V/rsinβ， =U/r
= - +
= - -
= - +
（2）
下降轨道上分析可得α是个定值，则V应始终为0，则软着力动力方程可简化为
=w, =U/r
= - (3)
=-F/(vege)
运用MATLAB软件对上述方程组（3）进行求解，结果如下， m =2400 - (F*t)/(Isp*g)
U =(pi*Isp*t)/900 - (pi*Isp)/2 + (Isp^2*M*pi*log(F*t - Isp*M))/(900*F) - (Isp^2*M*pi*log(450*F - Isp*M))/(900*F)
W =Isp/2 - (Isp*t)/450 + (225*um)/ryue^2 - (t*um)/ryue^2 + (Isp^2*M*log(225*F - Isp*M))/(450*F) - (Isp^2*M*log(F*t - Isp*M))/(450*F) + 57/2 r =
int(((pi*Isp)/2 - (pi*Isp*t)/900 - (Isp^2*M*pi*log(F*t - Isp*M))/(900*F) + (Isp^2*M*pi*log(450*F -Isp*M))/(900*F))/((Isp*t^2)/900-t*((Isp^2*M*log(-Isp*M))/(450*F)+1700)+ (t^2*um)/(2*ryue^2) -(Isp^3*M^2*log(F*t-Isp*M))/(450*F^2)-(Isp^2*M*t)/(450*F)+(Isp^3*M^2*log(-Isp*M))/(450*2) + (Isp^2*M*t*log(F*t - Isp*M))/(450*F) - 15000), t)
（注：：dbeta为?β ，deta为β）
程序见附件程序1.
通过以上求解，我们发现完全求解软着陆的动力学方程组(2)有一定难度，并且解得的解析解不利于数值的求取和问题分析，主要表现在两方面:(1)求出的解无法进行积分，导致后面不能求出某一时刻的速度和位置;(2)积分求出的解，带入不同的初始值，解的差异很大，且都不符合事实依据。但是以上求解的切向分速度U的结果比较符合事实和推理，且数值准确，容易运算分析，其图像如下图（3）所示，因此我们将动力学方法进行简化改进得方程组（3），进而求解。