2014A數(shù)學(xué)建模優(yōu)秀論文

1、2014 高教社杯全國大學(xué)生數(shù)學(xué)建模競賽承諾書我們仔細閱讀了全國大學(xué)生數(shù)學(xué)建模競賽章程和全國大學(xué)生數(shù)學(xué)建模競賽參賽規(guī)則（以下簡稱為“競賽章程和參賽規(guī)則”，可從全國大學(xué)生數(shù)學(xué)建模競賽網(wǎng)站下載）。我們完全明白，在競賽開始后參賽隊員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問題。我們知道，抄襲別人的成果是違反競賽章程和參賽規(guī)則的，如果引用別人的成果或其他公開的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻的表述方式在正文引用處和參考文獻中明確列出。我們鄭重承諾，嚴格遵守競賽章程和參賽規(guī)則，以保證競賽的公正、公平性。如有違反競賽章程和參賽規(guī)

2、則的行為，我們將受到嚴肅處理。我們授權(quán)全國大學(xué)生數(shù)學(xué)建模競賽組委會，可將我們的論文以任何形式進行公開展示（包括進行網(wǎng)上公示，在書籍、期刊和其他媒體進行正式或非正式發(fā)表等）。我們參賽選擇的題號是（從 A/B/C/D 中選擇一項填寫）： A 我們的報名參賽隊號為（8 位數(shù)字組成的編號）： 10009072所屬學(xué)校（請?zhí)顚懲暾娜?東南大學(xué) 參賽隊員 (打印并簽名) ：1. 吉張鶴軒2. 楊升3. 陳同廣 指導(dǎo)教師或指導(dǎo)教師組負責人(打印并簽名)： （論文紙質(zhì)版與電子版中的以上信息必須一致，只是電子版中無需簽名。以上內(nèi)容請仔細核對，提交后將不再允許做任何修改。如填寫錯誤，論文可能被取消評獎資格

3、。） 日期： 2014 年 09 月 15 日 賽區(qū)評閱編號（由賽區(qū)組委會評閱前進行編號）： 2014 高教社杯全國大學(xué)生數(shù)學(xué)建模競賽編 號 專 用 頁賽區(qū)評閱編號（由賽區(qū)組委會評閱前進行編號）： 賽區(qū)評閱記錄（可供賽區(qū)評閱時使用）： 評 閱人 評 分 備注 全國統(tǒng)一編號（由賽區(qū)組委會送交全國前編號）： 全國評閱編號（由全國組委會評閱前進行編號）： 嫦娥三號軟著陸軌道設(shè)計與控制策略摘要本題要求我們以嫦娥三號登月為背景，分析登月軌道參數(shù)，重點探討了登月過程最具難度的著陸軌道設(shè)計優(yōu)化，并對所使用的優(yōu)化方案進一步作出了誤差分析與靈敏度分析。對于第一問，由于正面求解條件有限，難以從已有的條件中得到近月

4、點和遠月點的位置以及準備軌道參數(shù)，因此巧妙的使用了逆推思路，通過已知條件求解主減速階段運動過程， 通過水平位移量反推近月點位置。首先，我們由已知的端點約束條件，結(jié)合所學(xué)物理知識， 并查閱相關(guān)資料，通過線性正切制導(dǎo)率等工程經(jīng)驗判斷減速過程中推力控制方案應(yīng)為推力大小始終為最大值，推力與速度反方向夾角也為恒量，由此建立微分方程模型。但在求解的過程中我們發(fā)現(xiàn)，正面求解難度十分大，于是對微分方程離散化，轉(zhuǎn)化為差分方程組，繼而通過計算機模擬行為，擬合出最接近解的軌跡，求得水平位移量為 385.21m，由此得到近月點位置為 19.51°W ，31.50°N ，距月面 15km，速度為 1

5、692.46m/s，平行月面指向預(yù)期落點方位； 遠月點位置為 160.49°E，31.50°S ，距月面 100km，速度為 1612.15m/s，方向與近日點反相平行。對于第二問，依據(jù)題意將其分為五個階段。對于前兩個階段，建立天體力學(xué)分析中常用的二體模型及力學(xué)方程組，并對約束條件進行歸一化處理。為了后續(xù)優(yōu)化，列出全部歸一化約束條件并建立目標函數(shù)，從而獲得非線性規(guī)劃模型。之后采用了序列化遺傳算法對其進行篩選，從而逼近該階段全局最優(yōu)解，最終此階段燃耗為 1055.39kg。第二階段同樣屬于復(fù)雜多變量優(yōu)化問題，同階段一建立優(yōu)化模型并通過遺傳算法求解該段全局最優(yōu)解，最終燃耗為26

6、.71kg。第三第四階段主要在于圖像處理與統(tǒng)計，第三階段，通過對高程圖進行 K 均值聚類分析，將圖中像素點分為安全點與危險點，在對地圖柵格化，并對方格內(nèi)點類型統(tǒng)計取整， 對方格二元化為安全格與危險格，在通過擴大尋找最大安全半徑，綜合考慮水平偏移量，建立合理的落點評價體系，最終找出最優(yōu)點坐標（1275,1000），燃耗為 86.97kg。第四階段同樣做聚類分析，并根據(jù)嫦娥三號實際體積選取合適的柵格大小，并對柵格通過最小二乘法對空間進行線性統(tǒng)計回歸，求出平均坡面與平均坡度，結(jié)合最大安全半徑建立最優(yōu)落點評價體系， 最終獲得最優(yōu)落點坐標為（88,56），燃耗為 20.68kg。第五階段最優(yōu)燃耗為 8.

