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

《数学建模全国赛A题嫦娥三号软着陆轨道设计与控制策略》由会员分享，可在线阅读，更多相关《数学建模全国赛A题嫦娥三号软着陆轨道设计与控制策略（18页珍藏版）》请在装配图网上搜索。

1、A题 嫦娥三号软着陆轨道设计与控制策略摘 要人类掌握航天技术之后，探测地外天体的首选目标就是月球。我国嫦娥三号于2013年12月6日成功登陆月球，嫦娥三号是我国首次地外天体软着陆任务，由于外太空的各种因素对探测器的影响很难进行人工干扰，为了保证登月探测器在月球表面平稳降落和应对外太空的影响，本文对探测器的软着陆过程的进行了深入的研究和设计。针对问题一，本文采用逆向推理和微元分析的思想方法，从着陆点进行倒推，将每段进行微分，分析受力和运动状态，在达到6个阶段状态要求的前提下，求解出探测器的水平位移为514.8km，通过坐标变换公式得出偏转角，最终确定近月点与远月点位置为：近月点：位置（19.51

2、W,27.08N）正上方15km处，速度为，方向为探测器俯仰姿态角83.17o远月点：位置（19.51E, 152.92S）正上方100km处，速度为，方向为远月点弧的切线方向。针对问题二，在轨道设计中，本文主要考虑粗避障与精避障阶段，为了避开月球表面的大型坑洞和障碍物，本文将附录中的两幅图像都分化为100100的小区域分别模拟着陆，最终利用优选法绘制出了安全区域和软着陆轨道的图像。着陆轨道的优化是一非线性、终端时间自由且带有控制约束的最优控制问题。本文利用着陆器质心动力学方程，对其进行归一化处理，采取直接求解法，将问题转化为目标函数为燃料最省的的优化问题，运用模拟退火算法求解，得出最小燃料消

3、耗为。关键词 嫦娥三号 软着陆 轨道优化 模拟退火算法一. 问题重述嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t，其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m。嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键

4、问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段，要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。根据上述的基本要求，请你们建立数学模型解决下面的问题：1. 确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。2. 确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。3. 对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。二. 问题分析嫦娥三号是我国首次地外天体软着陆任务，在世界上首次成功实现了利用机械视觉的地外天体软着陆自

5、主避障技术。为了保证登月探测器在月球表面平稳地降落并且有效应对外太空环境下各种因素造成的干扰，需要对着陆轨道的控制方案进行深入研究和认真的设计。针对上述问题，我们进行分析：Ø 对于着陆准备轨道的近月点和远月点位置的确定，首先我们根据着陆点及着陆轨道6个阶段的要求分析，由于在软着陆过程中的各个阶段的加速度不同，我们采用逆向思想，从探测器的着落点进行倒推，再结合运动学和微积分方法来对每段进行分析，从而可以得到每段的相关参数，确定探测器从近月点到着陆点的水平位移，进而得到近月点和远月点的位置；我们根据万有引力定律和探测器的运动，可知近月点和远月点的速度大少及方向。Ø 对于嫦娥三号

6、的着陆轨道的设计，我们通过探测器的质心动力学方程和基于遗传模拟退火，将是非线性月球软着陆轨道控制问题转换为以易于处理的优化问题，进行。对于6个阶段的最优控制策略，通过分析我们主要从降落避障和燃料消耗最两方面进行优化软着陆轨道，在主减速和快速调整阶段，因推力和探测器的姿态是可变的，故采用模拟退火方法，进行优化；在粗避障和精避障阶段，主要面临是快速确定降落点并达到避障目的，我们采用matlab工具，根据查询探测器下降对月面的要求对拍摄区域进行螺旋式搜索，既达到避障，也节省时间（即减少燃料消耗）。本文流程图如图1所示：图1 本文流程图三. 模型假设1) 假设月球为不旋圆球，引力场均匀；2) 假设在软

