數(shù)模競賽論文模板

1、嫦娥三號軟著陸軌道優(yōu)化模型摘要本文針對嫦娥三號軟著陸軌道最優(yōu)設(shè)計問題，確定了近、遠(yuǎn)月點的位置與速度，建立了嫦娥三號六個階段的最優(yōu)軌道控制模型，提出了相應(yīng)的最優(yōu)控制策略，最后做出了誤差分析和敏感性分析。針對問題一，本文建立空間直角坐標(biāo)系，在著陸準(zhǔn)備軌道平面內(nèi)建立動力學(xué)二階常微分方程模型，利用微元法的思想，求得近月點的經(jīng)緯度為（19.51W，31.29N），遠(yuǎn)月點的經(jīng)緯度為（160.49E，31.29S）。利用開普勒第二定律，得出嫦娥三號在近月點和遠(yuǎn)月點的速度大小分別為1692.2m/s,1614.4m/s，速度方向與橢圓切線方向相同。針對問題二，分別確定了嫦娥三號軟著陸六個階段軌道的最優(yōu)控制策略

2、。對于著陸準(zhǔn)備軌道，根據(jù)燃耗量最小的原則，借鑒霍曼轉(zhuǎn)移模型，得出嫦娥三號在此階段的軌道是月心為焦點，長半軸為1794.5km，短半軸為1794km的橢圓。對于主減速階段，根據(jù)動力學(xué)原理，建立軌跡優(yōu)化模型，用改進(jìn)的遺傳算法求解，得到該階段最低燃耗量為1060.71kg，軌道形狀為類拋物線。對于快速調(diào)整階段，將水平偏移量作為優(yōu)化目標(biāo)，建立微分方程模型，得到最小的水平偏移量276.3米。利用附件中的數(shù)字高程圖，分析得到各點的海拔。在粗避障階段，提出崎嶇度的概念，建立基于崎嶇度最小的水平軌道優(yōu)化模型和基于燃耗量最小的垂直軌道優(yōu)化模型，得到嫦娥三號在此階段的水平位移為234.31米，最小燃耗量為69.3

3、8千克。對于精避障階段，建立基于月面坡度最小、著陸器燃耗量最小的軌道優(yōu)化模型，解出嫦娥三號水平總位移為5米，最小燃耗量為14.29千克。在緩速下降和自由落體階段，利用動力學(xué)公式求解出運動時間為13秒。針對問題三，通過對著陸點和其它各關(guān)鍵點的位置進(jìn)行誤差分析，發(fā)現(xiàn)本文確定的著陸點與實際著陸點相差80千米，緯度相差2.14，偏差可以接受。最后依據(jù)主發(fā)動機作用力與運動反方向的夾角的變化對主減速階段和快速轉(zhuǎn)移階段進(jìn)行了敏感性分析。關(guān)鍵詞：優(yōu)化模型 微分方程 遺傳算法 MATLAB仿真34一、問題重述嫦娥三號于2013年12月2日1時30分成功發(fā)射，12月6日抵達(dá)月球軌道。嫦娥三號在著陸準(zhǔn)備軌道上的運行