7、09kg。最后還討論了簡單運動的局部最優(yōu)模型，簡化了后幾個階段的運動學(xué)分析與計算。最終綜合各段最優(yōu)解，獲得最優(yōu)著陸軌道與控制策略。對于第三問，首先總結(jié)了優(yōu)化模型中引入的一些誤差因素，并針對主要因素做了數(shù)值上的相對誤差分析，證明了誤差對于優(yōu)化方案并未產(chǎn)生很大影響。其次從初始變量和約束條件入手，分析了這些變量的波動對于結(jié)果產(chǎn)生的影響，最終發(fā)現(xiàn)角度控制向量的靈敏性較高， 而其他因素的靈敏性普遍處在較低水平，側(cè)面說明了優(yōu)化方案的對于初值的不敏感性與方案對于全局最優(yōu)解的逼近程度較高。關(guān)鍵詞：非線性規(guī)劃模型序列化遺傳算法K 均值聚類空間線性回歸二體模型一問題的提出1.1 背景介紹根據(jù)計劃，嫦娥三號將在北京

8、時間 12 月 14 號在月球表面實施軟著陸。嫦娥三號如何實現(xiàn)軟著陸以及能否成功成為外界關(guān)注焦點。目前，全球僅有美國、前蘇聯(lián)成功實施了 13 次無人月球表面軟著陸。北京時間 12 月 10 日晚，嫦娥三號已經(jīng)成功降軌進入預(yù)定的月面著陸準備軌道，這是嫦娥三號“落月”前最后一次軌道調(diào)整。在實施軟著陸之前，嫦娥三號還將在橢圓軌道上繼續(xù)飛行，做最后準備。嫦娥三號著陸地點選在較為平坦的虹灣區(qū)。但由于月球地形的不確定性，最終“落月”地點的選擇仍存在一定難度。在整個“落月”過程中，“動力下降”被業(yè)內(nèi)形容為最驚心動魄的環(huán)節(jié)。在這個階段，嫦娥三號要完全依靠自主導(dǎo)航控制，完成降低高度、確定著陸點、實施軟著陸等一系

9、列關(guān)鍵動作，人工干預(yù)的可能性幾乎為零。在距月面 100 米處時，嫦娥三號要進行短暫的懸停，掃描月面地形，避開障礙物，尋找著陸點。1.2 問題重述嫦娥三號在著陸準備軌道上的運行質(zhì)量為 2.4t，其安裝在下部的主減速發(fā)動機能夠產(chǎn)生的推力可調(diào)節(jié)，變化范圍為 1500N 到 7500N，其比沖為 2940m/s，可以滿足調(diào)整速度的控制要求。在四周安裝有姿態(tài)調(diào)整發(fā)動機，能夠自動通過多個發(fā)動機的脈沖組合實現(xiàn)各種姿態(tài)的調(diào)整控制。嫦娥三號的預(yù)定著陸點為 19.51W，44.12N，海拔為-2641m。嫦娥三號在高速飛行的情況下，要保證準確地在月球預(yù)定區(qū)域內(nèi)實現(xiàn)軟著陸，關(guān)鍵問題是著陸軌道與控制策略的設(shè)計。其著陸

10、軌道設(shè)計的基本要求：著陸準備軌道為近月點 15km， 遠月點 100km 的橢圓形軌道；著陸軌道為從近月點至著陸點，其軟著陸過程共分為 6 個階段，要求滿足每個階段在關(guān)鍵點所處的狀態(tài)，盡量減少軟著陸過程的燃料消耗。根據(jù)上述的基本要求，建立數(shù)學(xué)模型解決下面的問題：（1）確定著陸準備軌道近月點和遠月點的位置，以及嫦娥三號相應(yīng)速度的大小與方向。（2）確定嫦娥三號的著陸軌道和在 6 個階段的最優(yōu)控制策略。（3）對設(shè)計的著陸軌道和控制策略做相應(yīng)的誤差分析和敏感性分析。二問題的分析2.1 問題一由于題目已經(jīng)明確給出準備軌道的形狀參數(shù)，可以通過已有的物理知識與幾何關(guān)系明確計算出近月點與遠月點處速度大小和相對

11、于月面的速度方向。為了預(yù)期落點與著陸軌道在同一個平面內(nèi)，且準備軌道過月心，可以大致確定軌道所在平面有無數(shù)多個，無法確定近月點和遠月點的位置，因此需要其他條件來推測。本題現(xiàn)有的條件下，只有通過著陸軌道逆推近月點，并結(jié)合地月軌道制動的實際情況綜合考慮，才能得到。對于著陸過程一，由相關(guān)報道以及 NASA 在 1976 年提出的線性正切制導(dǎo)率1，得知主減速階段通常都是恒推力作用在軌道切線上，且嫦娥三號主減速階段實際也是保持著最大推力依照這一定律進行制導(dǎo)。我們由此出發(fā)，通過二體模型，結(jié)合已知條件，建立微分方程組，通過計算機模擬降落軌跡即可求出降落弧線距離，從而反推近月點，對稱得出遠月點。2.2 問題二由

12、問題一已經(jīng)得出近月點，即開始降落點位置，也知道每一階段的狀態(tài)，因此，降落軌道大致范圍基本確定，但六個過程的精確路徑是要通過策略優(yōu)化來控制的。由于燃料消耗表現(xiàn)在推力在時間上的積累量，即減小推力作用的沖量，即可優(yōu)化燃料消耗。對第一個過程，由于推力很大且歷時較長，因此燃料消耗主要體現(xiàn)在這一階段，對應(yīng)的， 優(yōu)化策略也應(yīng)重點體現(xiàn)，由于有二體模型，建立微分方程模型，并由初值條件以及階段限定條件，可以寫出非線性約束條件，本問題及轉(zhuǎn)化為軌道優(yōu)化中的非線性規(guī)劃問題，一般可通過成熟的 SQP 算法可以得到全局最優(yōu)解，但本題采用了更為常見也相對傳統(tǒng)的遺傳算法， 逼近全局最優(yōu)解，得出最優(yōu)方案。對于第二階段，僅僅是為了

