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

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

2、：模糊評(píng)判，聚類分析，流體交通量，排隊(duì)論，多元非線性回歸一、問(wèn)題重述嫦娥三號(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)通過(guò)多個(gè)發(fā)動(dòng)機(jī)的脈沖組合實(shí)現(xiàn)各種姿態(tài)的調(diào)整控制。嫦娥三號(hào)的預(yù)定著陸點(diǎn)為19.51W，44.12N，海拔為-2641m（見(jiàn)附件1）。嫦娥三號(hào)在高速飛行的情況下，要保證準(zhǔn)確地在月球預(yù)定區(qū)域內(nèi)實(shí)現(xiàn)

3、軟著陸，關(guān)鍵問(wèn)題是著陸軌道與控制策略的設(shè)計(jì)。其著陸軌道設(shè)計(jì)的基本要求：著陸準(zhǔn)備軌道為近月點(diǎn)15km，遠(yuǎn)月點(diǎn)100km的橢圓形軌道；著陸軌道為從近月點(diǎn)至著陸點(diǎn)，其軟著陸過(guò)程共分為6個(gè)階段（見(jiàn)附2），要求滿足每個(gè)階段在關(guān)鍵點(diǎn)所處的狀態(tài)；盡量減少軟著陸過(guò)程的燃料消耗。根據(jù)上述的基本要求，請(qǐng)你們建立數(shù)學(xué)模型解決下面的問(wèn)題：（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)的誤差分析和敏感性分析。二、問(wèn)題分析2.1問(wèn)題（1）的分析首先根據(jù)問(wèn)題的假設(shè)、題目中所提供的數(shù)據(jù)及圖片分析

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

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

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

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

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

