11一个飞行管理问题

1、在中国大学生数学建模竞赛（ Chi na un dergraduate mathematical con test in modeling,CUMCM）中，曾经出现过大量的优化建模赛题.本章从中选择了部分典 型赛题，举例分析其优化建模过程，说明如何应用LINDO/LINGO软件包求解这些 赛题.12.1个飞行管理问题12.1.1 问题描述1995年全国大学生数学建模竞赛中的 A题（“一个飞行管理问题”）.在约10000m高空的某边长为160km的正方形区域内，经常有若干架飞机作 水平飞行.区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行 飞行管理.当一架欲进入该区域的飞机到达区域边

2、缘时，记录其数据后，要立即 计算并判断是否会与区域内的飞机发生碰撞如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞.现假定条件如下：(1)不碰撞的标准为任意两架飞机的距离大于8km；(2)飞机飞行方向角调整的幅度不应超过30 °(3)所有飞机飞行速度均为800km/h；(4)进入该区域的飞机在到达该区域边缘时，在60km以上；与该区域内的飞机的距离应(5)最多需考虑6架飞机；(6)不必考虑飞机离开此区域后的状况.请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以 下数据进行计算（方向角误差不超过 0.01 °），要求飞机飞行方向角调

3、整的幅度 尽量小.设该区域4个顶点的坐标为（0, 0），（ 160, 0），（ 160，160），（0，160） 记录数据见表12-1.表12-1飞机位置和方向角记录数据飞机编号横坐标x纵坐标y方向角/（°）飞机编号横坐标x纵坐标y方向角/（°）11501402434145501592r 858523651301502303150155220.5新进入0052说明：方向角指飞行方向与x轴正向的夹角试根据实际应用背景对你的模型进行评价和推广.12.1.2 模型1及求解模型建立这个问题显然是一个优化问题.设第1架飞机在调整时的方向角为:（题目中 已经给出），调整后的方向角为 =

4、+蜒（=1,2,，6）.题目中就是要求飞机飞 行方向角调整的幅度尽量小，因此优化的目标函数可以是轨时.为了建立这个问题的优化模型，只须要明确约束条件就可以了.一个简单的约束是飞机飞行方向角调整的不应超过 30°,即II _30° .题中要求进入该区域的飞机在到达该区域边缘时，与该区域内飞机的距离应 在60km以上，这个条件是个初始条件，很容易验证目前所给数据是满足的，因 此本模型中可以不予考虑.剩下的关键是要满足题目中描述的任意两架飞机不碰 撞的要求，即任意两架位于该区域内的飞机的距离应大于8km.但这个问题的难点在于飞机是动态的，这个约束不好直接描述，为此我们首先需要描述

5、每架飞机 的飞行轨迹.记飞机飞行速率为v （=800km/h）,以当前时刻为0时刻.设第i架飞机在调 整时的位置坐标为（聲，yf）（已知条件），t时刻的位置坐标为（£,时），贝U-= +说cos壯歼=弓+谢血备如果要严格表示两架位于该区域内的飞机的距离应大于8km则需考虑每架(4)飞机在该区域内的飞行时间的长度.记 为第1架飞机飞出去与的时刻，即审=argmint>0:埠 +=0 或 160,或者即+龜韵L层=0或160.(5)记t时刻第1架飞机与第J架飞机的距离为Hj（ t ），并记& （t） =r*t）-64， 这时在该区域内飞机不相撞的约束条件就变成了希（t） =

6、 - . -640 （0t - 一）.其中=mi n ，.此外，经过计算，可以得到（t）=半罰脛幽皆肉曲曙其中（jf-僦sin血-y卜起和片）*-64(7)(8)t 9-6j=2vt,(9)(10)划=2卜(一二)：二_+ 二；二一,2 亠 2*凉椚+ -64.V i 1所以，一 （t）是一个关于 的二次函数，表示的是一条开口向上的抛物线9；-0；y当坷二陶/2，即t=-掩/4v：;j （记为-）时，函数長（t ）取最小值-再/4+乌. 注意到 （0） =q >0 （初始时刻不相撞）,如果-匕卜/20 （即 > 0），则此时约束条件（5） 定成立，所以对约束条件（5）只需考虑以下两

