2014高教社杯全国大学生数学建模竞赛
承诺书
我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛
参赛规则》(以下简称为“竞赛章程和参赛规则”,可从全国大学生数学建模竞赛网站
下载)。
我们完全明白,在竞赛开始后参赛队员不能以任何方式(包括电话、电子邮件、
网上咨询等)与队外的任何人(包括指导教师)研究、讨论与赛题有关的问题。
我们知道,抄袭别人的成果是违反竞赛章程和参赛规则的,如果引用别人的成果
或其他公开的资料(包括网上查到的资料),必须按照规定的参考文献的表述方式在正
文引用处和参考文献中明确列出。
我们郑重承诺,严格遵守竞赛章程和参赛规则,以保证竞赛的公正、公平性。如
有违反竞赛章程和参赛规则的行为,我们将受到严肃处理。
我们授权全国大学生数学建模竞赛组委会,可将我们的论文以任何形式进行公开
展示(包括进行网上公示,在书籍、期刊和其他媒体进行正式或非正式发表等)。
我们参赛选择的题号是(从A/B/C/D中选择一项填写): A 我们的报名参赛队号为(8位数字组成的编号): 19005007 所属学校(请填写完整的全名):
参赛队员 (打印并签名) :1.
2.
3.
指导教师或指导教师组负责人 (打印并签名):
(论文纸质版与电子版中的以上信息必须一致,只是电子版中无需签名。以上内容请仔细核对,提交后将不再允许做任何修改。如填写错误,论文可能被取消评奖资格。)
日期: 2014 年 9 月 14 日
2014高教社杯全国大学生数学建模竞赛
编号专用页
赛区评阅编号(由赛区组委会评阅前进行编号):
全国统一编号(由赛区组委会送交全国前编号):全国评阅编号(由全国组委会评阅前进行编号):
嫦娥三号软着陆轨道设计与控制策略
摘要
本文主要研究在嫦娥三号高速飞行的情况下,准确降落在月球预定区域的软着陆轨道设计与控制策略问题:
针对问题一,查资料可知嫦娥三号软着陆准备轨道是极月轨道,忽略月球自转公转的影响,则着陆准备轨道与软着陆轨道在一个平面上,近月点投影的经度为19.51W 。然后根据问题二得出的结果求出水平位移X ,再列出微分方程方程并求其数值解,最终得到近月点投影的纬度,进而根据近月点投影的经纬度算出远月点投影的经纬度,至此便得到近月点与远月点的位置;对于二者的速度,应用开普勒定律和机械能守恒定律可求出。最终得到近月点:59.07N , 19.51W , 距地高度15km ,31.69210m/s A V ≈?,方向为北极-虹湾-南极方向;远月点:59.07S , 160.49E ,距地高度100km ,
31.61410/B V m s ≈?,方向为南极-北极-虹湾方向。
针对问题二,在安全着陆的前提下,确定着陆轨道和最优控制策略,即转化为燃料消耗最小,最终转化为飞行时间tf 最短,经计算最优时间为633s ,消耗燃料的最小值为645.9Kg 。我们首先列出软着陆轨道动力学方程并做归一化处理,经过软着陆轨道的离散化,应用函数逼近方法拟合推力控制角β(t),拟合结果为
7342(t)7.4710 6.14100.026 1.72;x x x β--=-?+?-+。从而将轨道优化问题转化为参数优
化问题,最后利用MATLAB 随机寻优,进而得出轨道优化设计方案;。
针对问题三,我们着重对主减速阶段进行相应的误差分析和敏感性分析。建立月球软着陆主减速阶段的初始状态误差模型,敏感性大小S 是由f P -
里面的每一个元素除以i P -
里面对应的元素得到的值的集合,六个值分别对应v 、r 、θ、w 、m 、β,f P -
为主减速阶段末状态总误差向量,i
P -
为初始状态偏差向量。得到结果为( 3.862015.9200
0.1540
582.0000 1.0000
65.6890)S =--,由这六个值大小分析可
知w 最敏感,即w 有误差时会导致结果误差最大;其次是β,r ,v ,m ,θ。于是我们应该尽力避免w 初始值有误差的情况。
关键词:软着陆 动力学方程 微分方程 随机寻优 寻优变量
1 问题重述
嫦娥三号软着陆轨道设计与控制策略
嫦娥三号于2013年12月2日1时30分成功发射,12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t,其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力,其比冲(即单位质量的推进剂产生的推力)为2940m/s,可以满足调整速度的控制要求。在四周安装有姿态调整发动机,在给定主减速发动机的推力方向后,能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W,44.12N,海拔为-2641m。
嫦娥三号在高速飞行的情况下,要保证准确地在月球预定区域内实现软着陆,关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求:着陆准备轨道为近月点15km,远月点100km的椭圆形轨道;着陆轨道为从近月点至着陆点,其软着陆过程共分为6个阶段,要求满足每个阶段在关键点所处的状态;尽量减少软着陆过程的燃料消耗。
根据上述的基本要求,请你们建立数学模型解决下面的问题:
(1)确定着陆准备轨道近月点和远月点的位置,以及嫦娥三号相应速度的大小与方向。
(2)确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。
(3)对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析
2 模型假设
1、假设嫦娥3号主减速阶段完成时的速度为57m/s且垂直于地面,并且不存在水平速
度;
2、假设近月点投影到月面的地区与虹湾之间的月面弧线可以看作一条直线;
3、假设日月引力摄动、月球旋转和实际月球为非标准球导致引力不均匀对嫦娥3号的
运动无影响;
4、假设推力F=3000N的大小保持恒定;
5、假设收集到的数据准确无误。
3 符号说明
::v r r θβ着陆器沿半径方向的速度;月心距;
:极角;
w:角速度;m:质量;
:制动发动机的推力方向角
4 问题分析
2.1问题一:
月球自转和公转周期皆为27天7小时43分11.559秒,而嫦娥三号软着陆过程只有几百秒,所以我们可以忽略由于自转和公转造成软着陆轨道平面与着陆准备轨道平面的偏差,即可认为软着陆轨道平面与着陆准备轨道平面是在同一个平面上的。已知着陆准备轨道为极月轨道,则近月点经度就可确定,如果得到近月点投影与虹湾水平距离X ,就可算出两地夹角,那么近月点纬度便随之确定,而X 的值可由问题二得出的各物理量算得。远月点与近月点关于月心对称,由此可算出远月点位置。对于速度,使用开普勒定律和机械能守恒定律则很容易算出。 2.2问题二:
确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略,对于着陆轨道我们可以得出月心距随时间变化的图像,对于最优控制策略,我们可以得到其他参数(径向速度,角速度等)随时间变化的图像。优化的指标是在保证安全和可控条件下使得燃料消耗最少(也即着陆时间tf 最短),于是我们把问题简化为求解两点边值问题中的最优解。因此我们采用随机寻优的思想对待优参数进行优化,得到了各参数的最优解,再进行拟合函数。对于避障问题,我们可以对处理后的数据划分成网格,用网格体现所有数据的特征,在经过坡度pn 和平移指数A(i)判断平移方向和距离。 2.3问题三:
对设计的着陆轨道和控制策略做相应的误差分析和敏感性分析,易知在避障阶段及以后是无法进行误差和敏感性分析的,于是我们只需对主减速阶段分析。参考资料可知误差源包括初始条件误差和导航与控制传感器误差等,为简化模型,我们只分析初始条件误差。然后建立初始状态误差模型,通过将000000,,,,,v r w m θβ?????? 的值
(其值均考虑为典型误差值)分别与00000,,,,v r w m θ和0β(即初始状态着陆器沿r 方向的速度、月心距、极角、角速度、质量和制动发电机推力方向角)相加,然后得到非标准末状态和标准末状态的差值
f f f f f ,,,,,v r w m θβ
??????,再让它们与
000000,,,,,v r w m θβ??????对应相除得到S ,以S 中每个值的大小来描述其对应的变量
初始状态对该变量的末状态的影响即敏感性分析与误差分析。
5 模型的建立与求解
5.1问题一:求近月点与远月点的位置及嫦娥三号相应的速度大小和方向
5.1.1、求近月点远月点位置:查资料可知嫦娥三号着陆准备轨道是极月轨道【3】,
即通过月球南北极的轨道,其方向为北极-虹湾(着陆点)-南极,如图1所示:
则我们可知近月点的经度和虹湾一样为19.51W ,近月点到虹湾的水平距离X 为(,,v x x r v a 分别是切向
的速度、加速度和径向的速度,此处θ不是指极角):
00tan (t)f f
r x
r r t x x t x a a dv a dt v a dt X v dt
β?
=??
?
=
????
=???=?
?? 00((cot ))f f t t r dv X dt dt dt β=?? 解得54.51810X =?m
然后我们由余弦定理公式:222
cosA 2b c a bc
+-=,其中22222,b ,a A c R X θ====,
(虹湾(着陆点)海拔虽然为-2641m ,但与月球半径相比数值很小,所以我们可以忽略它,认为R 取月球平均半径1.737013×106m )可解得θ=14.95度,则得到近月点纬度Q=θ+44.12=59.07N ;
远月点与近月点连线过月心,则其纬度大小与近月点纬度大小相等,但在南半球,所以其纬度为59.07S ,其经度与近日点经度相差180度,则可以算得为180-19.51=160.49E 。
5.1.2、求近月点A 距地15km 和远月点B 距地100km 的速度VA 和VB :开普勒定律适用于宇宙中一切绕心的天体运动,自然嫦娥三号在着陆准备轨道上的绕月运动也
图1
自转方向
在此类,于是我们可以使用能量法[1]来计算V A 和V B :
如图2所示,由开普勒第一定律:每一个行
星都沿各自的椭圆轨道环绕太阳,而太阳则处在椭圆的一个焦点中。可知其中f 是月心即椭圆的
焦距,月球质量M=7.3477×1022kg ,月球平均半径r=1.737013×106m ,A 到f 的距离L A =(1737.013+15)×103=1.752013×106(忽略虹湾
与月平均半径的差2641m ),B 到f 的距离L B =1737.013+100=1.837013×106。在A,B 两点分别取?t →0,则嫦娥三号与月球的连线在这段时间内扫过的面积分别为:
11
t ,t 22
A A A
B B B S V L S V L ?=??=?
根据开普勒第二定律:在相等时间内,太阳和运动中的行星的连线(向量半径)所扫过的面积都是相等的。可得A B S S ?=?,带入得 A
B A
B
L V V L =(1)。 嫦娥三号运动的总机械能等于其动能和引力势能之和,故其经过A 、B 两点时的机械能为:
1(),(2)21(),(3)
2A A A
B B B
GMm E mV L GMm E mV L =
+-=+-,
(以无穷远处为引力势能零点, 11226.6710/kg G N m -=?)
由于嫦娥三号在着陆准备轨道上只受万有引力作用,所以遵循机械能守恒定律,即有:A B E E =,(4)
然后联立(1)(2)(3)(4)得到:
A B V V =
=
最后带入数据算得:
331.69210m/s, 1.61410/A B V V m s ≈?≈?
综上所述得出:近月点:59.07N , 19.51W , 距地高度15km ,31.69210m/s A V ≈?,
方向为北极-虹湾-南极方向;
远月点:59.07S , 160.49E ,距地高度100km ,31.61410/B V m s ≈?,方向为南极-北极-虹湾方向。
X
5.2问题二:求嫦娥三号的着陆轨道及六个阶段最优控制策略
嫦娥三号在15公里高度下降到月面只需要几百秒,所以月球引力非球项、日月引力摄动和月球旋转皆可忽略。嫦娥三号在着陆准备轨道近月点处时制动发动机点火,以抵消其的初始动能和势能,从而使它水平速度被抵消,最终以相对月面速度为零的垂直姿态落到月面。由于在距离月面3000m 处时嫦娥三号已经基本在目标点上方,所以可假设此处不存在水平速度,只有垂直向下的速度57m/s 。
5.2.1、我们假设着陆轨道在纵向平面内,并建立如图3所示的轨道平面坐标系,其中坐标原点在月心0,Y 轴指向着陆准备轨道近
月点,X 轴指向嫦娥三号在近月点开始运动的方
向。 软着陆的动力学方程描述如式子(1):
22sin cos 2e dv F u
rw dt
m r dr v dt d w dt
dw F vw dt mr r dm F dt v βθ
β?=
-+???=???=???=--??
?=-??
(1)【2】F 是制动发动机推力,其值恒定,e v 为比冲=2940m/s ,u 为月球引力常数=4.9028×1012m 3/s 2,β为制动发动机推力方向角,其定义为F 与当地水平方向的夹角。
由资料可知近日点曲率半径R A = 2
b a
(图一:a 为半长轴= (A B L L +)/2,b 2为半
短轴=a 2
-c 2
= 2
2()22A B A B B L L L L L ++??
-- ??
?)
【1】,则近月点的曲率半径也应为R A ,由向心力公式:22A A A
v m mw R R =可得A A v
w R ==9.431×10-4rad/s 。
月球软着陆轨道服从两点边值约束,近月点为起点,对于此处的我们认为近月点
轨道高度15公里是对海拔为-2641m 的月面而言,则r 0=L A +海拔=1.752013×106-2641=1.749372×106m ,
6+ 1.7374010+8m f r =?=月平均半径海拔距地高度3000,所以满足的初始条件为:
640000 1.74937210,0,w 9.431/r 0,1/0m r v m ad ra s d s θ-?==?== (2)
X
终端约束条件为软着陆,需满足: 【2】
6f , 1.73740810m ,w 0/57/=2f f f v rad s m ad s r r π
θ?===, (3)
待优化变量为制动发动机的推力方向角β(t )和终端时刻t f ,优化的性能指标为满足初始条件和终端约束条件的前提下,主减速阶段燃料消耗最少,即为:
'0min (t)dt f
t J m =? (4) 【2】
因为在发动机推力大小可调变的前提下得出的软着陆最优制动方案,实际上是一个推力大小恒定,推力与速度夹角β变化的制动过程,因此式子(4)描述的性能指标也相当于使得终端时刻f t 最短。 5.2.2、数值归一化处理:
为避免由于各变量的量级相差较大导致计算中有效数位的丢失,从而提高计算精度,我们需要对状态变量进行归一化处理:
令200m v ,m ,v F ,ref ref ref ref ref ref
ref ref ref ref
r r r m t r v =====, 【2】
则:
,,,
,,e ref ref ref
ref ref
v r F
v r v v F v r F m t m w t m t θθ-
----
---======== 【2】
经过归一化处理后,初始条件和终端约束条件改写为:
初始条件:---
0001,0,r v w w === (5)
终端条件:-
-
-
,0f f f f r v r v w r =
=
= (6)
动力学方程改写为:'22
''''sin 1cos 2e F v r w m r r v w F v w w m r r F m v βθβ-
---------
----
---
--
-
?
??=-+??
?=???=???=-??
??=???
(7)【2】
'
(t)d f
t J m t -
-
-
-
=? (8)
5.2.3、参数化方法:
由于优化变量β(t )的搜索空间是一个泛函空间,必须将控制变量参数化,本
文采用函数逼近的方法。函数逼近法假设推力方向角可用一个多项式(9)表示:
230123t t t βλλλλ+++(t )= (9) 【4】
如图4所示,函数逼近法首先将软着陆轨道做离散化,即将轨道分为N 段,然后
在每一段的节点处设一个推力方向
(i 0,1,2...
i β=,最后利用这些推力方向做多项式拟合,从而求出多项式系数i λ,
5.2.4、算法的求解:
传统的优化算法(如遗传算法,蚁群算法)虽然能够精确的计算出最优值,但由于它们对数据的寻优是逐个进行的,算法极其复杂繁琐,计算的速度也十分缓慢,因此我们寻找一种速度较快、优化性能较好的折中算法。对于繁琐的数据搜索,我们选
择每次搜索一组数据,这便大大减小了数据搜索的难度,同时我们也可以知道只要保
证优化次数足够大便能得到合理的优化结果,因此对传统算法的优化是合理的。 随机寻优步骤:
(1) 确定寻优空间。即根据实际需要确定待优化上下界,此题中为16个方向控制角和终端时刻的上下界,上界H=[0,10,20,30,40,50,60,70,80,90,90,90,700] ,下界L=[0,0,0,0,0,0,0,0,0,0,0,0,500] (2) 设置优化次数M ,此题中M=500
(3) 算法的应用实现。在寻优空间中随机产生一组随机数据X=L+(H-L )*rand,其中rand 是区间(0,1)之间的均匀随机数。
(4) 计算每组数据的适应度,比较适应度的大小,根据优化条件选择适应度的最值所对应的参数,模型中得到的最优参数为:
[0,4.5,6.0,16.0,25.0,25.0,25.0,25.0,25.0,86.0,86.0,86.0,633.0] 5.2.5、优化策略的求解:
通过对β(t )的拟合我们可以得到推力方向角β(t )与时间的关系,运用MATLAB 画图可得其图像。同理,我们也可以得到其他参数与时间的关系如下: 7342734243261137244
133122(t)7.4710 6.14100.026 1.72;
3.6610 1.85100.219 6.01;1.86100.1651
4.8 1.7510;
6.7310 4.0109.4810 1.3810;
w 2.2110 5.12107.4410r x x x v x x x r x x x x x x x x βθ-----------=-?+?-+=?+?-+=?-++?=?-?+?-?=?-?-?741931629.4410;1.9310 1.110 1.022400.0x m x x x ----+?=-?+?-+
(其他物理量的图
像见附页)
5.2.6、优化轨道的求解:
至于轨道的求解,我们认为求出月心距随时间变化的表达式即为轨道的求解根据上文的求解可知:
5.2.7、主减速阶段后的优化:
经过上一目的计算,我们可得出主减速阶段的最优控制策略,即各个待优参数(径
向速度,角速度,极角,推力方向角等)随时间变化的函数关系,同时也可以得到月
心距随时间变化的函数,即最优轨道。又有资料显示,嫦娥3号的燃料消耗很大部分
来自主减速阶段,因此我们可忽略快速调整阶段和避障阶段的燃料消耗。由附件2可知,嫦娥3号在主减速阶段完成后进入快速调整阶段,此后将近似沿直线垂直下降。
附件4和附件5给出了嫦娥3号在2400米和100米高度下拍色的虹湾的高程图,通MATLAB软件中imread函数可读入数据,但数据格式为uinnt8无法直接运算,故通过im2double函数将数据转化为浮点型,进而进行运算分析。
数据转化后我们分别得到高度矩阵a和b,下面我们以矩阵a为例进行分析。我们将2300*2300的矩阵a分为50*46个网格点,每个网格点含有46*50个数据。由于网格点足够密集,对任一网格点我们采用mean函数求其平均值,并用这一平均值代替该网格点的特征,最终2300*2300的矩阵a便可用50*46的矩阵x来表示,并且各网格点的相对位置不变,最后我们利用mesh函数绘出三维高度图如下:
由图可知,待降落区有明显的起伏,当拍摄图像显示待降落点为起伏点时,必须对降落轨道进行调整以达到避开高低和低洼,应立即启动姿态调整发动机对嫦娥3号进行平移,最终安全降落在预定地点。下面我们讨论如何进行平移 1. 计算平移平移点附件各个方向的坡度p1,p2,p3……pn。
2. 根据各坡度的大小进行加权获得各个方向的平移指数a(1),a(2),a(3)……a(n),坡
度p 越大,平移指数a 越小。
3. 由某一方向的综合平移指数Ai=a(i)-[a(1)+a(2)+……+a(i -1)+a(i+1)+……+a(n)]
来确定平移的方向,A(i)的最大值对应的方向即为平移的方向。
4. 对于平移距离我们采用在满足条件的情况下,平移距离最小的原则,因为这样可以
最大限度的节约燃料。 5.3问题三:对设计的着陆轨道和控制策略做相应的误差分析和敏感性分析 5.3.1、初始状态误差模型:
我们所得的轨道数据为标准,记嫦娥三号标准初始状态为0x -
,实际初始状态为0x ,则定义初始状态偏差i P -
为:0
i i P x x --
=-,i P -
的值均考虑为典型误差值(其值参考【5】
中表一),则可得:
66644
0000000000000
001.74937210 1.74947210109.010********,,102400431109.4411100i i v v r r x P x P x w w m m θθββ-----
--?????????
? ? ? ?? ? ? ? ?
? ? ? ??=====+= ? ? ? ?? ? ? ? ?
? ? ? ?
? ? ? ? ?
? ? ? ??????????????24101?? ? ? ? ? ? ? ? ?
??, 由0x 可求出实际末状态值f x ,嫦娥三号的标准末状态为f x -
,则终端总误差向量f
P -
为:f f f P x x -
-
=-。即为:
665718.38-38.6215921.73910-0.1541.4172,x ,=0.00058200.0005821018291839651.7.68967.3740810262f
f f
f f f f f f f f
v r x P x x w m πθπβ---???? ????? ? ? ? ?? ? ? ? ? ?
? ?====- ? ? ? ? ? ? ? ? ? ? ? ? ?????? ???
?则????????? 参考资料【5】可得敏感性大小S 是由f P -
里面的每一个元素除以i P -
里面对应的元
素得到的值的集合,为( 3.862015.92000.1540582.0000 1.000065.6890)S =--,它们分别表示的是:v 末状态对0v ?的的敏感性为-3.8620,r 末状态对0r ?的敏感性为15.9200,θ末状态对0θ?的敏感性为-0.1540,w 末状态对0w ?的敏感性为582.0000,m 末状态对0m ?的敏感性为1.0000,β末状态对0β?的敏感性为65.6890。通过对以上六个敏感性大小的对比可看出w 最敏感,即w 有误差时会导致结果误差最大;其次是
β,r ,v ,m ,θ。于是我们应该尽力避免w 初始值有误差的情况。
6.模型的评价与改进
优点:(1)题一使用物理和地理知识对近月点和远月点的位置及嫦娥三号对应的速度大小和方向直接列方程求解,方法简单,直白易懂;
(2)题二随机寻优的优点就是实现比较简单,它不像 GA 、ACO 等其他智能算法那样具有比较复杂的应用流程,也不需要设置太多的参数,因此随机寻优在应用操作上更加简便。计算结果表明了随机寻优能够获得较高的优化精度。
(3)题三建立的初始状态误差模型简单,直白,易懂。 缺点:(1)题二所建的随机寻优模型受代数影响较大,易出现较大波动。
(2)题三所建立的初始状态误差模型只讨论了变量末状态对自己这个变量初状态偏差的误差与敏感性分析,忽视了它对别的变量的误差与敏感性,存在缺陷。
参考文献:
【1】王建伟,李兴;《近日点和远日点速度的两种典型算法》;
https://www.docsj.com/doc/1613928790.html,/Periodical_wuljs201306027.aspx; 2014/9/12 【2】郭景录,付平;《登月软着陆轨道优化算法研究》;
https://www.docsj.com/doc/1613928790.html,/Article/CJFDTotal-JSJZ200912024.htm;2014/9/13 【3】国防科工局,中国科学院;《权威解析嫦娥三号任务全过程》;嫦娥三号专家解读_新闻中心_新浪网https://www.docsj.com/doc/1613928790.html,/c/z/changesanhao4/;2014/9/12; 【4】王劼,李俊峰,崔乃刚,刘暾;《登月飞行器软着陆轨道的遗传算法优化》;https://www.docsj.com/doc/1613928790.html,/p-125650348.html?qq-pf-to=pcqq.discussion;
【5】刘浩敏,冯军华,张泽旭;《月球软着陆制导律及其误差分析》;
https://www.docsj.com/doc/1613928790.html,/p-192729211.html;2014/9/14;
【6】赵静,但琪;《数学建模与数学实验(第三版)》;2008年1月1日。
附录:
题二:
程序:
Youhua.m
M=13;
H=[0 10 20 30 40 50 60 70 80 90 90 90 700]; %寻优上界
L=[0 0 0 0 0 0 0 0 0 0 0 0 500]; %寻优下界
N=500;
a=zeros(N,M);
for i=1:N %初始矩阵a,对a进行初次排序
a(i,:)=L+(H-L).*rand(1,13);
for k=2:M