2014高教社杯數(shù)學(xué)建模A題解法(共27頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上摘要本文針對嫦娥三號軟著陸軌道設(shè)計與控制策略的實際問題，以理論力學(xué)（萬有引力、開普勒定律、萬能守恒定律等）和衛(wèi)星力學(xué)知識為理論基礎(chǔ)，結(jié)合微分方程和微元法，借助MATLAB軟件解決了題目所要求解的問題。針對問題（1），在合理的假設(shè)基礎(chǔ)上，利用物理理論知識、解析幾何知識和微元法，分析并求解出近月點和遠月點的位置，即139.1097 。再運用能量守恒定律和相關(guān)數(shù)據(jù)，計算出速度v1（近月點的速度）=1750.78m/s,v2（遠月點的速度）=1669.77m/s，最后利用曲線的切線方程，代入點（近月點與遠月點）的坐標求值，計算出方向余弦即為相應(yīng)的速度方向。針對問題（2）關(guān)鍵詞

2、：模糊評判，聚類分析，流體交通量，排隊論，多元非線性回歸一、問題重述嫦娥三號于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（見附件1）。嫦娥三號在高速飛行的情況下，要保證準確地在月球預(yù)定區(qū)域內(nèi)實現(xiàn)

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

4、，可以知道嫦娥三號繞月球的軌道是由圓形軌道變?yōu)闄E圓形軌道，借助開普勒定律、能量守恒定律求解出近月點的速度。為了確定近月點和元月點的精確位置及相應(yīng)的速度方向，我們建立以赤道（月球的赤道）平面為xoy平面、月心為原點、月心與零度經(jīng)線和零度緯線交線的交點的連線為坐標軸的坐標系和赤道（月球的赤道）平面為xoy平面，為極軸（月球的極軸）為z軸建立空間直角坐標系，x軸與極坐標系的軸相重合。首先根據(jù)著陸點的經(jīng)度、緯度及月球的半徑求解出著陸點和近月點（帶參數(shù)a)的空間直角坐標。其次利用兩點間的距離公式，并借助MATLAB軟件求解出近月點與著陸點最短距離。從而計算出a（近月點的經(jīng)度）=。最后利用衛(wèi)星的軌跡是以月

5、心為其中一個焦點，以近月點與遠月點的距離為長軸的橢圓，從而求解出衛(wèi)星的軌跡方程，再運用隱函數(shù)求導(dǎo)的應(yīng)用的知識，求解出在近月點和遠月點的方向?qū)?shù)，進而求解近月點和遠月點方向余即為近月點和遠月點的速度的方向。2.2問題（2）的分析首先在根據(jù)題意，將嫦娥三號軟著陸問題，分為6個階段依次為主減速、快速調(diào)整、粗避障、精避障、緩慢下降、自由下降，我們先將6個階段分為4個階段，依次為第一階段（主減速和快速調(diào)整）、第二階段（粗避障）第三階段(精避障),第四階段（緩慢下降和自由下降）。其次在第一階段粗避障階段，嫦娥三號懸停在月球表面約2400米上方，對星下月表進行二維和三維成像，利用遺傳算法的思想，從圖像中先隨

6、機選取部分點，能直接從三維圖像中得知該點的海拔高度，再分別掃描這些點附近的地貌，找出一些地勢平坦的區(qū)域，我們用區(qū)域內(nèi)所有點與中心點海拔的均方差作為地勢判斷依據(jù)之一，保留這些坐標， 并進行重新組合，并改變某些坐標以便能獲得其他新區(qū)域的坐標，再次搜索地勢平坦的區(qū)域，重復(fù)進行多次搜索，直到?jīng)]有出現(xiàn)崎嶇地勢的時候， 我們將此時地勢最平坦的地方作為全局最優(yōu)降落地點三、模型假設(shè)1、不考慮空間飛行器上各點因燃料消耗而產(chǎn)生的位移；2、在對衛(wèi)星和空間飛行器進行軌道估計時，認為作用于其上的所有外力都通過其質(zhì)心；3、衛(wèi)星和空間飛行器的運動是在真空中進行的；4、衛(wèi)星只受重力影響，空間飛行器除自身推力外只受重力影響；5

7、、衛(wèi)星的觀測圖片及數(shù)據(jù)精準；6、四、變量與符號說明C0 一條車道的基本通行能力連續(xù)車流的車頭間距n 條車道的基本通行能力排隊長度車流量橫斷面通行能力系數(shù)車流量持續(xù)時間 L C y x1 x2 x3五、模型建立與求解5.1 問題(1)的分析、模型建立與求解5.1.1建模準備（1）開普勒定律開普勒第一定律開普勒第一定律開普勒第一定律，也稱橢圓定律：每一個行星都沿各自的橢圓軌道環(huán)繞太陽，而太陽則處在橢圓的一個焦點中。開普勒第二定律開普勒定律開普勒第二定律，也稱面積定律：在相等時間內(nèi)，太陽和運動著的行星的連線所掃過的面積都是相等的。 這一定律實際揭示了行星繞太陽公轉(zhuǎn)的角動量守恒。用公式表示為開普勒定律

