2023年高教社杯全国大学生数学建模竞赛(嫦娥3号软着陆)

2023高教社杯全国大学生数学建模竞赛2023高教社杯全国大学生数学建模竞赛承诺书我们认真阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛〔下载。〔包括、电子邮件、网上询问等〕与队外的任何人〔包括指导教师〕争论、争论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛章程和参赛规章的，假设引用别人的成果〔包括网上查到的资料文引用处和参考文献中明确列出。有违反竞赛章程和参赛规章的行为，我们将受到严峻处理。呈现〔包括进展网上公示，在书籍、期刊和其他媒体进展正式或非正式发表等。我们参赛选择的题号是〔从A/B/C/D中选择一项填写：A 我们的报名参赛队号为〔8位数字组成的编号19005007 所属学校〔请填写完整的全名：：1.2.3.指导教师或指导教师组负责人(打印并签名)：〔论文纸质版与电子版中的以上信息必需全都，只是电子版中无需签名。以上内容〕日期：2023 年9月14日2023高教社杯全国大学生数学建模竞赛2023高教社杯全国大学生数学建模竞赛赛区评阅编号〔由赛区组委会评阅前进展编号〕：评阅人评阅人评分备注全国统一编号〔由赛区组委会送交全国前编号〕：全国评阅编号〔由全国组委会评阅前进展编号〕：一次参赛获益终身一次参赛获益终身10嫦娥三号软着陆轨道设计与把握策略摘要轨道设计与把握策略问题：针对问题一，查资料可知嫦娥三号软着陆预备轨道是极月轨道，无视月球自转公19.51W。然后依据问题二得出的结果求出水平位移X，再列出微分方程方程并求其数值解，最终得到近月点投影的纬度，进而依据近月点投影的经纬度算出远月点投影的经纬度，至此便得到近月点与远月点的位置；对于二者的速度，应用开普勒定律和机械能守恒59.07N19.51W15kmVA

1.692103m/s，59.07S,160.49E100km，V 1.614103m/s，方向为南极-北极-虹湾方向。B针对问题二，在安全着陆的前提下，确定着陆轨道和最优把握策略，即转化为燃tf633s，消耗燃料的最645.9Kg。我们首先列出软着陆轨道动力学方程并做归一化处理，经过软着陆轨道的离散化，应用函数靠近方法拟合推力把握角 (t)，拟合结果为(t)7.47107x36.14104x20.026x1.72;。从而将轨道优化问题转化为参数优MATLAB针对问题三，我们着重对主减速阶段进展相应的误差分析和敏感性分析。建立月SPf

里面的每一个元素除以P里面对应的元素得到的值的集合，六个值分别对应v、r、θ、w、m、β，P 为以Pi fP主减速阶段末状态总误差向量，Pi

为初始状态偏差向量。得到结果为S(3.8620 15.9200 0.1540 582.0000 1.0000 65.6890) ,由这六个值大小分析可β，r，v，m，θ。于是我w关键词：软着陆动力学方程微分方程随机寻优寻优变量问题重述嫦娥三号软着陆轨道设计与把握策略130126抵达月球轨道。嫦〔即单位质量的推动剂产生的推力2940m/s，可以满足调整速度的把握要求。在四周安装有姿势调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿势的调整把握。嫦娥19.51W，44.12N，海拔为-2641m。嫦娥三号在高速飞行的状况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与把握策略的设计。其着陆轨道设计的根本要求：着陆预备轨道为15km，远月点100km6的燃料消耗。依据上述的根本要求，请你们建立数学模型解决下面的问题：方向。确定嫦娥三号的着陆轨道和在6个阶段的最优把握策略。对于你们设计的着陆轨道和把握策略做相应的误差分析和敏感性分析模型假设57m/s度；2、假设近月点投影到月面的地区与虹湾之间的月面弧线可以看作一条直线；3运动无影响；4F=3000N5、假设收集到的数据准确无误。符号说明r月心距;:极角;w:角速度;m:质量;:制动发动机的推力方向角问题分析：74311.559有几百秒，所以我们可以无视由于自转和公转造成软着陆轨道平面与着陆预备轨道平面的偏差，即可认为软着陆轨道平面与着陆预备轨道平面是在同一个平面上的。着陆预备轨道为极月轨道，则近月点经度就可确定，假设得到近月点投影与虹湾水平X使用开普勒定律和机械能守恒定律则很简洁算出。：〔径向速〕随时间变化的图像。优化的指标是在保证安全和可控条件下使得燃料消耗最少〔tfpnA(i)推断平移方向和距离。：初始条件误差。然后建立初始状态误差模型，通过将v,r,,w,m, 的值0 0 0 0 0 0〔其值均考虑为典型误差值〕分别与v,r,