13、是水平速度將為 0，且推力迅速減小，由上一階段的優(yōu)化結(jié)果，得出末速度水平分量，由于沖量可分解，則此階段分為水平方向和豎直方向分別優(yōu)化， 即分為了兩個變速直線運動模型，簡化了優(yōu)化模型，可以由這兩個局部最優(yōu)解加和得到該階段的全局最優(yōu)解。第三個階段水平速度初始為 0，經(jīng)過對月面成像分析，制定平坦度評價體系，選擇距離中心點最近且滿足平坦度要求的區(qū)域中心為粗調(diào)整目標降落點。由于這一段終點懸停，速度減為 0，因此可以對該段推力進行優(yōu)化，從而局部燃料最優(yōu)。第四階段懸停，精細成像并分析，同樣制定平坦度評價體系并選擇距離中心點最近且滿足平坦度要求的區(qū)域中心作為最終目標降落點，修正軌道。結(jié)束時水平速度依然為 0，

14、因此同樣存在優(yōu)化過程。第五階段與第六階段是減速至 0 然后自由落體的過程，針對減速階段也可考慮優(yōu)化，可經(jīng)過簡單討論得到結(jié)果。最后根據(jù)各個階段的最優(yōu)方案，模擬出嫦娥三號著陸軌道即可。2.3 問題三為了分析設(shè)計軌道和控制策略的誤差與敏感性，有必要制定誤差指標并考慮各部分誤差對于結(jié)果的最大影響，敏感性也同樣需要這一思路，各個階段的細微變化會對結(jié)果產(chǎn)生影響的衡量。另一方面，或許還有必要尋找參考物以顯示該方案的好壞。三模型假設(shè)1. 由于月球自轉(zhuǎn)速度為 27.3d，十分緩慢，而著陸過程僅有十幾分鐘，因此在本題中月球不考慮自轉(zhuǎn)；2. 由于月球扁率很小，可認為月球為球體，半徑以平均半徑為準，并且引力場分布均勻

15、；3. 由于側(cè)面姿態(tài)調(diào)整噴射裝置對燃料影響很小，為簡化模型，認為飛行器變換姿態(tài)的過程不消耗燃料；4. 認為飛行器變換姿態(tài)是瞬間完成的；5. 由于著陸時間較短，所以諸如月球引力非球項、日月引力攝動等影響因素均可忽略不計；6. 推力大小可瞬間改變；7. 燃料除供給推力外，無任何其他耗散方式；8. 月球空氣稀薄，不考慮任何摩擦力；9. 由于預(yù)定著陸點海拔-2641m，因此在著陸過程中所使用的高度均不應(yīng)是海拔高度， 而是相對于著陸點海拔的高度。四符號說明符號符號意義G引力常量M月球質(zhì)量m嫦娥三號質(zhì)量n比沖gm月球表面重力加速度F推力J燃耗P平坦度評價指標H綜合評價指標F平均坡度評價指標五模型建立與求解

16、5.1 問題一5.1.1 模型建立為了確定近月點和遠月點，正常的思路是求出橢圓軌道所在平面以及橢圓長軸在空間中的位置，但本題僅給出了軌道兩點的高度信息以及平均半徑，還有根據(jù)常識推斷出的預(yù)定著陸點在橢圓平面內(nèi)，除此之外并無其他信息，因此從現(xiàn)有條件準確判斷兩點的位置是不可能的。因此，本題應(yīng)采用逆推思路，由著陸軌道的第一個階段反推近月點。先對軌道參數(shù)進行分析，已知近月點高度 Hc=15km，遠月點高度 Hf=100km，月球平均半徑為 R=1737.013km，軌道為橢圓。如圖 5.1.1.1.圖 5.1.1.1 近月軌道示意圖（為了方便示意，本圖不符合比例）由開普勒定律，任何橢圓天體軌道的中心天體

17、一定在橢圓的一個焦點上。則由橢圓幾何性質(zhì)，可以近似得出方程：a+c=Hf+Ra-c=Hc+R聯(lián)系萬有引力定律與牛頓第二定律，可以列出嫦娥三號在近月點和遠月點的運動學(xué)方程：ìMmv2rïG 0 = mfï (a+ c)20ífïMmv2G 0 = m c îrcï(a- c)200其中 M = 7.3477´1022 kg 為月球質(zhì)量，G=6.672´10-11 為引力常量，m = 2.4t 為在近月軌道上飛行器的質(zhì)量， r f 與 rc 為遠月點和近月點的曲率半徑。由橢圓的幾何性質(zhì)，在長軸兩端點處的曲率半

18、徑分別為：rìïï cí2= b = (a - c)(1+ e) aa2ïr =ïî f聯(lián)立以上各式，帶入?yún)?shù)，得：= (a + c)(1- e)b近月點速度大小vc = 1692.46m / s ； 遠月點速度大小vf = 1614.15m / s 。由假設(shè) 2，月球視為球體，則近月點處速度方向方向應(yīng)平行于月面且為了著陸方便，方向指向預(yù)定著陸點所在方向，遠月點處同樣平行月面但速度方向與之相反。下面開始研究主減速過程，該過程由近月點開始，切向速度vq (0) = vc ，徑向速度vr 等于0，由切線正切制導(dǎo)率及相關(guān)報道可知，

19、該階段推力始終保持最大推力，且始終調(diào)整使推力與速度方向相反，即 F º 7500N ，F(xiàn) v ，最終應(yīng)大致達到預(yù)定著陸點目標上空，且豎直方向速度為 57m/s，高度共下降 12000m。由于有端點限制條件，可以通過物理知識建立微分方程組，即采用微分方程模型進行分析求解。由于飛行器繞月球表面飛行，且切向速度較大，在月心極坐標系下必須考慮向心力和科里奧利力，且涉及轉(zhuǎn)動變量，較為復(fù)雜。這里采用參考系轉(zhuǎn)化，將極坐標系轉(zhuǎn)化為非慣性參考系，并對該系中所有質(zhì)點提供一個向上的離心力，大小由水平速度與質(zhì)點到月心距離決定：v2F = m x dr這樣，問題就轉(zhuǎn)化到了平面直角坐標系下，并且可以對速度和受力

