12.1一個飛行管理問題

在中國大學生數學建模競賽（Chinaundergraduatemathematicalcontestinmodeling,CUMCM）中，曾經出現(xiàn)過大量的優(yōu)化建模賽題.本章從中選擇了部分典型賽題，舉例分析其優(yōu)化建模過程，說明如何應用LINDO/LINGO軟件包求解這些賽題.12.1一個飛行管理問題12.1.1問題描述1995年全國大學生數學建模競賽中的A題（“一個飛行管理問題”）.在約10000m高空的某邊長為160km的正方形區(qū)域內，經常有若干架飛機作水平飛行.區(qū)域內每架飛機的位置和速度向量均由計算機記錄其數據，以便進行飛行管理.當一架欲進入該區(qū)域的飛機到達區(qū)域邊緣時，記錄其數據后，要立即計算并判斷是否會與區(qū)域內的飛機發(fā)生碰撞.如果會碰撞，則應計算如何調整各架（包括新進入的）飛機飛行的方向角，以避免碰撞.現(xiàn)假定條件如下：不碰撞的標準為任意兩架飛機的距離大于8km；飛機飛行方向角調整的幅度不應超過30°；所有飛機飛行速度均為800km/h；進入該區(qū)域的飛機在到達該區(qū)域邊緣時，與該區(qū)域內的飛機的距離應在60km以上；最多需考慮6架飛機；不必考慮飛機離開此區(qū)域后的狀況.請你對這個避免碰撞的飛行管理問題建立數學模型，列出計算步驟，對以下數據進行計算（方向角誤差不超過0.01°），要求飛機飛行方向角調整的幅度盡量小.設該區(qū)域4個頂點的坐標為（0，0），（160，0），（160，160），（0，160）.記錄數據見表12-1.表12-1飛機位置和方向角記錄數據飛機編號橫坐標縱坐標方向角飛機編號橫坐標縱坐標方向角11501402434145501592858523651301502303150155220.5新進入0052說明：方向角指飛行方向與ｘ軸正向的夾角.試根據實際應用背景對你的模型進行評價和推廣.12.1.2模型1及求解模型建立這個問題顯然是一個優(yōu)化問題.設第ｉ架飛機在調整時的方向角為QUOTE（題目中已經給出），調整后的方向角為QUOTE=QUOTE+QUOTE（QUOTE=1,2,…，6）.題目中就是要求飛機飛行方向角調整的幅度盡量小，因此優(yōu)化的目標函數可以是QUOTE. （1）為了建立這個問題的優(yōu)化模型，只須要明確約束條件就可以了.一個簡單的約束是飛機飛行方向角調整的不應超過30°，即|QUOTE|QUOTE30°. (2)題中要求進入該區(qū)域的飛機在到達該區(qū)域邊緣時，與該區(qū)域內飛機的距離應在60km以上，這個條件是個初始條件，很容易驗證目前所給數據是滿足的，因此本模型中可以不予考慮.剩下的關鍵是要滿足題目中描述的任意兩架飛機不碰撞的要求，即任意兩架位于該區(qū)域內的飛機的距離應大于8km.但這個問題的難點在于飛機是動態(tài)的，這個約束不好直接描述，為此我們首先需要描述每架飛機的飛行軌跡.記飛機飛行速率為v（=800km/h）,以當前時刻為0時刻.設第ｉ架飛機在調整時的位置坐標為（QUOTE，QUOTE）（已知條件），t時刻的位置坐標為（QUOTE，QUOTE），則QUOTE=QUOTE+QUOTE，QUOTE=QUOTE+QUOTE.如果要嚴格表示兩架位于該區(qū)域內的飛機的距離應大于8km，則需考慮每架飛機在該區(qū)域內的飛行時間的長度.記QUOTE為第ｉ架飛機飛出去與的時刻，即QUOTE=QUOTEargmin{t>0:QUOTE+QUOTE=0或160， （4）或者QUOTE+QUOTE=0或160}. （5）記t時刻第ｉ架飛機與第ｊ架飛機的距離為QUOTE（t），并記QUOTE（t）=QUOTE-64，這時在該區(qū)域內飛機不相撞的約束條件就變成了QUOTE（t）=QUOTE-64QUOTE0（0QUOTEtQUOTE）.其中QUOTE=min{QUOTE，QUOTE}. (6)此外，經過計算，可以得到QUOTE（t）=QUOTE+QUOTE-64=QUOTE+QUOTE+QUOTE， （7）其中

