一個飛行管理問題數(shù)模競賽

1、 一個飛行管理問題 摘要 在某一空域里對飛機(jī)的飛行合理管理事關(guān)重大，比如乘客及機(jī)上工作人員生命財產(chǎn)安全和航空公司的運(yùn)作效益等。本文通過對飛機(jī)飛行管理問題的研究，得到了調(diào)整飛機(jī)架數(shù)較少同時調(diào)整幅度均最?。ㄆ椒胶妥钚。┑娘w行管理最優(yōu)安排的非線性模型，這樣既使得乘客所受影響達(dá)到最少，也便于飛機(jī)調(diào)整，還有利于飛機(jī)回到原來的航線，同時還在決策時間上對模型進(jìn)行了優(yōu)化和調(diào)整。本文不僅一般性地將不相撞的問題轉(zhuǎn)化為歐式距離控制，而且很巧妙的將不碰撞條件轉(zhuǎn)化成簡單的二次函數(shù)標(biāo)準(zhǔn)形式進(jìn)行含參討論，建立一個只含有轉(zhuǎn)向角變量的模型。并且再次很妙的具體化區(qū)域內(nèi)受控時間形成矩陣，大大得簡化運(yùn)算，節(jié)約了大量運(yùn)算的時間，便于管

2、理人員控制操作，從而確保飛機(jī)的安全。更重要的是最后結(jié)合實際縮短了搜索區(qū)間，并優(yōu)化算法，使得決策更加高效。最后的延時檢驗也充分體現(xiàn)了模型的可靠性。關(guān)鍵字：歐氏距離 約束轉(zhuǎn)化 縮短搜索區(qū)間 時間矩陣 延時檢驗16一、問題重述在約10000米的高空某邊長為160公里的正方形區(qū)域內(nèi)，經(jīng)常有若干架飛機(jī)作水平飛行。區(qū)域內(nèi)每架飛機(jī)的位置和速度向量均由計算機(jī)記錄其數(shù)據(jù)，以便進(jìn)行飛行管理。當(dāng)一駕欲進(jìn)入該區(qū)域的飛機(jī)到達(dá)區(qū)域邊緣時，記錄其數(shù)據(jù)后，要立即計算并判斷是否會與區(qū)域內(nèi)的飛機(jī)發(fā)生碰撞。如果會碰撞，則應(yīng)計算如何調(diào)整各架（包括新進(jìn)入的）飛機(jī)飛行的方向角，以避免碰撞?，F(xiàn)假定條件如下： 1）不碰撞的標(biāo)準(zhǔn)為任意兩架飛機(jī)

3、的距離大于8公里； 2）飛機(jī)飛行方向角調(diào)整的幅度不應(yīng)超過30度； 3）所有飛機(jī)飛行速度均為每小時800公里； 4）進(jìn)入該區(qū)域的飛機(jī)在到達(dá)區(qū)域邊緣時，與區(qū)域內(nèi)飛機(jī)的距離應(yīng)在60公里以上； 5）最多需考慮6架飛機(jī)； 6）不必考慮飛機(jī)離開此區(qū)域后的狀況。 請你對這個避免碰撞的飛機(jī)管理問題建立數(shù)學(xué)模型，列出計算步驟，對以下數(shù)據(jù)進(jìn)行計算（方向角誤差不超過0.01度），要求飛機(jī)飛行方向角調(diào)整的幅度盡量小。 設(shè)該區(qū)域內(nèi)4個頂點的坐標(biāo)為(0,0),(160,0),(160,160),(0,160)。 記錄數(shù)據(jù)為：飛機(jī)編號橫坐標(biāo)x縱坐標(biāo)y方向角（度）115040243285852363150155220.541