20、進行正交分解，如圖。圖 5.1.1.2 非慣性系下力學(xué)分析圖示由此可得微分方程模型：ì dvyGMv2ï= -+ x + a × sin bï dtr 2rïïdvx = a × cos b - vx vyïdtrïdr = vídtym0 - mtïïïa =Fïïïvy = tan bîvx其中 r 為飛行器到月心的距離， 為推力與水平方向的夾角，a 為推力產(chǎn)生的加速度， 整個過程中質(zhì)量是勻速減少的，減少系數(shù)由比沖與推力定

21、義，即 F =n m ，本題中n =2940N / kg 。至此，本問題已經(jīng)完全抽象為數(shù)學(xué)語言，并且有初始參數(shù)：vy = 0 ， vx = 1692.46m / s ， r = 1749.37km ， b =0理論上可求得當高度下降 12000m 時速度方向大小以及飛行距離，但這個方程組要想解出一個描述運動軌跡的函數(shù)是十分困難的，因此，需要另尋它路，逼近它的解。5.1.2 模型求解本題由方程組反推運動方程是十分困難的，但這類問題就像是解決 NP 難題一樣，無法由問題得到結(jié)論，但可以“猜測”結(jié)論從而驗證問題的正確性，繼而得到想要的數(shù)據(jù)?；仡櫛绢}，之所以無法解決是因為微分方程的連續(xù)性，這種連續(xù)性計

22、算機很難求解，但若是讓該組方程做某種近似，使之轉(zhuǎn)化為差分方程，從而利用差分方程的離散型導(dǎo)出差分方程組。輸入初始條件，使用計算機模擬其運動過程，最后打點畫出折線圖。由微元思想可知， 當時間被微分成足夠小的時候，模擬的運動軌跡可認為無限逼近真實運動軌跡。微元化后的差分方程依然可以有運動學(xué)規(guī)律得出：ì v= v- (a× cos b+ vxn vyn ) × tïxn+1ïïxnnnnrGMv2ïvyn+1ï= vyn + ( 2rn- an × sin bn- xn ) × t rnï&#

23、239;an =ïïím= mFmn- m × tïn+1nïrn+1 = rn - ynï12ïxn+1 = xn + vxn × t -ï2 an cos bn × tïïyn+1ïïî= yn + vyn× t - 1 a2n sin bn× t 2其中 為時間微元小量，由于減速段大約幾百秒，為了精確描繪運動曲線 一般取 0.1s。帶入初始條件：F = 7500N ， m0 = 2400kg ， b0 =0 ，

24、 x0 = 0 ， y0 = 0 ， vx0 = 1692.46m / s ， vy0 = 0即可通過迭代法，求出在 y=3000m 時的 x 的值，該值即為飛行器轉(zhuǎn)過的角度所對應(yīng)的月面弧長。但是仿真結(jié)果表明，當高度降至 3000m 時，速度并不為 57m/s，而且遠大于這一數(shù)值，顯然，是我們采用的控制方案不合題意。重新考慮這一過程，應(yīng)該是豎直方向分力不夠，導(dǎo)致豎直速度分量增加過快。之后在了解了嫦娥三號實際登月過程后，我們發(fā)現(xiàn)這一過程中推力并非與速度方向相反，而是始終與速度反方向保持一個角度 ，并向曲線內(nèi)側(cè)傾斜。于是我們通過改變 的大小，試圖輸出大量運動軌跡簇，這一過程由計算機模擬生成。通過計

25、算機仿真模擬，得到非慣性參考系下的運動軌跡簇。其中可以發(fā)現(xiàn)，在 取時，高度 3000m 處豎直方向速度為，在誤差允許的范圍內(nèi)可以接受，由此得出符合題意的主減速階段飛行軌跡圖，如圖 5.1.1.4.圖 5.1.1.4 主減速階段飛行軌跡圖（非慣性參考系）圖 5.1.1.5 水平與豎直速度變化示意圖圖 夾角變化示意圖得到弧長為，對應(yīng)轉(zhuǎn)過角度為，所用時間為。由于推力大小不變，由時間即可得出此階段燃耗。該階段內(nèi)水平方向與豎直方向速度變化如圖 5.1.1.5.可以看出，水平速度不斷下降，而豎直方向速度先增大后減小，合速度 57.12m/s，符合題意。水平總位移為 385.21km。雖然該階段末并不是準確

26、到達預(yù)定著陸點上空，需要經(jīng)過快速調(diào)整階段才能準確到達， 但題目表述中明確指出該處以基本到達目標上空，再加之段末水平速度分量相比初始速度已經(jīng)十分小，快速調(diào)整階段又很迅速，即快速調(diào)整階段對到目標點上空距離的影響十分小，相比主減速段弧長可以忽略，即認為主減速末段已經(jīng)到達目標上空。由此，我們結(jié)合弧長與轉(zhuǎn)角可以計算近月點所在范圍。由于該角度可在以預(yù)定落點為圓心的球面圓上任意選取，考慮到當天月球月面范圍，降落過程應(yīng)基本暴露在有光月面一側(cè)， 以便拍攝或信息收集，并且一般選擇軌道時為了方便計算控制，軌道平面會與落點所在經(jīng)線平面基本重合，即近月點經(jīng)度也為 19.51°W，緯度應(yīng)更加靠近赤道，從而保證降

27、落過程全程在光側(cè)面，最終確定緯度位置為 31.50°N 。由月球的球?qū)ΨQ性，遠月點緯度南北對調(diào)，經(jīng)度東西對調(diào)，角度互補。最終結(jié)論，近月點位置為 19.51°W ，31.50°N ，距月面 15km，速度為 1692.46m/s；遠月點位置為 160.49°E ，31.50°S ，距月面 100km，速度為 1612.15m/s。5.2 問題二本題為典型的多因素優(yōu)化問題，優(yōu)化目標是減少燃料損耗，并且通過圖像分析選擇最為平坦的地區(qū)著陸。由于燃料損耗直接與推力有關(guān)，即 F =n m ，其中n =2940m / s 為燃料比沖，也就有ò Fd