,w,m和

〔即初始状态着陆器沿r方向0 0 0 0 0 0的速度、月心距、极角、角速度、质量和制动发电机推力方向角〕相加，然后得到非v,r,

,,

f f f f f

S，以S中每个值的大小来描述其对应的变量0 0 0 0 0 0初始状态对该变量的末状态的影响即敏感性分析与误差分析。模型的建立与求解：求近月点与远月点的位置及嫦娥三号相应的速度大小和方向3即通过月球南北极的轨道，其方向为北极-虹湾〔着陆点〕1图1

则我们可知近月点的经度和虹湾一样为19.51W，虹X湾Nθ44.12NR19.51WS虹X湾Nθ44.12NR19.51WS远月点x

,a,vx

分别是切向θ不是指极角：tan

(t)ararxdva r

t t dvr dt

Xf(f(cot r)dt)dtv tx 0

x

0 0 dtXtf0

vdtxX4.518105m2 2 2 2 b2c2a然后我们由余弦定理公式：cosA ，其中A,b c R,a X ，2 2 2 2 2bc〔虹湾〔着陆点〕海拔虽然为-2641m，但与月球半径相比数值很小，所以我们可以无视Q=θ+44.12=59.07N；所以其纬度为59.07S，其经度与近日点经度相差180度，则可以算得为180-19.51=160.49E。在此类，于是我们可以使用能量法[1]来计算VA

VBS A

VtLA

,SB

1V2 B

XYA2YA2bfcaBB

如图2所示，由开普勒第确定律：每一个行M=7.3477×1022kg，月球平均半径r=1.737013106m,A到f的距离LA〔1737.013+15×103=1.752023×10〔无视虹湾与月平均半径的差2641mB到f的距离L=1737.013+100=1.837013×106A,BB别取t→0,则嫦娥三号与月球的连线在这段时间内扫过的面积分别为：依据开普勒其次定律：在相等时间内，太阳和运动中的行星的连线〔向量半径〕L所扫过的面积都是相等的。可得SA

SB

VB

VA〔1。ALBAB械能为：E 1mVA 2 A

GMm( L ),(2)E 1mV

A ,GMm( ),(3)B 2 B LB〔以无穷远处为引力势能零点,G6.671011Nm2/kg2〕三号在着陆预备轨道上只受万有引力作用，所以遵循机械能守恒定律，即有：E E〔4〕A B然后联立1〔2〔3〔4〕得到：2LGM2V BA

,V ,2L2LGML(LL)AB ABA A B最终带入数据算得：V 1.692103m/s,VA B

1.614103m/s59.07N19.51W15km，VA1.692103m/s，方向为北极-虹湾-南极方向;59.07S160.49E100kmB1.614103m/s方向为南极-北极-虹湾方向。：求嫦娥三号的着陆轨道及六个阶段最优把握策略57m/s。F3YβrθwO、F3YβrθwO月点，X向。软着陆的动力学方程描述如式子〔1：dv Fsin u dt

m r2

rw2drdtdt

vw 〔1【2】XdwFcos 2vwdt mr rdmFdt veF是制动发动机推力，其值恒定， v 为比冲=2940m/s，u为月球引力常数e=4.9028×112m3/s2βF由资料可知近日点曲率半径R=b2〔图一：a为半长轴=〔L 2A a A BL L L L 2短轴=a2-c2=( A B)2L A B【1，则近月点的曲率半径也应为R，由向2 B 2 Av2 v心力公式：m ARA

mw2RA

可得w A=9.431×10-4rad/s。RA月球软着陆轨道听从两点边值约束，近月点为起点，对于此处的我们认为近月点15-2641mr=L0 A=1.752023×106-2641=1.749372×106m ，r 月平均半径+海拔+距地高度30001.737408106m,所以满足的初始条件为：fv 0m/s,r0 0

1.749372106m,0