7、种可能情况：如果対处0且-一，只需要一 （t ）在右端点的函数值非负即可，即1| J(11)如果対试0且0，只需要1（ t ）的最小值（-）=-/4+0即可,即墉40.（12）实际上，约束（11）表示的是二（t）在右端点一的函数值非负，这个约束在（12）的条件下也是自然成立的，所以可以对约束（11）不再附加0且，:.的条件.于是，我们的模型就是；（13）(14)(15)ft (屬)_040.(当叮；0 且 0)(16)模型求解上面这是一个非线性规划模型，虽然是严格满足题目要求的模型，但得到的模型逻辑关系比较复杂，约束（16）是在一定条件下才成立的约束，而且其中- 的计算式（4）也是含有相当复杂

8、的关系式，使用 LINGC软件不太容易将模型很 方便地输入，因为逻辑处理不是 LINGO的优势所在.即使能想办法把把这个模型 输入到LING0也不一定能求出好的解（笔者尝试过,但LINGO运行时有时会出现 系统错误,可能是系统有漏洞，无法继续求解）.而且,在实时飞行调度中显然需要 快速求解,所以下面我们想办法简化模型这个模型麻烦之处就在于，要求严格表示两架位于该区域内的飞机距离应大 于8km,所以需要考虑每架飞机在该区域内的飞行时间，比较繁琐.注意到区域对角线的长度只有160 km任何一架飞机在所考虑的范围内停留的时间不会超过.=160 /800=0.2施0.283（h），因此这里我们简化一下

9、问题：不再单独考虑每架飞机在该区域停留的时间，而以最大时间（这里已经是一个常数）；代替之，此时所有一=.这实际上强化了问题的要求，即考虑了有些飞机可能已经飞出该区域,但仍不允许两架飞机的距离小于8km.这个简化的模型可以如下输入LINGC软件：MODEL:TITLE 飞行管理问题的非线性规划模型；SETS:Plane/1.6/: x0, y0, citaO, cita1,d_cita ;! citaO表示初始角度，cita1为调整后的角度,d_cita 为调整的角度；lin k（pla ne, pla ne）| &1 #LT# &2: b,c ;ENDSETSDATA:x0 y