28、t = n ò mdt該式左邊為推力對時間的積累量，即推力所做的沖量；右邊為燃料損耗對時間的積累， 即一定時間內(nèi)的燃料損耗總量。由于比沖為常數(shù)，因此推力沖量與燃耗成正比例關(guān)系，從而使優(yōu)化目標與力學(xué)量直接聯(lián)系，便于分析。另外，由動量定理mdv = Fdt ，在近似處理下，也可將燃耗與速度變化聯(lián)系在一起。5.2.1 過程一模型建立與求解由 5.1 計算過程可知，主減速段是著陸過程用時最長，燃耗最多的階段。該階段的主要任務(wù)是消除較大的初始水平速度, 因此推進劑消耗優(yōu)化是該階段優(yōu)化的主要設(shè)計目標，也是整個著陸軌道燃耗優(yōu)化的重中之重。嫦娥三號的主減速階段是從近月點開始下降到距離預(yù)定地點 3 k

29、m 處，這一過程的時間比較短，僅有幾百秒的時間，所以可以不考慮月球引力攝動。月球自轉(zhuǎn)速度比較小 , 也可忽略 。所以，可以僅僅對月球和嫦娥三號進行分析，將問題轉(zhuǎn)化為二體模型，如圖。以月心為原點，建立平面直角坐標系，設(shè)嫦娥三號的月心距為r ，極角為q ，角速度為w ，質(zhì)量為m ；又設(shè)v 為嫦娥三號沿r 方向上的速度， F 為主減速發(fā)動機的推力（在本題中為固定值），發(fā)動機推力與當?shù)厮骄€的夾角即為推力的方向j ，比沖為 ISP ；設(shè)月球的引力常數(shù)為 m 。圖 5.2.1.1 二體模型力學(xué)分析示意圖首先，從運動學(xué)角度進行分析。由動力學(xué)基本方程可得dr = tdv ，所以 dr = v dt同理，可以

30、得到角度變化的關(guān)系式為 dq = wdt圖 5.2.1.2 受力分析示意圖再從動力學(xué)角度進行分析。如圖，對嫦娥三號進行受力分析可知。在徑向上，受到萬有引力 FGM mmm=月 =以及主減速發(fā)動機在徑向上的分力Fsinj 。引r 2r 2由此，根據(jù)牛頓第二定律可以得到：=mmmv2r 2rFsin j = m dvdt聯(lián)立，可得 mm - Fsin j = mv -dv= w2m，利用公式vr 化簡該式可得：r 2rdtdv = F sin j - m + rw2dtmr 2再對切向上的受力進行分析，此處需要注意到我們選擇參考系月球自身的自轉(zhuǎn)，是一個非慣性參考系。因此，嫦娥三號除了受到主減速發(fā)動

31、機推動力的分量Fcosj 的作用，還需要計入大小為2vw 的科里奧利力，根據(jù)右手定則可以確定其方向為切向。由牛頓運動學(xué)第二定律可以得到F cosj + 2vw = -m dvy 。同樣，利用vdty= rw 化簡可以得到：F * cosj + 2* v *wdw =- m.dtr由于飛船在運行過程利用燃料反沖制動，因此飛船得而質(zhì)量是不斷變化的。考慮比沖的定義為“火箭發(fā)動機單位質(zhì)量推進劑產(chǎn)生的沖量”可得：dm =- F。dtISP綜上所述，聯(lián)立式，即可得到完整的描述嫦娥三號與月球這一二體運動模型的方程組：dr = v dtdv = F sin j - m + rw2dtmr 2dq = wdtF

32、 cosj + 2vwdw =- mdtrdm =- F dtISP對于嫦娥三號燃料使用指標的衡量，我們可以反映到其動量上去，即將動量作為衡量燃tf FI料消耗的性能指標。由質(zhì)量變化的推導(dǎo)，我們可以得到 J = òt0 SPdt 。在軌道優(yōu)化過程中 , 由于各狀態(tài)變量的量級相差較大 , 尋優(yōu)過程中可能會導(dǎo)致有效位數(shù)的丟失2。歸一化處理可以克服這一缺點，提高計算精度。另外，由于對軌道的優(yōu)化也要求優(yōu)化變量盡可能地保持在相同的量級，故作以下處理：令： rref= r0, mref = m0 .mrrefrrefmrv Fm v2 m則： r =, v =, v=, I= I, F =, F

33、=ref ref, m =,rrefvrefrefSPSPFrefrefrrefmrefr3refm r w = w, t =t tref, tref= refvref, q = q .那么，嫦娥三號的動力學(xué)方程可改寫為：dr = v dt d v = F sin j - 1 + rw2d tmr 2dq = wdt F cosj + 2vwdw =- mdtrd m =- F dtISP又由第一問可以得到飛行器的初始條件和終端條件分別為:vr 0 = 0v = 0vq 0= 1692.46m / srfq fv= 57m / sr0 = 1749.372km注：此處的終端條件中，徑向速度和切

34、向速度為理想情況下的結(jié)果，實際只能得到最優(yōu)化的值。使用歸一化條件可以得到初始條件和終端條件分別變?yōu)椋簐r 0 = 0vrf = 0vq 0= 1692.46 m / s vrefvq f= 57vrefm / s0r = 1749.372 km rref在考慮燃料最優(yōu)的情況，并要保障水平速度幾乎減到0m / s 的前提下，可以選擇推力較大的發(fā)動機完成軟著陸任務(wù)，即 F = 7500N 。由此可以得到優(yōu)化的目標函數(shù)為：tf FFIJ = òt0dt=t f SPISP其中，推力 F 和比沖 ISP都是定值，所以即要求時間最短的方案。最優(yōu)控制問題的數(shù)值方法主要包括直接法和間接法，對上述最

