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

一个飞行管理问题

摘要

在某一空域里对飞机的飞行合理管理事关重大，比如乘客及机上工作人员生命财产安全和航空公司的运作效益等。本文通过对飞机飞行管理问题的研究，得到了调整飞机架数较少同时调整幅度均最小（平方和最小）的飞行管理最优安排的非线性模型，这样既使得乘客所受影响达到最少，也便于飞机调整，还有利于飞机回到原来的航线，同时还在决策时间上对模型进行了优化和调整。

本文不仅一般性地将不相撞的问题转化为欧式距离控制，而且很巧妙的将不碰撞条件转化成简单的二次函数标准形式进行含参讨论，建立一个只含有转向角变量的模型。并且再次很妙的具体化区域内受控时间形成矩阵，大大得简化运算，节约了大量运算的时间，便于管理人员控制操作，从而确保飞机的安全。更重要的是最后结合实际缩短了搜索区间，并优化算法，使得决策更加高效。最后的延时检验也充分体现了模型的可靠性。

关键字：欧氏距离约束转化缩短搜索区间时间矩阵延时检验

一、问题重述

在约10000米的高空某边长为160公里的正方形区域内，经常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。当一驾欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞。现假定条件如下：

1）不碰撞的标准为任意两架飞机的距离大于8公里；

2）飞机飞行方向角调整的幅度不应超过30度；

3）所有飞机飞行速度均为每小时800公里；

4）进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60公里以上；

5）最多需考虑6架飞机；

6）不必考虑飞机离开此区域后的状况。

请你对这个避免碰撞的飞机管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过0.01度），要求飞机飞行方向角调整的幅度尽量小。

设该区域内4个顶点的坐标为(0,0),(160,0),(160,160),(0,160)。

记录数据为：

飞机编号横坐标x纵坐标y方向角（度）

115040243

28585236

3150155220.5

414550159

5130150230

新进入0052

注：方向角指飞行方向与x轴方向的夹角。

试根据实际应用背景对你的模型进行评价与推广。

二、问题分析

1初步分析

根据问题容易知道，这显然是一个优化问题，当两架飞机可能发生碰撞时，即在规定区域内某一时刻两架飞机之间的距离小于8公里，因此要调整飞行方向一定角度，保证任意两架飞机在区域内任意时刻，两者的距离均不小于8公里，避免相撞。考虑到调整角度应尽量小，可以简化飞行方向调整策略，降低调整难度，同时减轻机内乘客及工作人员的不适。此外由此初步确定了调整目标：所有六架飞机的飞行方向调整角度均尽量小。

2解决方案

由于所有飞机均处于1000米得高空作水平飞行，可将飞机飞行的空域视为二维平面xoy中的一个正方形区域,顶点为(0,0),(160,0),(160,160),(0,160)。于是可以引入时间变量后，确定每架飞机在任意时刻的坐标，列出任意两点的欧氏距离，令其恒大于8公里，则得出一个重要约束条件。再结合变化角度应小于30度，即可得出约束条件，然后运用LINGO软件编辑程序进行求解。

为提高决策效率，在反复试验中又可对约束条件进行调整。

三、条件假设

1.不碰撞的标准为任意两架飞机的距离在以后任何一个时间里大于8公里；

2.飞机飞行方向角调整的幅度不应超过30度；

3.所有飞机飞行速度均为每小时800公里；

4.进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60公里以上。即

在计算如何最优地调整各架（包括新进入的）飞机飞行的方向角时，飞行管理中心得出合理的最优调整措施，；

5.最多需考虑6架飞机。

6.此处忽略飞机在执行过程中所需耗费的时间，即假设从飞机管理中心发出的调整信

息飞机马上可以接收并执行，不存在滞后或延迟；

7.飞行管理中心在计算飞行调整信号和发出信号所需时间内，忽略各架飞机（包括刚

进入的飞机）调整航向前飞行数据的变化；

8.假定飞机在该区域内完全依赖飞行管理中心调度；

9.假设飞机在飞出区域之后，飞行员可以自觉调整飞行策略，回归原始航线，即飞行