8、開普勒第三定律開普勒定律開普勒第三定律，也稱調(diào)和定律：各個行星繞太陽公轉(zhuǎn)周期的平方和它們的橢圓軌道的半長軸的立方成正比。由這一定律不難導(dǎo)出：行星與太陽之間的引力與半徑的平方成反比。這是牛頓的萬有引力定a3律的一個重要基礎(chǔ)。用公式表示為2=K開普勒定律 T這里，是行星公轉(zhuǎn)軌道半長軸，是行星公轉(zhuǎn)周期，是常數(shù) 。（2）萬有引力萬有引力：任意兩個質(zhì)點有通過連心線方向上的力相互吸引。該引力大小與它們質(zhì)量的乘積成正比與它們距離的平方成反比，與兩物體的化學(xué)組成和其間介質(zhì)種類無關(guān)。即：M1M2， r2-11 其中M1，M2為兩物體的質(zhì)量，G=6.6710Nm.2kg.2（牛頓每平方米二次方千F=G克）5.1.

9、2 模型的建立根據(jù)以上的分析，建立以月球赤道平面為xOy平面,月心為原點O、Ox為月心與零度經(jīng)線和零度緯線交線的交點的連線，Oz為極軸（月球的極軸），Oy與Ox和Oz滿足右手標架，建立空間直角坐標系(如圖5-1所示）。圖5-1 衛(wèi)星繞月軌跡及軟著陸軌跡由于著陸點在球面上且近月點與遠月點是由月球的經(jīng)度、緯度及高度唯一確定，在此為了便于計算 將極坐標轉(zhuǎn)化為空間直角坐標，并代數(shù)題中相關(guān)數(shù)據(jù)，反解出經(jīng)度a。極坐標轉(zhuǎn)化為空間直角坐標x=rsinjcosq即：y=rsinjsinqz=rcosj （5.1.1）x=rsin(90-b）cos(-a)y=rsin(90-b）cos(-a) （5.1.2） z

10、=rcos(90-b）距離公式：d= （5.1.3） 其中：b為緯度；a為經(jīng)度；r為嫦娥三號距月心的距離；d為嫦娥三號距著陸點的距離；根據(jù)能量守恒、開普勒第二定律（面積定律），建立以下模型即： r1v1=r2v2 （5.1.4） 1122mv1+mgh=mv2+mgH22則近月點的速度，近月點的速度：v1= （5.1.5） v=2其中：m為衛(wèi)星的質(zhì)量，h1為海拔高度，h近月點距月球表面的距離； r1=h+r0+h1，r2=H+r0+h1，r0月球半徑，H遠月點距月球表面的距離， g月球重力加速度，v1 近月點的速度，v2 近月點的速度。5.1.3模型的求解5.1.3.1 近月點與遠月點的位置根

11、據(jù)題目所給數(shù)據(jù)以上分析，可知：b=0,h=15000m,r0=m,h1=-2641m將以上數(shù)據(jù)代入（5.1.1）式可得，著陸點及近月點的空間直角坐標分別為： x0=r0sin(90-b)cos(-a)=r0sin(90-19.51)cos(-44.12)y0=r0sin(90-b)sin(-a)=r0sin(90-19.51)sin(-44.12) （5.1.6） z=rcos(90-b)=r0cos(90-19.51)00x=rsin(90-b）cos(-a)=(r0+h)cosay=rsin(90-b）sin(-a)=-(r0+h)sina z=rcos(90-b）=0 （5.1.7）再將

12、（5.1.6）式和（5.1.7）式代入（5.1.3）式可得關(guān)于a與d（近月點和著陸點距離）的函數(shù)，？利用Mathematica 5.0編程求解可得：a=-139.1075.1.3.2近月點與遠月點的速度大小及方向近月點與遠月點的速度方向，即為相應(yīng)速度在x軸與y軸方向上的投影（如圖5-2所示）圖5-2 近月點與遠月點的速度方向示意圖由圖易知：5.2 模型二的建立5.2.1模型準備5.2.1.1系統(tǒng)模型1、著陸器的動力下降段一般從15km左右的軌道高度開始，下降到月球表面的時間比較短，在幾百秒范圍內(nèi)，所以可以不考慮月球引力攝動。月球自轉(zhuǎn)速度比較小，也可忽略。因此，可以利用二體模型描述系統(tǒng)的運動。建

13、立圖5-2所示的著陸坐標系，并假設(shè)著陸軌道在縱向平面內(nèi)，令月心為坐標原點，Oy指向動力下降段的開始制動點，Ox 指向著陸器的開始運動方向。則著陸器的質(zhì)心動力學(xué)方程可描述如下:r=vv=(F/m)siny-m/r2+rw2q=w w=-(F/m)cosy+2vw/r m=-F/ISP式中:r,q,w和m分別為著陸器的月心距、極角、角速度和質(zhì)量；v為著陸器沿r 方向上的速度；F為制動發(fā)動機的推力(固定的常值或0）；ISP為其比m為月球引力常數(shù)；y為發(fā)動機推力與當?shù)厮骄€的夾角即推力方向角。沖；圖5-3 月球軟著陸坐標系動力下降的初始條件由霍曼變軌后的橢圓軌道近月點確定，終端條件為著陸器在月面實現(xiàn)軟

14、著陸。令初始時刻t0=0，終端時刻tf不定，則相應(yīng)的初始條件為r0終端約束為rf=rL,vf=0,wf=0 =rL+h0,v0=0,w0=wo 式中:rL為月球半徑；h0為初始軌道高度；wo為軌道角速度。月球軟著陸的最優(yōu)軌道設(shè)計就是要在滿足上述初始條件和終端約束的前提下，調(diào)整推力大小和方向9使得著陸器實現(xiàn)燃料最優(yōu)軟著陸，即要求以下性能指標達最大。J=mdt 0tf5.2.1.2模型歸一化在軌道優(yōu)化過程中，由于各狀態(tài)變量的量級相差較大，尋優(yōu)過程中可能會導(dǎo)致有效位數(shù)的丟失。通過歸一化處理可以克服這一缺點9，提高。計算精度。令rref=r0,mtef=m0，則=r/rref,=v/vref,vref

15、=ISp=I72=F/Fref,Fref=mrefvref/rref,=m/mref,=t/tref,=rref/vref,=q。那么，著陸器的動力學(xué)方程可改為：=v22=(F/m)siny-m/+ =w=-(F/)cosy+2/=-F/ISP相應(yīng)的初始條件和終端約束變?yōu)椋?1,=0,=w 000/mf=r1/r0,vf=0,wf=0性能指標改寫為：=0第4期朱建豐等:基于自適應(yīng)模擬退火遺傳算法的月球軟著陸軌道優(yōu)化道優(yōu)化問題轉(zhuǎn)化為多參數(shù)優(yōu)化問題,再利用SQP方法求解。雖然避開了沒有明確物理意義的參數(shù)猜測,但是SQP的本質(zhì)仍然會使該方法遇到病態(tài)梯度、初始點敏感和局部收斂問題。曾國強6和徐敏7分別

16、用二進制和浮點數(shù)GA對著陸軌道進行了優(yōu)化,避免了初值猜測,得到的結(jié)果也比較滿意。但是,鑒于GA局部搜索能力較差的缺點,會使得GA的優(yōu)化精度不夠或優(yōu)化效率不高。相對而言,國外對月球軟著陸軌道的優(yōu)化問題研究比較少。GA最早是由Holland教授提出的8,它是一種隨機優(yōu)化方法,具有不依賴問題模型、適用面廣和魯棒性強的優(yōu)點,并已應(yīng)用在航天器的軌道優(yōu)化設(shè)計中1。然而,GA在實際應(yīng)用中存在收斂速度慢和早熟等問題,不具備“爬山”的能力。模擬退火算法(SAA)最早是由Kirkpatrick等提出的,它是一種啟發(fā)式隨機搜索算法,具有很強的局部搜索能力和“爬山”能力,但是SAA產(chǎn)生的新解不及GA豐富,對全局的了解

17、甚少,尋優(yōu)過程很慢。因此,可以將GA和SAA的優(yōu)點結(jié)合起來,揚長避短,構(gòu)成高效、魯棒的新算法。本文將GA和SAA有機地結(jié)合,形成自適應(yīng)模擬退火遺傳算法(ASAGA),并將其應(yīng)用到月球軟著陸的最優(yōu)軌道設(shè)計中。1系統(tǒng)模型著陸器的動力下降段一般從15 km左右的軌道高度開始,下降到月球表面的時間比較短,在幾百秒范圍內(nèi),所以可以不考慮月球引力攝動。月球自轉(zhuǎn)速度比較小,也可忽略。因此,可以利用二體模型描述系統(tǒng)的運動。建立圖1所示的著陸坐標系,并假設(shè)著陸軌道在縱向平面內(nèi),令月心O為坐標原點,Oy指向動力下降段的開始制動點,Ox指向著陸器的開始運動方向。則著陸器的質(zhì)心動力學(xué)方程可描述如下:r= vv=(F

18、/m)sin - /r2+ r 2= = -(F /m)cos+ 2v /rm= -F /ISP(1)式中:r,和m分別為著陸器的月心距、極角、角速度和質(zhì)量;v為著陸器沿r方向上的速度;F為制動發(fā)動機的推力(固定的常值或0);ISP為其比沖;為月球引力常數(shù);為發(fā)動機推力與當?shù)厮骄€的夾角即推力方向角。圖1月球軟著陸極坐標系Fig. 1Polar coordinate system of lunar soft landing動力下降的初始條件由霍曼變軌后的橢圓軌道近月點確定,終端條件為著陸器在月面實現(xiàn)軟著陸。令初始時刻t0= 0,終端時刻tf不定,則相應(yīng)的初始條件為r0= rL+ h0,v0=

19、0,0= o(2)終端約束為rf= rL,vf= 0,f= 0 (3)式中:rL為月球半徑;h0為初始軌道高度;o為軌道角速度。月球軟著陸的最優(yōu)軌道設(shè)計就是要在滿足上述初始條件和終端約束的前提下,調(diào)整推力大小和方向,使得著陸器實現(xiàn)燃料最優(yōu)軟著陸,即要求以下性能指標達最大。J=tf0mdt (4)2歸一化在軌道優(yōu)化過程中,由于各狀態(tài)變量的量級相差較大,尋優(yōu)過程中可能會導(dǎo)致有效位數(shù)的丟失。通過歸一化處理可以克服這一缺點9,提高計算精度。令rref= r0,mref= m0,則r= r /rref, v=v /vref,vref= /rref, ISP= ISPrref/, F= F /Fref,F

20、ref= mrefv2ref/rref, m= m /mref, = r3ref/,t= t /tref,tref= rref/vref,=。那么,著陸器的動力學(xué)方程可改寫為r= vv=( F / m)sin -1 /r2+r 2= = -( F / m)cos+ 2 v /rm= - F / ISP(5)相應(yīng)的初始條件和終端約束變航空學(xué)報第28卷4. 3算法的實現(xiàn)將提出的ASAGA應(yīng)用到月球軟著陸的軌道優(yōu)化中,具體的實現(xiàn)步驟如下:(1)設(shè)置初始參數(shù),包括種群規(guī)模M,最大遺傳代數(shù)Tmax,退火初始溫度T0,溫度下降系數(shù),最小新解接受次數(shù)Nmin,最大內(nèi)循環(huán)次數(shù)Cmax,軌道離散參數(shù)n。隨機產(chǎn)生

21、初始種群。(2)計算種群中各個個體的適應(yīng)度值,記錄最優(yōu)個體,并采用如下方法進行適應(yīng)度拉伸f= exp -(fmax- f)/Ti (22)式中f為拉伸后的適應(yīng)度值。采用輪盤選擇策略進行個體選擇,并將選擇的個體隨機兩兩配對,按照式(12)、式(13)、式(15)和式(16)進行交叉操作,再利用式(14)、式(17)和式(18)對個體的每個參數(shù)進行變異操作。(3)令新解接受次數(shù)na= 0,內(nèi)循環(huán)次數(shù)nc=0。對遺傳操作后的子代個體計算適應(yīng)度值,選擇前no個優(yōu)秀個體,并分別進行如下的模擬退火操作:計算(Ti),分別對個體的每個參數(shù)按式(20)進行解的變換,并按照Meteopolis準則來判斷是否接受

22、新解為新的當前解。如果接受,則na= na+ 1;反之,na不變。如果na Nmin,且nc Cmax,則nc= nc+ 1,并返回 ;否則,跳出循環(huán),并將群體中適應(yīng)度最差的no個個體替換成退火操作后的個體,進入步驟)。(4)刪除子代種群中的任意一個個體,并替換成步驟(2)記錄的最優(yōu)個體。(5)如果當前遺傳代數(shù)T Tmax,則按式(19)進行降溫,T= T+ 1,并返回步驟(2);否則結(jié)束整個優(yōu)化過程?？梢钥吹?在上述的ASAGA中,進行了適應(yīng)度拉伸。這樣處理的優(yōu)點是:在溫度高時(算法初期),適應(yīng)度相近的個體產(chǎn)生的后代概率相近,避免了早期個別好的個體的后代充斥整個種群,造成早熟;而當溫度不斷下

23、降后,拉伸作用加強,使適應(yīng)度相近的個體適應(yīng)度差異放大,從而使得優(yōu)秀個體更加突出。同時,算法采用了最優(yōu)個體保留策略,使得最優(yōu)個體不會被破壞,保證了算法的收斂性。只選擇前no個優(yōu)秀個體進行模擬退火操作,可以大大縮短優(yōu)化時間,也能夠使優(yōu)秀個體向最優(yōu)的方向進化。5仿真實例設(shè)定月球軟著陸的初始條件為r0= 1 753km,v0= 0,0= 9. 65 10 -4 rad/s,0= 0,m0=600 kg。終端約束為:rf= 1 738 km,vf= 0,f= 0。月球引力常數(shù)= 4 902. 75km3 /s2,制動發(fā)動機推力F= 3 450 N,比沖ISP= 300 9. 8 m /s。令軌道離散化參

24、數(shù)n= 9,即待優(yōu)化參數(shù)為10個推力方向角和1個終端時刻tf,共11個參數(shù)。由于尋優(yōu)邊界L和R決定了搜索空間的大小,而搜索空間小,則算法收斂速度快,反之則慢。因此需要估計L和R的值,加快收斂速度。根據(jù)齊奧爾科夫斯基公式和軟著陸初始條件,終端時刻可由下式估計tf=1 -exp(Vf-V0)/ISP(ISPm0/F)(23)式中:Vf和V0分別為著陸器的終端速度大小和初始速度大小。同時由經(jīng)驗可知,推力方向角是逐漸增大的,所以尋優(yōu)邊界可以確定為L= 0000000000500(24)R= 1 2 3 4 5 6 7 8 9 10 9 700(25)ASAGA的其他參數(shù)設(shè)定為:M= 20,Tmax=3

25、00,T0= 10 000,= 0. 97,Nmin= 10,Cmax= 10,no= 2,(Ti)= 0. 5 1 -(T -1)/Tmax,1=52= 5 max(1,10(2T/Tmax -1),Rmax= 10。由于GA與SAA都是隨機優(yōu)化算法,所以ASAGA也具有隨機性。因此為了驗證該算法的有效性,文中對以上的初值條件,進行10次數(shù)值仿真,結(jié)果全部收斂,說明該算法有很好的收斂性。10次結(jié)果的均值為:終端速度大小Vf= 0. 11 m/s,終端月心距rf= 1 737 999. 92 m,燃料消耗277. 677 5kg。其中最好的結(jié)果如圖3圖6所示。圖3推力方向角的時間歷程Fig.

26、3Time histories of direction angle of thrust朱建豐等:基于自適應(yīng)模擬退火遺傳算法的月球軟著陸軌道優(yōu)化式中:A,B和A,B分別為父代和子代的個體;,為(0,r)區(qū)間上的均勻分布隨機數(shù),r 1為交叉系數(shù),其大小可以控制交叉操作的變化范圍;L,R分別為尋優(yōu)參數(shù)的左右邊界。如果交叉后的子代超出尋優(yōu)邊界,則重復(fù)進行交叉操作,假如重復(fù)操作達到最大設(shè)定次數(shù)Rmax,子代仍然不能滿足邊界約束,則利用式(13)進行子代邊界修正。變異操作11采用如下形式C=C+ k(R -C)U(0,1)= 0C -k(C -L)U(0,1)= 1(14)式中:C和C分別為父代和子代的

27、個體;為(0,1)區(qū)間上的均勻分布隨機數(shù);k(0,1為變異系數(shù);U(0,1)為隨機產(chǎn)生的整數(shù)0或1。GA的交叉概率Pc與變異概率Pm對其性能影響很大。Pc越大,新個體產(chǎn)生的速度越快。然而,Pc過大會使得具有高適應(yīng)度的個體結(jié)構(gòu)很快被破壞;Pc過小會使搜索過程緩慢,以致停滯不前。如果變異概率Pm過小,就不易產(chǎn)生新的個體結(jié)構(gòu);Pm過大,GA就變成了純粹的隨機搜索算法。因此,必須采用自適應(yīng)的方法,讓交叉概率和變異概率隨著適應(yīng)度的變化而改變。同時,注意到交叉系數(shù)r與變異系數(shù)k的大小也會影響GA的性能,所以也采用自適應(yīng)策略。本文中的Pc,Pm,r和k按照如下公式進行自適應(yīng)調(diào)整:Pc=Pc1-Pc2(fc-

28、 favg)/(fmax-favg) fc favgPc1fc favg(15)r=r1+ r2(fmax- fc)/(fmax- favg) fc favgr1+ r2fc favg(16)Pm=Pm1-Pm2(fm-favg)/(fmax-favg) fm favgPm1fm favg(17)k=k1+ k2(fmax- fm)/(fmax- favg) fm favgk1+ k2fm favg(18)式中:fmax為群體中最大的適應(yīng)度值; favg為群體的平均適應(yīng)度值;fc為要交叉的兩個個體中較大的適應(yīng)度值; fm為要變異個體的適應(yīng)度值。Pc1= 0. 9,Pc2= 0. 2,r1= 0

29、. 5 -0. 25T/Tmax(T為當前遺傳代數(shù),Tmax為最大遺傳代數(shù)),r2= 0. 5,Pm1= 0. 01,Pm2= 0. 005,k1= 0. 5 -0. 4T/Tmax,k2= 0. 5。采取如上的自適應(yīng)策略后,使得適應(yīng)度高于群體平均適應(yīng)度的個體,對應(yīng)較小的r與k和較低的Pc與Pm。前者使該個體在原值附近小范圍內(nèi)進行交叉和變異,加快GA的收斂速度;后者使該個體得以保護進入下一代;而低于平均適應(yīng)度的個體,對應(yīng)較大的r與k和較高的Pc和Pm。前者使該個體進行交叉和變異的范圍變大,加強GA的全局搜索能力;后者使該個體很快被淘汰掉。同時,隨著優(yōu)化的進行,r和k的值逐漸減小,GA的局部搜索

30、能力也逐漸加強。自適應(yīng)策略使得GA在保持群體多樣性的同時,保證了算法的收斂性,也在一定程度上提高了局部搜索能力。4. 2模擬退火算法SAA是基于金屬退火的機理而建立起的一種隨機算法,它能夠通過控制溫度的變化過程來實現(xiàn)大范圍的粗略搜索與局部的精細搜索。本文采用指數(shù)降溫策略對溫度的變化進行控制,即Ti= T0()T-1 (19)式中:Ti為當前控制溫度;T0為初始設(shè)定溫度;為溫度下降系數(shù)。在SAA進行時,解的變換即新解的產(chǎn)生,是發(fā)生在當前解的鄰域內(nèi),采用如下公式進行解的變換:Y=Y+(R -Y)(Ti)U(0,1)= 0Y -(Y -L)(Ti)U(0,1)= 1(20)式中:Y和Y分別為當前解和

31、新解,為種群的個體;(Ti)為隨Ti減小而減小的擾動值,當Ti為T0時,(T0) 1,保證了最大擾動值不會使新解Y越界,當Ti0時,(Ti)0,滿足算法收斂的條件;為區(qū)間(0,1)上的均勻分布隨機數(shù)。狀態(tài)由Y變?yōu)閅的接受概率可由下面的Meteopolis準則來確定p=1 f(Y) f(Y)exp -(f(Y)-f(Y)/Ti f(Y) f(Y)(21)使用上述準則的優(yōu)點是:當新解更優(yōu)時,完全接受新解為新的當前解;而當新解為惡化解時,以概率p接受惡化解為新的當前解。這便使得SAA能夠避免陷入局部最優(yōu)。隨著優(yōu)化的進行,SAA的局部搜索能力也逐漸增強,確保算法有足夠的搜索精度。r0= 1, v0=

32、0, 0= or30/ (6)rf= rL/r0,vf= 0,f= 0 (7)性能指標改寫為J=tf0 m dt (8)3參數(shù)化方法因為月球軟著陸的軌道優(yōu)化搜索空間是一個泛函空間,不能直接用優(yōu)化算法處理,必須將控制變量參數(shù)化。參數(shù)化方法是軌道優(yōu)化與優(yōu)化算法的橋梁,其精度直接影響到最優(yōu)軌道的好壞程度。目前,主要的參數(shù)化方法有3種10:直接離散法、多段參數(shù)插值法和函數(shù)逼近法。直接離散法的精度最高,但計算最慢。多段參數(shù)插值法的計算最快,但精度最差。所以這兩種方法都不宜采用。函數(shù)逼近法的精度適中,計算速度比較快,但是多項式系數(shù)沒有明確的物理意義,很難進行初值猜測。因此,本文將函數(shù)逼近法進行改進。首先將

33、軌道離散化成許多小段,在各小段的節(jié)點設(shè)定待優(yōu)化的參數(shù),然后利用這些參數(shù)進行多項式擬合,從而得到整個軌道的控制曲線。將月球軟著陸軌道離散化,分割成n個小段,每段的節(jié)點設(shè)定一個推力方向角,如圖2所示。那么,可以令n+ 1個節(jié)點的推力方向角和終端時刻tf作為待優(yōu)化的參數(shù)。每個節(jié)點的時刻可以由下式得到ti= t0+ i(tf-t0)/n(i= 0,1, ,n)(9)這樣,就使得每個節(jié)點的推力方向角都有一個對應(yīng)的節(jié)點時刻。如果假設(shè)月球軟著陸的推力方向角可以表示成一個多項式,即(t)= 0+1t+2t2+3t3 (10)那么,可以將節(jié)點的推力方向角與對應(yīng)的節(jié)點時刻對上述多項式進行擬合,求得多項式的系數(shù)i(

34、i= 0,1,2,3),進而得到整個著陸軌道的推力方向角。圖2軌道離散化Fig. 2Trajectory discretization可以看到,采取如上改進的函數(shù)逼近法進行參數(shù)化,所選定的待優(yōu)化參數(shù)都具有明確的物理意義,從而避開對沒有物理含義的待優(yōu)化參數(shù)的初始猜測。同時,該參數(shù)化方法的本質(zhì)是函數(shù)逼近,因此它與函數(shù)逼近法的精度相當,計算速度也比較快。4自適應(yīng)模擬退火遺傳算法GA的局部搜索能力較差,容易陷入局部最優(yōu)解,但把握搜索過程總體的能力較強;而SAA具有較強的局部搜索能力,并能使搜索過程避免陷入局部最優(yōu)解,但SAA卻對整個搜索空間的狀況了解不多,不便于使搜索過程進入最有希望的搜索區(qū)域,從而使

35、得SAA的運算效率不高。ASAGA就是將SAA嵌入到GA的循環(huán)體中,互相取長補短,形成一種新的優(yōu)化算法。4. 1自適應(yīng)遺傳算法GA是模擬生物在自然環(huán)境中的遺傳和進化過程而形成的一種自適應(yīng)全局優(yōu)化隨機搜索算法。由于它對搜索空間不作任何假設(shè),它既不要求搜索空間是光滑的,也不要求它是處處可微的,因而它能解決很大一類問題,在飛行器軌道優(yōu)化中逐漸開始得到應(yīng)用,但是要將GA成功的應(yīng)用到軌道優(yōu)化中,還必須處理好如何將終端約束反映到GA的適應(yīng)度函數(shù)中?？紤]到軌道優(yōu)化的精度,采用動態(tài)罰函數(shù)法10對終端約束進行處理,則適應(yīng)度函數(shù)可表示成如下形式f(X)= -1(v(tf)-vf)2+(tf)r(tf)-frf)2

36、 1 /2 -2| r(tf)-rf|+ m0+ J(11)式中:X為種群的個體;v(tf),(tf)和r(tf)為終端時刻的狀態(tài)值;1,2為懲罰因子,隨著優(yōu)化的進行它們逐漸變大,使得對約束的限制也逐漸加強。在遺傳操作中,交叉操作是產(chǎn)生新個體的主要方法,它決定著GA的全局搜索能力;而變異操作只是產(chǎn)生新個體的輔助方法,它決定了GA的局部搜索能力。交叉與變異相互配合,共同完成對搜索空間的全局搜索和局部搜索。本文的交叉操作11采用如下方式:A=(1 -)A+BB=(1 -)B+A(12)如果A(B) R那么A(B)= R(13)第4期朱建豐等:基于自適應(yīng)模擬退火遺傳算法的月球軟著陸軌道優(yōu)化圖4徑向速

37、度的時間歷程Fig. 4Time histories of radial velocity圖5月心距的時間歷程Fig. 5Time histories of distance from lunar center to lander圖6角速度的時間歷程Fig. 6Time history of angle velocity從結(jié)果圖可以看出:雖然由ASAGA得到的最佳推力方向角與由Pontryagin極大值原理得到的最優(yōu)推力方向角存在較大的差別,但是這并不影響最后的著陸效果。著陸器的終端狀態(tài)為:Vf= 0. 06 m/s,rf= 1 738 000. 02 m。燃料消耗為:277. 622 4 k

38、g,比由Pontryagin極大值原理得到的最優(yōu)燃料消耗(277. 576 5 kg)僅多0. 045 9kg。以上數(shù)據(jù)說明了該算法的優(yōu)化精度很高。下面給出群體中最大適應(yīng)度值隨時間的變化曲線,如圖7所示。圖7說明:在優(yōu)化進程中,多次出現(xiàn)局部“早熟”現(xiàn)象,但都被ASAGA的搜索能力所克服,成功地跳出了局部最優(yōu)的誘惑。圖7最大適應(yīng)度的時間歷程Fig. 7Time history of maximum fitness8.參考文獻1陳剛,萬自明,徐敏,等.遺傳算法在航天器軌道優(yōu)化中的應(yīng)用J.彈道學(xué)報, 2006,18(1):1-5.Chen G, Wan Z M, Xu M, et al. Overv

39、iew of spacecrafttrajectory optimization using genetic algorithm J. Jour-nal of Ballistics, 2006,18(1):1-5.(in Chinese)2王大軼,李鐵壽,馬興瑞.月球最優(yōu)軟著陸兩點邊值問題的數(shù)值解法J.航天控制, 2000(3):44-49.Wang D Y, Li T S, Ma X R. Numerical solution of TPB-VP in optimal lunar soft landing J. Aerospace Control,2000(3):44-49.(in Chin

40、ese)3王劼,崔乃剛,劉暾,等.定常推力登月飛行器最優(yōu)軟著陸軌道研究J.高技術(shù)通訊, 2003,13(4):39-42.Wang J, Cui N G, Liu T, et al. Study on soft-landing trajec-tories of constant-thrust-amplitude lunar probe J. HighTechnology Letters, 2003,13(4):39-42.(in Chinese)4王劼,李俊峰,崔乃剛,等.登月飛行器軟著陸軌道的遺傳算法優(yōu)化 J.清華大學(xué)學(xué)報(自然科學(xué)版), 2003,43(8):1056-1059.Wang

41、J, Li J F,Cui N G, et al. Genetic algorithm optimi-zation of lunar probe soft-landing trajectories J. Journalof Tsinghua University(Sci and Tech), 2003,43(8):1056-1059.(in Chinese)9.附件附件1/求 m 階B 樣條：function y=N(m,x)s=factorial(m-1);sum=0;for k=0:mif(x-k)0 sum=sum+(-1)k*(nchoosek(n,k)*(x-k)(m-1); end

42、 y=s(-1)*sum/畫 m 階 B 樣條：function null=l(m)x=-1:0.01:m+1;a=size(x);new=;for k=1:size(x,2) new(k)=psi(m,x(k); end plot(x,new)/對 m 階 B 樣條求 n 階導(dǎo)數(shù)：function z=derive(m,n,x)sum=0; for k=0:m sum=sum+(-1)k*nchoosek(n,k)*N(m-n,x-k); end z=sum;/求 m 階 B 樣條小波:function g=psi(m,x)if m=1 if x0&x0.5&x1 g=-1;else g=0

43、; endelse s=2(-m+1);sum=0;for j=0:(2*m-2) sum=sum+(-1)j*N(2*m,j+1)*derive(2*m,m,2*x-j); end g=s*sum; end/畫 m 階樣條小波的圖像：function null=h(m)x=-1:0.01:2*m; a=size(x);new=;for k=1:size(x,2) new(k)=psi(m,x(k); end plot(x,new)附件2clc;clear;%繪照片圖Fujian_3 = imread(附件3 距2400m處的數(shù)字高程圖.tif);Fujian_4 = imread(附件4 距

44、月面100m處的數(shù)字高程圖.tif);x3=1:2300;y3=1:2300;Fujian_31=double(Fujian_3);Fujian_41=double(Fujian_4);xx3 yy3=meshgrid(x3,y3);x4=0.1:0.1:100;y4=0.1:0.1:100;xx4 yy4=meshgrid(x4,y4);figure(1): mesh(xx3,yy3,Fujian_31);figure(2): mesh(xx4,yy4,Fujian_41);附件3function P = mutation(P,Nb,Pm)%變異Nbb = length(Nb);for n

45、= 1:size(P,1) b2 = 0; for m = 1:Nbb if rand Pm b1 = b2 + 1; bi = b1 + mod(floor(rand*Nb(m),Nb(m); b2 = b2 + Nb(m); P(n,bi) = P(n,bi); end endendfunction xo,fo = Opt_Simu(f,x0,l,u,kmax,q,TolFun)% 模擬退火算法求函數(shù) f(x)的最小值點， 且 l = x = u% f為待求函數(shù)，x0為初值點，l，u分別為搜索區(qū)間的上下限，kmax為最大迭代次數(shù)% q為退火因子，TolFun為函數(shù)容許誤差%算法第一步根據(jù)輸

46、入變量數(shù)，將某些量設(shè)為缺省值if nargin 7 TolFun = 1e-8;endif nargin 6 q = 1;endif nargin 5 kmax = 100;end%算法第二步，求解一些基本變量N = length(x0); %自變量維數(shù)x = x0;fx = feval(f,x); %函數(shù)在初始點x0處的函數(shù)值xo = x;fo = fx;%算法第三步，進行迭代計算，找出近似全局最小點for k =0:kmax Ti = (k/kmax)q; mu = 10(Ti*100); % 計算mu dx = Mu_Inv(2*rand(size(x)-1,mu).*(u - l);%步長dx x1 = x + dx; %下一個估計點 x1 = (x1 l).*l +(l = x1).*(x1 = u).*x1 +(u x1).*u; %將x1限定在區(qū)間l,u上 fx1 = feval(f,x1); df = fx1- fx; if df 0|rand exp(-Ti*df/(abs(fx) +