2023年高教社杯全國大學(xué)生數(shù)學(xué)建模競賽(嫦娥3號軟著陸)

2023高教社杯全國大學(xué)生數(shù)學(xué)建模競賽2023高教社杯全國大學(xué)生數(shù)學(xué)建模競賽承諾書我們認(rèn)真閱讀了《全國大學(xué)生數(shù)學(xué)建模競賽章程》和《全國大學(xué)生數(shù)學(xué)建模競賽〔下載?！舶?、電子郵件、網(wǎng)上詢問等〕與隊外的任何人〔包括指導(dǎo)教師〕爭論、爭論與賽題有關(guān)的問題。我們知道，抄襲別人的成果是違反競賽章程和參賽規(guī)章的，假設(shè)引用別人的成果〔包括網(wǎng)上查到的資料文引用處和參考文獻(xiàn)中明確列出。有違反競賽章程和參賽規(guī)章的行為，我們將受到嚴(yán)峻處理。呈現(xiàn)〔包括進(jìn)展網(wǎng)上公示，在書籍、期刊和其他媒體進(jìn)展正式或非正式發(fā)表等。我們參賽選擇的題號是〔從A/B/C/D中選擇一項填寫：A 我們的報名參賽隊號為〔8位數(shù)字組成的編號19005007 所屬學(xué)?！舱?zhí)顚懲暾娜海?.2.3.指導(dǎo)教師或指導(dǎo)教師組負(fù)責(zé)人(打印并簽名)：〔論文紙質(zhì)版與電子版中的以上信息必需全都，只是電子版中無需簽名。以上內(nèi)容〕日期：2023 年9月14日2023高教社杯全國大學(xué)生數(shù)學(xué)建模競賽2023高教社杯全國大學(xué)生數(shù)學(xué)建模競賽賽區(qū)評閱編號〔由賽區(qū)組委會評閱前進(jìn)展編號〕：評閱人評閱人評分備注全國統(tǒng)一編號〔由賽區(qū)組委會送交全國前編號〕：全國評閱編號〔由全國組委會評閱前進(jìn)展編號〕：一次參賽獲益終身一次參賽獲益終身10嫦娥三號軟著陸軌道設(shè)計與把握策略摘要軌道設(shè)計與把握策略問題：針對問題一，查資料可知嫦娥三號軟著陸預(yù)備軌道是極月軌道，無視月球自轉(zhuǎn)公19.51W。然后依據(jù)問題二得出的結(jié)果求出水平位移X，再列出微分方程方程并求其數(shù)值解，最終得到近月點投影的緯度，進(jìn)而依據(jù)近月點投影的經(jīng)緯度算出遠(yuǎn)月點投影的經(jīng)緯度，至此便得到近月點與遠(yuǎn)月點的位置；對于二者的速度，應(yīng)用開普勒定律和機械能守恒59.07N19.51W15kmVA

1.692103m/s，59.07S,160.49E100km，V 1.614103m/s，方向為南極-北極-虹灣方向。B針對問題二，在安全著陸的前提下，確定著陸軌道和最優(yōu)把握策略，即轉(zhuǎn)化為燃tf633s，消耗燃料的最645.9Kg。我們首先列出軟著陸軌道動力學(xué)方程并做歸一化處理，經(jīng)過軟著陸軌道的離散化，應(yīng)用函數(shù)靠近方法擬合推力把握角 (t)，擬合結(jié)果為(t)7.47107x36.14104x20.026x1.72;。從而將軌道優(yōu)化問題轉(zhuǎn)化為參數(shù)優(yōu)MATLAB針對問題三，我們著重對主減速階段進(jìn)展相應(yīng)的誤差分析和敏感性分析。建立月SPf

里面的每一個元素除以P里面對應(yīng)的元素得到的值的集合，六個值分別對應(yīng)v、r、θ、w、m、β，P 為以Pi fP主減速階段末狀態(tài)總誤差向量，Pi

