12.1一个飞行管理问题

在中国大学生数学建模竞赛（Chinaundergraduatemathematicalcontestinmodeling,CUMCM）中，曾经出现过大量的优化建模赛题.本章从中选择了部分典型赛题，举例分析其优化建模过程，说明如何应用LINDO/LINGO软件包求解这些赛题.12.1一个飞行管理问题12.1.1问题描述1995年全国大学生数学建模竞赛中的A题（“一个飞行管理问题”）.在约10000m高空的某边长为160km的正方形区域内，经常有若干架飞机作水平飞行.区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理.当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞.如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞.现假定条件如下：不碰撞的标准为任意两架飞机的距离大于8km；飞机飞行方向角调整的幅度不应超过30°；所有飞机飞行速度均为800km/h；进入该区域的飞机在到达该区域边缘时，与该区域内的飞机的距离应在60km以上；最多需考虑6架飞机；不必考虑飞机离开此区域后的状况.请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过0.01°），要求飞机飞行方向角调整的幅度尽量小.设该区域4个顶点的坐标为（0，0），（160，0），（160，160），（0，160）.记录数据见表12-1.表12-1飞机位置和方向角记录数据飞机编号横坐标纵坐标方向角飞机编号横坐标纵坐标方向角11501402434145501592858523651301502303150155220.5新进入0052说明：方向角指飞行方向与ｘ轴正向的夹角.试根据实际应用背景对你的模型进行评价和推广.12.1.2模型1及求解模型建立这个问题显然是一个优化问题.设第ｉ架飞机在调整时的方向角为QUOTE（题目中已经给出），调整后的方向角为QUOTE=QUOTE+QUOTE（QUOTE=1,2,…，6）.题目中就是要求飞机飞行方向角调整的幅度尽量小，因此优化的目标函数可以是QUOTE. （1）为了建立这个问题的优化模型，只须要明确约束条件就可以了.一个简单的约束是飞机飞行方向角调整的不应超过30°，即|QUOTE|QUOTE30°. (2)题中要求进入该区域的飞机在到达该区域边缘时，与该区域内飞机的距离应在60km以上，这个条件是个初始条件，很容易验证目前所给数据是满足的，因此本模型中可以不予考虑.剩下的关键是要满足题目中描述的任意两架飞机不碰撞的要求，即任意两架位于该区域内的飞机的距离应大于8km.但这个问题的难点在于飞机是动态的，这个约束不好直接描述，为此我们首先需要描述每架飞机的飞行轨迹.记飞机飞行速率为v（=800km/h）,以当前时刻为0时刻.设第ｉ架飞机在调整时的位置坐标为（QUOTE，QUOTE）（已知条件），t时刻的位置坐标为（QUOTE，QUOTE），则QUOTE=QUOTE+QUOTE，QUOTE=QUOTE+QUOTE.如果要严格表示两架位于该区域内的飞机的距离应大于8km，则需考虑每架飞机在该区域内的飞行时间的长度.记QUOTE为第ｉ架飞机飞出去与的时刻，即QUOTE=QUOTEargmin{t>0:QUOTE+QUOTE=0或160， （4）或者QUOTE+QUOTE=0或160}. （5）记t时刻第ｉ架飞机与第ｊ架飞机的距离为QUOTE（t），并记QUOTE（t）=QUOTE-64，这时在该区域内飞机不相撞的约束条件就变成了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)模型求解上面这是一个非线性规划模型，虽然是严格满足题目要求的模型，但得到的模型逻辑关系比较复杂，约束（16）是在一定条件下才成立的约束，而且其中QUOTE的计算式（4）也是含有相当复杂的关系式，使用LINGO软件不太容易将模型很方便地输入，因为逻辑处理不是LINGO的优势所在.即使能想办法把把这个模型输入到LINGO,也不一定能求出好的解(笔者尝试过,但LINGO运行时有时会出现系统错误,可能是系统有漏洞,无法继续求解).而且,在实时飞行调度中显然需要快速求解,所以下面我们想办法简化模型.这个模型麻烦之处就在于,要求严格表示两架位于该区域内的飞机距离应大于8km,所以需要考虑每架飞机在该区域内的飞行时间,比较繁琐.注意到区域对角线的长度只有160QUOTEkm，任何一架飞机在所考虑的范围内停留的时间不会超过QUOTE=160QUOTE/800=0.2QUOTE0.283（h），因此这里我们简化一下问题：不再单独考虑每架飞机在该区域停留的时间QUOTE，而以最大时间（这里已经是一个常数）QUOTE代替之，此时所有QUOTE=QUOTE.这实际上强化了问题的要求,即考虑了有些飞机可能已经飞出该区域,但仍不允许两架飞机的距离小于8km.这个简化的模型可以如下输入LINGO软件：MODEL:TITLE飞行管理问题的非线性规划模型;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:飞行管理问题的非线性规划模型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这个结果得到的是一个局部极小点，调整角度较大.能找到更好的解吗?如果不用全局求解程序,通常很难得到稍大规模的非线性规划问题的全局最优解.所以我们启动LINGO全局求解程序求解这个模型,可以得到全局最优解如下:Globaloptimalsolutionfoundatiteration:93Objectivevalue:6.953944ModelTitle:飞行管理问题的非线性规划模型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.其实,使用全局求解程序,通常也不一定要等到得到全局最优解,而是观察求解状态窗口,看到一个较好的当前解(或当前最好解在较长时间内不发生变化)时,就可以终止程序,用当前最好的局部最优解作为最后的恶结果.列如,对于本列,LINGO求出全局最优解大约需要1min,而实际上5s内LINGO就得到了与全局最优解类似的解.此外,上面的模型还可以进一步简化,列如可以假设要求飞机永远不相撞,即认为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）；实际计算显示，此时得到的结果与前面计算的结果几乎没有差别.备注优化的目标函数除了QUOTE外，也可以设定为QUOTE或QUOTE等，用LINGO求解的过程是完全类似的，计算结果略有差异，这里就不再对两个目标函数具体计算了.甚至可以考虑让参与调整的飞机的数量尽量小，这种想法在实际中也不能说没有道理，但与题目中的要求不符，而且解题难度并没有减小，意义似乎不大.实际调度中，由于计算上面的调度方案需要时间，将调度信息告诉飞机驾驶员并作出调整方向角的操作也需要时间，因此如果考虑一定的反应滞后时间，应该是比较合理的.也就是说，如果反应时间是10s，则计算式中应采用飞机沿当前的方向角飞行10s以后的位置作为计算的基础.12.1.3模型2及求解从12.1.2节可以看出，求解模型1的非线性规划模型是比较困难的，输入后也很可能找不到好的解甚至出现错误.此外，演示版软件还会受到求解规模的限制，尤其可能无权使用全局求解程序.因此，如果能把这个问题简化成比较简单的规划模型，将是非常有价值的.模型建立如图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不相撞的充要条件是（实际上不只是在所考虑的区域内不相交，而是永远不会相交）|QUOTE|QUOTE. （19）如果调整前这个关系式成立，则不需要调整.否则，仍用QUOTE表示第i架飞机飞行方向角的调整量，并记由此引起的QUOTE的改变量为QUOTE.现在，问题的关键是如何弄清楚QUOTE如何随QUOTE和QUOTE变化.可以证明QUOTE=（QUOTE+QUOTE）/2. （20）下面利用复数的知识证明式（20）证明由题知.设改变前的速度分别为,改变方向后速度分别为，改变前相对速度，改变后相对速度，即与辐角相差.因此，可以得到如下的数学规划模型： （21） （22）

（23）这仍然是一个非线性规划模型.同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软件的优势所在.所以最好使用其他软件先计算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中计算.于是，该飞机的数学规划模型可如下输入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