0rad,w0

9.431104rad/s 〔2〕终端约束条件为软着陆，需满足： 【2】v 57m/s,rf

=1.737408106m，f

rad,w2

0rad/s 〔3〕待优化变量为制动发动机的推力方向角β〔t〕和终端时刻t ，优化的性能指标f为满足初始条件和终端约束条件的前提下，主减速阶段燃料消耗最少，即为：minJ

tfm(t)dt 〔4〕 【2】0〔4〕描述的性能指标也相当于使得终端时刻t 最短。f、数值归一化处理：度，我们需要对状态变量进展归一化处理：u m v 2 r令r

m,v

ref ref ,t

ref， 【2】ref

0

0

rref

ref

rref

ref

vref v r

r Fv ,r ,v

ref,F ,v则： refm

rref

e er3

u Freft

【2】m ,ww

ref,t ,mref

u t 0ref经过归一化处理后，初始条件和终端约束条件改写为：- - - r3r0

1,v0

0,w w0

0, 〔5〕u- r - v -rf

f,v r f0

f ,wuu/rref

0 〔6〕 sin 1 v””

rw2 m r2vw 〔7【2】w cos 2w”

F

mr rm”v ve Jtfm”(t)dt 〔8〕0、参数化方法：β〔t〕的搜寻空间是一个泛函空间，必需将把握变量参数化，本文承受函数靠近的方法。函数靠近法假设推力方向角可用一个多项式〔9〕表示：〔t〕tt2t3 〔9〕 【4】0 1 2 34Ni ,,.i

后利用这些推力方向做、算法的求解：

多项式拟合，从而求出多项式系数，i〔如遗传算法，蚁群算法〕虽然能够准确的计算出最优值，但由随机寻优步骤: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]设置优化次数M，此题中M=500X=L+〔H-L〕*rand,其中rand是区间(0,1)之间的均匀随机数。所对应的参数，模型中得到的最优参数为：[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]、优化策略的求解：〔t〔t与时间的关系，运用MATLAB画图可得其图像。同理，我们也可以得到其他参数与时间的关系如下：(t)7.47107x36.14104x20.026x1.72;v 3.66107x31.85104x20.219x6.01;rr1.86104x30.165x214.8x1.75106;6.731011x34.0107x29.48104x1.38104;w2.211013x35.121012x27.44107x9.44104;m1.931019x31.11016x21.02x2400.0〔 其 他 物 理 量 的 图 像 见 附 页 〕、优化轨道的求解：依据上文的求解可知：、主减速阶段后的优化：〔径向速度，角速度，极角，推力方向角等〕随时间变化的函数关系，同时也可以得到月323和附件5给出了嫦娥3号在2400米和100MATLAB软件中imread函数可读入数据，但数据格式为uinnt8无法直接运算，故通过im2double函数将数据转化为浮点型，进而进展运算分析。数据转化后我们分别得到高度矩阵a和b,下面我们以矩阵a为例进展分析。我们将2300*2300的矩阵a分为50*46个网格点，每个网格点含有46*50个数据。由于网格点足够密集，对任一网格点我们承受mean函数求其平均值，并用这一平均值代替该网格点的特2300*2300的矩阵a便可用50*46的矩阵x最终我们利用mesh函数绘出三维高度图如下：降落轨道进展调整以到达避开凹凸和低洼，应马上启动姿势调整发动机对嫦娥3号进展平移，最终安全降落在预定地点。下面我们争论如何进展平移计算平移平移点附件各个方向的坡度p1，p2，p3„„pn。依据各坡度的大小进展加权获得各个方向的平移指数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.对于平移距离我们承受在满足条件的状况下，平移距离最小的原则，由于这样可以最大限度的节约燃料。问题三：对设计的着陆轨道和把握策略做相应的误差分析和敏感性分析、初始状态误差模型：

x

x0，

x 中表一，则可得：

i为：i

i 0，

i的值均考虑为典型误差值〔其值参考【5】v 0

v

10

10 0

r 1.749372106 r

100

1.7494721060 0

1

1 x

0 ,P

0

,x

Px ，0 w

9.431104

i

106

i 0 9.441104 0 0 m 2400

m 100 0 0 0

2410 0 00

1

1 xx0

嫦娥三号的标准末状态为x 则终端总误差向量Pf f 为：P x x 。即为：f f f57 v f