9、2 模型的建立根據(jù)以上的分析，建立以月球赤道平面為xOy平面,月心為原點(diǎn)O、Ox為月心與零度經(jīng)線和零度緯線交線的交點(diǎn)的連線，Oz為極軸（月球的極軸），Oy與Ox和Oz滿足右手標(biāo)架，建立空間直角坐標(biāo)系(如圖5-1所示）。圖5-1 衛(wèi)星繞月軌跡及軟著陸軌跡由于著陸點(diǎn)在球面上且近月點(diǎn)與遠(yuǎn)月點(diǎn)是由月球的經(jīng)度、緯度及高度唯一確定，在此為了便于計(jì)算 將極坐標(biāo)轉(zhuǎn)化為空間直角坐標(biāo)，并代數(shù)題中相關(guān)數(shù)據(jù)，反解出經(jīng)度a。極坐標(biāo)轉(zhuǎn)化為空間直角坐標(biāo)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為嫦娥三號(hào)距月心的距離；d為嫦娥三號(hào)距著陸點(diǎn)的距離；根據(jù)能量守恒、開(kāi)普勒第二定律（面積定律），建立以下模型即： r1v1=r2v2 （5.1.4） 1122mv1+mgh=mv2+mgH22則近月點(diǎn)的速度，近月點(diǎn)的速度：v1= （5.1.5） v=2其中：m為衛(wèi)星的質(zhì)量，h1為海拔高度，h近月點(diǎn)距月球表面的距離； r1=h+r0+h1，r2=H+r0+h1，r0月球半徑，H遠(yuǎn)月點(diǎn)距月球表面的距離， g月球重力加速度，v1 近月點(diǎn)的速度，v2 近月點(diǎn)的速度。5.1.3模型的求解5.1.3.1 近月點(diǎn)與遠(yuǎn)月點(diǎn)的位置根

11、據(jù)題目所給數(shù)據(jù)以上分析，可知：b=0,h=15000m,r0=m,h1=-2641m將以上數(shù)據(jù)代入（5.1.1）式可得，著陸點(diǎn)及近月點(diǎn)的空間直角坐標(biāo)分別為： 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（近月點(diǎn)和著陸點(diǎn)距離）的函數(shù)，？利用Mathematica 5.0編程求解可得：a=-139.1075.1.3.2近月點(diǎn)與遠(yuǎn)月點(diǎn)的速度大小及方向近月點(diǎn)與遠(yuǎn)月點(diǎn)的速度方向，即為相應(yīng)速度在x軸與y軸方向上的投影（如圖5-2所示）圖5-2 近月點(diǎn)與遠(yuǎn)月點(diǎn)的速度方向示意圖由圖易知：5.2 模型二的建立5.2.1模型準(zhǔn)備5.2.1.1系統(tǒng)模型1、著陸器的動(dòng)力下降段一般從15km左右的軌道高度開(kāi)始，下降到月球表面的時(shí)間比較短，在幾百秒范圍內(nèi)，所以可以不考慮月球引力攝動(dòng)。月球自轉(zhuǎn)速度比較小，也可忽略。因此，可以利用二體模型描述系統(tǒng)的運(yùn)動(dòng)。建

13、立圖5-2所示的著陸坐標(biāo)系，并假設(shè)著陸軌道在縱向平面內(nèi)，令月心為坐標(biāo)原點(diǎn)，Oy指向動(dòng)力下降段的開(kāi)始制動(dòng)點(diǎn)，Ox 指向著陸器的開(kāi)始運(yùn)動(dòng)方向。則著陸器的質(zhì)心動(dòng)力學(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為制動(dòng)發(fā)動(dòng)機(jī)的推力(固定的常值或0）；ISP為其比m為月球引力常數(shù)；y為發(fā)動(dòng)機(jī)推力與當(dāng)?shù)厮骄€的夾角即推力方向角。沖；圖5-3 月球軟著陸坐標(biāo)系動(dòng)力下降的初始條件由霍曼變軌后的橢圓軌道近月點(diǎn)確定，終端條件為著陸器在月面實(shí)現(xiàn)軟

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

15、=ISp=I72=F/Fref,Fref=mrefvref/rref,=m/mref,=t/tref,=rref/vref,=q。那么，著陸器的動(dòng)力學(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性能指標(biāo)改寫為：=0第4期朱建豐等:基于自適應(yīng)模擬退火遺傳算法的月球軟著陸軌道優(yōu)化道優(yōu)化問(wèn)題轉(zhuǎn)化為多參數(shù)優(yōu)化問(wèn)題,再利用SQP方法求解。雖然避開(kāi)了沒(méi)有明確物理意義的參數(shù)猜測(cè),但是SQP的本質(zhì)仍然會(huì)使該方法遇到病態(tài)梯度、初始點(diǎn)敏感和局部收斂問(wèn)題。曾國(guó)強(qiáng)6和徐敏7分別

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

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

18、/m)sin - /r2+ r 2= = -(F /m)cos+ 2v /rm= -F /ISP(1)式中:r,和m分別為著陸器的月心距、極角、角速度和質(zhì)量;v為著陸器沿r方向上的速度;F為制動(dòng)發(fā)動(dòng)機(jī)的推力(固定的常值或0);ISP為其比沖;為月球引力常數(shù);為發(fā)動(dòng)機(jī)推力與當(dāng)?shù)厮骄€的夾角即推力方向角。圖1月球軟著陸極坐標(biāo)系Fig. 1Polar coordinate system of lunar soft landing動(dòng)力下降的初始條件由霍曼變軌后的橢圓軌道近月點(diǎn)確定,終端條件為著陸器在月面實(shí)現(xiàn)軟著陸。令初始時(shí)刻t0= 0,終端時(shí)刻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è)計(jì)就是要在滿足上述初始條件和終端約束的前提下,調(diào)整推力大小和方向,使得著陸器實(shí)現(xiàn)燃料最優(yōu)軟著陸,即要求以下性能指標(biāo)達(dá)最大。J=tf0mdt (4)2歸一化在軌道優(yōu)化過(guò)程中,由于各狀態(tài)變量的量級(jí)相差較大,尋優(yōu)過(guò)程中可能會(huì)導(dǎo)致有效位數(shù)的丟失。通過(guò)歸一化處理可以克服這一缺點(diǎn)9,提高計(jì)算精度。令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,=。那么,著陸器的動(dòng)力學(xué)方程可改寫為r= vv=( F / m)sin -1 /r2+r 2= = -( F / m)cos+ 2 v /rm= - F / ISP(5)相應(yīng)的初始條件和終端約束變航空學(xué)報(bào)第28卷4. 3算法的實(shí)現(xiàn)將提出的ASAGA應(yīng)用到月球軟著陸的軌道優(yōu)化中,具體的實(shí)現(xiàn)步驟如下:(1)設(shè)置初始參數(shù),包括種群規(guī)模M,最大遺傳代數(shù)Tmax,退火初始溫度T0,溫度下降系數(shù),最小新解接受次數(shù)Nmin,最大內(nèi)循環(huán)次數(shù)Cmax,軌道離散參數(shù)n。隨機(jī)產(chǎn)生

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

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