4、45501595130150230新進(jìn)入0052注：方向角指飛行方向與x軸方向的夾角。 試根據(jù)實際應(yīng)用背景對你的模型進(jìn)行評價與推廣。2、 問題分析1初步分析 根據(jù)問題容易知道，這顯然是一個優(yōu)化問題，當(dāng)兩架飛機(jī)可能發(fā)生碰撞時，即在規(guī)定區(qū)域內(nèi)某一時刻兩架飛機(jī)之間的距離小于8公里，因此要調(diào)整飛行方向一定角度，保證任意兩架飛機(jī)在區(qū)域內(nèi)任意時刻，兩者的距離均不小于8公里，避免相撞?？紤]到調(diào)整角度應(yīng)盡量小，可以簡化飛行方向調(diào)整策略，降低調(diào)整難度，同時減輕機(jī)內(nèi)乘客及工作人員的不適。此外由此初步確定了調(diào)整目標(biāo)：所有六架飛機(jī)的飛行方向調(diào)整角度均盡量小。 2解決方案 由于所有飛機(jī)均處于1000米得高空作水平飛行，

5、可將飛機(jī)飛行的空域視為二維平面xoy中的一個正方形區(qū)域,頂點為(0,0),(160,0),(160,160),(0,160)。于是可以引入時間變量后，確定每架飛機(jī)在任意時刻的坐標(biāo)，列出任意兩點的歐氏距離，令其恒大于8公里，則得出一個重要約束條件。再結(jié)合變化角度應(yīng)小于30度，即可得出約束條件，然后運(yùn)用LINGO軟件編輯程序進(jìn)行求解。為提高決策效率，在反復(fù)試驗中又可對約束條件進(jìn)行調(diào)整。三、條件假設(shè)1. 不碰撞的標(biāo)準(zhǔn)為任意兩架飛機(jī)的距離在以后任何一個時間里大于8公里；2. 飛機(jī)飛行方向角調(diào)整的幅度不應(yīng)超過30度； 3. 所有飛機(jī)飛行速度均為每小時800公里； 4. 進(jìn)入該區(qū)域的飛機(jī)在到達(dá)區(qū)域邊緣時，

6、與區(qū)域內(nèi)飛機(jī)的距離應(yīng)在60公里以上。即在計算如何最優(yōu)地調(diào)整各架（包括新進(jìn)入的）飛機(jī)飛行的方向角時，飛行管理中心得出合理的最優(yōu)調(diào)整措施，；5. 最多需考慮6架飛機(jī)。 6. 此處忽略飛機(jī)在執(zhí)行過程中所需耗費(fèi)的時間，即假設(shè)從飛機(jī)管理中心發(fā)出的調(diào)整信息飛機(jī)馬上可以接收并執(zhí)行，不存在滯后或延遲；7. 飛行管理中心在計算飛行調(diào)整信號和發(fā)出信號所需時間內(nèi)，忽略各架飛機(jī)（包括剛進(jìn)入的飛機(jī)）調(diào)整航向前飛行數(shù)據(jù)的變化；8. 假定飛機(jī)在該區(qū)域內(nèi)完全依賴飛行管理中心調(diào)度；9. 假設(shè)飛機(jī)在飛出區(qū)域之后，飛行員可以自覺調(diào)整飛行策略，回歸原始航線，即飛行管理中心不必考慮飛機(jī)離開此區(qū)域后的狀況。四、符號說明 符號含義第i架飛

7、機(jī)在初始時刻的橫坐標(biāo)第j架飛機(jī)在初始時刻的橫坐標(biāo)第i架飛機(jī)在t時刻的橫坐標(biāo)第j架飛機(jī)在t時刻的橫坐標(biāo)第i架飛機(jī)在初始時刻的縱坐標(biāo)第j架飛機(jī)在初始時刻的縱坐標(biāo)第i架飛機(jī)在t時刻的縱坐標(biāo)第j架飛機(jī)在t時刻的縱坐標(biāo)第i架飛機(jī)在初始時刻飛行方向與X軸正向的夾角第j架飛機(jī)在初始時刻飛行方向與X軸正向的夾角第i架飛機(jī)在t時刻飛行方向與X軸正向的夾角第j架飛機(jī)在t時刻飛行方向與X軸正向的夾角第i架飛機(jī)飛行方向角的調(diào)整幅度第j架飛機(jī)飛行方向角的調(diào)整幅度第i架飛機(jī)在規(guī)定區(qū)域內(nèi)可能飛行的最長時間第j架飛機(jī)在規(guī)定區(qū)域內(nèi)可能飛行的最長時間T是一個6*6矩陣；=min，V飛機(jī)的飛行速度飛機(jī)i與飛機(jī)j的歐氏距離 五、模型

