一个飞行管理问题数模竞赛

1、 一个飞行管理问题 摘要 在某一空域里对飞机的飞行合理管理事关重大，比如乘客及机上工作人员生命财产安全和航空公司的运作效益等。本文通过对飞机飞行管理问题的研究，得到了调整飞机架数较少同时调整幅度均最小（平方和最小）的飞行管理最优安排的非线性模型，这样既使得乘客所受影响达到最少，也便于飞机调整，还有利于飞机回到原来的航线，同时还在决策时间上对模型进行了优化和调整。本文不仅一般性地将不相撞的问题转化为欧式距离控制，而且很巧妙的将不碰撞条件转化成简单的二次函数标准形式进行含参讨论，建立一个只含有转向角变量的模型。并且再次很妙的具体化区域受控时间形成矩阵，大大得简化运算，节约了大量运算的时间，便于管理

2、人员控制操作，从而确保飞机的安全。更重要的是最后结合实际缩短了搜索区间，并优化算法，使得决策更加高效。最后的延时检验也充分体现了模型的可靠性。关键字：欧氏距离 约束转化 缩短搜索区间 时间矩阵 延时检验一、问题重述在约10000米的高空某边长为160公里的正方形区域，经常有若干架飞机作水平飞行。区域每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。当一驾欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞。现假定条件如下： 1）不碰撞的标准为任意两架飞机的距离大于8

3、公里； 2）飞机飞行方向角调整的幅度不应超过30度； 3）所有飞机飞行速度均为每小时800公里； 4）进入该区域的飞机在到达区域边缘时，与区域飞机的距离应在60公里以上； 5）最多需考虑6架飞机； 6）不必考虑飞机离开此区域后的状况。 请你对这个避免碰撞的飞机管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过0.01度），要求飞机飞行方向角调整的幅度尽量小。 设该区域4个顶点的坐标为(0,0),(160,0),(160,160),(0,160)。 记录数据为：飞机编号横坐标x纵坐标y方向角（度）115040243285852363150155220.54145501595

4、130150230新进入0052    注：方向角指飞行方向与x轴方向的夹角。 试根据实际应用背景对你的模型进行评价与推广。2、 问题分析1初步分析 根据问题容易知道，这显然是一个优化问题，当两架飞机可能发生碰撞时，即在规定区域某一时刻两架飞机之间的距离小于8公里，因此要调整飞行方向一定角度，保证任意两架飞机在区域任意时刻，两者的距离均不小于8公里，避免相撞。考虑到调整角度应尽量小，可以简化飞行方向调整策略，降低调整难度，同时减轻机乘客及工作人员的不适。此外由此初步确定了调整目标：所有六架飞机的飞行方向调整角度均尽量小。 2解决方案 由于所有飞机均处于1

5、000米得高空作水平飞行，可将飞机飞行的空域视为二维平面xoy中的一个正方形区域,顶点为(0,0),(160,0),(160,160),(0,160)。于是可以引入时间变量后，确定每架飞机在任意时刻的坐标，列出任意两点的欧氏距离，令其恒大于8公里，则得出一个重要约束条件。再结合变化角度应小于30度，即可得出约束条件，然后运用LINGO软件编辑程序进行求解。为提高决策效率，在反复试验中又可对约束条件进行调整。三、条件假设1. 不碰撞的标准为任意两架飞机的距离在以后任何一个时间里大于8公里；2. 飞机飞行方向角调整的幅度不应超过30度； 3. 所有飞机飞行速度均为每小时800公里； 4. 进入该区

6、域的飞机在到达区域边缘时，与区域飞机的距离应在60公里以上。即在计算如何最优地调整各架（包括新进入的）飞机飞行的方向角时，飞行管理中心得出合理的最优调整措施，；5. 最多需考虑6架飞机。 6. 此处忽略飞机在执行过程中所需耗费的时间，即假设从飞机管理中心发出的调整信息飞机马上可以接收并执行，不存在滞后或延迟；7. 飞行管理中心在计算飞行调整信号和发出信号所需时间，忽略各架飞机（包括刚进入的飞机）调整航向前飞行数据的变化；8. 假定飞机在该区域完全依赖飞行管理中心调度；9. 假设飞机在飞出区域之后，飞行员可以自觉调整飞行策略，回归原始航线，即飞行管理中心不必考虑飞机离开此区域后的状况。四、符号说

