12.1一個(gè)飛行管理問(wèn)題

在中國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽（Chinaundergraduatemathematicalcontestinmodeling,CUMCM）中，曾經(jīng)出現(xiàn)過(guò)大量的優(yōu)化建模賽題.本章從中選擇了部分典型賽題，舉例分析其優(yōu)化建模過(guò)程，說(shuō)明如何應(yīng)用LINDO/LINGO軟件包求解這些賽題.12.1一個(gè)飛行管理問(wèn)題12.1.1問(wèn)題描述1995年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽中的A題（“一個(gè)飛行管理問(wèn)題”）.在約10000m高空的某邊長(zhǎng)為160km的正方形區(qū)域內(nèi)，經(jīng)常有若干架飛機(jī)作水平飛行.區(qū)域內(nèi)每架飛機(jī)的位置和速度向量均由計(jì)算機(jī)記錄其數(shù)據(jù)，以便進(jìn)行飛行管理.當(dāng)一架欲進(jìn)入該區(qū)域的飛機(jī)到達(dá)區(qū)域邊緣時(shí)，記錄其數(shù)據(jù)后，要立即計(jì)算并判斷是否會(huì)與區(qū)域內(nèi)的飛機(jī)發(fā)生碰撞.如果會(huì)碰撞，則應(yīng)計(jì)算如何調(diào)整各架（包括新進(jìn)入的）飛機(jī)飛行的方向角，以避免碰撞.現(xiàn)假定條件如下：不碰撞的標(biāo)準(zhǔn)為任意兩架飛機(jī)的距離大于8km；飛機(jī)飛行方向角調(diào)整的幅度不應(yīng)超過(guò)30°；所有飛機(jī)飛行速度均為800km/h；進(jìn)入該區(qū)域的飛機(jī)在到達(dá)該區(qū)域邊緣時(shí)，與該區(qū)域內(nèi)的飛機(jī)的距離應(yīng)在60km以上；最多需考慮6架飛機(jī)；不必考慮飛機(jī)離開(kāi)此區(qū)域后的狀況.請(qǐng)你對(duì)這個(gè)避免碰撞的飛行管理問(wèn)題建立數(shù)學(xué)模型，列出計(jì)算步驟，對(duì)以下數(shù)據(jù)進(jìn)行計(jì)算（方向角誤差不超過(guò)0.01°），要求飛機(jī)飛行方向角調(diào)整的幅度盡量小.設(shè)該區(qū)域4個(gè)頂點(diǎn)的坐標(biāo)為（0，0），（160，0），（160，160），（0，160）.記錄數(shù)據(jù)見(jiàn)表12-1.表12-1飛機(jī)位置和方向角記錄數(shù)據(jù)飛機(jī)編號(hào)橫坐標(biāo)縱坐標(biāo)方向角飛機(jī)編號(hào)橫坐標(biāo)縱坐標(biāo)方向角11501402434145501592858523651301502303150155220.5新進(jìn)入0052說(shuō)明：方向角指飛行方向與ｘ軸正向的夾角.試根據(jù)實(shí)際應(yīng)用背景對(duì)你的模型進(jìn)行評(píng)價(jià)和推廣.12.1.2模型1及求解模型建立這個(gè)問(wèn)題顯然是一個(gè)優(yōu)化問(wèn)題.設(shè)第ｉ架飛機(jī)在調(diào)整時(shí)的方向角為QUOTE（題目中已經(jīng)給出），調(diào)整后的方向角為QUOTE=QUOTE+QUOTE（QUOTE=1,2,…，6）.題目中就是要求飛機(jī)飛行方向角調(diào)整的幅度盡量小，因此優(yōu)化的目標(biāo)函數(shù)可以是QUOTE. （1）為了建立這個(gè)問(wèn)題的優(yōu)化模型，只須要明確約束條件就可以了.一個(gè)簡(jiǎn)單的約束是飛機(jī)飛行方向角調(diào)整的不應(yīng)超過(guò)30°，即|QUOTE|QUOTE30°. (2)題中要求進(jìn)入該區(qū)域的飛機(jī)在到達(dá)該區(qū)域邊緣時(shí)，與該區(qū)域內(nèi)飛機(jī)的距離應(yīng)在60km以上，這個(gè)條件是個(gè)初始條件，很容易驗(yàn)證目前所給數(shù)據(jù)是滿足的，因此本模型中可以不予考慮.剩下的關(guān)鍵是要滿足題目中描述的任意兩架飛機(jī)不碰撞的要求，即任意兩架位于該區(qū)域內(nèi)的飛機(jī)的距離應(yīng)大于8km.但這個(gè)問(wèn)題的難點(diǎn)在于飛機(jī)是動(dòng)態(tài)的，這個(gè)約束不好直接描述，為此我們首先需要描述每架飛機(jī)的飛行軌跡.記飛機(jī)飛行速率為v（=800km/h）,以當(dāng)前時(shí)刻為0時(shí)刻.設(shè)第ｉ架飛機(jī)在調(diào)整時(shí)的位置坐標(biāo)為（QUOTE，QUOTE）（已知條件），t時(shí)刻的位置坐標(biāo)為（QUOTE，QUOTE），則QUOTE=QUOTE+QUOTE，QUOTE=QUOTE+QUOTE.如果要嚴(yán)格表示兩架位于該區(qū)域內(nèi)的飛機(jī)的距離應(yīng)大于8km，則需考慮每架飛機(jī)在該區(qū)域內(nèi)的飛行時(shí)間的長(zhǎng)度.記QUOTE為第ｉ架飛機(jī)飛出去與的時(shí)刻，即QUOTE=QUOTEargmin{t>0:QUOTE+QUOTE=0或160， （4）或者QUOTE+QUOTE=0或160}. （5）記t時(shí)刻第ｉ架飛機(jī)與第ｊ架飛機(jī)的距離為QUOTE（t），并記QUOTE（t）=QUOTE-64，這時(shí)在該區(qū)域內(nèi)飛機(jī)不相撞的約束條件就變成了QUOTE（t）=QUOTE-64QUOTE0（0QUOTEtQUOTE）.其中QUOTE=min{QUOTE，QUOTE}. (6)此外，經(jīng)過(guò)計(jì)算，可以得到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）是一個(gè)關(guān)于QUOTE的二次函數(shù)，表示的是一條開(kāi)口向上的拋物線.當(dāng)QUOTE=-QUOTE/2，即t=-QUOTE/4vQUOTE（記為QUOTE）時(shí)，函數(shù)QUOTE（t）取最小值-QUOTE/4+QUOTE.注意到QUOTE（0）=QUOTE0（初始時(shí)刻不相撞）,如果-QUOTE/2QUOTE0（即QUOTE0），則此時(shí)約束條件（5）一定成立，所以對(duì)約束條件（5）只需考慮以下兩種可能情況：如果QUOTE0且QUOTE，只需要QUOTE（t）在右端點(diǎn)QUOTE的函數(shù)值非負(fù)即可，即QUOTE（QUOTE）QUOTE0； （11）如果QUOTE0且0QUOTE，只需要QUOTE（t）的最小值QUOTE（QUOTE）=-QUOTE/4+QUOTE0即可，即QUOTE4QUOTE0. (12)實(shí)際上，約束（11）表示的是QUOTE（t）在右端點(diǎn)QUOTE的函數(shù)值非負(fù)，這個(gè)約束在（12）的條件下也是自然成立的，所以可以對(duì)約束（11）不再附加QUOTE0且QUOTE的條件.于是，我們的模型就是QUOTE; （13）S．t．|QUOTE|QUOTE30°， （14）QUOTE（QUOTE）QUOTE0 (15)QUOTE4QUOTE0.（當(dāng)QUOTE0且0QUOTE）. (16)模型求解上面這是一個(gè)非線性規(guī)劃模型，雖然是嚴(yán)格滿足題目要求的模型，但得到的模型邏輯關(guān)系比較復(fù)雜，約束（16）是在一定條件下才成立的約束，而且其中QUOTE的計(jì)算式（4）也是含有相當(dāng)復(fù)雜的關(guān)系式，使用LINGO軟件不太容易將模型很方便地輸入，因?yàn)檫壿嬏幚聿皇荓INGO的優(yōu)勢(shì)所在.即使能想辦法把把這個(gè)模型輸入到LINGO,也不一定能求出好的解(筆者嘗試過(guò),但LINGO運(yùn)行時(shí)有時(shí)會(huì)出現(xiàn)系統(tǒng)錯(cuò)誤,可能是系統(tǒng)有漏洞,無(wú)法繼續(xù)求解).而且,在實(shí)時(shí)飛行調(diào)度中顯然需要快速求解,所以下面我們想辦法簡(jiǎn)化模型.這個(gè)模型麻煩之處就在于,要求嚴(yán)格表示兩架位于該區(qū)域內(nèi)的飛機(jī)距離應(yīng)大于8km,所以需要考慮每架飛機(jī)在該區(qū)域內(nèi)的飛行時(shí)間,比較繁瑣.注意到區(qū)域?qū)蔷€的長(zhǎng)度只有160QUOTEkm，任何一架飛機(jī)在所考慮的范圍內(nèi)停留的時(shí)間不會(huì)超過(guò)QUOTE=160QUOTE/800=0.2QUOTE0.283（h），因此這里我們簡(jiǎn)化一下問(wèn)題：不再單獨(dú)考慮每架飛機(jī)在該區(qū)域停留的時(shí)間QUOTE，而以最大時(shí)間（這里已經(jīng)是一個(gè)常數(shù)）QUOTE代替之，此時(shí)所有QUOTE=QUOTE.這實(shí)際上強(qiáng)化了問(wèn)題的要求,即考慮了有些飛機(jī)可能已經(jīng)飛出該區(qū)域,但仍不允許兩架飛機(jī)的距離小于8km.這個(gè)簡(jiǎn)化的模型可以如下輸入LINGO軟件：MODEL:TITLE飛行管理問(wèn)題的非線性規(guī)劃模型;SETS:Plane/1..6/:x0,y0,cita0,cita1,d_cita;!cita0表示初始角度，cita1為調(diào)整后的角度,d_cita為調(diào)整的角度;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;);!避免碰撞的條件;!右端點(diǎn)非負(fù);@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);!最小點(diǎn)非負(fù);@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));!調(diào)整角度上下限，單位為角度;@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只接受弧度，所以程序中一律進(jìn)行了轉(zhuǎn)換.求解這個(gè)模型,得到Localoptimalsolutionfound.Objectivevalue：295.4937ModelTitle:飛行管理問(wèn)題的非線性規(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這個(gè)結(jié)果得到的是一個(gè)局部極小點(diǎn)，調(diào)整角度較大.能找到更好的解嗎?如果不用全局求解程序,通常很難得到稍大規(guī)模的非線性規(guī)劃問(wèn)題的全局最優(yōu)解.所以我們啟動(dòng)LINGO全局求解程序求解這個(gè)模型,可以得到全局最優(yōu)解如下:Globaloptimalsolutionfoundatiteration:93Objectivevalue:6.953944ModelTitle:飛行管理問(wèn)題的非線性規(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°的誤差要求下，需要調(diào)整第3，4，6三架飛機(jī)的角度，分別調(diào)整2.06°，-0.50°，1.57°.調(diào)整量的平方和為6.95.其實(shí),使用全局求解程序,通常也不一定要等到得到全局最優(yōu)解,而是觀察求解狀態(tài)窗口,看到一個(gè)較好的當(dāng)前解(或當(dāng)前最好解在較長(zhǎng)時(shí)間內(nèi)不發(fā)生變化)時(shí),就可以終止程序,用當(dāng)前最好的局部最優(yōu)解作為最后的惡結(jié)果.列如,對(duì)于本列,LINGO求出全局最優(yōu)解大約需要1min,而實(shí)際上5s內(nèi)LINGO就得到了與全局最優(yōu)解類似的解.此外,上面的模型還可以進(jìn)一步簡(jiǎn)化,列如可以假設(shè)要求飛機(jī)永遠(yuǎn)不相撞,即認(rèn)為QUOTE為無(wú)限大，這時(shí)顯然約束（15）也是多余的,而且約束(16)中只需要QUOTE0的條件就可以了.也就是說(shuō)，上面的程序中的對(duì)應(yīng)部分（約束[Right]和[Minimum]）可以改寫(xiě)為更簡(jiǎn)單的形式：！有端點(diǎn)非負(fù)，不再需要；！最小點(diǎn)非負(fù)，簡(jiǎn)化為以下形式；@for(link(i,j):@if(b（i，j）#lt#0，b（i，j）^2-4*c（i，j），-1）QUOTE0）；實(shí)際計(jì)算顯示，此時(shí)得到的結(jié)果與前面計(jì)算的結(jié)果幾乎沒(méi)有差別.備注優(yōu)化的目標(biāo)函數(shù)除了QUOTE外，也可以設(shè)定為QUOTE或QUOTE等，用LINGO求解的過(guò)程是完全類似的，計(jì)算結(jié)果略有差異，這里就不再對(duì)兩個(gè)目標(biāo)函數(shù)具體計(jì)算了.甚至可以考慮讓參與調(diào)整的飛機(jī)的數(shù)量盡量小，這種想法在實(shí)際中也不能說(shuō)沒(méi)有道理，但與題目中的要求不符，而且解題難度并沒(méi)有減小，意義似乎不大.實(shí)際調(diào)度中，由于計(jì)算上面的調(diào)度方案需要時(shí)間，將調(diào)度信息告訴飛機(jī)駕駛員并作出調(diào)整方向角的操作也需要時(shí)間，因此如果考慮一定的反應(yīng)滯后時(shí)間，應(yīng)該是比較合理的.也就是說(shuō)，如果反應(yīng)時(shí)間是10s，則計(jì)算式中應(yīng)采用飛機(jī)沿當(dāng)前的方向角飛行10s以后的位置作為計(jì)算的基礎(chǔ).12.1.3模型2及求解從12.1.2節(jié)可以看出，求解模型1的非線性規(guī)劃模型是比較困難的，輸入后也很可能找不到好的解甚至出現(xiàn)錯(cuò)誤.此外，演示版軟件還會(huì)受到求解規(guī)模的限制，尤其可能無(wú)權(quán)使用全局求解程序.因此，如果能把這個(gè)問(wèn)題簡(jiǎn)化成比較簡(jiǎn)單的規(guī)劃模型，將是非常有價(jià)值的.模型建立如圖12-1，把兩架飛機(jī)i，j分別看成半徑為4km的圓（圖中i，j為圓心），AB,CD為公切線，將AB和CD的夾角的一半稱為碰撞角.在調(diào)整時(shí)刻，第i架飛機(jī)與第j架飛機(jī)的碰撞角為QUOTE，則易知QUOTE=QUOTE）， （17）其中QUOTE為當(dāng)前這兩架飛機(jī)連線的長(zhǎng)度（距離）圖12-1第i架飛機(jī)與第j架飛機(jī)的碰撞角因?yàn)轱w機(jī)間的距離大于8km就不會(huì)相撞，所以這兩個(gè)圓隨著時(shí)間的推移不相交就可以了.為此，考慮第i架飛機(jī)相對(duì)于第j架飛機(jī)的相對(duì)速度（矢量，圖中記為QUOTE）是比較方便的，因?yàn)橄鄬?duì)速度的大小和方向在飛機(jī)飛行中會(huì)始終保持不變（除非調(diào)整飛行角度）.設(shè)QUOTE為調(diào)整前的相對(duì)速度QUOTE與這兩架飛機(jī)連線（從i指向j的矢量）的夾角（以連線矢量為基準(zhǔn)，逆時(shí)針?lè)较驗(yàn)檎?，順時(shí)針?lè)较驗(yàn)樨?fù)），則QUOTE=-QUOTE，具體來(lái)說(shuō)，QUOTE應(yīng)該如下計(jì)算：QUOTE=相對(duì)速度QUOTE的輻角-從i指向j的連線矢量的輻角=QUOTE-QUOTE，QUOTE）. （18）注意：標(biāo)準(zhǔn)的反正切函數(shù)的符號(hào)是QUOTE，返回主值；我們這里使用QUOTE表示一個(gè)特殊的返回象限輻角的反正切函數(shù)，即QUOTE返回向量（a，b）的-QUOTE到QUOTE之間輻角（或者返回0到2QUOTE之間的輻角也是可以的）.即使這樣，也還不能完全滿足要求，因?yàn)檫@樣得到的QUOTE取值位于-2QUOTE到2QUOTE之間，還需要將它轉(zhuǎn)換到-QUOTE到QUOTE之間才行（超過(guò)QUOTE時(shí)就減去2QUOTE，小于-QUOTE就加上2QUOTE）.從圖中可以看出（注意圖中的兩條輔助線QUOTEn//CD、QUOTEm//AB），兩架飛機(jī)i，j不相撞的充要條件是（實(shí)際上不只是在所考慮的區(qū)域內(nèi)不相交，而是永遠(yuǎn)不會(huì)相交）|QUOTE|QUOTE. （19）如果調(diào)整前這個(gè)關(guān)系式成立，則不需要調(diào)整.否則，仍用QUOTE表示第i架飛機(jī)飛行方向角的調(diào)整量，并記由此引起的QUOTE的改變量為QUOTE.現(xiàn)在，問(wèn)題的關(guān)鍵是如何弄清楚QUOTE如何隨QUOTE和QUOTE變化.可以證明QUOTE=（QUOTE+QUOTE）/2. （20）下面利用復(fù)數(shù)的知識(shí)證明式（20）證明由題知.設(shè)改變前的速度分別為,改變方向后速度分別為，改變前相對(duì)速度，改變后相對(duì)速度，即與輻角相差.因此，可以得到如下的數(shù)學(xué)規(guī)劃模型： （21） （22）

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

LINGO中并沒(méi)有能返回-QUOTE到QUOTE之間的輻角的反正切函數(shù)QUOTE.如果一定要用LINGO求QUOTE，就需要很仔細(xì)地利用LINGO中正常的@tan函數(shù)，通過(guò)判斷每個(gè)點(diǎn)的位置，來(lái)正確得到這種關(guān)系，這是很不方便的，不是LINGO軟件的優(yōu)勢(shì)所在.所以最好使用其他軟件先計(jì)算QUOTE以后直接輸入LINGO.這里假設(shè)已經(jīng)用其他方法（如MATLAB）計(jì)算得到了QUOTE的值，如表12-2所示（由于對(duì)稱性，只需要求出表中的一半元素的值）.表12-2其他方法計(jì)算得到了QUOTE的值單位：（°）ji1234561109.263642-128.25000024.179830173.06505114.4749342-88.871097-42.243563-92.3048479.000000312.476311-58.7862430.31080945.969234-3.52560651.9143836對(duì)于QUOTE，由式（17）知它的取值位于0到QUOTE/2之間，在反正切函數(shù)arcsin返回的角度的主值內(nèi)，用LINGO計(jì)算也不麻煩，所以我們直接在LINGO中計(jì)算.于是，該飛機(jī)的數(shù)學(xué)規(guī)劃模型可如下輸入LINGO求解：MODEL:!飛行管理模型;SETS:Plane/1..6/:x0,y0,d_cita;!d_cita為調(diào)整的角度;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!計(jì)算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計(jì)算結(jié)果如下（只顯示QUOTE和QUOTE的結(jié)果）：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這個(gè)結(jié)果與前面得到的結(jié)果幾乎是一樣的.注意上面顯示結(jié)果的最前面幾行，實(shí)際上是告訴我們LINGO對(duì)約束自動(dòng)進(jìn)行了線性化處理（“Linearizationcompon