8、建立與求解1問題簡化首先，如果對六架飛機(jī)在區(qū)域內(nèi)做實時監(jiān)控，再做多次調(diào)整，則每作一次航向調(diào)整都要進(jìn)行一次決策，這將使問題復(fù)雜化，總體計算量較大，同時實際問題中計算也要耗費(fèi)時間，效率大大降低，飛機(jī)控制的安全性必然會降低。并且對問題所給原始數(shù)據(jù)利用MATLAB軟件（程序見附錄1）作出原始航線圖，如圖1，可以粗略驗證一次調(diào)整可行，既可以避免相撞，又簡單易行。圖1結(jié)論一：我們認(rèn)為只做一次調(diào)整是優(yōu)于多次調(diào)整的。其次，分析飛機(jī)飛行軌跡，作如圖2，圖2某架飛機(jī)在點A的航向是由A到C，此時次飛機(jī)在點C可能遭碰撞，D點是一個安全位置，如果飛機(jī)在A點早做航向的調(diào)整，方向角變化角為 ；如果待飛機(jī)飛到B點時再做航向的

9、調(diào)整，則方向的變化角為 ，很顯然， 。由此我們得出結(jié)論，早調(diào)整航向優(yōu)于晚調(diào)整航向，因此我們可將問題簡化。結(jié)論二：在新進(jìn)入飛機(jī)到達(dá)區(qū)域邊緣時便對所有飛機(jī)做航向調(diào)整，避免碰撞。2模型的導(dǎo)出過程1）模型推導(dǎo)首先引入時間變量t,其中t的取值區(qū)間為（0，）；另外，為了提高程序運(yùn)行效率，簡化計算，我們引入一個矩陣來控制搜索循環(huán)，較合理的設(shè)定了變量的上下界，從而大大提高決策效率。如圖3所示，表示第i架飛機(jī)在規(guī)定區(qū)域內(nèi)可能飛行的最長時間，即該架飛機(jī)飛至區(qū)域邊界四個角上點所用時間的最大值，令=min，則所有的可以構(gòu)成一個N*N的矩陣。此矩陣可以運(yùn)用MATLAB軟件求解，程序見附錄2，計算結(jié)果見表1。圖3得出矩陣

10、各元素值如下表1表1T1T2T3T4T5T6T100.150260.2401170.2275030.2401170.240117T20.1502600.150260.150260.150260.15026T30.2401170.1502600.2275030.2481180.269621T40.2275030.150260.22750300.2275030.227503T50.2401170.150260.2481180.22750300.248118T60.2401170.150260.2696210.2275030.2481180將第i架飛機(jī)在t時刻的位置用坐標(biāo)（，）表示；其中=+Vt c

11、os（+）；=+Vt sin（+）；則歐式距離 =8（t）=-64= 0；綜上所述，可得初步的模型如下：Obj min S.t. （t）0 30（i=1，2，3，4，5，6）；2）模型優(yōu)化 在考慮LINGO軟件的程序編輯和運(yùn)行時，我們發(fā)現(xiàn)時間變量t不宜為連續(xù)變量，因此我們通過對函數(shù) （t）的化簡，對約束條件進(jìn)行轉(zhuǎn)化，消去變量t。過程如下：（t）=+ 其中=2Vt； =2-(-)+(-)；=+-64；很明顯模型的約束可以等價的轉(zhuǎn)化為一個一元二次函數(shù)在給定區(qū)間內(nèi)恒大于零的討論。表示的是一條開口向上的拋物線。分別從三個方面討論，如圖4所示的三種情況，其中對稱軸為-/2,設(shè)此時對應(yīng)的t值為。 0 0,