7、明 符号含义第i架飞机在初始时刻的横坐标第j架飞机在初始时刻的横坐标第i架飞机在t时刻的横坐标第j架飞机在t时刻的横坐标第i架飞机在初始时刻的纵坐标第j架飞机在初始时刻的纵坐标第i架飞机在t时刻的纵坐标第j架飞机在t时刻的纵坐标第i架飞机在初始时刻飞行方向与X轴正向的夹角第j架飞机在初始时刻飞行方向与X轴正向的夹角第i架飞机在t时刻飞行方向与X轴正向的夹角第j架飞机在t时刻飞行方向与X轴正向的夹角第i架飞机飞行方向角的调整幅度第j架飞机飞行方向角的调整幅度第i架飞机在规定区域可能飞行的最长时间第j架飞机在规定区域可能飞行的最长时间T是一个6*6矩阵；=min，V飞机的飞行速度飞机i与飞机j的欧

8、氏距离 五、模型建立与求解1问题简化首先，如果对六架飞机在区域做实时监控，再做多次调整，则每作一次航向调整都要进行一次决策，这将使问题复杂化，总体计算量较大，同时实际问题中计算也要耗费时间，效率大大降低，飞机控制的安全性必然会降低。并且对问题所给原始数据利用MATLAB软件（程序见附录1）作出原始航线图，如图1，可以粗略验证一次调整可行，既可以避免相撞，又简单易行。图1结论一：我们认为只做一次调整是优于多次调整的。其次，分析飞机飞行轨迹，作如图2，图2某架飞机在点A的航向是由A到C，此时次飞机在点C可能遭碰撞，D点是一个安全位置，如果飞机在A点早做航向的调整，方向角变化角为 ；如果待飞机飞到B

9、点时再做航向的调整，则方向的变化角为 ，很显然，> 。由此我们得出结论，早调整航向优于晚调整航向，因此我们可将问题简化。结论二：在新进入飞机到达区域边缘时便对所有飞机做航向调整，避免碰撞。2模型的导出过程1）模型推导首先引入时间变量t,其中t的取值区间为（0，）；另外，为了提高程序运行效率，简化计算，我们引入一个矩阵来控制搜索循环，较合理的设定了变量的上下界，从而大大提高决策效率。如图3所示，表示第i架飞机在规定区域可能飞行的最长时间，即该架飞机飞至区域边界四个角上点所用时间的最大值，令=min，则所有的可以构成一个N*N的矩阵。此矩阵可以运用MATLAB软件求解，程序见附录2，计算结果

10、见表1。图3得出矩阵各元素值如下表1表1T1T2T3T4T5T6T100.150260.2401170.2275030.2401170.240117T20.1502600.150260.150260.150260.15026T30.2401170.1502600.2275030.2481180.269621T40.2275030.150260.22750300.2275030.227503T50.2401170.150260.2481180.22750300.248118T60.2401170.150260.2696210.2275030.2481180将第i架飞机在t时刻的位置用坐标（，）表

11、示；其中=+Vt cos（+）；=+Vt sin（+）；则欧式距离 =8（t）=-64= 0；综上所述，可得初步的模型如下：Obj min S.t. （t）0 30°（i=1，2，3，4，5，6）；2）模型优化 在考虑LINGO软件的程序编辑和运行时，我们发现时间变量t不宜为连续变量，因此我们通过对函数 （t）的化简，对约束条件进行转化，消去变量t。过程如下：（t）=+ 其中=2Vt； =2-(-)+(-)；=+-64；很明显模型的约束可以等价的转化为一个一元二次函数在给定区间恒大于零的讨论。表示的是一条开口向上的抛物线。分别从三个方面讨论，如图4所示的三种情况，其中对称轴为-/2,

12、设此时对应的t值为。 >0 <0,0<< <0,> 图4当>0 ，可知，很明显，f(0)>0；当 <0,0<< ，则只要=；当 <0,>，只需最右端满足 （）0；则模型的约束条件可以转化如下，这个非线性模型仅是关于一个变量，即关于方向角的约束条件。S.t. （）0 (当且0<<) .接下来，通过查阅文献资料（见附录），了解了LINGO软件的运算性能和原理，综合前人的经验，我们将考虑对原目标函数 min 调整，在同样满足问题分析中的目标的情况下，即所有六架飞机的飞行方向调整角度均尽量小，对其适宜的优化。由于