4、質(zhì)量為2.4t，其安裝在下部的主減速發(fā)動機能夠產(chǎn)生1500N到7500N的可調(diào)節(jié)推力，其比沖（即單位質(zhì)量的推進(jìn)劑產(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。嫦娥三號在高速飛行的情況下，要保證準(zhǔn)確地在月球預(yù)定區(qū)域內(nèi)實現(xiàn)軟著陸，關(guān)鍵問題是著陸軌道與控制策略的設(shè)計。其著陸軌道設(shè)計的基本要求：著陸準(zhǔn)備軌道為近月點15km，遠(yuǎn)月點100km的橢圓形軌道；著陸軌道為從近月點至著陸點，其軟著陸過程共分為6個階段

5、，要求滿足每個階段在關(guān)鍵點所處的狀態(tài)；盡量減少軟著陸過程的燃料消耗。根據(jù)上述的基本要求，請你們建立數(shù)學(xué)模型解決下面的問題：（1）確定著陸準(zhǔn)備軌道近月點和遠(yuǎn)月點的位置，以及嫦娥三號相應(yīng)速度的大小與方向。（2）確定嫦娥三號的著陸軌道和在6個階段的最優(yōu)控制策略。（3）對于你們設(shè)計的著陸軌道和控制策略做相應(yīng)的誤差分析和敏感性分析。二、問題分析本題要求我們建立數(shù)學(xué)模型，研究嫦娥三號的著陸準(zhǔn)備軌道，著陸軌道與6個階段的最優(yōu)控制策略，并對建立的模型做相應(yīng)的誤差分析和敏感性分析。2.1問題一的分析問題一要求建立模型，確定著陸準(zhǔn)備軌道近月點和遠(yuǎn)月點的位置，并且求解出嫦娥三號在近、遠(yuǎn)月點的速度大小和方向。我們將著

6、陸過程看做發(fā)射的逆過程，將預(yù)定著陸點到其正上方3000米視為直線運動，將預(yù)定著陸點正上方3000米到近月點的軌跡視為類拋物線。以月心為坐標(biāo)原點，月心與近月點之間的連線為軸，著陸準(zhǔn)備軌道所在平面為平面，建立空間直角坐標(biāo)系。再建立平面直角坐標(biāo)系研究類平拋過程，建立微分方程模型，求解得到近月點的經(jīng)緯度。對近月點的位置進(jìn)行幾何分析，得到遠(yuǎn)月點的坐標(biāo)和經(jīng)緯度。嫦娥三號在近月點和遠(yuǎn)月點的運動方向均與橢圓軌道的長軸垂直，通過求得的近月點和遠(yuǎn)月點的具體位置，可以得到它在這兩個點的速度方向。利用面積速度定理和開普勒第二定律，求解出嫦娥三號在近月點和遠(yuǎn)月點的相應(yīng)速度大小。2.2問題二的分析問題二要求建立模型，確定

7、嫦娥三號的著陸軌道和在6個階段的最優(yōu)控制策略。首先，確定嫦娥三號的準(zhǔn)備著陸軌道，建立霍曼轉(zhuǎn)移模型，得出嫦娥三號的橢圓軌道。主減速階段是軟著陸過程用時最長，推進(jìn)力消耗燃料最多的階段。該段的主要任務(wù)是減速使嫦娥三號到達(dá)3km的速度為57米/秒。因此盡量減小燃料消耗是該階段的關(guān)鍵問題。利用標(biāo)稱軌道制導(dǎo)方法對嫦娥三號進(jìn)行最優(yōu)控制，建立含有約束條件的軌跡優(yōu)化問題模型，用改進(jìn)的遺傳算法來求解該優(yōu)化問題?？焖僬{(diào)整階段的主要任務(wù)是調(diào)整姿態(tài)。因此優(yōu)化調(diào)整姿態(tài)時水平移動距離最小。建立水平偏移量最小的優(yōu)化控制模型，并通過仿真求解出水平最小偏移量。粗避障階段，用MATLAB對附件中給的數(shù)字高程圖進(jìn)行處理，得到嫦娥三號

8、的在此過程的位移，將水平運動和豎直運動分開考慮，用崎嶇度函數(shù)表示水平方向的優(yōu)化模型，用燃耗量最小表示豎直方向的優(yōu)化模型。精避障階段最優(yōu)控制策略與粗避障的方法相似。在嫦娥三號緩速下降階段和自由落體階段，利用動力學(xué)可以求解緩速下降的推力和時間。2.3問題三的分析 問題三要求對問題二中設(shè)計的著陸軌道和控制策略做相應(yīng)的誤差分析和敏感性分析。對著陸點和其它各階段的關(guān)鍵點的位置做誤差分析，從而判斷模型的準(zhǔn)確性。在主減速階段和快速轉(zhuǎn)移階段中，分析主發(fā)動機的作用力與運動反方向之間的夾角的變化，判斷該夾角在兩種不同階段的敏感性。三、模型假設(shè)1、主減速階段姿態(tài)調(diào)整發(fā)動機不工作。2、嫦娥三號到達(dá)近月點時，速度方向與

9、該點在軌道的切線方向一致。3、假設(shè)嫦娥三號在類平拋運動中沿某一經(jīng)線飛行。4、假設(shè)類平拋運動的平面與橢圓軌道在同一平面。5、嫦娥三號運動的任意過程中，均不考慮向心力。4、 符號說明符號說明主發(fā)動機提供的減速動力嫦娥三號衛(wèi)星運行中的質(zhì)量主發(fā)動機提供的加速度嫦娥三號運行的速度主發(fā)動機提供的加速度與運行速度反方向的夾角嫦娥三號至月球中心的距離（注：其它未提及的符號在文中說明）五、模型建立與求解5.1問題一的模型建立與求解問題一要求確定著陸準(zhǔn)備軌道近月點和遠(yuǎn)月點的位置以及嫦娥三號相應(yīng)速度的大小與方向。我們認(rèn)為由以下步驟完成：步驟一：以月心為坐標(biāo)原點，月心與近月點之間的連線為軸，著陸準(zhǔn)備軌道所在平面為平面

10、，建立空間直角坐標(biāo)系。 步驟二：確定近月點的位置。將著陸過程看做發(fā)射的逆過程，將預(yù)定著陸點到其正上方3000米視為直線運動，將預(yù)定著陸點正上方3000米到近月點的軌跡視為類拋物線。分析這兩個過程，并結(jié)合模型假設(shè)，得到近月點的經(jīng)緯度。步驟三：對近月點的位置進(jìn)行幾何分析，確定出遠(yuǎn)月點的坐標(biāo)和經(jīng)緯度。步驟四：求解嫦娥三號在近月點及遠(yuǎn)月點的相應(yīng)速度大小和方向。5.1.1 近月點位置微分方程模型的建立首先確定近月點的位置。我們將軟著陸過程看做發(fā)射的逆過程，通過分析發(fā)射過程來推算出近月點的位置。為了簡化模型，我們假設(shè)從快速調(diào)整到自由落體5個階段嫦娥三號均只在豎直方向上運動。嫦娥三號的預(yù)定著陸點為月球表面上

11、（19.51W，44.12N）的點，因此設(shè)距離該點高為3000米的點為，近月點為點,近月點與月心連線與月面的交點為點。從近月點到點的軌跡為類拋物線，為了分析問題比較方便，以月心為坐標(biāo)原點，月心與近月點之間的連線為軸，著陸準(zhǔn)備軌道所在平面為平面，建立空間直角坐標(biāo)系。如圖1。 圖1 嫦娥三號軌道示意圖點的坐標(biāo)為（0,0,1752），求出點到點的類拋物線軌跡就可以確定點、點的相對位置和點的坐標(biāo)，從而確定出點的經(jīng)緯度，下面對間的軌跡進(jìn)行求解。 嫦娥三號在主減速階段受到月球引力和減速動力兩個力的影響，所以加速度可以分解為指向月心和沿運動軌跡切線兩個方向的加速度，分別記為和，嫦娥三號運動速度為。假設(shè)運動方

12、向與水平方向夾角為，重力加速度與水平方向夾角為，加速度與速度的反方向的夾角為。令的初值為0.1，通過遞增值來改變嫦娥三號在距離月面3km處的速度大小。當(dāng)嫦娥三號距離月面3km處的速度大小達(dá)到57m/s時，確定此時的角的大小。、和在運動的過程中不斷變化。在類拋物線運動中選取任意一點，進(jìn)行運動學(xué)分析，畫出示意圖如圖2，圖2中實線表示月面，虛線表示運動軌跡。圖2 主減速階段嫦娥三號在任意一點的運動示意圖對圖2進(jìn)行運動學(xué)分析并建立微分方程模型： （1）其中表示上述受力分析點的水平速度，表示其豎直速度。水平速度、豎直速度與合速度的關(guān)系為： （2）根據(jù)已知條件，速度分量和的初值、大小為： （3）根據(jù)牛頓第

13、二定律得到減速動力和沿運動軌跡加速度的關(guān)系： （4）在文獻(xiàn)1中給出了燃耗最優(yōu)著陸過程中推力奇異區(qū)間不存在的證明，即在某一時間段內(nèi)只能取1500N或7500N這兩個值。表示嫦娥三號在著陸準(zhǔn)備軌道上的運行質(zhì)量即為2.4噸，表示嫦娥三號運動過程中的質(zhì)量，是單位時間燃料消耗的公斤數(shù)，表示比沖，它是對一個推進(jìn)系統(tǒng)的燃燒效率的描述。首先定義末速度，在每一個階段結(jié)束時，嫦娥三號對應(yīng)的速度即為該階段的末速度。下面以主減速階段末速度大小無限接近于57m/s為目標(biāo)函數(shù)建立目標(biāo)規(guī)劃模型。 （5） （6）其中，和分別表示主減速階段末速度在軸和軸的分量，表示嫦娥三號運動的任一時刻，表示加速度和速度的反方向的夾角，通過仿

14、真得到最佳的角。表示時刻嫦娥三號距離月面的高度。5.1.2 近月點微分方程模型的求解（1）近月點的位置根據(jù)（1）式至（6）式，利用MATLAB（程序1，2）可以求解出近月點的經(jīng)緯度為（19.51W，31.29N）。最優(yōu)夾角為6.13。（2）遠(yuǎn)月點的位置根據(jù)近月點的位置確定遠(yuǎn)月點的位置。近月點高度為15公里、遠(yuǎn)月點高度為100公里，且嫦娥三號的著陸準(zhǔn)備軌道是以月心為焦點的橢圓。遠(yuǎn)月點距離月心的距離為1837km。遠(yuǎn)月點在圖1空間直角坐標(biāo)系中的坐標(biāo)為（0,0,1837）。遠(yuǎn)月點的經(jīng)緯度為（160.49E，31.29S）。（3）嫦娥三號在近月點和遠(yuǎn)月點相應(yīng)速度的大小和方向嫦娥三號在這兩點的運動方向均

15、為與長軸垂直的方向，并且在近月點的速度方向與類拋物線運動的初速度方向相同，遠(yuǎn)月點的速度方向與該方向相反。結(jié)合空間直角坐標(biāo)系，可以得到，沿近月點速度方向的單位方向向量為（1,0,0），即與軸正方向平行。沿遠(yuǎn)月點速度方向的單位方向向量為（-1,0,0），即與軸負(fù)方向平行。下面來確定嫦娥三號在近月點和遠(yuǎn)月點相應(yīng)速度大小。設(shè)月球的質(zhì)量為，嫦娥三號的著陸準(zhǔn)備軌道周期為，該橢圓軌道半長軸為，半短軸為，嫦娥三號運行時的線速度為，面積速度為，為嫦娥三號到月心間的距離。根據(jù)面積速度公式2： （7）面積速度公式是開普勒第二定律的定量形式，再根據(jù)開普勒第二定律，嫦娥三號在單位時間內(nèi)掃過的面積是常數(shù)，在周期內(nèi)，嫦娥三

16、號掃過的面積等于橢圓面積，所以有： （8）因而： （9）由解析幾何可知，其中是橢圓的半焦弦，代入（9）式得： （10）根據(jù)牛頓力學(xué)公式3： （11）聯(lián)立（10）和（11）式得到： （12）是萬有引力常量，其值為。代表月球的質(zhì)量，題目中已給出為。根據(jù)對著陸準(zhǔn)備橢圓軌道的分析計算，得到半焦弦為1793.5km，近月點，遠(yuǎn)月點的月心距分別為。將這些數(shù)值代入（12）式可以對應(yīng)求得近月點，遠(yuǎn)月點的速度大小分別為：5.2問題二的模型建立與求解問題二要求建立數(shù)學(xué)模型，確定嫦娥三號的著陸軌道和在軟著陸過程中的六個階段的最優(yōu)控制策略。根據(jù)模型假設(shè)，嫦娥三號做類拋物線運動時，它的運動軌跡和它在橢圓軌道上的運動軌跡

17、位于同一平面，即圖1中的平面上。在該坐標(biāo)系下分析6個階段的軌道和最優(yōu)控制策略。圖3 嫦娥三號六個階段運動示意圖5.2.1 著陸準(zhǔn)備軌道著陸準(zhǔn)備軌道的近月點是15km，遠(yuǎn)月點是100km。近月點在月心坐標(biāo)系的位置和軟著陸軌道形態(tài)共同決定了著陸點的位置。（1） 霍曼轉(zhuǎn)移模型的建立在進(jìn)入著陸準(zhǔn)備軌道之前，嫦娥三號運動的軌跡為以月心為圓心，遠(yuǎn)月點到月心的距離為半徑的圓。在太空動力學(xué)中，霍曼轉(zhuǎn)移軌道4是一種變換嫦娥三號軌道的方法，途中只需一次引擎推進(jìn)，相對地節(jié)省燃料。圖4 霍曼轉(zhuǎn)移軌道圖4表示嫦娥三號從高軌道也就是環(huán)月軌道進(jìn)入霍曼轉(zhuǎn)移軌道即橢圓軌道。嫦娥三號在原先軌道上瞬間減速后，進(jìn)入一個以月球為焦點的

18、橢圓形著陸準(zhǔn)備軌道，其中，近月點是15km，遠(yuǎn)月點是100km。霍曼轉(zhuǎn)移軌道也就是著陸準(zhǔn)備軌道上嫦娥三號的總機械能等于動能和重力勢能的和，也等于當(dāng)該軌道半徑為半長軸時的重力勢能的一半，即： （13）軌道能量守恒方程為： （14）其中，為嫦娥三號的速度，為月球的標(biāo)準(zhǔn)重力參數(shù)，為嫦娥三號至月球中心的距離，為嫦娥三號著陸準(zhǔn)備軌道的半長軸。假設(shè)速度改變瞬間完成，霍曼轉(zhuǎn)移所需的速度變化量為： （15）（2） 霍曼轉(zhuǎn)移模型的求解和分別對應(yīng)近月點，遠(yuǎn)月點的月心距。即，。在第一次減速過程中，損失的能量為： （16）求得軌道轉(zhuǎn)移過程中損失能量最小為23628J。如圖4可求得霍曼轉(zhuǎn)移軌道的軌跡方程為： （17）所

19、以，第一階段的軌道為以月心為焦點，半長軸為1794.5km，半短軸為1794km的橢圓，此時滿足消耗能量最小的最優(yōu)控制策略。5.2.2 主減速段主減速段的區(qū)間是距離月面15km到3km。該階段的主要任務(wù)是減速，實現(xiàn)到距離月面公里處嫦娥三號的速度降到57m/s。（1） 主減速段軌道優(yōu)化模型的建立由于嫦娥三號的預(yù)定著陸點海拔為-2641m，所以假定3km為嫦娥三號距離著陸點海拔的相對高度。而由問題一可知，嫦娥三號在主減速段移動的水平距離為388km，虹灣著陸區(qū)是長為356km、寬為91km的長方形區(qū)域，假設(shè)從15km開始進(jìn)入主減速段時，地面海拔高度為0。所以類拋物線運動在豎直方向運動了14641米

20、。下面利用標(biāo)稱軌道制導(dǎo)方法5對嫦娥三號主減速段進(jìn)行最優(yōu)控制并確定該階段的軌道。標(biāo)稱軌道制導(dǎo)方法是指按照一定的目標(biāo)及約束條件預(yù)先設(shè)計一條標(biāo)稱軌跡，將著陸器的位置和速度的測量信息與標(biāo)稱軌跡相比較，導(dǎo)引著陸器飛向著陸目標(biāo)點。標(biāo)稱軌道制導(dǎo)由于其在燃耗最優(yōu)、可實現(xiàn)定點著陸以及對軌道約束的滿足性等方面的優(yōu)勢，是軟著陸較為理想的制導(dǎo)方法。首先不考慮月球自轉(zhuǎn)，并且忽略其它天體的引力影響。在月心坐標(biāo)系中，著陸器的動力學(xué)方程為： （18）（18）式中：分別為嫦娥三號的位置和速度在各坐標(biāo)軸方向下的分量，為制動發(fā)動機的推力大小，為制動發(fā)動機推力在各坐標(biāo)軸方向下的分量，為嫦娥三號和月心間的距離，為嫦娥三號質(zhì)量，為萬有引

21、力常量，為月球的質(zhì)量，為發(fā)動機比沖。針對以上著陸器在坐標(biāo)系下的三維動力學(xué)模型，給出含有約束條件的定點軟著陸軌跡優(yōu)化問題模型。若使著陸過程燃耗量最小，可選擇著陸器終端質(zhì)量最大為目標(biāo)函數(shù)，由質(zhì)量方程可得： （19）則末質(zhì)量最大指標(biāo)等價于： （20）主減速段的軌跡需滿足的初始和終端狀態(tài)條件為： （21）（3） 主減速段軌道優(yōu)化模型的求解此問題以主減速段軌跡的初始狀態(tài)，終端狀態(tài)為條件，以推力的大小和方向為決策變量的非線性優(yōu)化問題。由于目標(biāo)函數(shù)復(fù)雜，決策變量過多，無法用傳統(tǒng)的非線性規(guī)劃方法求解，所以本問題應(yīng)用遺傳算法來求得最優(yōu)解。遺傳算法6是一種基于自然選擇原理和自然遺傳機制的搜索算法，它是模擬自然界中

22、的生命進(jìn)化機制，在人工系統(tǒng)中實現(xiàn)特定目標(biāo)的優(yōu)化。傳統(tǒng)遺傳算法是對單目標(biāo)函數(shù)的優(yōu)化算法，無法適應(yīng)本題的雙目標(biāo)優(yōu)化問題，所以本文對遺傳算法做出改進(jìn)，將目標(biāo)函數(shù)的適應(yīng)度值進(jìn)行加權(quán)求和，再應(yīng)用傳統(tǒng)遺傳算法進(jìn)行求解（程序3）。其實現(xiàn)方法如下：（1）確定推力大小和方向的取值范圍，將推力大小和方向按每60秒的間隔進(jìn)行編碼，形成染色體。（2）對每一解的兩個適應(yīng)度函數(shù)進(jìn)行加權(quán)求和，得到總適應(yīng)度函數(shù)，適應(yīng)度函數(shù)應(yīng)為非負(fù)函數(shù)。（3）確定進(jìn)化參數(shù)群體規(guī)模、交叉概率、變異概率、進(jìn)化終止條件。下面給出改進(jìn)遺傳算法的算法框圖：圖5 改進(jìn)遺傳算法的算法框圖（4） 結(jié)果分析將改進(jìn)遺傳算法得到的主減速段的最優(yōu)軌道的相關(guān)信息列表見

23、表1，并且我們定義主減速階段結(jié)束時嫦娥三號的速度為末速度。表1 主減速段最優(yōu)軌道相關(guān)參數(shù)燃耗量末速度大小末速度方向飛行時間飛行距離終點經(jīng)緯度1060.71kg56.06m/s74.66464s468.96km19.51W，46.76N 由表1可以得到結(jié)果：該運行軌道初始點的經(jīng)緯度為（19.51W，31.29N），終點的經(jīng)緯度為（19.51W，46.76N）。在該運行軌道下，飛行距離為468.96km，飛行時間為464秒。在距離月面3000米上空時，末速度的大小為56.06m/s，末速度與豎直向下方向夾角為74.66。最低燃耗量為1060.71kg。所以，該階段末態(tài)的嫦娥三號質(zhì)量為1.34噸。主

24、發(fā)動機的推力隨時間變化的曲線圖，如下所示 圖6 主發(fā)動機推力隨時間變化曲線從圖6中我們可以看出，主發(fā)動機在主減速階段0到60秒?yún)^(qū)間的推力為1500N,隨后，推力突然增大到7500N并在此階段一直維持這個值。發(fā)動機推力與速度的夾角，如下圖所示圖7 發(fā)動機推力與速度反方向的夾角變化曲線由圖7可以看出，發(fā)動機的推力與速度反方向的夾角隨時間的增長而增大，且變化過程比較平緩，值也比較小。圖8 嫦娥三號主減速階段運動軌跡的動態(tài)模擬圖圖8是我們的仿真結(jié)果，由圖可以看出主減速階段的軌跡的確是一個類拋物線。5.2.3 快速調(diào)整段快速調(diào)整段的主要任務(wù)是調(diào)整探測器姿態(tài)，要求嫦娥三號從距離月面3km到 2.4km處將

25、水平速度減為0m/s。（1）快速調(diào)整段優(yōu)化模型的建立根據(jù)主減速段的結(jié)果可知，當(dāng)嫦娥三號距離地面3000米時，速度大小為57m/s，速度方向與豎直向下方向的夾角為74.66 度。假設(shè)姿態(tài)調(diào)整發(fā)動機的合力方向始終與運動方向相垂直，主發(fā)動機的減速動力方向與運動方向的反方向夾角為只研究水平方向，受力圖如圖7。 圖9 快速調(diào)整階段嫦娥三號受力分析示意圖在快速調(diào)整階段中，主要優(yōu)化目標(biāo)為調(diào)整姿態(tài)需要移動的水平距離，約束條件為水平方向的速度大小、速度與豎直方向的夾角等。建立水平偏移量最小的控制優(yōu)化模型： (22) (23)對于豎直方向，動力學(xué)方程為： (24)我們可以求得快速調(diào)整段的末速度。（2）快速調(diào)整段優(yōu)

26、化模型的求解對以上過程進(jìn)行仿真，得到快速調(diào)整段最優(yōu)軌道的相關(guān)信息列表見表2。表2 快速調(diào)整段最優(yōu)軌道相關(guān)參數(shù)燃耗量末速度大小末速度方向飛行時間水平偏移距離42.34kg4.90m/s046s276.3m從表2看出，嫦娥三號在快速調(diào)整段，調(diào)整力的大小和方向，可以得到最小水平偏移距離為276.3米，燃耗量為42.34kg，此時嫦娥三號衛(wèi)星的質(zhì)量為1.30噸。末速度豎直向下，大小為4.9m/s，飛行時間為46秒。通過仿真，我們得到主發(fā)動機的控制方案：在1至42秒反推力為1500N，在42秒至46秒反推力為7500N，力與運動反方向夾角始終為6.02。5.2.4 粗避障段粗避障段的范圍是距離月面2.4

27、km到100m，此階段要求避開大的隕石坑，實現(xiàn)在設(shè)計著陸點上方100m處懸停，并初步確定落月地點。在快速調(diào)整段我們得到衛(wèi)星位于2.4km高度時，速度大小為4.9m/s，方向豎直向下。為了便于分析，我們把嫦娥三號在此階段的水平和豎直方向上的運動分開考慮。水平方向的運動由姿態(tài)發(fā)動機提供動力，而豎直方向的運動由主發(fā)動機和月球的引力共同控制。（1）粗避障軌道控制模型的建立對附件3中距2400米處的數(shù)字高程圖進(jìn)行分析處理，應(yīng)用MATLAB（程序10）讀取圖片并分析灰度值，以附件3中圖像的左上角為原點，兩個邊分別為軸和軸，灰度值為Z軸建立坐標(biāo)系，利用MATLAB畫出嫦娥三號正下方月面23002300m的三

28、維圖像，如圖10。x/m y/m 圖10 嫦娥三號正下方月面23002300m的三維圖像粗避障過程中，認(rèn)為嫦娥三號能自動選取一個的正方形平穩(wěn)區(qū)域，并且會水平移動到該平穩(wěn)區(qū)域中心的上空100米處。定義在某一個的正方形區(qū)域中，各個點高度的標(biāo)準(zhǔn)差為該區(qū)域的崎嶇度，則崎嶇度函數(shù)為水平軌道控制模型。 （25）其中，表示第個點的高度，表示上述的正方形區(qū)域中90000個點的平均高度。以坐標(biāo)為（150,150）的點為正方形區(qū)域中心，求出該區(qū)域的崎嶇度。然后將正方形區(qū)域中心沿軸和軸以10m為單位進(jìn)行平移，得到各區(qū)域的崎嶇度，畫出崎嶇度的三維圖像如圖11。圖11 粗避障中各區(qū)域崎嶇度三維圖像利用MATLAB可以求

29、解出崎嶇度最小的點的坐標(biāo)為（970,1300），所以，嫦娥三號在粗避障過程中，將從中心點（1150,1150）水平移動至（970,1300）點處。相比的正方形區(qū)域中心坐標(biāo)（1150,1150），坐標(biāo)為（970,1300）的點是粗避障過程中選取的最佳平移點。嫦娥三號從（1150,1150）點水平移動到（970,1300）點，由于無法確定姿態(tài)發(fā)動機提供的動力范圍，我們假設(shè)它豎直下落過程需要時間為，則在時間內(nèi)嫦娥三號水平移動至目標(biāo)點即可。下面對此階段豎直運行軌道進(jìn)行優(yōu)化和控制。假設(shè)在粗避障段嫦娥三號的初始質(zhì)量為，運行過程中質(zhì)量為，并且在此過程中，重力加速度恒為。以該階段消耗燃料最少為目標(biāo)函數(shù)，末速度

30、減為0，下降距離為2300米等關(guān)系為約束條件，建立豎直軌道控制模型。 （26） （27）在這個粗避障階段，初速度為4.9m/s，末速度為0。（2）粗避障軌道控制模型的求解通過MATLAB仿真得到粗避障的最優(yōu)軌道信息，見表3。表3 粗避障段最優(yōu)軌道相關(guān)參數(shù)燃耗量飛行時間水平偏移距離69.38kg97s234.31m在豎直方向上，1至87秒主發(fā)動機提供動力1500N，88秒至96秒主發(fā)動機提供動力7500N，燃耗量最小為69.38kg，在粗避障段結(jié)束后，嫦娥三號衛(wèi)星的質(zhì)量為1.23噸。5.2.5 精避障段精避障段的區(qū)間是距離月面100m到30m。經(jīng)過粗避障階段，嫦娥三號懸停于目標(biāo)上方100m，對星

31、下月面進(jìn)行二維和三維成像，利用三維數(shù)字高程圖分析著陸點附近區(qū)域100m范圍內(nèi)地勢情況，避開較大的隕石坑，確定最佳著陸地點，實現(xiàn)在著陸點上方30m處速度為水平速度為0的目標(biāo)。同樣，為了便于分析，我們把水平運動和豎直方向上的運動分開考慮，姿態(tài)發(fā)動機提供水平方向上運動的動力，主發(fā)動機和月球的引力共同控制豎直方向的運動。（1）精避障目標(biāo)規(guī)劃模型的建立對附件4中距離月面100米處的數(shù)字高程圖進(jìn)行分析處理，以附件4中圖像的左上角為原點，兩個邊分別為軸和軸，像素點的灰度值為軸，利用MATLAB畫出嫦娥三號正下方月面100100m的三維圖像，分辨率為0.1m/像素。yx圖12 嫦娥三號正下方月面100100m

32、的三維圖像應(yīng)用DEM提取坡度7計算模型，對該100100m的地表進(jìn)行坡度計算。地表任意一點坡度是指經(jīng)該點的切平面與水平面的夾角，在數(shù)值上等于過該點的地表微分單元的法矢量與軸的夾角，即： （28）式中，是方向的高程變化率，是方向的高程變化率。用MATLAB對該高程圖進(jìn)行差分計算（程序10），得到每個點的坡度值，可以發(fā)現(xiàn)在以嫦娥三號為中心的1515m范圍內(nèi)存在坡度值較低的區(qū)域，區(qū)域中心坐標(biāo)為（55，51），可以作為降落區(qū)域。下面利用與粗避障相似的方法，對該階段豎直運行軌道進(jìn)行優(yōu)化和控制，該階段的初速度為0，末速度不為0。利用（26）和（27）式的目標(biāo)規(guī)劃模型，利用MATLAB可以求得最優(yōu)軌道的相關(guān)

33、信息。（2）精避障目標(biāo)規(guī)劃模型的求解通過MATLAB仿真得到精避障的控制策略為：1至19秒主發(fā)動機提供1500N的力，19秒至20.5秒主發(fā)動機提供7500N的力。燃耗量最小為14.29kg。此階段末態(tài)的衛(wèi)星質(zhì)量為1.22噸。5.2.6 緩速下降和自由落體階段緩速下降階段的區(qū)間是距離月面30m到4m。該階段的主要任務(wù)控制著陸器在距離月面4m處的速度為0m/s，即實現(xiàn)在距離月面4m處相對月面靜止。在4m處自由落體，實現(xiàn)軟著陸。在緩速下降階段，燃料消耗已經(jīng)不是主要優(yōu)化目標(biāo)。該階段主要考慮嫦娥三號運行的穩(wěn)定性，實現(xiàn)在距離月面4m處相對月面靜止。當(dāng)主發(fā)動機提供的力為一個恒力時，嫦娥三號勻減速運動，速度

34、均勻變化，運動穩(wěn)定性最高。利用動力學(xué)公式： （29） 求得該階段的總時間為6.68秒，落地速度為4.3m/s，由于探測器的著陸緩沖機構(gòu)的存在，實現(xiàn)了軟著陸。5.3問題三的解答問題三要求對問題二中設(shè)計的著陸軌道和控制策略做相應(yīng)的誤差分析和敏感性分析。（1）誤差分析由于測量儀器，方法，外界條件的影響等因素的限制，使得計算值與真值存在誤差，通過研究誤差，可以驗證我們所用建立模型和算法的正確性。各階段誤差表如下表四 各階段誤差表項目計算值實際值相對誤差主減速完成時速度值56.0576180357-0.016533017粗避障階段速度為0時的位置1001000緩速下降階段速度為0時的位置5.840.45

35、落地速度4.333.60.202777778 除此之外，還計算了落地點坐標(biāo)的誤差，其計算值為（19.51W,46.76N），真實值為（19.51W,44.12N），可以計算出誤差為80km。對上述誤差進(jìn)行分析，發(fā)現(xiàn)總體誤差較小，但是個別相對誤差較大，說明本文模型還有待優(yōu)化。（2）敏感性分析在最優(yōu)化方法中經(jīng)常利用敏感度分析來研究原始數(shù)據(jù)不準(zhǔn)確或發(fā)生變化時最優(yōu)解的穩(wěn)定性。對主減速階段主發(fā)動機提供的力與運動反方向夾角進(jìn)行敏感度分析，如圖13所示。圖13 主減速階段中夾角與速度的關(guān)系由圖13可知，主減速階段中，嫦娥三號高度降為3000米時的速度大小受主發(fā)動機提供的力與運動反方向夾角影響大，說明夾角的敏

36、感性較高。對快速轉(zhuǎn)移階段的夾角進(jìn)行敏感度分析，如圖14所示。圖14 快速轉(zhuǎn)移階段中夾角與飛行時間的關(guān)系由圖15可知，快速轉(zhuǎn)移階段主發(fā)動機提供的力與運動反方向夾角對飛行時間幾乎沒有影響，說明它的敏感性差。七、模型評價 模型的優(yōu)點1、本文采用標(biāo)稱軌道制導(dǎo)方法對嫦娥三號主減速段進(jìn)行最優(yōu)控制，標(biāo)稱軌道制導(dǎo)由于其在燃耗最優(yōu)、可實現(xiàn)定點著陸以及對軌道約束的滿足性等方面的優(yōu)勢，是軟著陸較為理想的制導(dǎo)方法。2、本文處理存在復(fù)雜約束條件的著陸軌跡優(yōu)化問題時，建立的定點軟著陸軌跡優(yōu)化模型在燃耗最優(yōu)和收斂速度上有很大優(yōu)勢，并且具有一定的穩(wěn)定性。3、采用遺傳算法求解最優(yōu)化模型，它不是從單個解開始，而是從問題解的串集開

37、始搜索，覆蓋面大，利于全局擇優(yōu)。 模型的缺點1、霍曼轉(zhuǎn)移雖然用到的能量少，可是時間代價大，如果嫦娥三號執(zhí)行時間要求高的任務(wù)，我們建立的霍曼轉(zhuǎn)移模型就不適用了。2、在主減速階段建立的定點軟著陸軌跡優(yōu)化問題模型魯棒性差。八、模型推廣本文通過建立微分方程模型，求解近月點的位置。該模型具有很大的實際意義，通過查閱更多的參考數(shù)據(jù)和相關(guān)資料，對模型作出改進(jìn)，從而可以利用該模型對近月點位置作出更精確的預(yù)測，從而更好地預(yù)測出著陸點位置。這對探月具有深遠(yuǎn)意義。建立霍曼轉(zhuǎn)移模型確定準(zhǔn)備著陸軌道，此模型在太空動力學(xué)中有廣泛的應(yīng)用，很多太空船為節(jié)省燃料在變換軌道都采用此方法。本文通過利用標(biāo)稱軌道制導(dǎo)方法，建立定點軟著

38、陸軌跡優(yōu)化模型，確定主減速階段的軌跡優(yōu)化策略。此模型適合處理存在復(fù)雜約束條件的各種著陸軌跡優(yōu)化問題。此外論文中大量使用了目標(biāo)規(guī)劃衡量決策的優(yōu)劣，這種模型在實際生活中有著廣泛的應(yīng)用，可以利用該模型處理生產(chǎn)計劃、環(huán)境保護(hù)、土地利用等各領(lǐng)域的決策問題。九、參考文獻(xiàn)1 林曉輝，于文進(jìn).基于凸優(yōu)化理論的含約束月球定點著陸軌道優(yōu)化J.宇航學(xué)報,2013,34(7):901-9082 吳樹鵬.利用幾何方法研究開普勒第二定律J.哈爾濱師范大學(xué)自然科學(xué)學(xué)報,2000,16(2):56-593楊文熊.現(xiàn)代牛頓力學(xué)M.上海:上海交通大學(xué)出版社，2011.4霍曼轉(zhuǎn)移軌道 ， http:/zh.wikipedia.or

39、g/wiki/霍曼轉(zhuǎn)移軌道 2014,9,145 梁棟，劉良棟，何英姿.月球精確軟著陸最優(yōu)標(biāo)稱軌跡在軌制導(dǎo)方法J.中國空間科學(xué)技術(shù),2011,6:27-356張文修，梁怡.遺傳算法的數(shù)學(xué)基礎(chǔ)M.西安:西安交通大學(xué)出版社，2003.7賈敦新，湯國安，王春，賈旖旎.DEM數(shù)據(jù)誤差與地形描述誤差對坡度精度的影響J.地球信息科學(xué)學(xué)報，2009，11（1）：43-49十、附錄程序所用軟件：MATLAB 程序1.1： Q1_try.mmine = inf;for gama = 0.100:0.001:0.120 tryQ1_paracurve; if abs(min(v)-57) 1737000 + 36

40、0 x(i) = x(i-1) + vx(i-1); y(i) = y(i-1) - vy(i-1); vx(i) = vx(i-1) + ax(i-1); vy(i) = vy(i-1) + ay(i-1); v(i) = sqrt(vx(i)2 + vy(i)2); theta(i) = atan(y(i)/(-x(i); alpha(i) = atan(vy(i)/(-vx(i); m(i) = m(i-1) - F/2940; a(i) = F/m(i); ax(i) = a(i) * cos(alpha(i) + gama) + g * cos(theta(i); ay(i) = -

41、a(i) * sin(alpha(i) + gama) + g * sin(theta(i); i = i + 1; if i 800 %防止陷入死循環(huán) break; endend clear d1 d2 d 43411184 t xc yc 程序2：Q1_paracurve.mclear F R a alpha ax ay i m t theta v vx vy x xc y yc% 初始化各變量初值gama = 0.107; %推力和軌跡的夾角consumption = 0;F = 7500;x(1) = 0;y(1) = 1752000;vx(1) = - 1692.2;vy(1) =

42、0;v(1) = sqrt(vx(1)2 + vy(1)2);m(1) = 2400;a(1) = F/m(1);g = 1.62;theta(1) = atan(y(1)/(x(1);%重力加速度方向和x軸夾角alpha(1) = atan(vy(1)/(-vx(1);%速度方向和x軸夾角ax(1) = a(1) * cos(alpha(1) + g * cos(theta(1);ay(1) = -a(1) * sin(alpha(1) + g * sin(theta(1);i = 2;% 開始迭代while sqrt(x(i-1)2 + y(i-1)2) 1737000 + 360 x(i

43、) = x(i-1) + vx(i-1); y(i) = y(i-1) - vy(i-1); vx(i) = vx(i-1) + ax(i-1); vy(i) = vy(i-1) + ay(i-1); v(i) = sqrt(vx(i)2 + vy(i)2); theta(i) = atan(y(i)/(-x(i); alpha(i) = atan(vy(i)/(-vx(i); m(i) = m(i-1) - F/2940; a(i) = F/m(i); ax(i) = a(i) * cos(alpha(i) + gama) + g * cos(theta(i); ay(i) = -a(i)

44、* sin(alpha(i) + gama) + g * sin(theta(i); i = i + 1; consumption = consumption + F/2940; if i 800 %防止程序陷入死循環(huán) break; endend% 畫出月球表面圖plot(x,y,-.)R = 1737000;t = 0 : pi/360 : 2*pi;xc = R * cos(t);yc = R * sin(t);hold onplot (xc,yc,r)xlabel(/m)ylabel(/m)axis equal% 求得軌跡長度和經(jīng)緯度s = 0;for cy = 1:i-2 d1 = x

45、(cy+1) - x(cy); d2 = y(cy+1) - y(cy); d = sqrt(d12 + d22); s = s + d;enddelta_degree = s / R / 2 / pi * 360;clear d1 d2 d 43411184 t xc yc 程序3.1：g_main.m % g_main.mmaxgen = 2000; % 迭代次數(shù)sizetop = 1000; % 種群數(shù)量lenchrom = ones(1,24); % the length of each variablep_cross = 0.70;p_mutation = 0.02; % 初始化種群F = round(rand(sizetop,12);F(F = 0) = 1500;F(F = 1) = 7500;chromB = round(84 + 50 *