管理中心不必考虑飞机离开此区域后的状况。

四、符号说明

符号

含义

0i X 第i 架飞机在初始时刻的横坐标 0j X

第j 架飞机在初始时刻的横坐标 t i X 第i 架飞机在t 时刻的横坐标 t j X

第j 架飞机在t 时刻的横坐标 0i Y 第i 架飞机在初始时刻的纵坐标 0j Y

第j 架飞机在初始时刻的纵坐标 t i Y 第i 架飞机在t 时刻的纵坐标 t j Y

第j 架飞机在t 时刻的纵坐标

0i θ 第i 架飞机在初始时刻飞行方向与X 轴正向的夹角

0j θ

第j 架飞机在初始时刻飞行方向与X 轴正向的夹角

i θ

第i 架飞机在t 时刻飞行方向与X 轴正向的夹角 j θ

第j 架飞机在t 时刻飞行方向与X 轴正向的夹角 i θ?

第i 架飞机飞行方向角的调整幅度 j θ?

第j 架飞机飞行方向角的调整幅度

i T

第i 架飞机在规定区域内可能飞行的最长时间 j T 第j 架飞机在规定区域内可能飞行的最长时间 ij T T 是一个6*6矩阵；ij T =min{i T ，j T } V

飞机的飞行速度

ij d

飞机i 与飞机j 的欧氏距离

五、模型建立与求解

1问题简化

首先，如果对六架飞机在区域内做实时监控，再做多次调整，则每作一次航向调整都要进行一次决策，这将使问题复杂化，总体计算量较大，同时实际问题中计算也要耗费时间，效率大大降低，飞机控制的安全性必然会降低。并且对问题所给原始数据利用MATLAB软件（程序见附录1）作出原始航线图，如图1，可以粗略验证一次调整可行，既可以避免相撞，又简单易行。

图1

结论一：我们认为只做一次调整是优于多次调整的。

其次，分析飞机飞行轨迹，作如图2，

图2

某架飞机在点A的航向是由A到C，此时次飞机在点C可能遭碰撞，D点是一个安全位置，如果飞机在A点早做航向的调整，方向角变化角为α；如果待飞机飞到B 点时再做航向的调整，则方向的变化角为β，很显然，β> α。由此我们得出结论，早调整航向优于晚调整航向，因此我们可将问题简化。

结论二：在新进入飞机到达区域边缘时便对所有飞机做航向调整，避免碰撞。

2模型的导出过程

1）模型推导

首先引入时间变量t,其中t的取值区间为（0，

T）；另外，为了提高程序运行效率，

ij

简化计算，我们引入一个矩阵来控制搜索循环，较合理的设定了变量的上下界，从而大大提高决策效率。

T表示第i架飞机在规定区域内可能飞行的最长时间，即该架飞机飞至区如图3所示，

i

T=min{i T，j T}，则所有的ij T可以构成一个域边界四个角上点所用时间的最大值，令

ij

N*N的矩阵。此矩阵可以运用MATLAB软件求解，程序见附录2，计算结果见表1。

图3

得出矩阵各元素值如下表1

表1

ij T

T1 T2 T3 T4 T5 T6 T1 0 0.15026 0.240117 0.227503 0.240117 0.240117 T2 0.15026 0 0.15026 0.15026 0.15026 0.15026 T3 0.240117 0.15026 0 0.227503 0.248118 0.269621 T4 0.227503 0.15026 0.227503 0 0.227503

0.227503 T5 0.240117 0.15026 0.248118 0.227503 0

0.248118

T6 0.240117 0.15026 0.269621 0.227503 0.248118

将第i 架飞机在t 时刻的位置用坐标（t i X ，t i Y ）表示； 其中

t i X =0i X +Vt cos （0i θ+i θ?）； t i Y =0i Y +Vt sin （0i θ+i θ?）