7、着陆过程中，不考虑地球，太阳对探测器的引力作用；3) 由于月球的形状变率太小，忽略对软着陆轨道设计的影响；四. 模型建立与求解4.1. 模型一：着陆准备轨道的设计与求解4.1.1 模型准备软着陆：是指月球着陆器经地月转移到达月球附近后, 在制动系统的作用下以很小的速度近乎垂直地降落到月面上, 以保证宇航员的安全和试验设备的完好。主减速：该段中，着陆器距离月面相对较高，且着陆器走过的月面距离比较长，由主减速发动机产生1500N到7500N的可调节推力，主要任务是减速制动。快速调整：利用16台小型姿态调整发动机进行快速调整探测器姿态，主要任务是快速衔接主减速和后续的接近段，快速姿态机动到接近段入口

8、姿态，发动机推力同步减到低推力水平。粗避障：为了在较大着陆范围内剔除明显危及着陆安全的大尺度障碍，为精避障提供较好的安全点选择区域，避免近距离精避障无可避的风险。考虑到探测器运动速度较大，要求成像快，计算快。其飞行轨迹要保证成像敏感器能够持续观测预定着陆区，以及接近段飞行轨迹为满足特定姿态和下降轨迹要求接近的目标着陆点轨迹。精避障：为了在粗避障选取比较安全区域内进行精确的障碍检测，识别剔除危及安全的小尺度障碍，确保了落点安全。缓速下降：为了保证着陆月面的速度和姿态控制精确度，缓慢下降段要以较小的设定的速度匀速垂直下降，消除水平速度和加速度，保持着陆器水平位置。自由落体：关闭发动机和推力器，着陆

9、器自由下降到月面。4.1.2 模型的建立对月球软着陆主减速段、快速调整段、粗避障段、精避障段、缓速下降段和自由落体段的飞行动力学模型进行了研究, 同时基于动力学模型对各阶段制导进行了计算；在接近段（主减速和快速调整），考虑到探测器的速度大，7500N主发动机羽流带来的不可见区域为半锥角约的锥体，而成像敏感器视场为，为了避免主发动机羽流对成像敏感器的影响，且使成像敏感器的视线距离尽可能短，取成像敏感器视线偏置角。为了保证在接近段成像敏感器视场能够观测到着陆区，确定采用下降轨迹接近于水平夹角1 的直线下降方式逐步接近着陆区。我们首先对探测器处于主减速段情况进行受力分析如图2；图2 主减速阶段的受力

10、分析接着结合运动学进行分析，将其近视为抛体运动；水平方向： 垂直方向： 式中：为水平方向的力；F为主减速发动机提供的推力；m为探测器的质量，水平速度；垂直速度。图3 主减速阶段速度与时间关系图从而可得出水平位移。在快速调整过程中，通过图3的受力分析，相对于主减速阶段增加了对探测器的姿势调整，则该阶段的加速度在不断变化，通过受力分析如图3所示；图4 快速调整阶段的受力分析该过程中，将从逐步变化到，所以运用微积分思想，把的变化范围进行等份，每等份记为。即第1次细分时夹角变为，第2次细分时夹角变为那第次细分时夹角变为直到第n次变为。当足够小时，那么被细分出来的线段就相当于是一条倾斜的线段，探测器在在

11、这段线段中受力不变。再对第段线段中的运动进行分析： 从而得到水平位移 在避障段，为了保证最后的避障落点精度和节省推进剂，着陆器精确避障和下降同时进行。嫦娥三号在距离月面2.4km处对正下方月面 的范围进行拍照，在降落可以用螺旋式搜索进行障碍的筛选，从上述已求得的水平距离为909.9053km,而在水平偏移少于2300m，对结果影响小，故忽略。在精避障段对正下方月面 的范围进行拍照，与精避障方法相同，因而将其简化成垂直向下运动。对其受力分析如图4；图5 避障阶段的受力分析综合对各个阶段的分析，我们得到最终的水平位移为我们建立月心系坐标，对探测器的软着陆轨道进行进一步分析，从上述中已知探测器从近月