35、優(yōu)控制問題，使用一種直接方法參數(shù)化方法進行求解3。即作如下的假設(shè)：推力方向角q 可以表示成一個多項式的形式：q = å a t6iii =0由此，對于所描述的問題可由一個有約束條件的非線性規(guī)劃問題描述。由于確定了嫦娥三號的推力，那么在確定變量為 a0 到a6 七個變量后，根據(jù)所建立的二體模型方程即可確定運行的軌跡，所以此處需要優(yōu)化的參數(shù)為 a0 到 a6 七個變量，問題的約束條件為歸一化后的f動力學(xué)方程以及初值條件和終端條件。由此需要求出燃料最優(yōu)化方案，即性能指標 J = F tISP的最小值。對以上轉(zhuǎn)化出來的問題，我們利用浮點數(shù)編碼的遺傳算法（FGA）進行求解。其步驟如下：1. 將

36、 7 個取值范圍給定的優(yōu)化參量按一定的浮點數(shù)編碼排列在一起成為一個個體， 隨機產(chǎn)生 2000 個這樣的個體作為初始種群；2. 通過編寫的適應(yīng)函數(shù)計算每一個個體的性能指標；3. 使用輪盤賭法作為選擇算子并對這 N 個個體進行排序；4. 選擇出若干個性能指標取值較優(yōu)的個體保留，并將其遺傳到下一代；5. 將個體隨機兩兩配對，按照既定的交叉概率進行交叉操作；6. 對每一個個體中的每一個參數(shù)，按照既定的變異概率進行編譯操作；7. 若滿足收斂條件則輸出最優(yōu)解，否則繼續(xù)進行編碼、評價、選擇、交叉和變異等操作。通過對約束條件的分析，可以選擇交叉概率為 0.6，變異概率為 0.05。-2通過計算機仿真模擬，由遺

37、傳算法規(guī)劃出的角度隨時間控制函數(shù)參數(shù)向量的取值為：(a0a6 ) = (7.58´107.01´10-54.30 ´10-6-8.12 ´10-86.03´10-10-1.60 ´10-121.40 ´10-15 )運動軌跡見圖 5.2.1.2.圖 5.2.1.3 運動軌跡圖圖 5.2.1.4 速度對比圖圖 5.2.1.5 角度變化圖當下降 12000m 時速度 v=57.12m/s， 水平方向與垂直方向分別為 vx = 48.63m / s ，vy = 29.96m / s ，符合題意。整 個 過 程 耗 時 413.2s

38、 ， 由 于 推 力 F 恒為 7500N ， 因 此 該 階 段 燃 料 消 耗 為J = Ft = 1055.39kg ，剩余質(zhì)量 1344.61kg，水平位移為 385.21km。1n5.2.2 過程二模型建立與求解快速調(diào)整階段主要任務(wù)有兩個：一是在主減速階段完成的基礎(chǔ)上，將嫦娥三號水平方向的速度減到0m / s ；二是快速調(diào)整嫦娥三號的姿態(tài)，使推力方向和所收到的引力方向基本一致。在該階段，考慮讓 F 的方向逐漸調(diào)整到豎直方向，并且在調(diào)整的過程中（即下落的600m 距離內(nèi)），將水平方向的速度減少到0m / s 。圖 5.2.1.6 受力分析示意圖由此，首先進行運動學(xué)分析。在水平方向上，速

39、度最終減小為有0m / s ，所以有dvx = adtx在豎直方向上，由加速度的定義，有這樣的關(guān)系：dvy = adty然后，在進行動力學(xué)分析。在水平方向上，根據(jù)牛頓第二定律可以得到：F cosj = max在豎直方向上，不妨設(shè)月球的重力加速度為 gm2，其值由黃金代換公式GM = gm R 可以確定。那么同理可以得到：mgm - F sinj = may至此，受力分析完畢。而另一方面，考慮到減速過程中燃料的消耗對質(zhì)量產(chǎn)生影響， 有：dm =- FdtISP綜合，式可以得到描述該階段運行的方程：F cosj = m dvxdtmgm- F sin j = m dvydtdm =- F dtIS

40、P由主減速階段的終值以及快速調(diào)整階段的要求可以得到該方程初始條件和終端條件分別為：vx0 = 48.63m / svxf= 0m / svy 0 = 29.96m / sm0 = 1344.61kgj f = 90ooq f = 90而力 F 的取值范圍為1500N <= F <= 7500N 。最優(yōu)控制問題的數(shù)值方法主要包括直接法和間接法，對上述最優(yōu)控制問題，使用一種直接方法參數(shù)化方法進行求解3。即作如下的假設(shè)：推力 F 可以表示成一個多項式的形式：6t fF = å a ti 。另外，在快速調(diào)整階段，燃料的性能指標同樣為 J =F dt 。ii =0òt0

41、SPI綜上所述，即可將求解快速調(diào)整階段的問題轉(zhuǎn)換成這樣的一個非線性規(guī)劃問題：需要優(yōu)化的參數(shù)為a0 到a6 以及j 、q 九個變量，問題的約束條件為方程組以及初值條件和終端條件。由此需要求出燃料最優(yōu)化方案，即性能指標 J 的最小值。同樣使用浮點數(shù)編碼的遺傳算法來解決該非線性規(guī)劃問題，由于具體步驟與主減速階段基本一致，故此處不再贅述。在這樣的算法下，求得的最優(yōu)解為：vx = 0m / svy = 29.96m / sj = 90oq = 90oF = 3928.115Nt = 20.027s J2 = 26.71kg5.2.3 過程三模型建立與求解粗避障段的范圍是距離月面 2.4km 到 100m

