建模案例-飞行管理问题

建模案例—飞行管理问题第1页/共62页一个飞行管理问题

1.1问题描述1995年全国大学生数学建模竞赛中的A题（“一个飞行管理问题”）。在约10000米高空的某边长为160km的正方形区域内，经常有若干架飞机做水平飞行，区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理，当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞。

第2页/共62页现假定条件如下：不碰撞的标准为任意两架飞机的距离大于8km；飞机飞行方向角调整幅度不应超过30°；所有飞机飞行速度均为800km／h；进入该区域的飞机在到达该区域边缘时，与区域内飞机的距离应在60km以上；最多需考虑6架飞机；不必考虑飞机离开此区域后的情况。第3页/共62页请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过0.01°），要求飞机飞行方向角调整的幅度尽量小。该区域四个定点的坐标为（0，0）、（160，0）、（160，160）、（0，160）。记录数据见表2－1。

第4页/共62页表2－1飞机位置和方向角记录数据飞机编号横坐标纵坐标方向角（°）飞机编号横坐标x纵坐标y方向角／（°）飞机编号横坐标纵坐标方向角（°）飞机编号横坐标x纵坐标y方向角／（°）11501402434145501592858523651301502303155155220．5新进入0052说明：方向角指飞行方向与x轴正向的夹角。试根据实际应用背景对你的模型进行评价和推广

第5页/共62页*对问题仔细阅读,首先抓住题目中的关键词“管理”进行联想.抓住诸如“碰撞”、“调整”、“避免碰撞”、“立即”、“判断”等等词语.

*联系解决问题的方案,不加约束继续联想，再将关键词搭配起来.1.问题的前期分析第6页/共62页VIIIIIIIVIV160km160km飞行位置示意图第7页/共62页立即判断

碰撞

条件

实时

算法

避免碰撞

调整方向角

实时

幅度尽量小相对距离优化问题优化算法优化调整方案第8页/共62页问题的初步理解和想法

飞行管理问题是优化问题,在调整方向角的幅度尽量小的同时，还必须注意调整方案及算法的实时性.第9页/共62页2.问题探究（1）优化问题的目标函数为何？方向角调整的尽量小方向角如何表示方向角的概念是什么方向角的平方和目标函数第10页/共62页任意两架飞机的距离大于8公里；飞机飞行方向角调整的幅度不应超过30°；（2）优化问题的约束条件为何？两点间距离表示方法判断避免碰撞的依据把飞机视为点总结：目标函数和约束条件中都含有方向角。第11页/共62页（3）分析6架飞机目前碰撞情况描点作图；分析飞机在飞行区域的时间；判断5架飞机此刻的情况；判断新进入飞机与其他5架飞机碰撞的情况。第12页/共62页（4）求解方法特殊到一般：先考虑2架，然考虑3架~优化问题为非线性规划问题，编程求解。总结：初等算法和高等算法都可。创新之处：算法。第13页/共62页模型一及求解模型建立这个问题显然是一个优化问题。设第i架飞机在调整时的方向角为（题目中已给出），调整后的方向为，题目中就是要求飞机飞行方向角调整的幅度尽量小，因此有化的目的函数可以是：

（1）第14页/共62页为了建立这个问题的优化模型，只需要明确约束条件就可以了。一个简单的约束是飞机飞行方向角调整的幅度不应超过30°，即

（2）

第15页/共62页题目中要求进入该区域的飞机在到达该区域边缘时，与区域内的飞机的距离应在60km以上。这个条件是个初始条件，很容易验证目前所给的数据是满足的，因此本模型中可以不予考虑。剩下的关键是要满足题目中描述的任意两架位于该区域内的飞机的距离应该大于8km。但这个问题的难点在于飞机是动态的，这个约束不好直接描述，为此我们首先需要描述每架飞机的飞行轨迹。第16页/共62页记飞机飞行速率为（800km／h），以当前时刻为0时刻。设第架飞机在调整时的位置坐标为（已知条件），时刻的位置坐标为，则（3）第17页/共62页如果要严格表示两架位于该区域内的飞机的距离应大于8km，则需要考虑每架飞机在区域内的飞行时间的长度。记Ti