為初始狀態(tài)偏差向量。得到結(jié)果為S(3.8620 15.9200 0.1540 582.0000 1.0000 65.6890) ,由這六個值大小分析可β，r，v，m，θ。于是我w關(guān)鍵詞：軟著陸動力學(xué)方程微分方程隨機尋優(yōu)尋優(yōu)變量問題重述嫦娥三號軟著陸軌道設(shè)計與把握策略130126抵達(dá)月球軌道。嫦〔即單位質(zhì)量的推動劑產(chǎn)生的推力2940m/s，可以滿足調(diào)整速度的把握要求。在四周安裝有姿勢調(diào)整發(fā)動機，在給定主減速發(fā)動機的推力方向后，能夠自動通過多個發(fā)動機的脈沖組合實現(xiàn)各種姿勢的調(diào)整把握。嫦娥19.51W，44.12N，海拔為-2641m。嫦娥三號在高速飛行的狀況下，要保證準(zhǔn)確地在月球預(yù)定區(qū)域內(nèi)實現(xiàn)軟著陸，關(guān)鍵問題是著陸軌道與把握策略的設(shè)計。其著陸軌道設(shè)計的根本要求：著陸預(yù)備軌道為15km，遠(yuǎn)月點100km6的燃料消耗。依據(jù)上述的根本要求，請你們建立數(shù)學(xué)模型解決下面的問題：方向。確定嫦娥三號的著陸軌道和在6個階段的最優(yōu)把握策略。對于你們設(shè)計的著陸軌道和把握策略做相應(yīng)的誤差分析和敏感性分析模型假設(shè)57m/s度；2、假設(shè)近月點投影到月面的地區(qū)與虹灣之間的月面弧線可以看作一條直線；3運動無影響；4F=3000N5、假設(shè)收集到的數(shù)據(jù)準(zhǔn)確無誤。符號說明r月心距;:極角;w:角速度;m:質(zhì)量;:制動發(fā)動機的推力方向角問題分析：74311.559有幾百秒，所以我們可以無視由于自轉(zhuǎn)和公轉(zhuǎn)造成軟著陸軌道平面與著陸預(yù)備軌道平面的偏差，即可認(rèn)為軟著陸軌道平面與著陸預(yù)備軌道平面是在同一個平面上的。著陸預(yù)備軌道為極月軌道，則近月點經(jīng)度就可確定，假設(shè)得到近月點投影與虹灣水平X使用開普勒定律和機械能守恒定律則很簡潔算出。：〔徑向速〕隨時間變化的圖像。優(yōu)化的指標(biāo)是在保證安全和可控條件下使得燃料消耗最少〔tfpnA(i)推斷平移方向和距離。：初始條件誤差。然后建立初始狀態(tài)誤差模型，通過將v,r,,w,m, 的值0 0 0 0 0 0〔其值均考慮為典型誤差值〕分別與v,r,

,w,m和

〔即初始狀態(tài)著陸器沿r方向0 0 0 0 0 0的速度、月心距、極角、角速度、質(zhì)量和制動發(fā)電機推力方向角〕相加，然后得到非v,r,

,,

f f f f f

S，以S中每個值的大小來描述其對應(yīng)的變量0 0 0 0 0 0初始狀態(tài)對該變量的末狀態(tài)的影響即敏感性分析與誤差分析。模型的建立與求解：求近月點與遠(yuǎn)月點的位置及嫦娥三號相應(yīng)的速度大小和方向3即通過月球南北極的軌道，其方向為北極-虹灣〔著陸點〕1圖1

則我們可知近月點的經(jīng)度和虹灣一樣為19.51W，虹X灣Nθ44.12NR19.51WS虹X灣Nθ44.12NR19.51WS遠(yuǎn)月點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，但與月球半徑相比數(shù)值很小，所以我們可以無視Q=θ+44.12=59.07N；所以其緯度為59.07S，其經(jīng)度與近日點經(jīng)度相差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,則嫦娥三號與月球的連線在這段時間內(nèi)掃過的面積分別為：依據(jù)開普勒其次定律：在相等時間內(nèi)，太陽和運動中的行星的連線〔向量半徑〕L所掃過的面積都是相等的?？傻肧A

SB

VB

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

GMm( L ),(2)E 1mV

A ,GMm( ),(3)B 2 B LB〔以無窮遠(yuǎn)處為引力勢能零點,G6.671011Nm2/kg2〕三號在著陸預(yù)備軌道上只受萬有引力作用，所以遵循機械能守恒定律，即有：E E〔4〕A B然后聯(lián)立1〔2〔3〔4〕得到：2LGM2V BA

,V ,2L2LGML(LL)AB ABA A B最終帶入數(shù)據(jù)算得：V 1.692103m/s,VA B

1.614103m/s59.07N19.51W15km，VA1.692103m/s，方向為北極-虹灣-南極方向;59.07S160.49E100kmB1.614103m/s方向為南極-北極-虹灣方向。：求嫦娥三號的著陸軌道及六個階段最優(yōu)把握策略57m/s。F3YβrθwO、F3YβrθwO月點，X向。軟著陸的動力學(xué)方程描述如式子〔1：dv Fsin u dt

m r2

rw2drdtdt

vw 〔1【2】XdwFcos 2vwdt mr rdmFdt veF是制動發(fā)動機推力，其值恒定， v 為比沖=2940m/s，u為月球引力常數(shù)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，則近月點的曲率半徑也應(yīng)為R，由向2 B 2 Av2 v心力公式：m ARA

mw2RA

可得w A=9.431×10-4rad/s。RA月球軟著陸軌道聽從兩點邊值約束，近月點為起點，對于此處的我們認(rèn)為近月點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〕待優(yōu)化變量為制動發(fā)動機的推力方向角β〔t〕和終端時刻t ，優(yōu)化的性能指標(biāo)f為滿足初始條件和終端約束條件的前提下，主減速階段燃料消耗最少，即為：minJ