12、0 圖4當(dāng)0 ，可知，很明顯，f(0)0；當(dāng) 0,0 ，則只要=；當(dāng) ，只需最右端滿足 （）0；則模型的約束條件可以轉(zhuǎn)化如下，這個非線性模型僅是關(guān)于一個變量，即關(guān)于方向角的約束條件。S.t. （）0 (當(dāng)且0) .接下來，通過查閱文獻(xiàn)資料（見附錄），了解了LINGO軟件的運(yùn)算性能和原理，綜合前人的經(jīng)驗，我們將考慮對原目標(biāo)函數(shù) min 調(diào)整，在同樣滿足問題分析中的目標(biāo)的情況下，即所有六架飛機(jī)的飛行方向調(diào)整角度均盡量小，對其適宜的優(yōu)化。由于min 是不光滑的，我們考慮用一個光滑的目標(biāo)函數(shù)來替代，即min 。這樣在運(yùn)用LINGO軟件進(jìn)行求解時，運(yùn)行速度和求解精度上都將大大提高。綜上所述，最終的模型為

13、：Obj min S.t. （）0 (當(dāng)且0) .3） 模型的求解 我們針對最終的模型，帶入N*N矩陣，運(yùn)用LINGO軟件編輯程序（見附錄3），求解模型。合理編輯算法并使用全局求解程序運(yùn)行，運(yùn)行結(jié)果見表2，發(fā)現(xiàn)運(yùn)行時間已經(jīng)達(dá)到預(yù)想的反應(yīng)速度（多次試驗平均時間3s）.表2飛機(jī)123456角度改變0.0000000.0000002.062449-0.49543730.0000001.5670116.954677運(yùn)行結(jié)果如上表，可得出飛行調(diào)整策略，即只改變其中兩架飛機(jī)，飛機(jī)1,2,5保持航向，飛機(jī)3調(diào)整+2.062449度，飛機(jī)4調(diào)整-0.4954373，飛機(jī)6調(diào)整+1.567011度。= 6.95

14、4677，綜合整體調(diào)整量較小已達(dá)到最優(yōu)，結(jié)果較理想。 為達(dá)到更加理想的決策速度（實際上這是符合現(xiàn)實需求的），我們考慮加快程序的運(yùn)行速度，這樣機(jī)上的乘客和工作人員的安全將得到更好的保證。我們考慮是否可以通過對模型中的飛行角度進(jìn)行預(yù)估，對函數(shù)中的飛行角度進(jìn)行更好的約束預(yù)估，縮小搜索區(qū)間，以優(yōu)化算法，提高運(yùn)行速度。理由一： 同時多次實踐發(fā)現(xiàn)，最優(yōu)解總是落在10度以內(nèi)，甚至更小。我們分析各架飛機(jī)的方位圖，作出假設(shè) ，如右圖，區(qū)域內(nèi)的飛機(jī)做可以忽略不計的微調(diào)，只對新進(jìn)入的飛機(jī)6進(jìn)行航向的調(diào)整，此時飛機(jī)6要做出的最優(yōu)調(diào)整角度應(yīng)該是小于原來的搜索區(qū)間最大值30度。理由二： 我們分析即使兩架飛機(jī)相向飛行，根據(jù)

15、題目條件4）可知進(jìn)入該區(qū)域的飛機(jī)在到達(dá)區(qū)域邊緣時，與區(qū)域內(nèi)飛機(jī)的距離應(yīng)在60公里以上，將此數(shù)據(jù)帶入，我們可做一定的估算。以飛機(jī)j為參考，飛機(jī)i正對飛機(jī)j的飛行，調(diào)整一定角度到兩架飛機(jī)恰好不相擦，如圖5圖5其中為飛機(jī)i調(diào)整后對j的相對飛行速度，角是調(diào)整角，則有：=8km/60km;因此，我們最后在編程時將搜索區(qū)間縮小在（，），如此再將新的程序?qū)υ瓟?shù)據(jù)進(jìn)行計算，最后時間成功縮短至1.5s。計算結(jié)果如下表3 表3飛機(jī)123456角度改變0.0000000.0000002.062456-0.49542980.0000001.5670046.954677 六、模型檢驗與分析模型在實際應(yīng)用中應(yīng)該是要將決策