为第架飞机飞出区域的时间，即

（4）

第18页/共62页记时刻第架飞机与第架飞机的距离为，并记，这时在区域内飞机不相撞的约束条件就变成了

（5）

其中（6）此外，经过计算可以得到第19页/共62页

（7）

（8）（9）

（10）

第20页/共62页所以是一个关于t的二次函数，表示的是一条开口向上的抛物线。当即（记为）时，函数取最小值。注意到（初始时刻不相撞），如果（即）则此时约束条件（5）一定成立，所以如果且，只要在右端点的函数值非负即可，即第21页/共62页

（11）如果且，只需要求最小值即可，即（12）实际上，约束（11）表示的是在右端点的函数值非负，这个约束在（12）的条件下也是自然成立的，所以可以是对约束（11）不再附加且的条件。

第22页/共62页于是我们的模型就是（13）

（14）

（15）（16）

第23页/共62页模型求解上面这是一个非线性规划模型，虽然是严格满足题目要求的模型，但得到的模型逻辑关系比较复杂，约束（16）是在一定条件下才成立的约束，而且其中的计算式（4）也含有相当复杂的关系式，使用LINGO软件不太容易将模型很方便的输入，因为逻辑处理不是LINGO的优势所在。即使想办法把这个模型输入到LINGO，也不一定能求出好的解（笔者尝试过，但是LINGO运行时有时会出现系统内部错误，可能是系统有问题，无法继续求解）。而且，在实时飞行调度中显然需要快速求解，所以下面我们想办法简化模型。第24页/共62页这个模型麻烦之处就在于，要求严格表示两架飞机的飞行距离应大于8km，所以需要考虑每架飞机在区域内的飞行时间的长度，比较繁琐。注意到区域对角线的长度只有，任何一架飞机在所考虑的区域内停留的时间不会超过。因此这里我们简化一下问题；

第25页/共62页不再单独考虑每架飞机在区域内停留的时间，而是以最大时间（这是已经是一个常数）代替之，此时所有，这实际上强化了问题的要求，即考虑了有些飞机可能已经飞出区域，但仍不允许两架飞机的距离小于8km。

第26页/共62页程序：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:第27页/共62页1501402438585236150155220.5145501591301502300052x0y0cita0=第28页/共62页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]