12、点到着陆点的水平距离，来确定探测器在软着陆过程的偏角值，从而确定近月点的位置。 图6 探测器下降轨道分段示意图通过上图分析，我们采用三角形的余弦定律；可求得 通过资料得知，嫦娥三号在变月轨道是月球的卫星，而卫星分为赤道卫星和极地卫星。从嫦娥三号的软着陆点来看，其轨道为极地卫星，从而可知近月点的位置为（19.51W,27.08N）正上方15km处，则远月点的位置为（19.51E, 152.92S）正上方100km处。我们利用牛顿的万有引力定律分别对近月点和远月点的速度大小进行求解；近月点：得; 根据对近月点出受力分析，可知速度方向为探测器俯仰姿态角83.17o远月点：我们通过角动量守恒定律进行求

13、解；根据质点系的角动量守恒定律：当质点系所受到的外力对某参考点的力矩的矢量和为零时，则质点系对该参考点的总角动量不随时间变化。得; 方向远月点弧的切线方向。4.2. 模型二：软着陆轨道的设计及6个阶段的最优控制策略4.2.1 软着陆轨道的设计建立模型时考虑月球表面没有大气层，且软着陆过程的时间较短，在几百米范围内，可以不考虑月球引力摄动，且月球自转速度比较小，也可以忽略。我们可以利用二体模型描述系统的运动。图7 月球软着陆极坐标系如图6所建立的月心极坐标系，在假设着陆轨道在纵向平面内，月心o为坐标原点，oy指向动力下降段的开始制动点，ox指向探测器的开始运动方向。则探测器在月球引力场中的动力学

14、方程2，对其进行归一化处理，以避免各状态变量的量级相差较大，寻找过程中可能会导致有效位数的丢失。则可表示为：表1 符号说明符号意义符号意义 探测器的月心距F制动发动机的推力（固定值或零）探测器的极角其比冲探测器的角速度探测器沿r方向上的速度 探测器的质量发动机推力与当地水平线的夹角，即推力方向角动力下降的初始条件由霍曼变轨后的椭圆轨道近月点确定,终端条件为着陆器在月面实现软着陆。令初始时刻,终端时刻不定,则相应的初始条件为终端约束为式中：为月球半径；为初始轨道高度；为轨道角速度。我们要对月球软着陆的最优轨道设计就要在满足上述初始条件和终端约束的条件下，调整推力大小和方向，使得探测器实现燃料最优

15、软着陆，即要求以下性能指标达到最大。可得 式中：T为发动机的推力；发动机的比冲；为月球的重力加速度。当发动机的类型已经确定时，发动机的比冲和推力幅值都已确定，即发动机的燃料消耗率已经确定。此时性能指标变为飞行时间极小，即：所以在我们在设计软着陆轨道时，将发动机燃料最优3转变为飞行时间最短问题，进行优化得到最优轨道。4.2.2 6个阶段的最优控制策略月球软着陆轨道优化问题是一非线性，终端时间自由且带有控制约束的最优控制问题。我们在进行轨道优化时需要满足减少软着陆过程的燃料消耗及降落的避障。在进行轨道优化时，最关键的地方就是对粗避障和精避障阶段路径的设计，根据题目给出的附件作出探测图像，分别对粗避

16、障和精避障阶段进行讨论；粗避障阶段：此阶段主要目的是避开大的陨石坑，实现在设计着陆点上方100m处悬停，并初步确定落月地点，已知在2400米处的图像，我们对其进行像素处理见如下图：图8 （附录1）2400km月面拍摄根据查资料知月球撞击坑内斜度一般为25° 50°, 平均35°, 坑外一般为3o8°, 平均5°. 撞击坑图像特征可归纳为:1) 在太阳照射的阳面将出现亮区域; 2) 在未照到的阴面将出现阴影区域;3) 暗区域的外边缘呈现圆弧根据以上特征，从图中可以看出，着陆器在2400米时所拍摄的陆地存在多个大型的坑洞，坑洞平均海拔约140160