18.38 -38.62 1.737408106 rf

1.739106 1592

1.417 -0.154 x

f

,x

,则P x

x = f w

f

f f

0.000582 f

10 mf

1829

183967.26

65.689 f 2 参考资料【5SPf

里面的每一个元素除以Pi

里面对应的元素得到的值的集合，为S(3.8620 15.9200 0.1540 582.0000 1.0000 它们分别表示的是：v末状态对v0

的的敏感性为-3.8620，r末状态对r0

的敏感性为15.9200，θ末状态对

的敏感性为-0.1540，w末状态对w0

582.0000，m末状态对m0

1.0000，β0

六个敏感性大小的比照可看出w最敏感，即wβ，r，v，m，θ。于是我们应当尽力避开w初始值有误差的状况。6.模型的评价与改进优点〔1〕大小和方向直接列方程求解，方法简洁，直白易懂；题二随机寻优的优点就是实现比较简洁，它不像GA、ACO等其他智能算作上更加简便。计算结果说明白随机寻优能够获得较高的优化精度。题三建立的初始状态误差模型简洁，直白，易懂。缺点〔1〕题二所建的随机寻优模型受代数影响较大，易消灭较大波动。题三所建立的初始状态误差模型只争论了变量末状态对自己这个变量初状态偏差的误差与敏感性分析，无视了它对别的变量的误差与敏感性，存在缺陷。参考文献：【1“://d.g.wanfangdata.cn/Periodical_wuljs202306027.aspx%3B“://d.g.wanfangdata.cn/Periodical_wuljs202306027.aspx;2023/9/12【2“://cnki.cn/Article/CJFDTotal-JSJZ202312024.htm“://cnki.cn/Article/CJFDTotal-JSJZ202312024.htm；2023/9/13【3读_闻中心_“://news.sina.cn/c/z/changesanhao4/“://news.sina.cn/c/z/changesanhao4/；2023/9/12;【4“://docin/p-125650348.html?qq-pf-to=pcqq.discussion“://docin/p-125650348.html?qq-pf-to=pcqq.discussion；【5“://docin/p-192729211.html“://docin/p-192729211.html；2023/9/14;【6〔第三版202311附录：题二：程序：Youhua.mYouhua.mM=13;H=[01020304050607080909090700]; %寻优上界L=[000000000000500]; N=500;a=zeros(N,M);fori=1:N a(i,:)=L+(H-L).*rand(1,13);fork=2:Mifa(i,k)<a(i,k-1)a(i,k)=a(i,k-1);end

end

endfori=1:N fitness(i)=fun(a(i,:));end[p,q]=sort(fitness); kkkk=a(q(1),:)fitness(q(1));x=a(q(1),:);aa=x(1:12);tf=x(13);t=0:tf/11:tf;b=polyfit(t,aa,3);beta=poly2sym(b) figure(1)plot(t,polyval(b,t))xlabel(”时间”)ylabel(”推力方向角”)[T,Y]=ode45(”weifen”,[0:0.5:tf],[0,1749372,0,0.0009431,2400],[],b); 微分方程数值解q1=Y(:,1);q2=Y(:,2);q3=Y(:,3);q4=Y(:,4);q5=Y(:,5);figure(2)plot(T,q1,”b”)xlabel(”时间”)figure(3)plot(T,q2,”r”)ylabel(”月心距”)figure(4)plot(T,q3,”k”)ylabel(”极角”)figure(5)plot(T,q4)ylabel(”角速度”)figure(6)plot(T,q5,”g”)ylabel(”质量”)Weifen.mWeifen.mfunctiondy=weifen(t,y,flag,b)dy=zeros(5,1);M=3000;dy(1)=M*sin((b(1)*t^3+b(2)*t^2+b(3)*t^1+b(4))*pi/180)/y(5)-(4.9*10^12)/(y(2))^2+y(2)*y(4)^2;dy(2)=y(1);dy(3)=y(4);dy(4)=-M*cos((b(1)*t^3+b(2)*t^2+b(3)*t^1+b(4))*pi/180)/(y(5)*y(2))-2*y(1)*y(4)/y(2);dy(5)=-M/2940;fun.mfun.mfunctiony=fun(x)tf=x(13);a=x(1:12);t=