16、時間縮短趨近于零，即時間短到可以忽略不計，這樣才能給飛機(jī)以足夠的時間做機(jī)械性調(diào)整，如果可以把程序運(yùn)行時間進(jìn)行補(bǔ)償，即對飛機(jī)在決策延遲延遲時間后調(diào)接到調(diào)整信息之前的位置預(yù)測，將其帶入程序運(yùn)算，得出一個較精確的解，這個決策應(yīng)該是更加可靠的，可以用來檢驗我們的模型。現(xiàn)在我們將預(yù)測后數(shù)據(jù)如表3（計算程序見附錄5）表4飛機(jī)編號橫坐標(biāo)x縱坐標(biāo)y方向角（度）1149.798239.6040243284.751584.63152363149.6620154.7114220.54144.585150.1593159512937143149.6595230新進(jìn)入0.27360.350252得到新的矩陣如下表5表5

17、T1T2T3T4T5T6T100.1497150.240230.226970.240230.24023T20.14971500.1497150.1497150.1497150.149715T30.240230.14971500.226970.2475620.269068T40.226970.1497150.2269700.226970.22697T50.240230.1497150.2475620.2269700.247562T60.240230.1497150.2690680.226970.2475620將新矩陣數(shù)據(jù)代入程序，計算出此時的調(diào)整方案，并與調(diào)整兩秒之間的計算結(jié)果通過作差和誤差計算

18、作比較，如下表6表6項目角度飛機(jī)延時0秒延時兩秒差值誤差10.000000 0.0000000.000000020.000000 0.0000000.000000032.062456 2.0707610.0083050.401%4-0.4954298 -0.4970576-0.00162780.327%50.000000 0.0000000.000000061.567004 1.5736810.0066770.424%6.9546777.0115880.0569110.812%很顯然誤差較小，在保證決策精確度的情況下，該誤差可以忽略不計，在實際問題中意味著我們的模型精確度已經(jīng)達(dá)到，并且運(yùn)行時間

19、較好的控制在在允許范圍之內(nèi)，也就是說決策效率依然很高，完全滿足要求。7、 模型的評價與推廣1）縱觀整個模型，我們的模型具有一般性，用兩架飛機(jī)的歐式距離的限制來表示不碰撞條件，對于三維甚至多維的情況可以推廣。2）該模型對1000米的高空平面內(nèi)的飛機(jī)飛行進(jìn)行管理，是滿足實際的。實際應(yīng)用中，為避免兩機(jī)碰撞，管理中心還可以控制飛機(jī)調(diào)整在統(tǒng)一空域的不同高度，這樣調(diào)整難度大大降低，此模型考慮處于同一高度平面內(nèi)的6架飛機(jī)已經(jīng)足夠，已具有實際意義。3)經(jīng)優(yōu)化模型及其算法，決策較準(zhǔn)確，運(yùn)行較快速，已經(jīng)將決策效率提高，滿足實際要求。4）模型也有一些不足，比如在計算組成的N*N矩陣是使用的是MATLAB軟件，后面的

20、模型優(yōu)化運(yùn)算卻使用LINGO軟件，如果可以深入研究軟件，將兩步運(yùn)算連接在統(tǒng)一軟件中運(yùn)算，效率將更理想。5)模型計算思路也有待進(jìn)一步改善，比如，在實際中，可以將在模型檢驗與分析環(huán)節(jié)的思路進(jìn)一步深入，不僅考慮決策時間的延遲，甚至加入飛機(jī)進(jìn)行機(jī)械調(diào)整所需的經(jīng)驗平均時間。這樣將更加精準(zhǔn)化。6）該模型其實可能應(yīng)該有更加一般的解法，這是物理學(xué)里面一個運(yùn)動的問題，我們可以將在模型求解環(huán)節(jié)對縮短搜索區(qū)間的證明里的思路進(jìn)一步深化，即是否可以繼續(xù)采用相對運(yùn)動原理優(yōu)化模型，列出新的不碰撞控制條件。八、參考文獻(xiàn)1顧其行 ，國際航空運(yùn)輸管理, 北京 ：知識出版社，19872朱道元，數(shù)學(xué)建模精品案例，南京：東南大學(xué)出版社