第29页/共62页(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第30页/共62页运行这个程序，结果得到的是一个局部极小点，调整角度较大。能找到更好的解吗？如果不用全局求解程序，通常很难得到稍大规模的非线性规划问题的全局最优解。所以我们启动LINGO全局求解程序求解这个模型（过程省略），可以得到全局最优解。可以看出，在0。01度的误差要求下，需要调整第3、4、6三架飞机的角度，分别调整2.06度，－0.5度,1.57度，调整量的平方和为6.95。其实，使用全局变量求解程序，通常也不一定要等到全局最优解，而是观察求解状态窗口，看到一个较好的当前解（或当前最好解在较长时间内不发生变化）时，就可以中止程序，用当前最好的局部最优解作为最后结果。例如对于本例，LINGO求出最优解大约需要一分，而实际上5秒内LINGO得到了与全局最优解类似的解。第31页/共62页此外，上面的模型还可以进一步简化，例如可以假设要求飞机永远不相撞，即认为为无穷大，这时显然约束条件（15）是多余的，而且约束（16）中只需要的条件就可以了。也就是说上面程序中的对应部分（约束和可以改写为更简单的形式：！右端点非负，不再需要；！最小点非负，需化为以下形式：实际计算结果显示，此时得到的结果与前面计算的结果几乎没有差别。第32页/共62页备注：优化的目标函数除了

外，也可以假定为

或

等，用LINGO求解的过程的是完全类似的，计算结果略有差异，这里就不再对这个目标函数进行计算了，甚至可以考虑让参与调整的飞机的数量尽量少，这种想法在实际中也不能说没有道理，但与题目的要求不符，而且解题难度并没有减小，意义似乎不大。在实际调度中，由于计算上面的调度方案，需要时间将调度信息告知飞机驾驶员并做出调整方向角的操作也需要时间，因此如果考虑一定的反应滞后时间应该是比较合理的。也就是说，如果反应时间是10秒，则计算中应采用飞机沿当前方向角飞行10秒以后的位置作为计算的基础。

第33页/共62页问题分析

根据题目的条件，可将飞机飞行的空域视为二维平面xoy中的一个正方形，顶点在(0，0),(160，0),(160，160),(0，160）。各架飞机的飞行方向角为飞行方向与x轴正向夹角（转角）。根据两飞机不碰撞的标准为二者距离大于8km，可将每架飞机视为一个以飞机为圆心、以4为半径的圆状物体（每架飞机在空域中的状态由圆心的位置矢量和飞行速度矢量确定）。这样两架飞机是否碰撞就化为两圆在运动中是否相交的问题。两圆是否相交只要讨论它们的相对运动即可。CABDnimilijαijαijαijβijvijij第34页/共62页建模时补充假定条件:

1.飞机在所定区域内作直线飞行，不偏离航向；

2.飞机管理系统内不发生意外，如发动机失灵，或其他意外原因迫使飞机改变航向；4.飞机管理系统发出的指令应被飞机立即执行，即认为转向是瞬间完成的（忽略飞机转向的影响，即转弯半径和转弯时间的影响）；3.飞机进入区域边缘时，立即作出计算，每架飞机按照计算后的指示立即作方向角改变；5.每架飞机在在整个过程中指点改变一次方向6.新飞机进入区域时，已在区域内部的飞机的飞行方向已调整合适，不会碰撞；7.对每架飞机方向角的相同调整量的满意程度是一样的。建模方案二：第35页/共62页模型的建立（1）圆状模型ABDnimilijαijαijαijβijvijij采用相对速度作为研究对象符号说明：i,j——第i，第j架飞机的圆心；αij——第i，第j架飞机的碰撞角，αij=αji;vij——第i架飞机相对第j架飞机的相对飞行速度

;lij——第i，第j架飞机的圆心距θi——第i架飞机的飞行方向与x轴正向夹角（逆时针为正）xi——第i架飞机的位置矢量vi——第i架飞机的的速度矢量βij——第i飞机对第j架飞机的相对速度与两架飞机圆心连线的夹角（逆时针为正）不碰撞|βij|>αij第36页/共62页(2)由圆状模型导出的方程讨论βij的改变量与第i第j两架飞机飞行方向角改变量Δθi,Δθj的关系由题目条件知|vi|=A=800，可用复数表示速度设第i，j飞机飞行方向改变前的速度分别为改变后的速度分别为改变前后相对速度分别为两者之商的幅角就是Δβij第37页/共62页定理：对第i，第j两架飞机，其相对速度方向βij的改变量Δβij等于两飞机飞行方向角改变量之和的一半，即第38页/共62页模型目标函数：Minθ其中θ为各飞机方向角调整量的最大值或为约束条件：调整方向角时不能超过300：调整飞行方向后飞机不能碰撞：模型为第39页/共62页化为线性规划模型由于Δθi可正可负，为使各变量均非负，引入新变量：模型化为第40页/共62页模型求解αij的计算model:sets:plane/1..6/:x0,y0;link(plane,plane):alpha,sin2;endsets@for(link(i,j)|i#ne#j:sin2(i,j)=64/((x0(i)-x0(j))^2+(y0(i)-y0(j))^2););@for(link(i,j)|i#ne#j:(@sin(alpha*3.14159265/180.0))^2=sin2;);data:x0=150,85,150,145,130,0;y0=140,85,155,50,150,0;enddataend第41页/共62页ALPHA(1,1)1.234568ALPHA(1,2)5.391190ALPHA(1,3)752.2310ALPHA(1,4)5.091816ALPHA(1,5)2000.963ALPHA(1,6)2.234507ALPHA(2,1)5.391190ALPHA(2,2)1.234568ALPHA(2,3)4.804024ALPHA(2,4)6.613460ALPHA(2,5)5.807866ALPHA(2,6)3.815925ALPHA(3,1)752.2310ALPHA(3,2)4.804024ALPHA(3,3)1.234568ALPHA(3,4)4.364672ALPHA(3,5)1102.834ALPHA(3,6)2.125539ALPHA(4,1)5.091816ALPHA(4,2)6.613460ALPHA(4,3)4.364672ALPHA(4,4)1.234568ALPHA(4,5)4.537692ALPHA(4,6)2.989819ALPHA(5,1)2000.963ALPHA(5,2)5.807866ALPHA(5,3)1102.834ALPHA(5,4)4.537692ALPHA(5,5)1.234568ALPHA(5,6)2.309841ALPHA(6,1)2.234507ALPHA(6,2)3.815925ALPHA(6,3)2.125539ALPHA(6,4)2.989819ALPHA(6,5)2.309841ALPHA(6,6)1.234568第42页/共62页αijJ=123456i=10.0000005.39119032.2309535.09181620.9633612.23450725.3911900.0000004.8040246.6134605.8078663.815925332.2309534.8040240.0000004.36467222.8336542.12553945.0918166.6134604.3646720.0000004.5376922.989819520.9633615.80786622.8336544.5376920.0000002.30984162.2345073.8159252.1255392.9898192.3098410.000000整理可得αij的值（单位角度）也可以用MATLAB计算αij的值第43页/共62页x=[150,85,150,145,130,0];y=[140,85,155,50,150,0];k=length(x);alpha=zeros(k);fori=1:kforj=1:kifi==jalpha(i,j)=0;elsealpha(i,j)=(180/3.14159265)*asin(8/sqrt((x(i)-x(j))^2+(y(i)-y(j))^2));endendendalpha计算αij的值程序为计算结果为第44页/共62页

05.3911902372235.391190237223032.2309526723314.8040239337975.0918164485506.61346048987220.9633608931285.8078662434212.2345067369953.81592477539932.2309526723315.0918164485504.8040239337986.61346048987204.3646718991114.364671899111022.8336542040094.5376924624022.1255388575512.98981913904520.9633608931282.2345067369955.8078662434213.81592477539922.8336542040092.1255388575514.5376924624032.98981913904502.3098413654052.3098413654050βij的计算:a=[150,85,150,145,130,0];b=[140,85,155,50,150,0];x=a+b*i;c=[243,236,220.5,159,230,52]*pi/180;v=exp(i*c);k=length(a);fori=1:k

forj=1:kbeita(i,j)=(angle(v(i)-v(j))-angle(x(j)-x(i)))*180/pi;

endendbeita用matlab程序编写第45页/共62页beita=0109.2636-128.250024.1798-186.934914.4749109.26360-88.8711-42.2436-92.30489.0000231.7500271.1289012.4763301.21380.310824.1798-42.243612.476305.9692-3.5256173.0651267.6952-58.78625.969201.914414.47499.00000.3108-3.52561.91440运算结果最优解的计算用LINGO求解程序如下第46页/共62页model:sets:plane/1..6/:cita;link(plane,plane):alpha,beta;endsetsmin=@sum(plane:@abs(cita));@for(plane(i):@bnd(-30,cita(i),30););@for(link(i,j)|i#ne#j:@abs(beta(i,j)+0.5*cita(i)+0.5*cita(j))>alpha(i,j););第47页/共62页data:alpha=0.000000,5.391190,32.230953,5.091816,20.963361,2.234507,5.391190,0.000000,4.8040024,6.813460,5.807866,3.815925,32.230953,4.804024,0.000000,4.364672,22.833654,2.125539,5.091816,6.613460,4.363673,0.000000,4.537692,2.989819,20.963361,5.807866,22.833654,4.537692,0.000000,2.309841,2.234507,3.815925,2.125539,2.989819,2.309841,0.000000;beta=0.000000109.263642-128.25000024.179830173.06505113.474934109.2636420.000000-88.871096-42.243563-92.3048479.000000-128.250000-88.870960.00000012.476311-58.7862430.31080924.179830-42.24356312.4763110.0000005.969234-3.525606174.065051-92.304846-58.7862445.9692340.0000001.91438314.4749349.0000000.310809-3.5256061.9133830.000000;enddataend用MATLAB计算编程如下第48页/共62页function[f,g]=plane(x)alph=[0.000000,5.391190,32.230953,5.091816,20.963361,2.234507,5.391190,0.000000,4.8040024,6.813460,5.807866,3.815925,32.230953,4.804024,0.000000,4.364672,22.833654,2.125539,5.091816,6.613460,4.363673,0.000000,4.537692,2.989819,20.963361,5.807866,22.833654,4.537692,0.000000,2.309841,2.234507,3.815925,2.125539,2.989819,2.309841,0.000000];bet=[0.000000109.263642-128.25000024.179830173.06505113.474934109.2636420.000000-88.871096-42.243563-92.3048479.000000-128.250000-88.870960.00000012.476311-58.7862430.31080924.179830-42.24356312.4763110.0000005.969234-3.525606174.065051-92.304846-58.7862445.9692340.0000001.91438314.4749349.0000000.310809-3.5256061.9133830.000000];f=abs(x(1))+abs(x(2))+abs(x(3))+abs(x(4))+abs(x(5))+abs(x(6));g(1)=alpha(1,2)-abs(beta(1,2)+0.5*x(1)+0.5*x(2));g(2)=alpha(1,3)-abs(beta(1,3)+0.5*x(1)+0.5*x(3));g(3)=alpha(1,4)-abs(beta(1,4)+0.5*x(1)+0.5*x(4));g(4)=alpha(1,5)-abs(beta(1,5)+0.5*x(1)+0.5*x(5));第49页/共62页g(5)=alpha(1,6)-abs(beta(1,6)+0.5*x(1)+0.5*x(6));g(6)=alpha(2,3)-abs(beta(2,3)+0.5*x(2)+0.5*x(3));g(7)=alpha(2,4)-abs(beta(2,4)+0.5*x(2)+0.5*x(4));g(8)=alpha(2,5)-abs(beta(2,5)+0.5*x(2)+0.5*x(5));g(9)=alpha(2,6)-abs(beta(2,6)+0.5*x(2)+0.5*x(6));g(10)=alpha(3,4)-abs(beta(3,4)+0.5*x(3)+0.5*x(4));g(11)=alpha(3,5)-abs(beta(3,5)+0.5*x(3)+0.5*x(5));g(12)=alpha(3,6)-abs(beta(3,6)+0.5*x(3)+0.5*x(6));g(13)=alpha(4,5)-abs(beta(4,5)+0.5*x(4)+0.5*x(5));g(14)=alpha(4,6)-abs(beta(4,6)+0.5*x(4)+0.5*x(6));g(15)=alpha(5,6)-abs(beta(5,6)+0.5*x(5)+0.5*x(6));执行程序x0=[0,0,0,0,0,0];v1=-30*ones(1,6);v2=30*ones(1,6);opt=[];x=constr('plane',x0,opt,v1,v2)第50页/共62页结果：x=-0.00000576637983-0.000005766379832.58794980234726-0.000012434879850.000036204730951.04151019765274最优解：第51页/共62页模型检验各飞行方向按此方案调整后，系统各架飞机均满足|βij+(ΔθI+Δθj)/2|>αij,结果是正确的。模型的评价与推广（1）此模型采用圆状模型分析碰撞问题是合理的，同时采用相对速度作为判别标准，既体现了碰撞的本质（相对运动），又简化了模型的计算；（2）建模中用了适当的简化，将一个复杂的非线性规划问题简化为线性规划问题，既求到合理的解，又提高了运算速度，这对解决高速飞行的飞机碰撞问题是十分重要的。此模型对题目所提供的例子计算得出的结果是令人满意的。（3）由对称性知模型中的约束个数是（n是飞机数），所有约束条件数是，计算量增加不大。第52页/共62页建模方案二：1.问题分析目的：不碰撞手段：调整飞行方向角要求：调整的幅度尽量小求解思路：(1)找出不碰撞的条件,可用matlab画出各架飞机在t=0时刻的位置和飞行方向(2)求调整幅度的极小值