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

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

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

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

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

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

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

7、區(qū)域后的狀況。四、符號說明 符號含義第i架飛機在初始時刻的橫坐標第j架飛機在初始時刻的橫坐標第i架飛機在t時刻的橫坐標第j架飛機在t時刻的橫坐標第i架飛機在初始時刻的縱坐標第j架飛機在初始時刻的縱坐標第i架飛機在t時刻的縱坐標第j架飛機在t時刻的縱坐標第i架飛機在初始時刻飛行方向與X軸正向的夾角第j架飛機在初始時刻飛行方向與X軸正向的夾角第i架飛機在t時刻飛行方向與X軸正向的夾角第j架飛機在t時刻飛行方向與X軸正向的夾角第i架飛機飛行方向角的調(diào)整幅度第j架飛機飛行方向角的調(diào)整幅度第i架飛機在規(guī)定區(qū)域內(nèi)可能飛行的最長時間第j架飛機在規(guī)定區(qū)域內(nèi)可能飛行的最長時間T是一個6*6矩陣；=min，V飛機

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

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

10、軟件求解，程序見附錄2，計算結(jié)果見表1。圖3得出矩陣各元素值如下表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架

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

12、示的三種情況，其中對稱軸為-/2,設此時對應的t值為。 >0 <0,0<< <0,> 圖4當>0 ，可知，很明顯，f(0)>0；當 <0,0<< ，則只要=；當 <0,>，只需最右端滿足 （）0；則模型的約束條件可以轉(zhuǎn)化如下，這個非線性模型僅是關于一個變量，即關于方向角的約束條件。S.t. （）0 (當且0<<) .接下來，通過查閱文獻資料（見附錄），了解了LINGO軟件的運算性能和原理，綜合前人的經(jīng)驗，我們將考慮對原目標函數(shù) min 調(diào)整，在同樣滿足問題分析中的目標的情況下，即所有六架飛機的飛行方向調(diào)整

13、角度均盡量小，對其適宜的優(yōu)化。由于min 是不光滑的，我們考慮用一個光滑的目標函數(shù)來替代，即min 。這樣在運用LINGO軟件進行求解時，運行速度和求解精度上都將大大提高。綜上所述，最終的模型為：Obj min S.t. （）0 (當且0<<) .3） 模型的求解 我們針對最終的模型，帶入N*N矩陣，運用LINGO軟件編輯程序（見附錄3），求解模型。合理編輯算法并使用全局求解程序運行，運行結(jié)果見表2，發(fā)現(xiàn)運行時間已經(jīng)達到預想的反應速度（多次試驗平均時間3s）.表2飛機123456角度改變0.0000000.0000002.062449-0.49543730.0000001.5670

14、116.954677運行結(jié)果如上表，可得出飛行調(diào)整策略，即只改變其中兩架飛機，飛機1,2,5保持航向，飛機3調(diào)整+2.062449度，飛機4調(diào)整-0.4954373，飛機6調(diào)整+1.567011度。= 6.954677，綜合整體調(diào)整量較小已達到最優(yōu)，結(jié)果較理想。 為達到更加理想的決策速度（實際上這是符合現(xiàn)實需求的），我們考慮加快程序的運行速度，這樣機上的乘客和工作人員的安全將得到更好的保證。我們考慮是否可以通過對模型中的飛行角度進行預估，對函數(shù)中的飛行角度進行更好的約束預估，縮小搜索區(qū)間，以優(yōu)化算法，提高運行速度。理由一： 同時多次實踐發(fā)現(xiàn)，最優(yōu)解總是落在10度以內(nèi)，甚至更小。我們分析各架飛機

15、的方位圖，作出假設 ，如右圖，區(qū)域內(nèi)的飛機做可以忽略不計的微調(diào)，只對新進入的飛機6進行航向的調(diào)整，此時飛機6要做出的最優(yōu)調(diào)整角度應該是小于原來的搜索區(qū)間最大值30度。理由二： 我們分析即使兩架飛機相向飛行，根據(jù)題目條件4）可知進入該區(qū)域的飛機在到達區(qū)域邊緣時，與區(qū)域內(nèi)飛機的距離應在60公里以上，將此數(shù)據(jù)帶入，我們可做一定的估算。以飛機j為參考，飛機i正對飛機j的飛行，調(diào)整一定角度到兩架飛機恰好不相擦，如圖5圖5其中為飛機i調(diào)整后對j的相對飛行速度，角是調(diào)整角，則有：=8km/60km;<因此，我們最后在編程時將搜索區(qū)間縮小在（，），如此再將新的程序?qū)υ瓟?shù)據(jù)進行計算，最后時間成功縮短至1.

16、5s。計算結(jié)果如下表3 表3飛機123456角度改變0.0000000.0000002.062456-0.49542980.0000001.5670046.954677 六、模型檢驗與分析模型在實際應用中應該是要將決策時間縮短趨近于零，即時間短到可以忽略不計，這樣才能給飛機以足夠的時間做機械性調(diào)整，如果可以把程序運行時間進行補償，即對飛機在決策延遲延遲時間后調(diào)接到調(diào)整信息之前的位置預測，將其帶入程序運算，得出一個較精確的解，這個決策應該是更加可靠的，可以用來檢驗我們的模型?，F(xiàn)在我們將預測后數(shù)據(jù)如表3（計算程序見附錄5）表4飛機編號橫坐標x縱坐標y方向角（度）1149.798239.604024

