12.1一个飞行管理问题

在中国大学生数学建模竞赛（China undergraduate mathematical contest in modeling,CUMCM）中，曾经出现过大量的优化建模赛题.本章从中选择了部分典型赛题，举例分析其优化建模过程，说明如何应用LINDO/LINGO软件包求解这些赛题.12.1 一个飞行管理问题12.1.1 问题描述1995年全国大学生数学建模竞赛中的A题（“一个飞行管理问题”）.在约10000m高空的某边长为160km的正方形区域内，经常有若干架飞机作水平飞行.区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理.当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞.如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞.现假定条件如下：（1）不碰撞的标准为任意两架飞机的距离大于8km；（2）飞机飞行方向角调整的幅度不应超过30°；（3）所有飞机飞行速度均为800km/h；（4）进入该区域的飞机在到达该区域边缘时，与该区域内的飞机的距离应在60km以上；（5）最多需考虑6架飞机；（6）不必考虑飞机离开此区域后的状况.请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过0.01°），要求飞机飞行方向角调整的幅度尽量小.设该区域4个顶点的坐标为（0，0），（160，0），（160，160），（0，160）.记录数据见表12-1.试根据实际应用背景对你的模型进行评价和推广.12.1.2 模型1及求解模型建立这个问题显然是一个优化问题.设第ｉ架飞机在调整时的方向角为（题目中已经给出），调整后的方向角为=+（=1,2,…，6）.题目中就是要求飞机飞行方向角调整的幅度尽量小，因此优化的目标函数可以是. （1）为了建立这个问题的优化模型，只须要明确约束条件就可以了.一个简单的约束是飞机飞行方向角调整的不应超过30°，即||30°. (2)题中要求进入该区域的飞机在到达该区域边缘时，与该区域内飞机的距离应在60km以上，这个条件是个初始条件，很容易验证目前所给数据是满足的，因此本模型中可以不予考虑.剩下的关键是要满足题目中描述的任意两架飞机不碰撞的要求，即任意两架位于该区域内的飞机的距离应大于8km.但这个问题的难点在于飞机是动态的，这个约束不好直接描述，为此我们首先需要描述每架飞机的飞行轨迹.记飞机飞行速率为v（=800km/h）,以当前时刻为0时刻.设第ｉ架飞机在调整时的位置坐标为（，）（已知条件），t时刻的位置坐标为（，），则=+，=+.如果要严格表示两架位于该区域内的飞机的距离应大于8km，则需考虑每架飞机在该区域内的飞行时间的长度.记为第ｉ架飞机飞出去与的时刻，即=argmin{t>0:+=0或160，（4）或者+=0或160}. （5）记t时刻第ｉ架飞机与第ｊ架飞机的距离为（t），并记（t）=-64，这时在该区域内飞机不相撞的约束条件就变成了（t）=-640（0t）.其中=min{，}. (6)此外，经过计算，可以得到（t）=+-64=++，（7）其中=2vt，（8）=2[-(-)+, （9）=+-64. （10）所以，（t）是一个关于的二次函数，表示的是一条开口向上的抛物线.当=-/2，即t=-/4v（记为）时，函数（t）取最小值-/4+.注意到（0）=0（初始时刻不相撞）,如果-/20（即0），则此时约束条件（5）一定成立，所以对约束条件（5）只需考虑以下两种可能情况：如果0且，只需要（t）在右端点的函数值非负即可，即（）0；（11）如果0且0，只需要（t）的最小值（）=-/4+0即可，即40. (12)实际上，约束（11）表示的是（t）在右端点的函数值非负，这个约束在（12）的条件下也是自然成立的，所以可以对约束（11）不再附加0且的条件.于是，我们的模型就是; （13）S．t．||30°，（14）（）0 (15)40.（当0且0）. (16)模型求解上面这是一个非线性规划模型，虽然是严格满足题目要求的模型，但得到的模型逻辑关系比较复杂，约束（16）是在一定条件下才成立的约束，而且其中的计算式（4）也是含有相当复杂的关系式，使用LINGO软件不太容易将模型很方便地输入，因为逻辑处理不是LINGO的优势所在.即使能想办法把把这个模型输入到LINGO,也不一定能求出好的解(笔者尝试过,但LINGO运行时有时会出现系统错误,可能是系统有漏洞,无法继续求解).而且,在实时飞行调度中显然需要快速求解,所以下面我们想办法简化模型.这个模型麻烦之处就在于,要求严格表示两架位于该区域内的飞机距离应大于8km,所以需要考虑每架飞机在该区域内的飞行时间,比较繁琐.注意到区域对角线的长度只有160km，任何一架飞机在所考虑的范围内停留的时间不会超过=160/800=0.20.283（h），因此这里我们简化一下问题：不再单独考虑每架飞机在该区域停留的时间，而以最大时间（这里已经是一个常数）代替之，此时所有=.这实际上强化了问题的要求,即考虑了有些飞机可能已经飞出该区域,但仍不允许两架飞机的距离小于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:;max_cita = 30;T_max = 0.283;V=800;ENDDATAINIT:d_cita = 0 0 0 0 0 0;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只接受弧度，所以程序中一律进行了转换.求解这个模型,得到Local optimal solution found.Objective value： 295.4937Model Title: 飞行管理问题的非线性规划模型Variable Value Reduced CostD_CITA(1) -10.5.980 0.000000D_CITA(2) 0.000000 0.000000D_CITA(3) 0.000000 0.000000D_CITA(4) 6.515425 0.000000D_CITA(5) 10.00681 0.000000D_CITA(6) 6.505425 0.000000这个结果得到的是一个局部极小点，调整角度较大.能找到更好的解吗?如果不用全局求解程序,通常很难得到稍大规模的非线性规划问题的全局最优解.所以我们启动LINGO全局求解程序求解这个模型,可以得到全局最优解如下: Global optimal solution found at iteration: 93Objective value: 6.953944Model Title: 飞行管理问题的非线性规划模型Variable ValueD_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(这里只给出的值)可以看到，在0.01°的误差要求下，需要调整第3，4，6三架飞机的角度，分别调整2.06°，-0.50°，1.57°.调整量的平方和为6.95.其实,使用全局求解程序,通常也不一定要等到得到全局最优解,而是观察求解状态窗口,看到一个较好的当前解(或当前最好解在较长时间内不发生变化)时,就可以终止程序,用当前最好的局部最优解作为最后的恶结果.列如,对于本列,LINGO求出全局最优解大约需要1min,而实际上5s内LINGO就得到了与全局最优解类似的解.此外,上面的模型还可以进一步简化,列如可以假设要求飞机永远不相撞,即认为为无限大，这时显然约束（15）也是多余的,而且约束(16)中只需要0的条件就可以了.也就是说，上面的程序中的对应部分（约束[Right]和[Minimum]）可以改写为更简单的形式：！有端点非负，不再需要；！最小点非负，简化为以下形式；@for(link(i,j): @if(b（i，j）#lt#0，b（i，j）^2-4* c（i，j），-1）0）；实际计算显示，此时得到的结果与前面计算的结果几乎没有差别.备注优化的目标函数除了外，也可以设定为或等，用LINGO求解的过程是完全类似的，计算结果略有差异，这里就不再对两个目标函数具体计算了.甚至可以考虑让参与调整的飞机的数量尽量小，这种想法在实际中也不能说没有道理，但与题目中的要求不符，而且解题难度并没有减小，意义似乎不大.实际调度中，由于计算上面的调度方案需要时间，将调度信息告诉飞机驾驶员并作出调整方向角的操作也需要时间，因此如果考虑一定的反应滞后时间，应该是比较合理的.也就是说，如果反应时间是10s，则计算式中应采用飞机沿当前的方向角飞行10s以后的位置作为计算的基础.12.1.3 模型2及求解从12.1.2节可以看出，求解模型1的非线性规划模型是比较困难的，输入后也很可能找不到好的解甚至出现错误.此外，演示版软件还会受到求解规模的限制，尤其可能无权使用全局求解程序.因此，如果能把这个问题简化成比较简单的规划模型，将是非常有价值的.模型建立如图12-1，把两架飞机i，j分别看成半径为4km的圆（图中i，j为圆心），AB,CD为公切线，将AB和CD的夹角的一半称为碰撞角.在调整时刻，第i架飞机与第j架飞机的碰撞角为，则易知=），（17）其中为当前这两架飞机连线的长度（距离）图12-1 第i架飞机与第j架飞机的碰撞角因为飞机间的距离大于8km就不会相撞，所以这两个圆随着时间的推移不相交就可以了.为此，考虑第i架飞机相对于第j架飞机的相对速度（矢量，图中记为）是比较方便的，因为相对速度的大小和方向在飞机飞行中会始终保持不变（除非调整飞行角度）.设为调整前的相对速度与这两架飞机连线（从i指向j的矢量）的夹角（以连线矢量为基准，逆时针方向为正，顺时针方向为负），则=-，具体来说，应该如下计算：=相对速度的辐角-从i指向j的连线矢量的辐角=-，）. （18）注意：标准的反正切函数的符号是，返回主值；我们这里使用表示一个特殊的返回象限辐角的反正切函数，即返回向量（a ，b ）的-到之间辐角（或者返回0到2之间的辐角也是可以的）.即使这样，也还不能完全满足要求，因为这样得到的取值位于-2到2之间，还需要将它转换到-到之间才行（超过时就减去2，小于-就加上2）.从图中可以看出（注意图中的两条辅助线n//CD 、m//AB ），两架飞机i ，j 不相撞的充要条件是（实际上不只是在所考虑的区域内不相交，而是永远不会相交）||.（19）如果调整前这个关系式成立，则不需要调整.否则，仍用表示第i 架飞机飞行方向角的调整量，并记由此引起的的改变量为.现在，问题的关键是如何弄清楚如何随和变化.可以证明=（+）/2.（20）下面利用复数的知识证明式（20）证明 由题知800i v Km A ==.设改变前的速度分别为'',ji i i i j v Ae v Ae θθ==,改变方向后速度分别为()2i ii i v Aeθθ+∆=，()2j ji j v Aeθθ+∆=改变前相对速度''j ii i ij i j v v v A e e θθ⎡⎤=-=-⎣⎦，改变后相对速度''j ii i ij i j v v v A e e θθ⎡⎤=-=-⎣⎦()()'22j jii i i ij i j v v v A e eθθθθ+∆+∆⎡⎤=-=-⎢⎥⎣⎦，()()'j j i i jii i iji i ijA e evv A e e θθθθθθ+∆+∆⎡⎤-⎢⎥⎣⎦=⎡⎤-⎣⎦()()()()sin sin cos sin cos sin cos sin i i i ijjjji i j ji i i i θθθθθθθθθθθθ+∆++∆-+∆-+∆=+--()2sin sin cos 2222sin sin cos 222i i j j i i j j i i j j i j i j i j i i θθθθθθθθθθθθθθθθθθ+∆-+∆+∆-+∆+∆++∆⎛⎫- ⎪⎝⎭=-++⎛⎫- ⎪⎝⎭()2sin2sin2i j iijji i jeθθθθθθθθ∆+∆⎛⎫⎪ ⎪⎝⎭+∆-+∆=-即'ij v 与ij v 辐角相差2i jθθ∆+∆.因此，可以得到如下的数学规划模型：（21）()1.. ,2ij i j ij s t βθθα+∆+∆≥（22）30,1,6i i θ∆≤=…,（23）这仍然是一个非线性规划模型.同一样，这个模型中的+()12i j θθ∆+∆的取值也需要转换到-到之间才合理.通常情况下调整量很小，即（+）很小，因此只需要位于-到之间就差不多了（除非很接近-和，下面的表12-2显示本题并非这种情况）. 模型求解为了编写LINGO 程序求解式（21）（23），必须解决如何用式（18）求的问题，因为21min ;nii θ=∑,1,6,,i j i j =≠…,LINGO中并没有能返回-到之间的辐角的反正切函数.如果一定要用LINGO 求，就需要很仔细地利用LINGO中正常的@tan函数，通过判断每个点的位置，来正确得到这种关系，这是很不方便的，不是LINGO软件的优势所在.所以最好使用其他软件先计算以后直接输入LINGO.这里假设已经用其他方法（如MATLAB ）计算得到了的值，如表12-2所示（由于对称性，只需要求出表中的一半元素的值）.表12-2 其他方法计算得到了的值单位：（°）对于，由式（17）知它的取值位于0到/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:x0 y0 =150 14085 85150 155145 50130 1500 0 ;beta= 109.263642 -128.250000 24.179830 173.065051 14.474934-88.871097 -42.243563 -92.304847 9.00000012.476311 -58.786243 0.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计算结果如下（只显示和的结果）：Linearization components added:Constraints: 60Variables: 60Integers: 15Local optimal solution found at iteration: 575Objective value: 6.954676Variable Value Reduced Cost D_CITA(1) -0.2622117E-07 -0.1776357E-07 D_CITA(2) -0.249.247E-07 0.000000 D_CITA(3) 2.062448 0.000000 D_CITA(4) -0.4954375 0.000000 D_CITA(5) -0.2482437E-07 0.000000 D_CITA(6) 1.567011 0.000000 ALPHA（1，2） 5.391190 0.000000 ALPHA（1，3） 32.23095 0.000000 ALPHA（1，4） 5.091816 0.000000 ALPHA（1，5） 20.96336 0.000000 ALPHA（1，6） 2.234507 0.000000 ALPHA（2，3） 4.804024 0.000000 ALPHA（2，4） 6.613460 0.000000 ALPHA（2，5） 5.807866 0.000000ALPHA（2，6） 3.815925 0.000000ALPHA（3，4） 4.364672 0.000000ALPHA（3，5） 22.83365 0.000000ALPHA（3，6） 2.125539 0.000000ALPHA（4，5） 4.537692 0.000000ALPHA（4，6） 2.989819 0.000000ALPHA（5，6） 2.309841 0.000000这个结果与前面得到的结果几乎是一样的.注意上面显示结果的最前面几行，实际上是告诉我们LINGO对约束自动进行了线性化处理（“Linearization components added“），这是通过加入15个整数变量做到的（“Integers：15”）.可见，对一些可以线性化的约束或目标（如这里的约束是变量线性函数的绝对值的形式的情行），LINGO具有自动线性化的功能，以便找到更好的解.应此，这时我们可以不用自己亲自对模型进行线性化（有时这是一件很困难的事情）.事实上，我们看到此时不使用全局求解程序，就很容易得到了很好的解（不过由于目标还是非线性的,所以LINGO仍然只是报告了局部极小点）.备注如果目标函数也采用绝对值和的形式，即|,则LINGO就能够自动实现整个模型线性化了.这只需将上面LINGO程序中的目标函数改写为：min=@sum（plane：@abs（d_cita））;求解得到的显示为Linearization components added:Constraints: 84Variables: 84Integers: 21Global optimal solution found.Objective value: 3.629460Extended solver steps 7Variable Value Reduced CostD_CITA(1) 0.000000 0.000000D_CITA(2) 0.000000 0.000000D_CITA(3) 2.557886 0.000000D_CITA(4) 0.000000 0.000000D_CITA(5) 0.000000 0.000000D_CITA(6) 1.071574 0.000000由此可知最优解为 2.56°， 1.07°（其他调整角度为0）.此时LINGO线性化时引入了21个整数解变量.由于转化后完全是（整数）线性规划模型，因此直接就可以得到全局最优解（不需要使用全局最优求解程序）.需要指出的是，这个模型中的和+（+）的取值需要转换到-到之间才合理，对于一般情形的飞机初始位置，可能会有出现错误的时候，所以最好对最后求到的解进行一次可行性检验.T = 0 0.1503 0.2401 0.2275 0.2401 0.24010.1503 0 0.1503 0.1503 0.1503 0.15030.2401 0.1503 0 0.2275 0.2481 0.26960.2275 0.1503 0.2275 0 0.2275 0.22750.2401 0.1503 0.2481 0.2275 0 0.24810.2401 0.1503 0.2696 0.2275 0.2481 0 ;。