21、，19993程極泰，最優(yōu)設(shè)計的數(shù)學(xué)方法，北京：國防工業(yè)出版社，19944張懷中，運(yùn)籌學(xué)基礎(chǔ)教程，上海：上海人民出版社，20105韓中庚，實用運(yùn)籌學(xué)，北京：清華大學(xué)出版社 ， 20076謝金星、薛毅，優(yōu)化建模與LINDO/LINGO軟件，北京：清華大學(xué)出版社，20057曹華林，湯志高，金平，Lingo在飛行管理中的應(yīng)用，科技信息(科學(xué)教研) 2007年17期 九、附錄1% 函數(shù)picture做飛機(jī)航線預(yù)測圖function picture(X,Y,a)% 輸入X:飛機(jī)原始橫坐標(biāo)向量% 輸入Y:飛機(jī)原始縱坐標(biāo)向量% 輸入a:飛機(jī)原始方向角向量figureaxis(0,160,0,160)for i

22、=1:length(X) hold on x=0:0.1:160; y=tan(a(i)*pi/180)*(x-X(i)+Y(i); if a(i)=180 plot(x,y,-r) %紅色虛線為向上飛的飛機(jī)航線 else plot(x,y,-b) %藍(lán)色直線為向下飛的飛機(jī)航線 endendhold onfor i=1:length(X) plot(X(i),Y(i),*k) %黑色的“*”號為飛機(jī)初始位置end2% 函數(shù)time計算任意兩架飛機(jī)中最先飛出區(qū)域的一架所用最少時間矩陣function A=time(X,Y,Z,W)% 輸入X:飛機(jī)原始橫坐標(biāo)向量% 輸入Y:飛機(jī)原始縱坐標(biāo)向量% 輸

23、入Z：區(qū)域四角橫坐標(biāo)向量% 輸入W：區(qū)域四角縱坐標(biāo)向量% 計算每架飛機(jī)到區(qū)域四角的距離矩陣Sfor i=1:length(X) for j=1:length(Z) A(i,j)=sqrt(X(i)-Z(j)*(X(i)-Z(j)+(Y(i)-W(j)*(Y(i)-W(j); endend% 計算每架飛機(jī)距區(qū)域四角最長距離向量ss=zeros(1,length(X);for i=1:length(X) for j=1:length(Z) if s(i)=t(j) T(i,j)=t(j); else T(i,j)=t(i); end end endend% 輸出：任意兩架飛機(jī)中最先飛出區(qū)域的一架所

24、用最少時間矩陣TT3MODEL:TITLE 飛行管理問題;SETS:Plane/1.6/: x0, y0, cita0, cita1, d_cita;link(plane, plane)|&1 #LT# &2: b,c,z;link1(plane, plane):T;ENDSETSDATA:!輸入飛機(jī)初始坐標(biāo)和角度（角度制）;x0 y0 cita0 =150402438585236150155220.5145501591301502300052; max_cita = 30;!每兩架飛機(jī)飛出區(qū)域所用最短時間矩陣，數(shù)據(jù)來自time.m文件 ;T = 0 0.1503 0.2401 0.2275

25、0.2401 0.2401 0.1503 0 0.1503 0.1503 0.1503 0.1503 0.2401 0.1503 0 0.2275 0.2481 0.2696 0.2275 0.1503 0.2275 0 0.2275 0.2275 0.2401 0.1503 0.2481 0.2275 0 0.2481 0.2401 0.1503 0.2696 0.2275 0.2481 0 ;V=800;ENDDATAINIT:d_cita = 0 0 0 0 0 0;ENDINITfor(plane(i): cita1(i) - cita0(i) = d_cita(i);!模型中二次函數(shù)系數(shù);for(link(i,j): b(i,j) = -2*(x0(i) -x0(j)*sin (cita1(i)+cita1(j)*3.14159265/360) +2*(y0(i) -y0(j)*cos (cita1(i)+cita1(j)*3