10、0 cita0 =1501402438585236150155220.5145501591301502300052max_cita = 30;T_max = 0.283;V=800;ENDDATAINIT :d_cita = 0 0 0 0 0 0;ENDINIT0 = d_cita);for ( plane: cita1 - citafor ( 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.1415926

11、5/360 );c(i,j) = (xO(i) -xO(j) A2 + (yO(i) -yO(j) A2 - 64;);! 避免碰撞的条件 ;! 右端点非负 ;for(link(i,j): Right( 2*V*T_max* sin(cita1(i)-cita1(j)+ b(i,j)*( 2*V*T_max* sin(cita1(i)-cita1(j)*3.14159265/360 )A 2*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(c

12、ita1(i)-cita1(j)#and#-b(i,j)/4/V/sin(cita1(i)-cita1(j)#lt # T_max ,b(i,j)A2-4*c(i,j),-1) < 0);!for(link (i,j): if(b(i,j)#lt#0, b(i,j)A2-4*c(i,j), -1) < 0for (link:free (b);! 调整角度上下限，单位为角度 ;for ( plane:bnd( - max_cita, d_ cita, max_cita ) ;obj MIN = SUM(plane: ( d_cita)A 2);!obj MIN = SUM(plane

13、: abs( d_cita);END*3.14159265/360 ) # gt #0*3.14159265/360 );注意：上面的模型中的方向角单位一律用角度,程序中一律进行了转换 .求解这个模型 , 得到Local optimal solution found.Objective valueModel Title:LINGO只接受弧度,所以：295.4937飞行管理问题的非线性规划模型这个结果得到的不用全局求解程序 , 通常很难得到稍大规模的非线性规划问题的全局最优解 以我们启动LINGO全局求解程序求解这个模型,可以得到全局最优解如下：个局部极小点,调整角度较大. 能找到更好的解吗?如

14、果.所VariableValueReduced CostD_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.000000Global optimal solution found at iteration： 93Objective value:6.953944Model Title:飞行管理问题的非线性规划模型Variable ValueD_CITA(1) 0

15、.2719480E-02D_CITA(2) 0.5613433E-02D_CITA (3)2.059140D_CITA -0.4985421D_CITA(5)-0.5407837E-03D_CITA(6)1.570129(这里只给出磁*的值)可以看到，在0.01。的误差要求下，需要调整第 3, 4, 6三架飞机的角度， 分别调整2.06 °，-0.50 °，1.57 ° .调整量的平方和为6.95.其实,使用全局求解程序，通常也不一定要等到得到全局最优解，而是观察求 解状态窗口，看到一个较好的当前解(或当前最好解在较长时间内不发生变化)时, 就可以终止程序，用当前

16、最好的局部最优解作为最后的恶结果.列如,对于本 列,LINGO求出全局最优解大约需要1min,而实际上5s内LINGO就得到了与全局 最优解类似的解.此外,上面的模型还可以进一步简化，列如可以假设要求飞机永远不相撞，即 认为二为无限大，这时显然约束(15)也是多余的，而且约束(16)中只需要0的条件就可以了 .也就是说，上面的程序中的对应部分(约束Right和 Minimum)可以改写为更简单的形式：!有端点非负，不再需要；!最小点非负，简化为以下形式；for(link(i,j):if(b( i，j)#lt#0，b( i，j)A2-4* c( i，j)，-1)0)；实际计算显示，此时得到的结果

17、与前面计算的结果几乎没有差别.备注优化的目标函数除了 .U計外，也可以设定为鳶：.：閃或贡喺1彌讣遞寸等，用LINGO求解的过程是完全类似的，计算结果略有差异，这里就不再对两个目标函数具体计算了 .甚至可以考虑让参与调整的飞机的数量尽量 小，这种想法在实际中也不能说没有道理，但与题目中的要求不符，而且解题难度并没有减小，意义似乎不大.实际调度中，由于计算上面的调度方案需要时间，将调度信息告诉飞机驾驶 员并作出调整方向角的操作也需要时间，因此如果考虑一定的反应滞后时间，应该是比较合理的.也就是说，如果反应时间是10s，则计算式中应采用飞机沿当 前的方向角飞行10s以后的位置作为计算的基础.12.

18、1.3 模型2及求解从12.1.2节可以看出，求解模型1的非线性规划模型是比较困难的，输入 后也很可能找不到好的解甚至出现错误.此外，演示版软件还会受到求解规模的 限制，尤其可能无权使用全局求解程序.因此，如果能把这个问题简化成比较简单的规划模型，将是非常有价值的模型建立如图12-1，把两架飞机i ，j分别看成半径为4km的圆（图中i , j为圆心）， AB,CD为公切线，将AB和CD的夹角的一半称为碰撞角.在调整时刻，第i架飞 机与第j架飞机的碰撞角为ijj，则易知-Qu1其中.为当前这两架飞机连线的长度（距离）图12-1第i架飞机与第j架飞机的碰撞角因为飞机间的距离大于8km就不会相撞，所

19、以这两个圆随着时间的推移不相 交就可以了 为此，考虑第i架飞机相对于第j架飞机的相对速度（矢量，图中记为対j）是比较方便的，因为相对速度的大小和方向在飞机飞行中会始终保持不 变（除非调整飞行角度）.设为调整前的相对速度 与这两架飞机连线（从i 指向j的矢量）的夹角（以连线矢量为基准，逆时针方向为正，顺时针方向为负）， 则朝=-£:，具体来说，应该如下计算:月=相对速度逹的辐角-从i指向j的连线矢量的辐角(18)arctan； l注意：标准的反正切函数的符号是，返回主值；我们这里使用：-表示一个特殊的返回象限辐角的反正切函数，即 以比匕飞返回向量（a, b）的-i: 到 之间辐角（或者

20、返回0到2之间的辐角也是可以的）.即使这样，也还不能 完全满足要求，因为这样得到的取值位于-2到2之间，还需要将它转换到-到一之间才行（超过一时就减去2，小于-一就加上2一）.从图中可以看出（注意图中的两条辅助线jn/CD、m/AB）,两架飞机i ，j 不相撞的充要条件是（实际上不只是在所考虑的区域内不相交， 而是永远不会相 交）I為 I .（ 19）如果调整前这个关系式成立，则不需要调整否则，仍用城匕表示第i架飞机 飞行方向角的调整量，并记由此引起的 的改变量为邮沪现在，问题的关键是如 何弄清楚:巒如何随團和盘j变化可以证明巒=（幽+ :）/2.（20）下面利用复数的知识证明式（20）证明

21、由题知y =800Km = A设改变前的速度分别为 v； = AeHvj= A® 賂七i.cos 1 i sin 1 cos '. i sin ijjq (巧 +A0j)r 色 _8j +A&jq +ieeej)2sinsini cos 一2(2 2 丿日j e. 了 色 +e. q +e.、j，改变万向后速度分别为v； = Ae '（日钿）v2 =Ae® 颁）改变前相对速度vij v' _ v j 二 A ” e"寸 - e | 二二 1 ?改变后相对速度Vij = Vi - Vj = A 聲-e旧Vj = Vi2sin . si

22、n - i cos - v2 = A 严2+朋）-e® 呵）1VjA)eT6 +AOj )1VijA严-e卍sin 二 -n i j 亠 i sin i - i cos V j - . '： t i sin V , - . : V ,222sin 宀 ide2.战一日j sin 2即vj与Vj辐角相差，"j因此，可以得到如下的数学规划模型:min-i2(21)s.tPij十1(心日i + A日j 工5i, j1,6(22)（23）这仍然是一个非线性规划模型| - 30,1.同 一样，这个模型中的+- .-: j的取值也需要转换到-到之间才合理通常情况下调整量很小，即

23、（蛰+）很小，因 此只需要 位于-到之间就差不多了（除非 很接近-和，下面的表12-2显 示本题并非这种情况）.模型求解为了编写LINGO程序求解式（21）和（23），必须解决如何用式（18）求陥的 问题，因为LINGO中并没有能返回-一到一之间的辐角的反正切函数.如果一定要用LINGO求，就需要很仔细地利用LINGO中正常的tan函数，通过判断每个点的 位置，来正确得到这种关系，这是很不方便的，不是 LINGO软件的优势所在.所 以最好使用其他软件先计算 以后直接输入LINGO这里假设已经用其他方法（如MATLAB计算得到了 的值，如表12-2所示（由于对称性，只需要求出表中的 一半元素的值

24、）表12-2 其他方法计算得到了 ：的值单位：（°）jji1234561109.263642-128.25000024.179830173.06505114.4749342-88.871097-42.243563-92.3048479.000000312.476311-58.7862430.31080945.969234-3.52560651.9143836对于陶，由式(17)知它的取值位于 0到/2之间，在反正切函数arcsin返回的角度的主值内，用LINGO计算也不麻烦，所以我们直接在LINGO中计算.于是，该飞机的数学规划模型可如下输入LINGO求解：MODEL:!飞行管理模型

25、；SETS:Plane/1.6/: x0, y0,d_cita ;! d_cita为调整的角度；lin k(pla ne, pla ne)| &1 #LT# &2: alpha, beta ;ENDSETSDATA:x0 y0 =150 1408585150 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.

26、5256061.914383;ENDDATA!计算 alpha;FORLINK(I,J):SIN(alpha*3.14159265/180.0)=8 / ( (X0(I)-X0(J)A2 +(Y0(I)-Y0(J)A2 )A.5 );cita(J)for (link(i,j):bs(beta (i,j)+0.5* d_cita(l)+0.5*d> alpha(I,J););for (link:bnd(0,alpha,90);for (plane:bnd(-30,d_cita,30) );!min=sum(pla ne: sqr(d_cita);min=sur(pla ne:ab$d_ci

27、ta);END计算结果如下(只显示喊*和附的结果)：Lin earizati on comp onents added:Con strai nts:60Variables:60In tegers:15Local optimal soluti on found at iterati on: 575Objective value:6.954676VariableValueReduced CostD_CITA(1)-0.2622117E-07-0.1776357E-07D_CITA (2)-0.249.247E-070.000000D_CITA (3)2.0624480.000000D_CITA -

28、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.000000ALPHA(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