13、min 是不光滑的，我们考虑用一个光滑的目标函数来替代，即min 。这样在运用LINGO软件进行求解时，运行速度和求解精度上都将大大提高。综上所述，最终的模型为：Obj min S.t. （）0 (当且0<<) .3） 模型的求解 我们针对最终的模型，带入N*N矩阵，运用LINGO软件编辑程序（见附录3），求解模型。合理编辑算法并使用全局求解程序运行，运行结果见表2，发现运行时间已经达到预想的反应速度（多次试验平均时间3s）.表2飞机123456角度改变0.0000000.0000002.062449-0.49543730.0000001.5670116.954677运行结果如上表

14、，可得出飞行调整策略，即只改变其中两架飞机，飞机1,2,5保持航向，飞机3调整+2.062449度，飞机4调整-0.4954373，飞机6调整+1.567011度。= 6.954677，综合整体调整量较小已达到最优，结果较理想。 为达到更加理想的决策速度（实际上这是符合现实需求的），我们考虑加快程序的运行速度，这样机上的乘客和工作人员的安全将得到更好的保证。我们考虑是否可以通过对模型中的飞行角度进行预估，对函数中的飞行角度进行更好的约束预估，缩小搜索区间，以优化算法，提高运行速度。理由一： 同时多次实践发现，最优解总是落在10度以，甚至更小。我们分析各架飞机的方位图，作出假设 ，如右图，区域的

15、飞机做可以忽略不计的微调，只对新进入的飞机6进行航向的调整，此时飞机6要做出的最优调整角度应该是小于原来的搜索区间最大值30度。理由二： 我们分析即使两架飞机相向飞行，根据题目条件4）可知进入该区域的飞机在到达区域边缘时，与区域飞机的距离应在60公里以上，将此数据带入，我们可做一定的估算。以飞机j为参考，飞机i正对飞机j的飞行，调整一定角度到两架飞机恰好不相擦，如图5图5其中为飞机i调整后对j的相对飞行速度，角是调整角，则有：=8km/60km;<因此，我们最后在编程时将搜索区间缩小在（，），如此再将新的程序对原数据进行计算，最后时间成功缩短至1.5s。计算结果如下表3 表3飞机1234

16、56角度改变0.0000000.0000002.062456-0.49542980.0000001.5670046.954677 六、模型检验与分析模型在实际应用中应该是要将决策时间缩短趋近于零，即时间短到可以忽略不计，这样才能给飞机以足够的时间做机械性调整，如果可以把程序运行时间进行补偿，即对飞机在决策延迟延迟时间后调接到调整信息之前的位置预测，将其带入程序运算，得出一个较精确的解，这个决策应该是更加可靠的，可以用来检验我们的模型。现在我们将预测后数据如表3（计算程序见附录5）表4飞机编号横坐标x纵坐标y方向角（度）1149.798239.6040243284.751584.63152363

17、149.6620154.7114220.54144.585150.1593159512937143149.6595230新进入0.27360.350252得到新的矩阵如下表5表5T1T2T3T4T5T6T100.1497150.240230.226970.240230.24023T20.14971500.1497150.1497150.1497150.149715T30.240230.14971500.226970.2475620.269068T40.226970.1497150.2269700.226970.22697T50.240230.1497150.2475620.2269700.24

18、7562T60.240230.1497150.2690680.226970.2475620将新矩阵数据代入程序，计算出此时的调整方案，并与调整两秒之间的计算结果通过作差和误差计算作比较，如下表6表6项目角度飞机延时0秒延时两秒差值误差10.000000 0.0000000.000000020.000000 0.0000000.000000032.062456 2.0707610.0083050.401%4-0.4954298 -0.4970576-0.00162780.327%50.000000 0.0000000.000000061.567004 1.5736810.0066770.424%

19、6.9546777.0115880.0569110.812%很显然误差较小，在保证决策精确度的情况下，该误差可以忽略不计，在实际问题中意味着我们的模型精确度已经达到，并且运行时间较好的控制在在允许围之，也就是说决策效率依然很高，完全满足要求。7、 模型的评价与推广1）纵观整个模型，我们的模型具有一般性，用两架飞机的欧式距离的限制来表示不碰撞条件，对于三维甚至多维的情况可以推广。2）该模型对1000米的高空平面的飞机飞行进行管理，是满足实际的。实际应用中，为避免两机碰撞，管理中心还可以控制飞机调整在统一空域的不同高度，这样调整难度大大降低，此模型考虑处于同一高度平面的6架飞机已经足够，已具有实际