23、降后,拉伸作用加強(qiáng),使適應(yīng)度相近的個(gè)體適應(yīng)度差異放大,從而使得優(yōu)秀個(gè)體更加突出。同時(shí),算法采用了最優(yōu)個(gè)體保留策略,使得最優(yōu)個(gè)體不會(huì)被破壞,保證了算法的收斂性。只選擇前no個(gè)優(yōu)秀個(gè)體進(jìn)行模擬退火操作,可以大大縮短優(yōu)化時(shí)間,也能夠使優(yōu)秀個(gè)體向最優(yōu)的方向進(jìn)化。5仿真實(shí)例設(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,制動(dòng)發(fā)動(dòng)機(jī)推力F= 3 450 N,比沖ISP= 300 9. 8 m /s。令軌道離散化參

24、數(shù)n= 9,即待優(yōu)化參數(shù)為10個(gè)推力方向角和1個(gè)終端時(shí)刻tf,共11個(gè)參數(shù)。由于尋優(yōu)邊界L和R決定了搜索空間的大小,而搜索空間小,則算法收斂速度快,反之則慢。因此需要估計(jì)L和R的值,加快收斂速度。根據(jù)齊奧爾科夫斯基公式和軟著陸初始條件,終端時(shí)刻可由下式估計(jì)tf=1 -exp(Vf-V0)/ISP(ISPm0/F)(23)式中:Vf和V0分別為著陸器的終端速度大小和初始速度大小。同時(shí)由經(jīng)驗(yàn)可知,推力方向角是逐漸增大的,所以尋優(yōu)邊界可以確定為L(zhǎng)= 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都是隨機(jī)優(yōu)化算法,所以ASAGA也具有隨機(jī)性。因此為了驗(yàn)證該算法的有效性,文中對(duì)以上的初值條件,進(jìn)行10次數(shù)值仿真,結(jié)果全部收斂,說(shuō)明該算法有很好的收斂性。10次結(jié)果的均值為:終端速度大小Vf= 0. 11 m/s,終端月心距rf= 1 737 999. 92 m,燃料消耗277. 677 5kg。其中最好的結(jié)果如圖3圖6所示。圖3推力方向角的時(shí)間歷程Fig.

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

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

29、. 5 -0. 25T/Tmax(T為當(dāng)前遺傳代數(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)度的個(gè)體,對(duì)應(yīng)較小的r與k和較低的Pc與Pm。前者使該個(gè)體在原值附近小范圍內(nèi)進(jìn)行交叉和變異,加快GA的收斂速度;后者使該個(gè)體得以保護(hù)進(jìn)入下一代;而低于平均適應(yīng)度的個(gè)體,對(duì)應(yīng)較大的r與k和較高的Pc和Pm。前者使該個(gè)體進(jìn)行交叉和變異的范圍變大,加強(qiáng)GA的全局搜索能力;后者使該個(gè)體很快被淘汰掉。同時(shí),隨著優(yōu)化的進(jìn)行,r和k的值逐漸減小,GA的局部搜索

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

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

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

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

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

35、得SAA的運(yùn)算效率不高。ASAGA就是將SAA嵌入到GA的循環(huán)體中,互相取長(zhǎng)補(bǔ)短,形成一種新的優(yōu)化算法。4. 1自適應(yīng)遺傳算法GA是模擬生物在自然環(huán)境中的遺傳和進(jìn)化過(guò)程而形成的一種自適應(yīng)全局優(yōu)化隨機(jī)搜索算法。由于它對(duì)搜索空間不作任何假設(shè),它既不要求搜索空間是光滑的,也不要求它是處處可微的,因而它能解決很大一類問(wèn)題,在飛行器軌道優(yōu)化中逐漸開(kāi)始得到應(yīng)用,但是要將GA成功的應(yīng)用到軌道優(yōu)化中,還必須處理好如何將終端約束反映到GA的適應(yīng)度函數(shù)中?？紤]到軌道優(yōu)化的精度,采用動(dòng)態(tài)罰函數(shù)法10對(duì)終端約束進(jìn)行處理,則適應(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為種群的個(gè)體;v(tf),(tf)和r(tf)為終端時(shí)刻的狀態(tài)值;1,2為懲罰因子,隨著優(yōu)化的進(jìn)行它們逐漸變大,使得對(duì)約束的限制也逐漸加強(qiáng)。在遺傳操作中,交叉操作是產(chǎn)生新個(gè)體的主要方法,它決定著GA的全局搜索能力;而變異操作只是產(chǎn)生新個(gè)體的輔助方法,它決定了GA的局部搜索能力。交叉與變異相互配合,共同完成對(duì)搜索空間的全局搜索和局部搜索。本文的交叉操作11采用如下方式:A=(1 -)A+BB=(1 -)B+A(12)如果A(B) R那么A(B)= R(13)第4期朱建豐等:基于自適應(yīng)模擬退火遺傳算法的月球軟著陸軌道優(yōu)化圖4徑向速

37、度的時(shí)間歷程Fig. 4Time histories of radial velocity圖5月心距的時(shí)間歷程Fig. 5Time histories of distance from lunar center to lander圖6角速度的時(shí)間歷程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ù)說(shuō)明了該算法的優(yōu)化精度很高。下面給出群體中最大適應(yīng)度值隨時(shí)間的變化曲線,如圖7所示。圖7說(shuō)明:在優(yōu)化進(jìn)程中,多次出現(xiàn)局部“早熟”現(xiàn)象,但都被ASAGA的搜索能力所克服,成功地跳出了局部最優(yōu)的誘惑。圖7最大適應(yīng)度的時(shí)間歷程Fig. 7Time history of maximum fitness8.參考文獻(xiàn)1陳剛,萬(wàn)自明,徐敏,等.遺傳算法在航天器軌道優(yōu)化中的應(yīng)用J.彈道學(xué)報(bào), 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)軟著陸兩點(diǎn)邊值問(wèn)題的數(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é)報(bào)(自然科學(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)/對(duì) 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)的最小值點(diǎn)， 且 l = x = u% f為待求函數(shù)，x0為初值點(diǎn)，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ù)在初始點(diǎn)x0處的函數(shù)值xo = x;fo = fx;%算法第三步，進(jìn)行迭代計(jì)算，找出近似全局最小點(diǎn)for k =0:kmax Ti = (k/kmax)q; mu = 10(Ti*100); % 計(jì)算mu dx = Mu_Inv(2*rand(size(x)-1,mu).*(u - l);%步長(zhǎng)dx x1 = x + dx; %下一個(gè)估計(jì)點(diǎn) 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) +