17、3284.751584.63152363149.6620154.7114220.54144.585150.1593159512937143149.6595230新進入0.27360.350252得到新的矩陣如下表5表5T1T2T3T4T5T6T100.1497150.240230.226970.240230.24023T20.14971500.1497150.1497150.1497150.149715T30.240230.14971500.226970.2475620.269068T40.226970.1497150.2269700.226970.22697T50.240230.149715

18、0.2475620.2269700.247562T60.240230.1497150.2690680.226970.2475620將新矩陣數(shù)據(jù)代入程序，計算出此時的調(diào)整方案，并與調(diào)整兩秒之間的計算結(jié)果通過作差和誤差計算作比較，如下表6表6項目角度飛機延時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.

19、5736810.0066770.424%6.9546777.0115880.0569110.812%很顯然誤差較小，在保證決策精確度的情況下，該誤差可以忽略不計，在實際問題中意味著我們的模型精確度已經(jīng)達到，并且運行時間較好的控制在在允許范圍之內(nèi)，也就是說決策效率依然很高，完全滿足要求。7、 模型的評價與推廣1）縱觀整個模型，我們的模型具有一般性，用兩架飛機的歐式距離的限制來表示不碰撞條件，對于三維甚至多維的情況可以推廣。2）該模型對1000米的高空平面內(nèi)的飛機飛行進行管理，是滿足實際的。實際應用中，為避免兩機碰撞，管理中心還可以控制飛機調(diào)整在統(tǒng)一空域的不同高度，這樣調(diào)整難度大大降低，此模型考慮

20、處于同一高度平面內(nèi)的6架飛機已經(jīng)足夠，已具有實際意義。3)經(jīng)優(yōu)化模型及其算法，決策較準確，運行較快速，已經(jīng)將決策效率提高，滿足實際要求。4）模型也有一些不足，比如在計算組成的N*N矩陣是使用的是MATLAB軟件，后面的模型優(yōu)化運算卻使用LINGO軟件，如果可以深入研究軟件，將兩步運算連接在統(tǒng)一軟件中運算，效率將更理想。5)模型計算思路也有待進一步改善，比如，在實際中，可以將在模型檢驗與分析環(huán)節(jié)的思路進一步深入，不僅考慮決策時間的延遲，甚至加入飛機進行機械調(diào)整所需的經(jīng)驗平均時間。這樣將更加精準化。6）該模型其實可能應該有更加一般的解法，這是物理學里面一個運動的問題，我們可以將在模型求解環(huán)節(jié)對縮短

21、搜索區(qū)間的證明里的思路進一步深化，即是否可以繼續(xù)采用相對運動原理優(yōu)化模型，列出新的不碰撞控制條件。八、參考文獻1顧其行 ，國際航空運輸管理, 北京 ：知識出版社，19872朱道元，數(shù)學建模精品案例，南京：東南大學出版社，19993程極泰，最優(yōu)設計的數(shù)學方法，北京：國防工業(yè)出版社，19944張懷中，運籌學基礎教程，上海：上海人民出版社，20105韓中庚，實用運籌學，北京：清華大學出版社 ， 20076謝金星、薛毅，優(yōu)化建模與LINDO/LINGO軟件，北京：清華大學出版社，20057曹華林，湯志高，金平，Lingo在飛行管理中的應用，科技信息(科學教研) 2007年17期 九、附錄1% 函數(shù)pi

22、cture做飛機航線預測圖function picture(X,Y,a)% 輸入X:飛機原始橫坐標向量% 輸入Y:飛機原始縱坐標向量% 輸入a:飛機原始方向角向量figureaxis(0,160,0,160)for i=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') %紅色虛線為向上飛的飛機航線 else plot(x,y,'-b') %藍色直線為向下飛的飛機航線 endendhold onfor i=1:length

23、(X) plot(X(i),Y(i),'*k') %黑色的“*”號為飛機初始位置end2% 函數(shù)time計算任意兩架飛機中最先飛出區(qū)域的一架所用最少時間矩陣function A=time(X,Y,Z,W)% 輸入X:飛機原始橫坐標向量% 輸入Y:飛機原始縱坐標向量% 輸入Z：區(qū)域四角橫坐標向量% 輸入W：區(qū)域四角縱坐標向量% 計算每架飛機到區(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% 計算每架飛機距區(qū)域四角

24、最長距離向量ss=zeros(1,length(X);for i=1:length(X) for j=1:length(Z) if s(i)<=A(i,j) s(i)=A(i,j); end endend% 計算任意兩架飛機飛出區(qū)域所用最少時間矩陣Tt=s./800T=zeros(length(t),length(t);for i=1:length(t) for j=1:length(t) if t(i)=t(j) T(i,j)=0; else if t(i)>=t(j) T(i,j)=t(j); else T(i,j)=t(i); end end endend% 輸出：任意兩架飛

25、機中最先飛出區(qū)域的一架所用最少時間矩陣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:!輸入飛機初始坐標和角度（角度制）;x0 y0 cita0 =150402438585236150155220.5145501591301502300052; max_cita = 30;!每兩架飛機飛出區(qū)域所用最短時間矩陣，數(shù)據(jù)來自time.m文件 ;T = 0 0

26、.1503 0.2401 0.2275 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)*