17、之间，而适合着陆的地区海拔应在90110米之间，为避开坑洞，作出等高线图（见附录2）如下：图9 月面的等高图观察等高线图，为实现“避坑”的目的，应当避开等高线密集的地方（等高线稠密意味着高度变化大），将等高线稠密的地方设为危险区，着陆区从（1200，1200）处对各个区域进行螺旋式搜索。图10 满足降落的区域图图中（见附录3）有灰格的地方为危险区域，从中心点开始采用螺旋式搜索，但考虑到燃料最省的原则，将搜索到的第一个安全区域作为着陆点，由此选出的安全区域（见附录4）如下图：图11 搜索的第一个安全区域分析图像可以发现，月球表面多倾斜角，适合月球表面着陆高度应在90到110之间，此段区域内月球表

18、面斜度约为为10o左右，最高点与最低点高度差小于30米，根据坡度和安全半径的加权判断，适合选取为安全着陆区。精避障阶段：目的分析三维数字高程图，避开较大的陨石坑，确定最佳着陆地点观察着陆器在100米拍摄的照片可以发现，图中存在一个大型坑洞，多个小型坑洞，为了成功避开这些坑洞，我们将区域划为为多个小区域，采用螺旋式搜索寻找最适着陆区域：图12 分析等高图图中(见附录5)存在4个大型明显的坑洞，着陆点应当避开这些危险区域。将适合的区域分别绘图可以观察到无障碍物且适合降落的安全区域，最后选取离中点最近且无坑洞存在区域作出下图：图13 选取最佳降落区域从图片可以发现，倾斜仰角在正切值在75/100=0

19、.75以内，即小于等于30度左右，适合降落，且不存在大型障碍物。此片区域地势为扇形分布，扇形区域内存在大量的平缓区域，适合降落。为节省能源，选取与点（500,500）欧氏距离最近的点（590,300）附近为安全着落点。图14最佳降落点通过粗避障和精避障对月球表面的扫描避障，既可以达到降落避障的目的，还能减少在粗避障和精避障过程中燃料消耗。在主减速和快速调整过程，由于在推力和加速度是可变的，但有约束条件即燃料消耗最优，基于上文的分析，我们采用模拟退火方法进行优化；模拟退火是基于金属退火的机理而建立起的一种随机算法，它能够通过控制温度的变化过程来实现大范围的粗略搜索与局部的精细搜索。本文采用指数降

20、温策略对温度的变化进行控制，即式中: 为当前控制温度;为初始设定温度;为温度下降系数。首先将轨道离散化成许多小段,在各小段的节点设定待优化的参数,然后利用这些参数进行多项式拟合,从而得到整个轨道的控制曲线。将月球软着陆轨道离散化,分割成n个小段,每段的节点设定一个推力方向角,如图2所示。那么,可以令n+1个节点的推力方向角和终端时刻 作为待优化的参数。每个节点的时刻可以由下式得到：这样,就使得每个节点的推力方向角都有一个对应的节点时刻。如果假设月球软着陆的推力方向角可以表示成一个多项式,即那么,可以将节点的推力方向角与对应的节点时刻对上述多项式进行拟合,求得多项式的系数 ,进而得到整个着陆轨道

21、的推力方向角 。图15 轨道离散化可以看到,采取如上改进的函数逼近法进行参数化,所选定的待优化参数都具有明确的物理意义,从而避开对没有物理含义的待优化参数的初始猜测。同时,该参数化方法的本质是函数逼近,因此它与函数逼近法的精度相当,计算速度也比较快。为了使得到的M最小，那就不断的随机的从1500-7500N中选取适合的f。最后得到最优解f值（见附录2）。4.2.3 求解结果为了避免程序陷入局部最优解，我们通过模拟退火多次运行，得到燃料消耗量与时间的关系曲线如图16；图16 燃料消耗与时间的关系曲线从而求得最优燃料消耗J为468.25kg。五. 模型的改进在进行轨道优化中，模拟退火算法把握搜索过