QUOTE=2vtQUOTE， （8）QUOTE=2[-(QUOTE-QUOTE)QUOTE+QUOTE, （9）QUOTE=QUOTE+QUOTE-64. （10）所以，QUOTE（t）是一個關于QUOTE的二次函數，表示的是一條開口向上的拋物線.當QUOTE=-QUOTE/2，即t=-QUOTE/4vQUOTE（記為QUOTE）時，函數QUOTE（t）取最小值-QUOTE/4+QUOTE.注意到QUOTE（0）=QUOTE0（初始時刻不相撞）,如果-QUOTE/2QUOTE0（即QUOTE0），則此時約束條件（5）一定成立，所以對約束條件（5）只需考慮以下兩種可能情況：如果QUOTE0且QUOTE，只需要QUOTE（t）在右端點QUOTE的函數值非負即可，即QUOTE（QUOTE）QUOTE0； （11）如果QUOTE0且0QUOTE，只需要QUOTE（t）的最小值QUOTE（QUOTE）=-QUOTE/4+QUOTE0即可，即QUOTE4QUOTE0. (12)實際上，約束（11）表示的是QUOTE（t）在右端點QUOTE的函數值非負，這個約束在（12）的條件下也是自然成立的，所以可以對約束（11）不再附加QUOTE0且QUOTE的條件.于是，我們的模型就是QUOTE; （13）S．t．|QUOTE|QUOTE30°， （14）QUOTE（QUOTE）QUOTE0 (15)QUOTE4QUOTE0.（當QUOTE0且0QUOTE）. (16)模型求解上面這是一個非線性規(guī)劃模型，雖然是嚴格滿足題目要求的模型，但得到的模型邏輯關系比較復雜，約束（16）是在一定條件下才成立的約束，而且其中QUOTE的計算式（4）也是含有相當復雜的關系式，使用LINGO軟件不太容易將模型很方便地輸入，因為邏輯處理不是LINGO的優(yōu)勢所在.即使能想辦法把把這個模型輸入到LINGO,也不一定能求出好的解(筆者嘗試過,但LINGO運行時有時會出現(xiàn)系統(tǒng)錯誤,可能是系統(tǒng)有漏洞,無法繼續(xù)求解).而且,在實時飛行調度中顯然需要快速求解,所以下面我們想辦法簡化模型.這個模型麻煩之處就在于,要求嚴格表示兩架位于該區(qū)域內的飛機距離應大于8km,所以需要考慮每架飛機在該區(qū)域內的飛行時間,比較繁瑣.注意到區(qū)域對角線的長度只有160QUOTEkm，任何一架飛機在所考慮的范圍內停留的時間不會超過QUOTE=160QUOTE/800=0.2QUOTE0.283（h），因此這里我們簡化一下問題：不再單獨考慮每架飛機在該區(qū)域停留的時間QUOTE，而以最大時間（這里已經是一個常數）QUOTE代替之，此時所有QUOTE=QUOTE.這實際上強化了問題的要求,即考慮了有些飛機可能已經飛出該區(qū)域,但仍不允許兩架飛機的距離小于8km.這個簡化的模型可以如下輸入LINGO軟件：MODEL:TITLE飛行管理問題的非線性規(guī)劃模型;SETS:Plane/1..6/:x0,y0,cita0,cita1,d_cita;!cita0表示初始角度，cita1為調整后的角度,d_cita為調整的角度;link(plane,plane)|&1#LT#&2:b,c;ENDSETSDATA:x0y0cita0=1501402438585236150155220.5145501591301502300052;max_cita=30;T_max=0.283;V=800;ENDDATAINIT:d_cita=000000;ENDINIT@for(plane:cita1-cita0=d_cita);@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.14159265/360);c(i,j)=(x0(i)-x0(j))^2+(y0(i)-y0(j))^2-64;);!避免碰撞的條件;!右端點非負;@for(link(i,j):[Right](2*V*T_max*@sin((cita1(i)-cita1(j))*3.14159265/360))^2+b(i,j)*(2*V*T_max*@sin((cita1(i)-cita1(j))*3.14159265/360))+c(i,j)>0);!最小點非負;@for(link(i,j):[Minimum]@if(b(i,j)#lt#0#and#-b(i,j)/4/V/@sin((cita1(i)-cita1(j))*3.14159265/360)#gt#0#and#-b(i,j)/4/V/@sin((cita1(i)-cita1(j))*3.14159265/360)#lt#T_max,b(i,j)^2-4*c(i,j),-1)<0);!@for(link(i,j):@if(b(i,j)#lt#0,b(i,j)^2-4*c(i,j),-1)<0);@for(link:@free(b));!調整角度上下限，單位為角度;@for(plane:@bnd(-max_cita,d_cita,max_cita));[obj]MIN=@SUM(plane:(d_cita)^2);![obj]MIN=@SUM(plane:@abs(d_cita));END注意：上面的模型中的方向角單位一律用角度，而LINGO只接受弧度，所以程序中一律進行了轉換.求解這個模型,得到Localoptimalsolutionfound.Objectivevalue：295.4937ModelTitle:飛行管理問題的非線性規(guī)劃模型VariableValueReducedCostD_CITA(1)-10.5.9800.000000D_CITA(2)0.0000000.000000D_CITA(3)0.0000000.000000D_CITA(4)6.5154250.000000D_CITA(5)10.006810.000000D_CITA(6)6.5054250.000000這個結果得到的是一個局部極小點，調整角度較大.能找到更好的解嗎?如果不用全局求解程序,通常很難得到稍大規(guī)模的非線性規(guī)劃問題的全局最優(yōu)解.所以我們啟動LINGO全局求解程序求解這個模型,可以得到全局最優(yōu)解如下:Globaloptimalsolutionfoundatiteration:93Objectivevalue:6.953944ModelTitle:飛行管理問題的非線性規(guī)劃模型VariableValueD_CITA(1)0.2719480E-02D_CITA(2)0.5613433E-02D_CITA(3)2.059140D_CITA(4)-0.4985421D_CITA(5)-0.5407837E-03D_CITA(6)1.570129(這里只給出QUOTE的值)可以看到，在0.01°的誤差要求下，需要調整第3，4，6三架飛機的角度，分別調整2.06°，-0.50°，1.57°.調整量的平方和為6.95.其實,使用全局求解程序,通常也不一定要等到得到全局最優(yōu)解,而是觀察求解狀態(tài)窗口,看到一個較好的當前解(或當前最好解在較長時間內不發(fā)生變化)時,就可以終止程序,用當前最好的局部最優(yōu)解作為最后的惡結果.列如,對于本列,LINGO求出全局最優(yōu)解大約需要1min,而實際上5s內LINGO就得到了與全局最優(yōu)解類似的解.此外,上面的模型還可以進一步簡化,列如可以假設要求飛機永遠不相撞,即認為QUOTE為無限大，這時顯然約束（15）也是多余的,而且約束(16)中只需要QUOTE0的條件就可以了.也就是說，上面的程序中的對應部分（約束[Right]和[Minimum]）可以改寫為更簡單的形式：！有端點非負，不再需要；！最小點非負，簡化為以下形式；@for(link(i,j):@if(b（i，j）#lt#0，b（i，j）^2-4*c（i，j），-1）QUOTE0）；實際計算顯示，此時得到的結果與前面計算的結果幾乎沒有差別.備注優(yōu)化的目標函數除了QUOTE外，也可以設定為QUOTE或QUOTE等，用LINGO求解的過程是完全類似的，計算結果略有差異，這里就不再對兩個目標函數具體計算了.甚至可以考慮讓參與調整的飛機的數量盡量小，這種想法在實際中也不能說沒有道理，但與題目中的要求不符，而且解題難度并沒有減小，意義似乎不大.實際調度中，由于計算上面的調度方案需要時間，將調度信息告訴飛機駕駛員并作出調整方向角的操作也需要時間，因此如果考慮一定的反應滯后時間，應該是比較合理的.也就是說，如果反應時間是10s，則計算式中應采用飛機沿當前的方向角飛行10s以后的位置作為計算的基礎.12.1.3模型2及求解從12.1.2節(jié)可以看出，求解模型1的非線性規(guī)劃模型是比較困難的，輸入后也很可能找不到好的解甚至出現(xiàn)錯誤.此外，演示版軟件還會受到求解規(guī)模的限制，尤其可能無權使用全局求解程序.因此，如果能把這個問題簡化成比較簡單的規(guī)劃模型，將是非常有價值的.模型建立如圖12-1，把兩架飛機i，j分別看成半徑為4km的圓（圖中i，j為圓心），AB,CD為公切線，將AB和CD的夾角的一半稱為碰撞角.在調整時刻，第i架飛機與第j架飛機的碰撞角為QUOTE，則易知QUOTE=QUOTE）， （17）其中QUOTE為當前這兩架飛機連線的長度（距離）圖12-1第i架飛機與第j架飛機的碰撞角因為飛機間的距離大于8km就不會相撞，所以這兩個圓隨著時間的推移不相交就可以了.為此，考慮第i架飛機相對于第j架飛機的相對速度（矢量，圖中記為QUOTE）是比較方便的，因為相對速度的大小和方向在飛機飛行中會始終保持不變（除非調整飛行角度）.設QUOTE為調整前的相對速度QUOTE與這兩架飛機連線（從i指向j的矢量）的夾角（以連線矢量為基準，逆時針方向為正，順時針方向為負），則QUOTE=-QUOTE，具體來說，QUOTE應該如下計算：QUOTE=相對速度QUOTE的輻角-從i指向j的連線矢量的輻角=QUOTE-QUOTE，QUOTE）. （18）注意：標準的反正切函數的符號是QUOTE，返回主值；我們這里使用QUOTE表示一個特殊的返回象限輻角的反正切函數，即QUOTE返回向量（a，b）的-QUOTE到QUOTE之間輻角（或者返回0到2QUOTE之間的輻角也是可以的）.即使這樣，也還不能完全滿足要求，因為這樣得到的QUOTE取值位于-2QUOTE到2QUOTE之間，還需要將它轉換到-QUOTE到QUOTE之間才行（超過QUOTE時就減去2QUOTE，小于-QUOTE就加上2QUOTE）.從圖中可以看出（注意圖中的兩條輔助線QUOTEn//CD、QUOTEm//AB），兩架飛機i，j不相撞的充要條件是（實際上不只是在所考慮的區(qū)域內不相交，而是永遠不會相交）|QUOTE|QUOTE. （19）如果調整前這個關系式成立，則不需要調整.否則，仍用QUOTE表示第i架飛機飛行方向角的調整量，并記由此引起的QUOTE的改變量為QUOTE.現(xiàn)在，問題的關鍵是如何弄清楚QUOTE如何隨QUOTE和QUOTE變化.可以證明QUOTE=（QUOTE+QUOTE）/2. （20）下面利用復數的知識證明式（20）證明由題知.設改變前的速度分別為,改變方向后速度分別為，改變前相對速度，改變后相對速度，即與輻角相差.因此，可以得到如下的數學規(guī)劃模型： （21） （22）

（23）這仍然是一個非線性規(guī)劃模型.同QUOTE一樣，這個模型中的QUOTE+的取值也需要轉換到-QUOTE到QUOTE之間才合理.通常情況下調整量很小，即（QUOTE+QUOTE）很小，因此只需要QUOTE位于-QUOTE到QUOTE之間就差不多了（除非QUOTE很接近-QUOTE和QUOTE，下面的表12-2顯示本題并非這種情況）.模型求解為了編寫LINGO程序求解式（21）QUOTE（23），必須解決如何用式（18）求QUOTE的問題，因為

LINGO中并沒有能返回-QUOTE到QUOTE之間的輻角的反正切函數QUOTE.如果一定要用LINGO求QUOTE，就需要很仔細地利用LINGO中正常的@tan函數，通過判斷每個點的位置，來正確得到這種關系，這是很不方便的，不是LINGO軟件的優(yōu)勢所在.所以最好使用其他軟件先計算QUOTE以后直接輸入LINGO.這里假設已經用其他方法（如MATLAB）計算得到了QUOTE的值，如表12-2所示（由于對稱性，只需要求出表中的一半元素的值）.表12-2其他方法計算得到了QUOTE的值單位：（°）ji1234561109.263642-128.25000024.179830173.06505114.4749342-88.871097-42.243563-92.3048479.000000312.476311-58.7862430.31080945.969234-3.52560651.9143836對于QUOTE，由式（17）知它的取值位于0到QUOTE/2之間，在反正切函數arcsin返回的角度的主值內，用LINGO計算也不麻煩，所以我們直接在LINGO中計算.于是，該飛機的數學規(guī)劃模型可如下輸入LINGO求解：MODEL:!飛行管理模型;SETS:Plane/1..6/:x0,y0,d_cita;!d_cita為調整的角度;link(plane,plane)|&1#LT#&2:alpha,beta;ENDSETSDATA:x0y0= 150 140 85 85 150 155 145 50 130 150 00;beta=109.263642-128.25000024.179830173.06505114.474934-88.871097-42.243563-92.3048479.00000012.476311-58.7862430.3108095.969234-3.5256061.914383;ENDDATA!計算alpha;@FOR(LINK(I,J):@SIN(alpha*3.14159265/180.0)=8/((X0(I)-X0(J))^2+(Y0(I)-Y0(J))^2)^.5);@for(link(i,j):@abs(beta(i,j)+0.5*d_cita(I)+0.5*d_cita(J))>alpha(I,J););@for(link:@bnd(0,alpha,90));@for(plane:@bnd(-30,d_cita,30));!min=@sum(plane:@sqr(d_cita));min=@sum(plane:@abs(d_cita));END計算結果如下（只顯示QUOTE和QUOTE的結果）：Linearizationcomponentsadded:Constraints:60Variables:60Integers:15Localoptimalsolutionfoundatiteration:575Objectivevalue:6.954676VariableValueReducedCostD_CITA(1)-0.2622117E-07-0.1776357E-07D_CITA(2)-0.249.247E-070.000000D_CITA(3)2.0624480.000000D_CITA(4)-0.49543750.000000D_CITA(5)-0.2482437E-070.000000D_CITA(6)1.5670110.000000ALPHA（1，2）5.3911900.000000ALPHA（1，3）32.230950.000000ALPHA（1，4）5.0918160.000000ALPHA（1，5）20.963360.000000ALPHA（1，6）2.2345070.000000ALPHA（2，3）4.8040240.000000ALPHA（2，4）6.6134600.000000ALPHA（2，5）5.8078660.000000

ALPHA（2，6）3.8159250.000000ALPHA（3，4）4.3646720.000000ALPHA（3，5）22.833650.000000ALPHA（3，6）2.1255390.000000ALPHA（4，5）4.5376920.000000ALPHA（4，6）2.9898190.000000ALPHA（5，6）2.3098410.000000這個結果與前面得到的結果幾乎是一樣的.注意上面顯示結果的最前面幾行，實際上是告訴我們LINGO對約束自動進行了線性化處理（“Linearizationcompon