42、 區(qū)間，其主要是要求避開大的隕石坑，實現(xiàn)在設(shè)計著陸點上方 100m 處懸停，并初步確定落月地點。嫦娥三號在距離月面 2.4km 處對正下方月面 2300×2300m 的范圍進行拍照，獲得數(shù)字高程，嫦娥三號在月面的垂直投影位于預(yù)定著陸區(qū)域的中心位置。該高程圖的水平分辨率是 1m/像素，其數(shù)值的單位是 1m。由附件 3 所提供的灰度圖，獲取各像素點高程，統(tǒng)計各高程數(shù)據(jù)出現(xiàn)的個數(shù)，繪制高程數(shù)量分布曲線圖，如圖 。圖 5.2.3.1 高程數(shù)量分布曲線圖由高程圖大致可以發(fā)現(xiàn)，在高程約為 95m 時有最高點，120m 以上或 80m 以下數(shù)量明顯下降，據(jù)觀測可以認為平坦的區(qū)域的高程大致在 80m

43、 到 120m 之間，于是應(yīng)大致將高程劃分為三個組，低于 80m 的區(qū)域較暗，高于 120m 的區(qū)域較亮，輸出的圖像上大致可以分辨出較大的隕石坑，但是小隕石坑仍然無法與平地區(qū)分，因此需要從統(tǒng)計學(xué)角度合理劃分。通過對全部高程數(shù)據(jù)進行 K 均值聚類分析，并設(shè)置聚落數(shù)，經(jīng)過嘗試發(fā)現(xiàn)當聚落數(shù)為 5時，可以較為清楚的區(qū)分出隕石坑與平地，如圖 。圖 5.2.3.2 聚落色階圖聚類中心值分別為 92.4595，113.8753，138.4552，168.0559，44.6447，平坦區(qū)域聚類編號為 1。圖中可以發(fā)現(xiàn)，對于平坦度造成影響的障礙地形大致分為三類：過高障礙區(qū)域，過低障礙區(qū)域，條紋障礙區(qū)域。通過排除

44、 1 以外的聚類編號，可以準確排除前兩種障礙區(qū)域，但是第三種條紋障礙區(qū)域無法去除，地圖上會留下數(shù)個環(huán)形地帶。至此為止，圖像分析基本完成，但由于此處降落點的調(diào)整同時受制于落點平坦度與落點偏移量，即落點到圖像中心的直線距離，因為這一距離會導(dǎo)致飛行器水平?jīng)_量的積累，燃耗增加，偏離最優(yōu)解。因此，有必要制定數(shù)字評價手段，通過地圖中現(xiàn)有參數(shù)，對平坦度以及偏移量做出數(shù)字化的評價，綜合權(quán)重分配即可得出落點的綜合評價指標。下面先制定平坦度評價指標 P 。由圖像聚類結(jié)果，可以粗略將圖中像素點分為安全點與危險點。但考慮到 2300×2300 點數(shù)量實在太大，為了減少數(shù)據(jù)處理量，對地圖進行柵格化， 每 5m

45、×5m，即 25 像素點作為一個柵格單元，根據(jù)格內(nèi)安全點與危險點數(shù)量的多少定義該格為安全格或危險格，這樣，地圖便被分為了 460×460 格的二元化柵格網(wǎng)絡(luò)。以任意一個安全格為中心，搜索其外圍 8 個方格，確認其中無危險格后，繼續(xù)擴大一周考察范圍，如此循環(huán)，直到出現(xiàn)危險格停止，記錄無危險格的最大安全半徑 RS （以格數(shù)衡量，若第一圈就出現(xiàn)危險格，則認為半徑為 0 格）。 RS 越大，則認為該格所在區(qū)域越安全，P 越高。由于平坦度僅與該因素有關(guān)，因此，若要將 RS 轉(zhuǎn)化為實際距離長度，正比例系數(shù)應(yīng)取 5m，即 P = 5RS 。條紋障礙由于 P 非常小，可以有效排除。接下來考

46、慮落點對中心的偏移量的評估指標。由于最終優(yōu)化目標是燃耗最小，因此有必要先確定粗調(diào)整階段下落規(guī)則。由于該階段推力不算很大，時間較短，可在分析階段忽略質(zhì)量變化。假設(shè)水平偏移量為 X，下降高度 H3 =3200m 已知，在忽略質(zhì)量變化的前提下，水平方向運動可以與豎直方向運動在滿足 F 2 £ F 2 + F 2 £ F 2的約束下完全獨立，為了進一步簡minxymax化水平方向運動優(yōu)化過程，認為豎直方向為勻速下降過程，這樣，下落時間變?yōu)槎ㄖ怠?由于動量定理我們知道，水平方向動量的變化就是該階段水平運動消耗的燃料，m 一定， 則速度vx 變化越小，沖量越小，燃耗越小，若在下面的v

47、- t 圖像 5.2.3.3 中，可以很明確的知道，速度曲線 2 要比速度曲線 1 燃耗小，由此，可以建立這樣的局部優(yōu)化策略。圖 5.2.3.3 速度曲線對比圖通過該階段的始末速度限制可知，初始水平速度為 0，豎直速度vy0 由上一階段得出；段末懸停，速度為 0。由此可見，要節(jié)省水平方向的燃耗，應(yīng)最快的增大vx ，接著勻速運動， 最后對應(yīng)最快的減小vx 至 0，從而完成該偏移過程。因此在水平加速與減速階段所受推力應(yīng)有 Fx =F2max- F 2 =ma ，同時，豎直方向勻減速yv2H運動，應(yīng)有 F 為常數(shù)，水平方向位移與速度關(guān)系為 X = v t - x ，同時t = 3 ，沖量變化即可yy

48、xav表示燃耗變化，有J = 2mvx ，基本認為 J µ v 且聯(lián)立前面各式可以發(fā)現(xiàn)v 與 X 有以下關(guān)系：3n3xxHv2x 3 v vy- x = X a即vx 與 X 有過原點的二次曲線遞增關(guān)系，X 在 500m 以內(nèi)較為平緩，可認為是線性關(guān)系，但隨著 X 繼續(xù)增大， vx 將明顯增大，付出的燃耗代價也會加速上升。因此，對于水平偏移量的評估權(quán)重應(yīng)為非線性加權(quán)。設(shè)偏移量評估指標為l ，則為了盡量淡化邊緣區(qū)域的影響，應(yīng)以偏移中心的距離 X 為自變量建立指數(shù)遞減模型，即l = A - BepX ，其中 A 決定了中心點指標大小，B 決定了指標衰減幅度，p 決定了指標靠近邊緣時的衰減