20、意义。3)经优化模型及其算法，决策较准确，运行较快速，已经将决策效率提高，满足实际要求。4）模型也有一些不足，比如在计算组成的N*N矩阵是使用的是MATLAB软件，后面的模型优化运算却使用LINGO软件，如果可以深入研究软件，将两步运算连接在统一软件中运算，效率将更理想。5)模型计算思路也有待进一步改善，比如，在实际中，可以将在模型检验与分析环节的思路进一步深入，不仅考虑决策时间的延迟，甚至加入飞机进行机械调整所需的经验平均时间。这样将更加精准化。6）该模型其实可能应该有更加一般的解法，这是物理学里面一个运动的问题，我们可以将在模型求解环节对缩短搜索区间的证明里的思路进一步深化，即是否可以继续

21、采用相对运动原理优化模型，列出新的不碰撞控制条件。八、参考文献1顾其行 ，国际航空运输管理, ：知识，19872朱道元，数学建模精品案例，：东南大学，19993程极泰，最优设计的数学方法，：国防工业，19944怀中，运筹学基础教程，：人民，20105中庚，实用运筹学，：清华大学 ， 20076金星、薛毅，优化建模与LINDO/LINGO软件，：清华大学，20057华林，汤志高，金平，Lingo在飞行管理中的应用，科技信息(科学教研) 2007年17期 九、附录1% 函数picture做飞机航线预测图function picture(X,Y,a)% 输入X:飞机原始横坐标向量% 输入Y:飞机原始

22、纵坐标向量% 输入a:飞机原始方向角向量figureaxis(0,160,0,160)for i=1:length(X) hold on x=0:0.1:160; y=tan(a(i)*pi/180)*(x-X(i)+Y(i); if a(i)<=180 plot(x,y,'-r') %红色虚线为向上飞的飞机航线 else plot(x,y,'-b') %蓝色直线为向下飞的飞机航线 endendhold onfor i=1:length(X) plot(X(i),Y(i),'*k') %黑色的“*”号为飞机初始位置end2% 函数time计

23、算任意两架飞机中最先飞出区域的一架所用最少时间矩阵function A=time(X,Y,Z,W)% 输入X:飞机原始横坐标向量% 输入Y:飞机原始纵坐标向量% 输入Z：区域四角横坐标向量% 输入W：区域四角纵坐标向量% 计算每架飞机到区域四角的距离矩阵Sfor i=1:length(X) for j=1:length(Z) A(i,j)=sqrt(X(i)-Z(j)*(X(i)-Z(j)+(Y(i)-W(j)*(Y(i)-W(j); endend% 计算每架飞机距区域四角最长距离向量ss=zeros(1,length(X);for i=1:length(X) for j=1:length(Z

24、) if s(i)<=A(i,j) s(i)=A(i,j); end endend% 计算任意两架飞机飞出区域所用最少时间矩阵Tt=s./800T=zeros(length(t),length(t);for i=1:length(t) for j=1:length(t) if t(i)=t(j) T(i,j)=0; else if t(i)>=t(j) T(i,j)=t(j); else T(i,j)=t(i); end end endend% 输出：任意两架飞机中最先飞出区域的一架所用最少时间矩阵TT3MODEL:TITLE 飞行管理问题;SETS:Plane/1.6/: x0,

25、 y0, cita0, cita1, d_cita;link(plane, plane)|&1 #LT# &2: b,c,z;link1(plane, plane):T;ENDSETSDATA:!输入飞机初始坐标和角度（角度制）;x0 y0 cita0 =150402438585236150155220.5145501591301502300052; max_cita = 30;!每两架飞机飞出区域所用最短时间矩阵，数据来自time.m文件 ;T = 0 0.1503 0.2401 0.2275 0.2401 0.2401 0.1503 0 0.1503 0.1503 0.15

26、03 0.1503 0.2401 0.1503 0 0.2275 0.2481 0.2696 0.2275 0.1503 0.2275 0 0.2275 0.2275 0.2401 0.1503 0.2481 0.2275 0 0.2481 0.2401 0.1503 0.2696 0.2275 0.2481 0 ;V=800;ENDDATAINIT:d_cita = 0 0 0 0 0 0;ENDINITfor(plane(i): cita1(i) - cita0(i) = d_cita(i);!模型中二次函数系数;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 (c