tfm(t)dt 〔4〕 【2】0〔4〕描述的性能指標(biāo)也相當(dāng)于使得終端時刻t 最短。f、數(shù)值歸一化處理：度，我們需要對狀態(tài)變量進(jìn)展歸一化處理：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經(jīng)過歸一化處理后，初始條件和終端約束條件改寫為：- - - 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、參數(shù)化方法：β〔t〕的搜尋空間是一個泛函空間，必需將把握變量參數(shù)化，本文承受函數(shù)靠近的方法。函數(shù)靠近法假設(shè)推力方向角可用一個多項式〔9〕表示：〔t〕tt2t3 〔9〕 【4】0 1 2 34Ni ,,.i

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

多項式擬合，從而求出多項式系數(shù)，i〔如遺傳算法，蟻群算法〕雖然能夠準(zhǔn)確的計算出最優(yōu)值，但由隨機尋優(yōu)步驟: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]設(shè)置優(yōu)化次數(shù)M，此題中M=500X=L+〔H-L〕*rand,其中rand是區(qū)間(0,1)之間的均勻隨機數(shù)。所對應(yīng)的參數(shù)，模型中得到的最優(yōu)參數(shù)為：[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]、優(yōu)化策略的求解：〔t〔t與時間的關(guān)系，運用MATLAB畫圖可得其圖像。同理，我們也可以得到其他參數(shù)與時間的關(guān)系如下：(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〔 其 他 物 理 量 的 圖 像 見 附 頁 〕、優(yōu)化軌道的求解：依據(jù)上文的求解可知：、主減速階段后的優(yōu)化：〔徑向速度，角速度，極角，推力方向角等〕隨時間變化的函數(shù)關(guān)系，同時也可以得到月323和附件5給出了嫦娥3號在2400米和100MATLAB軟件中imread函數(shù)可讀入數(shù)據(jù)，但數(shù)據(jù)格式為uinnt8無法直接運算，故通過im2double函數(shù)將數(shù)據(jù)轉(zhuǎn)化為浮點型，進(jìn)而進(jìn)展運算分析。數(shù)據(jù)轉(zhuǎn)化后我們分別得到高度矩陣a和b,下面我們以矩陣a為例進(jìn)展分析。我們將2300*2300的矩陣a分為50*46個網(wǎng)格點，每個網(wǎng)格點含有46*50個數(shù)據(jù)。由于網(wǎng)格點足夠密集，對任一網(wǎng)格點我們承受mean函數(shù)求其平均值，并用這一平均值代替該網(wǎng)格點的特2300*2300的矩陣a便可用50*46的矩陣x最終我們利用mesh函數(shù)繪出三維高度圖如下：降落軌道進(jìn)展調(diào)整以到達(dá)避開凹凸和低洼，應(yīng)馬上啟動姿勢調(diào)整發(fā)動機對嫦娥3號進(jìn)展平移，最終安全降落在預(yù)定地點。下面我們爭論如何進(jìn)展平移計算平移平移點附件各個方向的坡度p1，p2，p3??pn。依據(jù)各坡度的大小進(jìn)展加權(quán)獲得各個方向的平移指數(shù)a(1),a(2),a(3)??a(n),坡度p越大，平移指數(shù)a越小。3.由某一方向的綜合平移指數(shù)Ai=a(i)-[a(1)+a(2)+??+a(i-1)+a(i+1)+??+a(n)]來確定平移的方向，A(i)的最大值對應(yīng)的方向即為平移的方向。4.對于平移距離我們承受在滿足條件的狀況下，平移距離最小的原則，由于這樣可以最大限度的節(jié)約燃料。問題三：對設(shè)計的著陸軌道和把握策略做相應(yīng)的誤差分析和敏感性分析、初始狀態(tài)誤差模型：

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

嫦娥三號的標(biāo)準(zhǔn)末狀態(tài)為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

里面對應(yīng)的元素得到的值的集合，為S(3.8620 15.9200 0.1540 582.0000 1.0000 它們分別表示的是：v末狀態(tài)對v0

的的敏感性為-3.8620，r末狀態(tài)對r0

的敏感性為15.9200，θ末狀態(tài)對

的敏感性為-0.1540，w末狀態(tài)對w0

582.0000，m末狀態(tài)對m0

1.0000，β0

六個敏感性大小的比照可看出w最敏感，即wβ，r，v，m，θ。于是我們應(yīng)當(dāng)盡力避開w初始值有誤差的狀況。6.模型的評價與改進(jìn)優(yōu)點〔1〕大小和方向直接列方程求解，方法簡潔，直白易懂；題二隨機尋優(yōu)的優(yōu)點就是實現(xiàn)比較簡潔，它不像GA、ACO等其他智能算作上更加簡便。計算結(jié)果說明白隨機尋優(yōu)能夠獲得較高的優(yōu)化精度。題三建立的初始狀態(tài)誤差模型簡潔，直白，易懂。缺點〔1〕題二所建的隨機尋優(yōu)模型受代數(shù)影響較大，易消滅較大波動。題三所建立的初始狀態(tài)誤差模型只爭論了變量末狀態(tài)對自己這個變量初狀態(tài)偏差的誤差與敏感性分析，無視了它對別的變量的誤差與敏感性，存在缺陷。參考文獻(xiàn)：【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]; %尋優(yōu)上界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); 微分方程數(shù)值解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(”質(zhì)量”)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=