49、速度。經(jīng)過與平坦度評價指標的具體值相比較， 最終決定，A 取 200，B 取 2.4，p 取 0.0027， l = 200 - 2.4e0.0027 X ，方程曲線如圖 5.2.3.4.圖 5.2.3.4 偏移量評價函數(shù)曲線圖最終制定綜合評估指標 H = P + l ，作為一個方格對于降落點優(yōu)先級的考量，指標評估值越高，則該方格作為最佳著陸點的優(yōu)先級越高。對柵格化地圖遍歷并計算后得到如圖 5.2.3.5 的結(jié)果。圖 5.2.3.5 綜合評估地圖可以發(fā)現(xiàn)，圖中靠近中央有一峰值明顯比周圍峰值高，參照圖 5.2.3.3 中坐標，該最優(yōu)方格的坐標為（1275m，1000m），水平偏移量 X 為 19

50、5.26m，由之前公式可以粗略計算出水平方向最大速度vx = 1.87m/s，用時t3 = 102.51s，燃耗 J3 » 86.97kg。5.2.4 過程四模型建立與求解精細避障段的區(qū)間是距離月面 100m 到 30m。要求嫦娥三號懸停在距離月面 100m 處， 對著陸點附近區(qū)域 100m 范圍內(nèi)拍攝圖像，并獲得三維數(shù)字高程圖。分析三維數(shù)字高程圖， 避開較大的隕石坑，確定最佳著陸地點，實現(xiàn)在著陸點上方 30m 處水平方向速度為 0m/s。該數(shù)字高程的水平分辨率為 0.1m/像素，高度數(shù)值的單位是 0.1m。此階段由于地圖僅有 100m×100m 的大小，水平偏移量最大也只

51、是 71m，對燃耗的影響實在太小，作為最終綜合評估指標中所占的權(quán)值實在太小，可忽略不計，但由于要精確降落，對于地圖平坦度的分析所要考慮的因素不僅有安全半徑大小 RS ，而且還要同時考慮到安全格內(nèi)的平均坡度j 。由于在過程三中我們已經(jīng)詳細討論了最大安全半徑的獲得方法，這里不再贅述。本次柵格化考慮到像素為 0.1m，而飛行器本身體積較大，為了綜合考慮，應(yīng)將單元格大小取為 5m×5m，即 50×50.聚類分析結(jié)果如圖 5.2.4.1.圖 5.2.4.1 聚類色階圖本次五類聚落中心值按照編號由小到大分別為 41.2624，111.6007，90.2679，157.2982， 圖中

52、3 表示平坦地區(qū)，即安全格。最大安全半徑分析計算，結(jié)果如圖 5.2.4.2.圖 5.2.4.2 最大安全半徑色階圖下面重點探討如何得到柵格內(nèi)的平均坡度。首先，需要通過一定方法擬合出每一個方格內(nèi)平均坡面，再由該坡面法向量方向計算出該坡面法向量與豎直方向的夾角（弧度制）作為衡量該坡面平均坡度的標準。對圖像做如此處理后得到的坡度分析如圖 5.2.4.3.圖 5.2.4.3 平均坡度色階圖下面制定坡度評價指標F 。由于代表平均坡度的夾角越小越好，并且應(yīng)突出小坡度的比重，因此，依然考慮非線性加權(quán)法，本次加權(quán)可直接對弧度進行統(tǒng)計歸一，找出上下界，并將下界作為 0，上界作為 1，然后進行數(shù)值翻轉(zhuǎn)，即可得到弧

53、度由小到大對應(yīng)指標由 0 到 1的坡度評價指標，即F = 1- j Î (0,1) ，然后將其作為優(yōu)化系數(shù)與平坦度評價指標相耦jmax - jmin合即可完成非線性加權(quán)，最終所得到的綜合評價指標 H = F × P 作為著陸點優(yōu)先級的評估指標，得到的結(jié)果如圖 5.2.4.4.圖 5.2.4.4 綜合評價指標色階圖可以看到，圖中有兩部分區(qū)域明顯高于其他區(qū)域，最終選擇指標最大值，以圖 5.2.4.2 中坐標系參考坐標位于（88m，56m）的方格作為最終著陸點，水平偏移 X = 38.47m。為了獲得此階段最優(yōu)策略，各方向應(yīng)按照先盡快加速至勻速在盡快減速，力學(xué)方程和限制條件如下：

54、ìv2ïv t - x = Xxïxaïïívy t -ï2 (a2v)y= Hy - gmïFï(a2 + a2 ) = max 4ïîxym目標函數(shù)為 J4 = 2mvx + mvy + mgmt 要盡可能小。在該動態(tài)規(guī)劃下可以得出最佳數(shù)值分配為：豎直緩速下降所用時間t = 25.65s，水平最大速度為vx = 1.44m/s，燃耗 J4 = 20.68kg，段末質(zhì)量m4 = 1210.25kg。5.2.5 過程五模型建立與求解緩速下降階段的區(qū)間是距離月面 30m 到 4m。該階段的主要任務(wù)是控制著陸器在距離月面 4m 處的速度為 0m/s，即實現(xiàn)在距離月面 4m 處相對月面靜止，之后關(guān)閉發(fā)動機，使嫦娥三號自由落體到精確落點。由上一個階段末狀態(tài)可知，飛行器水平速度分量已減為 0，則顯然本次運動為一個減速運動，為了保障安全性，應(yīng)將速度和加速度減至 0，則該階段應(yīng)設(shè)計為勻減速直線運動。由于此階段時間較短，推力不大，質(zhì)量減少（幾十 kg）相比飛行器本身質(zhì)量（1t）來說可以忽略不計，則由動量定理，設(shè)此階段初始狀態(tài)豎直速度分量為v5 ，質(zhì)量為m5 ，則ò Fdt = m5