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