- 相關(guān)推薦
數(shù)學建模國賽A題優(yōu)秀論文
嫦娥三號軟著路軌道設(shè)計與控制策略
摘要
本文主要為分階段研究嫦娥三號的軟著陸軌道設(shè)計與最優(yōu)控制策略。
建立模型一確定近月點和遠月點的位置，以及嫦娥三號速度大小與方向。首
先以月球中心為坐標原點建立空間坐標系，根據(jù)計算的作用力可知地球影響較小，故忽略不計。然后將嫦娥三號軟著陸看作拋物線的運動過程，計算在最大推力下的減速運動，求得月面偏移距離為，由此計算出偏移角度為15.25°。從而得出近月點和遠月點的經(jīng)緯度分別為（34.76°W，44.12°N）和（34.76°E，44.12°S）。最后在軟著陸的橢圓軌道上，由動力勢能和重力勢能的變化，計算出嫦娥三號在遠月點和近月點的速度分別為，沿軌道切線方向。
建立模型二和模型三確定著陸軌道和在6個階段的最優(yōu)控制策略。模型二主
要對主減速階段和快速調(diào)整階段進行初步分析。模型三分六個階段確定軌道和最優(yōu)控制策略，主減速階段建立目標函數(shù)燃料，假設(shè)推力最大,將最優(yōu)燃耗軟著陸問題轉(zhuǎn)化為最短時間控制問題，然后采用擬牛頓法和四階Admas預(yù)測-校正得到；快速調(diào)整階段采用重力轉(zhuǎn)彎制導，在假設(shè)條件下對嫦娥三號進行受力分析，得到嫦娥三號的動力學模型，然后通過開關(guān)控制得到燃耗最優(yōu)控制，并畫出仿真圖；粗避障階段采用多項式制導，通過初始狀態(tài)和末端狀態(tài)反解多項式系數(shù)進而求取標稱軌跡，然后將避障區(qū)域網(wǎng)格化，比較網(wǎng)格內(nèi)的方差大小確定最優(yōu)區(qū)域范圍；精避障階段需在滿足本文提出的避障原則式下搜索全局最優(yōu)解，以網(wǎng)格區(qū)域總體得分作為目標函數(shù)，得到最優(yōu)區(qū)域為坐標附近，并以螺旋搜索法搜索安全半徑的個數(shù)。其余階段僅對其做簡單物理分析后繪制出六個階段的著陸軌道。
建立模型四做相應(yīng)的誤差分析和敏感性分析。首先以模型二為基礎(chǔ)進行誤差
分析，當主減速階段的推力、初始質(zhì)量變化時，計算嫦娥三號質(zhì)量和燃料消耗速率的變化趨勢。再以模型三為基礎(chǔ)進行分析，對初始高度變化前后主減速階段的的偏角和和著陸軌道進行對比分析并計算誤差。然后進行敏感性分析，主要利用蒙特卡洛分析著陸軌道的粗避障階段和精避障階段月面不同地形高度，對嫦娥三號降落時所需調(diào)整概率大小的影響，接著分析嫦娥三號著陸占地面積大小對著陸調(diào)整概率的影響。
關(guān)鍵字：拋物線、燃料、擬牛頓法、Admas、網(wǎng)格化、蒙特卡洛模擬
1.問題重述
嫦娥三號于2013年12月2日1時30分成功發(fā)射，12月6日抵達月球軌道。
嫦娥三號在著陸準備軌道上的運行質(zhì)量為2.4t，其安裝在下部的主減速發(fā)動機能夠產(chǎn)生1500N到7500N的可調(diào)節(jié)推力，其比沖（即單位質(zhì)量的推進劑產(chǎn)生的推
力）為2940m/s，可以滿足調(diào)整速度的控制要求。在四周安裝有姿態(tài)調(diào)整發(fā)動機，在給定主減速發(fā)動機的推力方向后，能夠自動通過多個發(fā)動機的脈沖組合實現(xiàn)各
種姿態(tài)的調(diào)整控制。嫦娥三號的預(yù)定著陸點為19.51W，44.12N，海拔為-2641m。
嫦娥三號在高速飛行的情況下，要保證準確地在月球預(yù)定區(qū)域內(nèi)實現(xiàn)軟著陸，關(guān)鍵問題是著陸軌道與控制策略的設(shè)計。其著陸軌道設(shè)計的基本要求www.stanzs.com：著陸準備軌道為近月點15km，遠月點100km的橢圓形軌道；著陸軌道為從近月點至著陸點，其軟著陸過程共分為6個階段，要求滿足每個階段在關(guān)鍵點所處的狀態(tài)；盡量減少軟著陸過程的燃料消耗。
根據(jù)上述的基本要求，請你們建立數(shù)學模型解決下面的問題：
（1）確定著陸準備軌道近月點和遠月點的位置，以及嫦娥三號相應(yīng)速度的大小與方向。
（2）確定嫦娥三號的著陸軌道和在6個階段的最優(yōu)控制策略。
（3）對于你們設(shè)計的著陸軌道和控制策略做相應(yīng)的誤差分析和敏感性分析。
2.問題的分析
本文所研究的問題一主要為基礎(chǔ)計算和物理知識,首先我們需要根據(jù)預(yù)定的
著陸點的經(jīng)緯度確定軌道，然后通過拋物線的運動計算出在月球著陸時的水平路程，然后計算出偏移角度，據(jù)此確定近月點的經(jīng)緯度，而嫦娥三號的著陸軌道為過月球中心點的橢圓軌道，所以遠月點的經(jīng)緯度和近月點對稱，則可以由近月點計算出遠月點的經(jīng)緯度。最后因為在著陸軌道上衛(wèi)星的能量守恒，則可以通過勢能和動能的轉(zhuǎn)換來計算嫦娥三號的速度和方向。
本文所研究的問題二主要為過程的最優(yōu)控制和建立嫦娥三號軟著陸軌道。因
為嫦娥三號的軟著陸主要分為六個階段，所以此問應(yīng)分為六個階段來求解。主減速階段采用燃料最優(yōu)制導律來分析，建立著陸坐標系，將最優(yōu)燃耗軟著陸問題轉(zhuǎn)化為最短時間控制問題，然后得到目標函數(shù)；快速調(diào)整階段采用重力轉(zhuǎn)彎制導，對嫦娥三號進行受力分析，得到嫦娥三號的動力學模型，然后計算出燃耗最優(yōu)控制，并畫出仿真圖；粗避障階段采用多項式制導，首先列出加速度、速度、位移的多項式，然后通過初始狀態(tài)和末端狀態(tài)反解多項式系數(shù)進而求取標稱軌跡；精避障階段首先設(shè)定嫦娥三號的體型大小，然后處理數(shù)據(jù)的數(shù)量級不同，最后在整個降落區(qū)域的范圍內(nèi)搜索最優(yōu)著陸點；由于在緩速下降和自由落體階段中，發(fā)動機已經(jīng)關(guān)閉，故僅對其做簡單物理分析。最后通過整個分析得出總的著陸軌道。
本文所研究的問題三主要為著陸軌道和控制策略做誤差分析和敏感度分析，
需要對問題二所設(shè)計的著陸軌道和控制策略中的發(fā)動機推力、初始速度、初始高度進行誤差分析。然后進行敏感度分析，即對著陸軌道的粗避障階段和精避障階段月面不同地形高度對嫦娥三號降落時所需調(diào)整概率大小的影響，最后分析嫦娥三號著陸占地面積大小對著陸調(diào)整概率的影響。
3.模型的假設(shè)
假設(shè)一：嫦娥三號與月球均不受其他行星及衛(wèi)星的影響
假設(shè)二：不考慮月球繞地及其他星球的公轉(zhuǎn)和月球的自轉(zhuǎn)
假設(shè)三：將月球近似的看做標準球體
假設(shè)四：嫦娥衛(wèi)星的燃料消耗主要是在著陸的主減速階段
假設(shè)五：軟著陸的四、五、六階段著陸軌跡基本在同一平面內(nèi)
4.符號與公式的約定和說明
： G=為引力常量，m、M分別為兩物體質(zhì)量，R為兩物體距離，為兩
物體間的作用力
： 為物體質(zhì)量，為物體在作用下產(chǎn)生的加速度
： 軟著陸起始速度
： 加速度
： 平拋產(chǎn)生的距離
： 物體的動能（
： 物體的重力勢能（
： 嫦娥三號的推力
： 偏好系數(shù)
： 降落地點總體得分
： 第段離散段的平均加速度
由于本文使用參數(shù)和公式較多，其他公式和符號在具體模型中再做說明。
5.模型的建立與求解
5.1模型一的建立
5.1.1模型的假設(shè)
由萬有引力公式計算，再由牛頓第二定律計算地球和月球在近月點和遠月點處的重力加速度。
三號與月球影響很小，故可忽略不計。所以本模型只考慮月球?qū)︽隙鹑柕挠绊憽?/p>
5.1.2模型的分析
根據(jù)附件2給出的軟著陸過程示意圖，即嫦娥三號將在近月點15公里處以拋物線下降，相對速度從每秒1.7公里逐漸降為零。整個過程大概需要750秒，我們將其看作勻減速運動過程。利用matlab繪制嫦娥三號繞月飛行的三維動態(tài)圖，更直觀的反應(yīng)嫦娥三號的環(huán)月飛行，如圖3（源程序見附錄）：
圖2 嫦娥三號繞月軌道坐標圖 圖3 嫦娥三號環(huán)月飛行
同時由附件二所給的嫦娥三號著陸區(qū)域和著陸點示意圖可知，只要保證嫦娥三號的著陸區(qū)域在虹灣著陸區(qū)，則認為著陸成功。
為保證嫦娥三號以最大概率降落到精準的著陸點和虹灣著陸區(qū)，經(jīng)分析后得出，選擇以北緯44.12°作為軟著陸的繞月軌道。在這種確定緯度的繞月軌道中，月球?qū)︽隙鹑柕娜f有引力，可以分解為兩個方向。一個是繞月的向心力，一個是與繞行面相切的力，則選擇最終狀態(tài)為繞赤道運行更為準確。故根據(jù)實際分析，嫦娥三號的繞月平面應(yīng)與南北極軸重合。
圖4 嫦娥三號繞月飛行軌道分析
5.1.3模型的建立與計算
據(jù)了解，嫦娥三號主發(fā)動機是目前中國航天器上最大推力的發(fā)動機，能夠產(chǎn)生從1500牛到7500牛的可調(diào)節(jié)推力，故可根據(jù)推力范圍求取嫦娥三號的加速度范圍。并用最大的加速度計算平拋產(chǎn)生的距離。
主減速段看作平拋運動：
起始速度
加速度的取值范圍
平拋產(chǎn)生的距離 （
圖5 嫦娥三號拋物示意圖
由上圖，并結(jié)合計算所得的拋物距離，得到準備著陸的點與軟著陸點相差15.25°，即可算出近月點的經(jīng)緯度，同時根據(jù)對稱性，又可求得遠月點的經(jīng)緯度。
由附件所給條件可知距離月球表面15km時，速度的大小為，則此速度看作近月點速度，在穩(wěn)定的軌道下，從近月點到遠月點可看作重力勢能和動能相互轉(zhuǎn)換的過程，而遠月點距離地球表面為100km，可以計算重力勢能的變化，即可算出遠月點的速度：
(1)
根據(jù)以上公式可得出近月點與遠月點的速度（速度方向沿軌道切線方向），連同經(jīng)緯度，如下表所示：
表6 近月點、遠月點位置與速度
5.2模型二的建立
5.2.1模型的分析
本模型主要對主減速階段和快速調(diào)整階段進行初步分析
首先分析嫦娥三號在此階段的的受力情況，假設(shè)受力與豎直方向的夾角為：
圖7主減速階段受力分析圖 圖8 不考慮質(zhì)量變化時的受力分析
利用動量守恒定律可得：
(2)
(3)
由題目和附件可知，嫦娥三號在運行過程中有燃料的消耗，本模型分為兩種情況考慮，一種為考慮質(zhì)量變化，另一種為不考慮質(zhì)量變化。由于主減速階段燃料消耗很大，故作為質(zhì)量變化考慮；而快速調(diào)整階段速度很小，質(zhì)量變化很小，故作為質(zhì)量不變考慮。
考慮質(zhì)量變化（主減速階段），推力大小
此階段的燃料的消耗量為
不考慮質(zhì)量變化（快速調(diào)整階段）：由于值較小，可以通過姿態(tài)調(diào)整發(fā)動機進行微調(diào),假設(shè)此階段質(zhì)量的變化較小，則可以假設(shè)質(zhì)量基本保持不變。
通過受力分析，可得到以下分析式：
最后得到燃料消耗為
(4)
5.1.2模型的建立
建立目標規(guī)劃函數(shù)，計算最少的燃料消耗。由分析階段的計算可以得出總?cè)剂舷牧浚?/p>
(5)
由表達式可以畫出總?cè)剂舷牧颗c質(zhì)量和時間的關(guān)系
圖9 總?cè)剂舷牧颗c時間的關(guān)系
由圖可以看出，嫦娥三號的質(zhì)量隨時間遞增而減少，而燃料的消耗隨著時間遞增而增加。
5.3模型三的建立
本模型為分階段深入分析嫦娥三號的著陸軌道和在6個階段的最優(yōu)控制策略。 5.3.1主減速階段制導控制律（燃料最優(yōu)率制導[2]）
? 模型的準備
擬牛頓法是求解非線性優(yōu)化問題最有效的方法之一。擬牛頓法只要求每一步迭代時知道目標函數(shù)的梯度。通過測量梯度的變化，構(gòu)造一個目標函數(shù)的模
型使之足產(chǎn)生超線性收斂性。構(gòu)造目標函數(shù)在當前迭代的二次模型和割線公式
預(yù)估—校正算法的方法包括三步四階Adams外插法和三步四階Adams內(nèi)插法為了保證計算得精度，本文采用內(nèi)插法
? 模型的分析與建立 嫦娥三號主減速階段從距離月球表面15km開始，由初
速度為 開始主減速。建立二維模型描述嫦娥三號在此階段的運動。令月心O為坐標原點, y 指向動力下降段的開始制動點, x 向著陸器的開始運動方向，見下圖：
圖10 著陸坐標系
由坐標系可建立嫦娥三號的質(zhì)心動力學方程，描述如方程組（6）：
（6）
式中: ,,和 分別為嫦娥三號的月心距、極角、角速度和質(zhì)量;
為嫦娥三號沿方向上的速度;
為制動發(fā)動機的推力(固定的常值或0) ;
為其比沖;
為月球引力常數(shù);
為發(fā)動機推力與當?shù)厮骄�的夾角即推力方向角。
動力下降的初始條件由霍曼變軌后的橢圓軌道的近月點確定,終端條件為嫦娥三號在月面實現(xiàn)軟著陸。令初始時刻,終端時刻 不定,則此過程的約束條件可以表示為方程組（7）：
（7）
? 對的求解 月球軟著陸的最優(yōu)軌道設(shè)計就是要在滿足上述初始條件和終端
約束條件的前提下, 調(diào)整推力大小和方向,使得嫦娥三號實現(xiàn)燃料最優(yōu)軟著陸,則設(shè)燃料最優(yōu)目標函數(shù)為表達式（8）：
（8）
在無奇異情況下,推力應(yīng)為開關(guān)控制。要么以最大推力工作,要么以最小推力工作。但為了簡化問題,采用常值推力假設(shè)，即認為制動發(fā)動機一直以最大推力工作。這一方法一方面有利于優(yōu)化，另一方面可降低發(fā)動機復(fù)雜性。采用常值推力假設(shè)后,月球最優(yōu)燃耗軟著陸問題轉(zhuǎn)化為最短時間控制問題,即尋找實現(xiàn)軟著陸的最短時間，求解步驟如下:
： 確定一終端時間,滿足條件
： 求解無約束最優(yōu)控制問題狀態(tài)方程式,終端時間為，性能指標為：
（9）
其中下標表示在時刻的取值。
： 根據(jù)終端能量特性修正，然后返回，直到。
終端時刻的初始值估計，由于軟著陸時著陸器能量為零,可知推力作用主要是抵消能量，將該能量等效為動能,則可推出等效速度為
假設(shè)采用脈沖推力模式,將該速度抵消需要消耗的燃料量為
而對于實際的有限推力模式，與相對應(yīng)的時間為
（10）
式中為發(fā)動機燃料秒流量
最終得計算結(jié)果為：
因脈沖推力比有限推力消耗的燃料量少，所以使得該計算結(jié)果偏小。 ? 目標函數(shù)的求解 第二階段垂直方向上的減速最大值為
由文獻可知,為使衛(wèi)星在第六階段自由落體，則快速調(diào)整階段的速度范圍為：
假設(shè)主減速階段衛(wèi)星以一定角度提供向上的推動力，則等效速度為
由于值較小，故可以忽略不計。
此問題為終端時間固定型無約束最優(yōu)控制問題,本模型將其轉(zhuǎn)化為非線性規(guī)劃問題,然后借助于擬牛頓法和四階Admas 預(yù)測-校正積分格式快速求解。為保證優(yōu)化精度,轉(zhuǎn)化方法采用計算量稍大但精度較高的直接離散化方法。
直接離散化方法將整個最優(yōu)控制過程分成若干個時間段,時間段之間的端點稱為節(jié)點;選擇節(jié)點處的控制變量作為未知參數(shù),通過插值得到整個最優(yōu)控制過程的控制變量積分狀態(tài)方程;根據(jù)這些控制變量積分狀態(tài)方程形成目標函數(shù),得到一個無約束數(shù)學規(guī)劃問題。具體如下：
(1) 將整個飛行時間分為N 個時間段,形成N+ 1 個時間節(jié)點 ( i = 0 ,1 , ?,
N) ,取時刻的控制量 為優(yōu)化變量,共有N + 1 個變量；
(2) 整個飛行過程的控制量可以通過在各時間節(jié)點處線性插值得到；
(3) 采用擬牛頓法和四階Admas預(yù)測-校正積分，得到從到 積分狀態(tài)方程(6)
和目標函數(shù)(9)。
圖11 偏角和垂直速度隨時間變化的趨勢
5.3.2快速調(diào)整段制導律（重力轉(zhuǎn)彎制導[4]）
? 模型的分析 由于在最終著陸段中，嫦娥三號的距月面距離只有 2 千米左
右，遠遠小于月球的半徑 1738 千米，因此在建模時可以忽略月球的曲率，將月面近似看為水平面；且考慮到在最終著陸段中嫦娥三號的切向速度只有幾十米每秒，設(shè)切向速度給嫦娥三號所帶來的離心加速度為，月球半徑為。因為嫦娥三號的切向速度為，則計算切向速度給嫦娥三號所帶來的離心加速度公式為：
因此可以忽略嫦娥三號的離心加速度，只考慮重力加速度。
? 模型的建立 假設(shè)嫦娥三號的下降軌跡在一個平面內(nèi)，設(shè)制動發(fā)動機的比沖
為，秒耗量為，嫦娥三號的垂直高度為，切向速度為，質(zhì)量為，制動發(fā)動機的推力方向與垂直方向夾角為。在以上假設(shè)條件下，我們對嫦娥三號進行受力分析，可以得到嫦娥三號的動力學模型為：
（12）
? 模型的最優(yōu)解 為了使嫦娥三號在最終著陸段中的燃料消耗達到最小，則設(shè)
嫦娥三號軟著陸燃料消耗為：
（13）
對于重力轉(zhuǎn)彎制導法下的軟著陸模型，推力的燃耗最優(yōu)控制是開關(guān)控制，而且開關(guān)次數(shù)最多不會超過 1 次。要實現(xiàn)嫦娥三號的終端狀態(tài)約束，嫦娥三號只能先進行自由落體，直到開關(guān)切換函數(shù)為 0 時，制動發(fā)動機工作，嫦娥三號進行制動減速，直至在到達月面時減速為 0,仿真圖如下所示：
圖12 快速調(diào)整階段運動狀態(tài)
5.3.3粗避障段制導律（多項式制導[5] ）
? 模型的分析 嫦娥三號軟著陸粗避障階段持續(xù)時間較短，所以需要設(shè)計有效
的制導律使探測器能在有限的時間內(nèi)跟蹤上標稱軌跡，外部環(huán)境的干擾是影響著陸精度的主要因素。所以，本模型首先給出了多項式，然后通過初始狀態(tài)和末端狀態(tài)反解多項式系數(shù)進而求取標稱軌跡，然后設(shè)計終端滑模制導律跟蹤標稱軌跡。
? 模型的建立 多項式形式的標稱軌跡規(guī)劃一般假設(shè)系統(tǒng)狀態(tài)變量為多項式，
基于邊界條件和著陸時間解相關(guān)系數(shù)。對于嫦娥三號粗避障階段，首先可以將著陸器的加速度表示為二次多項式的形式：
(14)
其中，和分別為待定常數(shù)矢量。對式（14）等式兩邊積分可以得到嫦娥三號的速度矢量和位置矢量的表達式為：
+ (15) + (16)
給定著陸時間和初末端狀態(tài)的情況下，可以解出：
=
? 模型的計算和分析 生成標稱軌道的仿真參數(shù)為著陸器在著陸點平移坐標
下的初始位置矢量 ,初始速度矢量，著陸時間為，將參數(shù)代入到式（17）可得常矢量為：
基于光學圖像的粗障礙檢測就是利用月球巖石和坑的圖像特征識別大障礙, 確定安全區(qū)域。根據(jù)巖石和坑的特征，本文選取避障原則如下式：
圖13 粗避障階段的等高線
將此區(qū)域圖片看做的矩陣，進一步分割為個的矩陣。根據(jù)組成地面高度的矩陣，利用var函數(shù)求解計算每一個矩陣的方差。方差的大小代表地面的平坦程度。
圖14 粗避障階段最優(yōu)著陸點
圖中白色區(qū)域為方差最小點，即為不考慮避障階段速度增量的值時，需要搜尋的最優(yōu)著陸區(qū)域。 5.3.4精避障階段
精避障階段，推力和姿態(tài)發(fā)動機的比沖較小且時間短，不將比沖燃料消耗計算在內(nèi)。為了在整個降落區(qū)域的范圍內(nèi)搜索最優(yōu)著陸點，將圖片區(qū)域網(wǎng)格化處理為的矩陣，選擇最優(yōu)區(qū)域的準則為總高度和總平坦度值的大小。
用Min-Max標準化方法消除數(shù)量級的不同
設(shè)置偏好系數(shù)表示區(qū)域總高度對降落點得分的影響，表示區(qū)域總平坦度對降落點得分的影響。則降落地點總體得分。
圖15 精避障階段的等高線
對著陸所占用的不同區(qū)域下的計算，得出結(jié)論在占用區(qū)域面積時，最優(yōu)點為的附近區(qū)域。
表28 精避障階段最優(yōu)降落點
根據(jù)需要著陸的大小，對整個各個區(qū)域進行搜索滿足的點，即為可選擇的降落點
5.3.5軌道的確定
上文對著陸軌道的六個階段進行分析，主減速階段嫦娥衛(wèi)星的速度和質(zhì)量變化最大，對軌道的計算也最為重要。對于緩速下降和自由落體階段，由于發(fā)動機已經(jīng)關(guān)閉，則對于最優(yōu)控制和軌道設(shè)計不必過多分析。通過前面四個階段的分析和自由落體的規(guī)律，得出最終的著陸軌道如下圖：
圖16 最終著陸軌道的設(shè)計
5.4誤差分析與敏感度分析
主要對模型三設(shè)計的著陸軌道和控制策略做相應(yīng)的誤差分析和敏感性分析。 5.4.1誤差分析
本模型主要分析發(fā)動機推力誤差、初始速度誤差、初始高度誤差等。 發(fā)動機推力誤差：主要分析為主減速階段推力變化和嫦娥三號初始質(zhì)量變化對嫦娥三號質(zhì)量和燃料消耗的影響。
首先設(shè)定嫦娥三號的推力為最大推力7500N,然后將分別乘以1.1、0.9，觀察的變化對嫦娥三號質(zhì)量和燃料消耗的影響，如下圖：
圖17 推力改變時的誤差分析
由圖可以看出，嫦娥三號的推力變化會引起嫦娥衛(wèi)星的質(zhì)量和燃料消耗的變化，推力越大，質(zhì)量改變越小，燃料消耗越少。
由題目所給條件可知嫦娥三號的初始質(zhì)量為 =2400kg，然后將嫦娥三號初始質(zhì)量乘以1.1和0.9，觀察此時嫦娥三號的質(zhì)量和燃料消耗的變化。
圖18嫦娥三號質(zhì)量改變時的誤差分析
由圖可知，嫦娥三號的初始質(zhì)量的變化會引起嫦娥三號的質(zhì)量和燃料消耗的變化，初始質(zhì)量越大嫦娥三號的質(zhì)量變化越大，燃料消耗的越多。
對主減速階段的初始速度和初始高度進行誤差分析，嫦娥三號的預(yù)定著陸點海拔為-2641m，則將主減速階段的高度設(shè)置為15Km至17.641Km之間。將其與原有狀態(tài)下的運動狀態(tài)相互比較。僅考慮切向速度變化，根據(jù)燃料最優(yōu)制導模型的計算方法，利用四階龍格-庫塔公式和擬牛頓法將主減速的30個階段嫦娥三號偏角的變化與原變化進行比較，如下圖：
圖19 偏角的變化
上圖藍線表示原的變化，綠線為改變切向速度時的變化，紅線為兩者的誤差，可以看出前期原偏角大于改變后的偏角，后期則相反。誤差也隨著時間變少。由誤差計算公式 ，計算偏角總誤差為-9.49% 。
根據(jù)已求得的偏角的的值，將主減速段運動路徑分割為30個階段，并將軌道離散化
圖20 初始高度變化時軌道的變化
圖21 初始高度變化的軌道離散化
圖18的紅線為原高度時軌道變化，粉紅線為改變原高度時的軌道變化。由誤差公式可得，在主減速階段的誤差為,誤差率為 。 已知嫦娥三號的初始比沖量為2940，將其分別乘以0.9、1.1，即改變比沖量，觀察嫦娥三號質(zhì)量和燃料消耗的變化。
圖22 比沖量變化時軌道特性的變化
由圖可以看出，比沖量的值越大，嫦娥三號的質(zhì)量變化越大，燃料消耗越大。 5.4.2 敏感度分析
粗避障段：粗避障段的范圍是距離月面2.4km到100m區(qū)間，其主要是要求避開大的隕石坑，實現(xiàn)在設(shè)計著陸點上方100m處懸停，并初步確定落月地點。 將附件所給圖片網(wǎng)格化為2300×2300的矩陣，本文根據(jù)處的月球高度，得到避障原則： （19）
使用matlab軟件并采用用蒙特卡洛的方法進行1000次仿真（源程序見附錄），模擬分析月面不同地形高度對嫦娥三號降落時所需調(diào)整概率大小的影響。
圖23 粗避障階段地形 圖24 粗避障階段不同降落高度所需調(diào)整
精避障段：精細避障段的區(qū)間是距離月面100m到30m。要求嫦娥三號懸停在距離月面100m處，對著陸點附近區(qū)域100m范圍內(nèi)拍攝圖像，并獲得三維數(shù)字高程圖。分析三維數(shù)字高程圖，避開較大的隕石坑，確定最佳著陸地點，實現(xiàn)在著陸點上方30m處水平方向速度為0m/s。 與粗避障一樣，在滿足同樣的避障原則下，分析月面不同降落地形高度對嫦娥三號降落時所需調(diào)整概率大小的影響。
圖25 精避障階段地形 圖26 精避障階段不同降落高度所需調(diào)整
在精確避障階段的避障原則下，為了研究嫦娥三號在降落時占地面積大小對軌道降落的敏感度的影響，選擇和這兩個數(shù)據(jù)作為嫦娥三號降落時的占地面積。并用這兩個數(shù)據(jù)在matlab軟件中做50次模擬比較（源程序見附錄）。
圖27 兩種不同著陸占地面積著陸比較
由圖可以看出，的占地面積的非調(diào)整降落次數(shù)高于。，由此推測著陸占地面積越大，可直接順利著陸的概率越小。
分別選用六個不同的占地面積，對其進行1000次模擬，計算出1000次模擬中無需調(diào)整即可順利著陸的次數(shù)，如下表：
表29 不同占地面積無需調(diào)整即可順利著陸的次數(shù)
的概率越大。
6.模型的評價與改進方向
6.1模型的評價
6.1.1模型一的評價
模型一分別以著陸點的經(jīng)度和緯度作為準備著陸軌道，選取經(jīng)度不變的軌道處于穩(wěn)定狀態(tài)，不需要產(chǎn)生推力，此種軌道保證了燃料消耗的最優(yōu)。選取緯度不變的軌道則保證了當平拋距離較大且難以精確確定時以最大概率降落在著陸區(qū)域內(nèi)。
6.1.2 模型三的評價
在模型三中建立的主減速階段燃料最優(yōu)目標規(guī)劃函數(shù)，利用時間逼近法快速求解月球最優(yōu)軟著陸問題。對于終端時間固定型最優(yōu)控制問題,將其直接離散化為非線性規(guī)劃問題,采用擬牛頓法和四階Admas預(yù)測-校正積分方法快速求解。此方法優(yōu)化精度較高,收斂速度快，比近年較為流行的智能算法(如遺傳算法等)，減少了計算量且更符合實際需求和精度要求。
對嫦娥三號軟著陸的其他階段也分別建立了動力學模型。并且分段建立了最優(yōu)目標函數(shù)。確定了著陸軌道。在精避障階段，綜合考慮了著陸位置的總高度和總平坦度，對不同數(shù)量級的數(shù)據(jù)標準化，設(shè)置偏好系數(shù)后對所有點進行全局搜索，得出了最優(yōu)降落策略和最優(yōu)降落點。 6.1.3 模型四的評價
模型四對誤差進行了多方面的分析，包括對的最大推力、初始速度的變化的軌道特性分析、進一步對主減速階段的偏角的趨勢分析、計算誤差和靈敏度。 6.1.4 模型的不足與改進方向
模型的不足：由于軌道的復(fù)雜多變性，本文簡化了模型的數(shù)學推導，將高度變化引起的軌道路徑長度變化忽略，只重點考慮和計算了主減速階段的軌道特性，造成了設(shè)計的軌道系統(tǒng)的誤差。且未對軌道路徑做出明顯的全局優(yōu)化。 改進方向：
1. 將著陸軌道的六個階段燃料量作為規(guī)劃函數(shù)，將自適應(yīng)遺傳算法與模擬退火算法相結(jié)合,形成一種自適應(yīng)模擬退火遺傳算法，增強軌道路徑設(shè)計的整體搜索能力。
2. 根據(jù)月球巖石和坑的特征, 設(shè)計了粗障礙識別和安全著陸區(qū)選取算法: 1) 圖像直方圖分析; 2)K 均值聚類; 3) 過亮障礙識別; 4) 過暗障礙識別; 5)紋理障礙識別; 6) 采用螺旋搜索算法確定每個單元格的安全半徑(圖3); 7) 根據(jù)安全半徑, 選取候選安全著陸點; 8) 評估候選安全著陸點避障所需的速度增量; 9) 根據(jù)安全半徑和速度增量評價值, 綜合確定安全著陸點.
7.參考文獻
[1] 張德豐，MATLAB數(shù)值分析，北京：機械工業(yè)出版社，2012。
[2] 趙吉松，谷良賢，高原，月球軟著陸軌道的時間逼近法快速優(yōu)化設(shè)計[J]，宇航學報，第29卷第5期：1-5，2008. 9。 [3] 朱建豐，徐世杰，基于自適應(yīng)模擬退火遺傳算法的月球軟著陸軌道優(yōu)化[J]，航空學報，第28卷第4期：2-3，2007.7。 [4] 于彥波，火星探測器動力下降段制導律研究[D]，哈爾濱，哈爾濱工業(yè)大學， 2013。
[5] 張仲滿，月球軟著陸的制導算法研究[D]，哈爾濱，哈爾濱工業(yè)大學， 2009。 [6]田青，常微分方程初值問題數(shù)值解的實現(xiàn)與分析，http://www.doc88.com/p-682406421063.html,2014.09.13。
[7]張洪華，梁俊等，嫦娥三號自主避障軟著陸控制技術(shù)[J]，中國科學：科學計算，第44卷第6期：2-4,2014。
附錄
蒙特卡洛分析不同降落地形高度選取調(diào)整概率 I2=imread ('data.tif'); p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y); pp2=double(p2); c=zeros(200,1); d=zeros(200,1); for j=1:1:200 for i=1:1000
aa=floor(2299*rand(1)+1);
bb=floor(2299*rand(1)+1); if pp2(aa,bb)
d(j,1)=d(j,1)+1; end end end
plot(1:200,c/1000,1:200,d/1000);
title('不同降落地形高度選取調(diào)整概率'); xlabel('地形高度');
ylabel('調(diào)整與非調(diào)整概率'); legend('不需調(diào)整','需要調(diào)整');
精避障階段50次模擬分析嫦娥安好占地面積大小對著陸的影響： d1=zeros(50,1); for j=1:50 for i=1:1000
aa=floor(979*rand(1)+11); bb=floor(979*rand(1)+11); aa1=aa-7; aa2=aa+7; bb1=bb-7; bb2=bb+7; c=0;
for cc=aa1:aa2 for dd=bb1:bb2
if pp2(cc,dd)>=50&pp2(cc,dd)
if c==15*15;
d(j,1)=d(j,1)+1; end end
end end end
plot(1:50,d)
快速調(diào)整階段運動狀態(tài) v1=-90;
a=(7500-1200*1.6)/1200; d=0;
for t=0.1:0.1:14
d=d+1;
x(1,d)=3000+v1*t+0.5*a*t^2; end
plot(0.1:0.1:14,x')
title('快速調(diào)整段制導律下運動狀態(tài)'); xlabel('時間/t'); ylabel('高度/m');
總?cè)剂舷牧颗c時間的關(guān)系 F=7500; m=2400; ve=2940; v=22; g=1.63; a=3.125;
m1=zeros(1,487); for i=1:487 m1=F/ve;
m=m-m1-m*g/ve; A(1,i)=m; F=m*a; end
subplot(1,2,1) plot(1:487,A) title('速度'); xlabel('時間/t');
ylabel('嫦娥衛(wèi)星質(zhì)量/Kg'); subplot(1,2,2)
plot(1:486,A(2:487)-A(1:486)) title('能料消耗的變化'); xlabel('時間/t'); ylabel('質(zhì)量/Kg');
快速調(diào)整階段最優(yōu)降落點
I2=imread ('data1.tif'); p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y); pp2=double(p2);
a=input('請輸入嫦娥號的大小a='); for aa=a+1:1:1000-a for bb=a+1:1:1000-a
AA=pp2(aa-a+1:aa+a,bb-a+1:bb+a); A(aa-a,bb-a)=sum(sum(AA));
if A(aa-a,bb-a)==2389
aa
bb
end
end
end
初始高度變化時軌道的變化
H1=15000;
H2=15000+2641;
X=0;
X2=0;
i=0:1:510;
I=0:1:600;
X=1700.*i-0.5*3.125*(i.^2);
X2=1700.*I-0.5*2.878*(I.^2);
Y=H1-22.*i;
Y7=min(H2-22*1.1*I);
Y5=H2-22.*I*1.1;
i1=510:559;
I1=600:640;
X1=57*(i1-510)+460590-0.5*2.*((i1-510).^2);
Y1=min(Y)-22.*(i1-510);
Y8=Y7-22*1.1.*(I1-600);
X3=460982:0.1:460990;
Y3=0:2700/80:2700;
Y_1=0:100:3500;
X_1=linspace(501967,500720,36)
plot(X,Y,'r',X2,Y5,'m',X1,Y1,'k',X3,Y3,'g',X_1,Y_1); title('初始高度變化時軌道的變化');
xlabel('和近月點的距離');
ylabel('距地面高度/m');
text(250000,11000,'主減速階段');
text(420000,3500,'快速調(diào)整階段');
text(380000,2000,'粗避障階段');
text(380000,1000,'其他階段');
網(wǎng)格化矩陣方差
I2=imread ('data.tif');
p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y);
pp2=double(p2);
A=zeros(100,100);
for m=1:22
for n=1:22
for i=(1:100)+m*100
for j=(1:100)+m*100
A((i-100)*m,(j-100)*n)=pp2(i,j);
c=var(A(:))
end
end
end
end
搜索可選擇的著陸點
I2=imread ('data1.tif');
p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y);
pp2=double(p2);
a=zeros(1000,1000);
for i=1:1000
for j=1:1000
if pp2(i,j)>=60&pp2(i,j)
a(i,j)=1;
else
a(i,j)=0;
end
end
end
bb=input('輸入bb=');
for m1=bb+1:999-bb
for m2=bb+1:999-bb
if sum(sum(a(m1-bb:m1+bb,m2-bb:m2+bb)))==(2*bb+1)^2 m1
m2
else
continue;
end
end
end
【數(shù)學建模國賽A題優(yōu)秀論文】相關(guān)文章：
數(shù)學建模論文07-02
數(shù)學建模論文模板01-25
數(shù)學建模范文03-13
數(shù)學建模論文模板07-22
數(shù)學建模在數(shù)學中的應(yīng)用04-29
數(shù)學建模范文（優(yōu)）03-13
數(shù)學建模策略的教學原則05-02
數(shù)學建模教學例談04-29
數(shù)學建模的學習心得11-07