； 则欧式距离 ij d =22)()(t j t i t j t i Y Y X X -+-≥8

ij f （t ）=2

ij

d -64 =-?++-?++2

0000))]cos(())(tcos [(j j j i i i Vt X V X θθθθ

64))]sin(())sin([(2

0000-?++-?++j j j i i i Vt Y Vt Y θθθθ ≥0；

综上所述，可得初步的模型如下： Obj min

∑=?6

i i

θ

S.t. ij f （t ）≥0

i θ?≤30°（i=1，2，3，4，5，6）；

2）模型优化

在考虑LINGO 软件的程序编辑和运行时，我们发现时间变量t 不宜为连续变量，因此我们通过对函数 ij f （t ）的化简，对约束条件进行转化，消去变量t 。 过程如下：

ij f （t ）=2

ij z +ij b ij z +ij c

其中

ij z =2Vt 错误！未找到引用源。

； ij b =2[-(错误！未找到引用源。-错误！未找到引用源。)错误！未找到引用源。+(0i Y 错误！未找到引用源。2

cos

j

i θθ+] ； ij c =错误！未找到引用源。+错误！未找到引用源。-64；

很明显模型的约束可以等价的转化为一个一元二次函数在给定区间内恒大于零的讨论。表示的是一条开口向上的抛物线。分别从三个方面讨论，如图4所示的三种情况，

其中对称轴为-ij b /2,设此时对应的t 值为?

ij t 。

ij b >0 ij b <0,0

ij t >ij T

图4

当ij b >0 ，可知，很明显，f(0)>0；当 ij b <0,0

≤-ij ij c b ；

当 ij b <0,?

ij t >ij T ，只需最右端满足 ij f （ij T ）≥0；

则模型的约束条件可以转化如下，这个非线性模型仅是关于一个变量，即关于方向角θ的约束条件。

S.t. ij f （ij T ）≥0

042

≤-ij ij c b (当0

接下来，通过查阅文献资料（见附录），了解了LINGO 软件的运算性能和原理，综合前人的经验，我们将考虑对原目标函数 min

∑=?6

i i

θ

调整，在同样满足问题分析中的

目标的情况下，即所有六架飞机的飞行方向调整角度均尽量小，对其适宜的优化。由于min

∑=?6

i i

θ

是不光滑的，我们考虑用一个光滑的目标函数来替代，即min

∑=?6

2

)

(i i

θ。

这样在运用LINGO 软件进行求解时，运行速度和求解精度上都将大大提高。 综上所述，最终的模型为： Obj min

∑=?6

2

)

(i i

θ

S.t. ij f （ij T ）≥0

042

≤-ij ij c b (当0

3）模型的求解

我们针对最终的模型，带入N*N 矩阵，运用LINGO 软件编辑程序（见附录3），求解模型。合理编辑算法并使用全局求解程序运行，运行结果见表2，发现运行时间已经达到预想的反应速度（多次试验平均时间3s ）.

表2

运行结果如上表，可得出飞行调整策略，即只改变其中两架飞机，飞机1,2,5保持航向，飞机3调整+2.062449度，飞机4调整-0.4954373，飞机6调整+1.567011度。∑?2)(i θ= 6.954677，综合整体调整量较小已达到最优，结果较理想。

为达到更加理想的决策速度（实际上这是符合现实需求的），我们考虑加快程序的运行速度，这样机上的乘客和工作人员的安全将得到更好的保证。我们考虑是否可以通过

飞机

1 2

3

4

5

6

∑?2

)

(i

θ

角度改变 0.00000

0.000000

2.062449

-0.4954373

0.000000

1.567011

6.954677

对模型中的飞行角度进行预估，对函数中的飞行角度进行更好的约束预估，缩小搜索区间，以优化算法，提高运行速度。

理由一：

同时多次实践发现，最优解总是落在10度以内，甚至更小。我们分析各架飞机的方位图，作出假设 ，如右图，区域内的飞机做可以忽略不计的微调，只对新进入的飞机6进行航向的调整，此时飞机6要做出的最优调整角度应该是小于原来的搜索区间最大值30度。 理由二：

我们分析即使两架飞机相向飞行，根据题目条件4）可知进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60公里以上，将此数据带入，我们可做一定的估算。以飞机j 为参考，飞机i 正对飞机j 的飞行，调整一定角度到两架飞机恰好不相擦，如图5

图5

其中ij V 为飞机i 调整后对j 的相对飞行速度，角γ是调整角，则有：

γsin =8km/60km; γ≈ 62.67< 10

因此，我们最后在编程时将搜索区间缩小在（ 10-， 10），如此再将新的程序对原数据进行计算，最后时间成功缩短至1.5s 。计算结果如下表3

表3

六、模型检验与分析

飞机

1

2

3

4

5

6

∑?2

)

(i

θ

角度改变

0.00000

0.000000

2.062456

-0.4954298

0.000000

1.567004

6.954677

模型在实际应用中应该是要将决策时间缩短趋近于零，即时间短到可以忽略不计，这样才能给飞机以足够的时间做机械性调整，如果可以把程序运行时间进行补偿，即对飞机在决策延迟延迟时间后调接到调整信息之前的位置预测，将其带入程序运算，得出一个较精确的解，这个决策应该是更加可靠的，可以用来检验我们的模型。现在我们将预测后数据如表3（计算程序见附录5）

表4

飞机编号横坐标x纵坐标y方向角（度）

1149.798239.6040243

284.751584.6315236

3149.6620154.7114220.5

4144.585150.1593159

512937143149.6595230

新进入0.27360.350252

得到新的矩阵如下表5

表5

T T1 T2 T3 T4 T5 T6 ij

T1 0 0.149715 0.24023 0.22697 0.24023 0.24023

T2 0.149715 0 0.149715 0.149715 0.149715 0.149715

T3 0.24023 0.149715 0 0.22697 0.247562 0.269068

T4 0.22697 0.149715 0.22697 0 0.22697 0.22697

T5 0.24023 0.149715 0.247562 0.22697 0 0.247562

T6 0.24023 0.149715 0.269068 0.22697 0.247562 0

将新矩阵数据代入程序，计算出此时的调整方案，并与调整两秒之间的计算结果通过作差和误差计算作比较，如下表6

表6

延时0秒延时两秒差值误差

1 0.000000 0.0000000.000000 0

2 0.000000 0.0000000.000000 0

3 2.062456 2.0707610.008305 0.401%

很显然误差较小，在保证决策精确度的情况下，该误差可以忽略不计，在实际问题中意味着我们的模型精确度已经达到，并且运行时间较好的控制在在允许范围之内，也就是说决策效率依然很高，完全满足要求。

七、模型的评价与推广

1）纵观整个模型，我们的模型具有一般性，用两架飞机的欧式距离的限制来表示不碰撞条件，对于三维甚至多维的情况可以推广。

2）该模型对1000米的高空平面内的飞机飞行进行管理，是满足实际的。实际应用中，为避免两机碰撞，管理中心还可以控制飞机调整在统一空域的不同高度，这样调整难度大大降低，此模型考虑处于同一高度平面内的6架飞机已经足够，已具有实际意义。

3)经优化模型及其算法，决策较准确，运行较快速，已经将决策效率提高，满足实际要求。

4）模型也有一些不足，比如在计算ij T 组成的N*N 矩阵是使用的是MATLAB 软件，后面的模型优化运算却使用LINGO 软件，如果可以深入研究软件，将两步运算连接在统一软件中运算，效率将更理想。

5)模型计算思路也有待进一步改善，比如，在实际中，可以将在模型检验与分析环节的思路进一步深入，不仅考虑决策时间的延迟，甚至加入飞机进行机械调整所需的经验平均时间。这样将更加精准化。

6）该模型其实可能应该有更加一般的解法，这是物理学里面一个运动的问题，我们可以将在模型求解环节对缩短搜索区间的证明里的思路进一步深化，即是否可以继续采用相对运动原理优化模型，列出新的不碰撞控制条件。

八、参考文献

[1]顾其行 ，《国际航空运输管理》, 北京 ：知识出版社，1987 [2]朱道元，《数学建模精品案例》，南京：东南大学出版社，1999 [3]程极泰，《最优设计的数学方法》，北京：国防工业出版社，1994 [4]张怀中，《运筹学基础教程》，上海：上海人民出版社，2010 [5]韩中庚，《实用运筹学》，北京：清华大学出版社 ， 2007 [6]谢金星、薛毅，《优化建模与LINDO/LINGO 软件》，北京：清华大学出版社，2005 [7]曹华林，汤志高，金平，《Lingo 在飞行管理中的应用》，《科技信息(科学教研)》 2007年17期

4 -0.4954298 -0.4970576 -0.0016278

0.327% 5 0.000000 0.000000 0.000000 0

6

1.567004 1.573681

0.006677 0.424% ∑?2)

(i

θ

6.954677

7.011588

0.056911

0.812%

九、附录

[1]

% 函数picture做飞机航线预测图

function picture(X,Y,a)

% 输入X:飞机原始横坐标向量

% 输入Y:飞机原始纵坐标向量

% 输入a:飞机原始方向角向量

figure

axis([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') %蓝色直线为向下飞的飞机航线

end

end

hold on

for i=1:length(X)

plot(X(i),Y(i),'*k') %黑色的“*”号为飞机初始位置

end

[2]

% 函数time计算任意两架飞机中最先飞出区域的一架所用最少时间矩阵function A=time(X,Y,Z,W)

% 输入X:飞机原始横坐标向量

% 输入Y:飞机原始纵坐标向量

% 输入Z：区域四角横坐标向量

% 输入W：区域四角纵坐标向量

% 计算每架飞机到区域四角的距离矩阵S

for 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)));

end

end

% 计算每架飞机距区域四角最长距离向量s

s=zeros(1,length(X));

for i=1:length(X)

for j=1:length(Z)

if s(i)<=A(i,j)

s(i)=A(i,j);

end

end

end

% 计算任意两架飞机飞出区域所用最少时间矩阵T

t=s./800

T=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

end

end

% 输出：任意两架飞机中最先飞出区域的一架所用最少时间矩阵T T

[3]

MODEL:

TITLE 飞行管理问题;

SETS:

Plane/1..6/: x0, y0, cita0, cita1, d_cita;

link(plane, plane)|&1 #LT# &2: b,c,z;

link1(plane, plane):T;

ENDSETS

DATA:

!输入飞机初始坐标和角度（角度制）;

x0 y0 cita0 =

150 40 243

85 85 236

150 155 220.5

145 50 159

130 150 230

0 0 52

;

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.1503 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;

ENDDATA

INIT:

d_cita = 0 0 0 0 0 0;

ENDINIT

@for(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 ((cita1(i)+cita1(j))*3.14159265/360);

c(i,j) = (x0(i) -x0(j)) ^2 + (y0(i) -y0(j)) ^2 - 64;

z(i,j)=2*V*T(i,j)*@sin((cita1(i)-cita1(j))*3.14159265/360);

);

!函数值约束;

@for(link(i,j): z(i,j)^2+ b(i,j)*z(i,j)+c(i,j)> 0);

@for(link(i,j): @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(i,j),

b(i,j)^2-4*c(i,j),-1) < 0);

@for(link: @free(b));

@for(link: @free(z));

!角度搜索范围;

@for(plane: @bnd( - max_cita, d_cita, max_cita));

!目标函数：角度改变量平方和最小;

[obj] MIN = @SUM(plane(i):(d_cita(i)^2));

END

[5]

% 函数forecast 预测区域中的飞机在t时刻的位置

function forecast(X,Y,a,t)

% 输入X:飞机原始横坐标向量

% 输入Y:飞机原始纵坐标向量

% 输入a:飞机原始方向角（角度）向量

% 输入t：预测t时间（秒）后的坐标

x=zeros(1,length(X)); % 飞机t时刻横坐标向量初始化

y=zeros(1,length(X)); % 飞机t时刻纵坐标向量初始化

% 计算第i架飞机在时间t时的横坐标、纵坐标

for i=1:length(X)

x(i)=X(i)+800*t/3600*cos(a(i)*pi/180);

y(i)=Y(i)+800*t/3600*sin(a(i)*pi/180);

end

%输出：飞机在t时刻的横坐标、纵坐标向量

x

y