22、程总体的能力较强，但局部搜索能力差，容易陷入局部最优解。遗传算法和粒子群算法具有较强的局部搜索能力，并能使搜索过程避免陷入局部最优解，但此类算法却对整个搜索空间的状况了解不多，不便于使搜索过程进入最有希望的搜索区域从而使得算法的运算效率不高。因此，如果将模拟退火与粒子群算法，遗传算法等具有较强的局部搜索能力的算法相结合，就能解决多类大型优化问题，适合于科研，工程，优化等多方面，有很强的推广性。六. 模型评价6.1. 模型的优点1. 本文设计软着陆轨道思想循序渐进，将空间复杂的曲线运动进行分解成简单运动；2. 采用微元思想，将曲线运动进行分割一小段一小段直线进行计算；3. 将非线性控制问题转化为

23、易于处理的优化问题；6.2. 模型的缺点1. 轨道设计时，由于各个阶段的约束条件无法完全考虑到，所以模拟出效果不是很好；2. 在粗避障和精避障中，探测器的移动空间变动，我们只考虑到纵向切面，对近月点的位置确定有一点影响。七. 参考文献1张洪华,梁俊,黄翔宇,赵宇,王立,关轶峰,程铭,李骥,王鹏基,于洁,袁利. 嫦娥三号自主避障软着陆控制技术J. 中国科学:技术科学,2014,06:559-568.2朱峰建,徐世杰.基于自适应模拟退火遗传算法的月球软着陆轨道优化J.航空学报,2007,04:806-812.3Meditch J S. On the problem of optimal thrus

24、t programming for a lunar soft landingJ. IEEE Transactionson Automatic Control, 1964, AC-9(4): 477-484.附录附录1a=imread('附件3 距2400m处的数字高程图.tif');a=double(a);x=1:2300;y=1:2300;meshc(x,y',a);colorbar;grid onCM=0.5,0.5,0.5;CM=repmat(CM,64,1);colormap(hsv);hold on% x=800 850;% y=1200 1300;% z=1

25、00 115;% plot3(x,y,z,'-r');xlabel('x'),ylabel('y'),zlabel('z');figure(2),contour(x,y,a,'k');grid on% shading flat附录2clear all;clc;F=1500:0.1:7500;k=length(F);m=2400;g=1.63;i=2940;s=514.8*(103);e=0.110;at=0.999;T=1;l=10000;r=0.1;v1=;z=;v=1700;h=0;n=1;for n=1:l

26、c=randperm(k); c1=c(1); f=F(1,c1); while h<=12000 t=(v*sin(r)*m)/(m*g-f*sin(r); h=(m*g-f*sin(r)*t2/(2*m); v1(1,n)=t; z(1,n)=h; h=sum(z,2); n=n+1; r=r+0.1; end vb=sum(v1,2); mi=f/i; M=mi*t; df=M-m; T=T*at; if df<0 M=M+df; elseif exp(-df/T)>=rand M=M+df; end if T<e break endendM=mi*t附录3cle

27、ar all;clc;vt=57;v0=1700;s=12000;m=2400;g=1.5;r=0.48;a=(vt*sin(r)2;a1=a*m;c=2.4194 55.6452 55.6452 147.5806 252.4194 344.3548 494.3548 594.3548;a2=a1/(2*s); a3=m*g-a2; f=a3/sin(r); b=(v0)2; b1=b-a; b2=b1*m; b3=b2/(2*f*cos(r); d=1.0e+003 *1.6816 1.3745 1.3026 0.9132 0.5553 0.2816 0.0974 0.0447; t=0:0.1:620;x=1700-(2.7719*t);y=0.057*t;x1=x.2;y1=y.2;v=sqrt(x1+y1);% plot(t,v);plot(c,d);18