飞行管理 数学建模
B 题:飞行管理问题
摘要:
飞行管理问题是一个既现实又重要的课题,本文利用偏转角度尽可能的小建立两个非线性规划模型。
模型一:时间模型。考虑到各架飞机的偏转角有正有负,在此模型中,对于各架飞机调整选取各个偏转角的绝对值的和作为目标函数,要求任意两架飞机任意时刻的距离大于8公里,则可以求出任意两架飞机的距离ij d 。此时,两架飞机距离ij d 是时间t 与各个飞机偏转角i θ?的函数,编写程序时将t 离散化,且t 有最大值0.2828s (沿对角线飞过的时间),这样可得到表1-1的结果:
表1-1
模型二:闭塞区域模型。在两架飞机中,将其中一架看成“静止”,另一架相对于它而运动。而以“静止”飞机为圆心,km 8为半径的圆形区域构成该飞机的闭塞区域,任意一架飞机的方向角均不能在此区域内,则为不相撞。为此,本文用复变函数的知识表示各架飞机的速度,从而算出相对速度,再求出相对位移,以相对速度与相对位移的夹角大于每两架飞机的临界夹角来刻画不相撞。目标函数为每架飞机偏转角的平方和。利用计算机编程得到表1-2的结果:
表1-2
对于上述两个非线性规划,在理论方面,本文利用SUTM 内点法(障碍函数法)进行算法描述,在操作方面,分别利用lingo 语言与MATLAB 语言直接编写程序进行计算
关键词:非线性规划、复变函数、SUMT 内点法、闭塞区域、禁飞角
一、问题重述
1.背景知识
与其他交通工具相比,飞机以其速度快、安全舒适等特点在交通领域占据了绝对地位。而近年来飞机事故的频繁发生也预示着飞机存在一定的安全隐患。经调查造成飞机相撞事故的原因主要是人、飞机(设备)、环境,而人的因素是事故中通常起主体作用的因素,直接影响事故的发生和结局。飞机事故的发生难以预测且死亡率极高,所以航空安全机制的健全,航空人员素质的提高已变得刻不容缓。
2.问题重述
在约10000米的高空某边长为160公里的正方形区域内,经常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据,以便进行飞行管理。当一驾欲进入该区域的飞机到达区域边缘时,记录其数据后,要立即计算并判断是否会与区域内的飞机发生碰撞。如果会碰撞,则应计算如何调整各架(包括新进入的)飞机飞行的方向角,以避免碰撞。现假定条件如下:
1)不碰撞的标准为任意两架飞机的距离大于8公里;
2)飞机飞行方向角调整的幅度不应超过30度;
3)所有飞机飞行速度均为每小时800公里;
4)进入该区域的飞机在到达区域边缘时,与区域内飞机的距离应在60公里以上;
5)最多需考虑6架飞机;
6)不必考虑飞机离开此区域后的状况。
请你对这个避免碰撞的飞机管理问题建立数学模型,列出计算步骤,对以下数据进行计算(方向角误差不超过0.01度),要求飞机飞行方向角调整的幅度尽量小。
设该区域内4个顶点的坐标为(0,0),(160,0),(160,160),(0,160)。
记录数据为:
二、 基本假设
1.飞机在规定区域内飞行,不会发生任何自然灾害;
2.飞行人员能够探知飞行区域内所有飞机的飞行轨迹,并能够及时作出相应调整;
3.飞机在规定区域内最多只调整一次,且在第六架飞机刚进入区域时调整;
4.飞机在调整过程中是瞬间完成,飞行轨迹可看成直线;
5.不碰撞的标准是任何两架飞机的距离大于km 8;
6.飞机飞行方向的调整幅度不应超过?30; 7.所有飞机调整角度总和最小为最优方案; 8.暂时不考虑飞机离开此区域的情况; 9.飞机在调整之后按固定速度飞行; 10.飞机在飞行过程中可看为质点。
三、 参数说明
1.()i i y x ,表示第i 架飞机的初始位置;
2.i θ表示第i 架飞机初始飞行角度;
3.i θ?表示第i 架飞机的调整角度;
4.ij s 表示第i 架飞机和第j 架飞机的距离; 5.ij d 表示第i 架飞机和第j 架飞机距离的平方;
5.i v 表示第i 架飞机的飞行速度;
6.t 表示飞机的飞行时间;
7.i t 表示第i 架飞机飞出飞行区域的时间; 10.;架飞机的飞行角度架飞机相对于第表示第i j ij θ 11.架飞机的速度矢量;表示第k Z k
12.飞机的飞行矢量架飞机指向第表示第j k f kj ; 13.架飞机的禁飞角的大小
架飞机对第表示第j k kj α。 四、 问题分析
图4-1:飞机初始状态
从距离约束条件来说,在初始状态下,我们可以通过调整飞行方向,目的使飞机在飞行区域内满足任意两架飞机的距离都大于8公里。调整角度我们可以假设已知,将上述初始角替换成调整角,通过计算,我们知道存在调整方法能够使飞机安全飞出飞行区域,且调整方式多样。若我们假设调整角度对于每架飞机来说,其满意度是一样的,那在这一系列调整方式中,必定存在调整角度之和达到最小的调整方式。在飞行过程中,初始位置和飞行方向对飞机的飞行路径有决定性作用,而在速度已知的情况下,时间的确定就可以直接决定飞机的具体位置。由此,我们想到,如果取适当的一段时间为步长,确定检测时间点,若这些时间点各架飞机都能满足不相撞的条件,那我们可以近似认为,飞机在这些飞行方向上是符合安全飞行条件的,之后只需找出调整最小的方式即可。
从角度约束条件来说,将飞机看做质点,在飞行区域内,以飞机为圆心,4公里为半径的圆形区域内是不能存在其他飞机的。即每架飞机都有其自身的闭塞区域,而且是动态的闭塞区域。在边长为160公里的正方形二维平面内,这6个闭塞区域在运动过程中是不允许存在重合部分的。飞机在飞行区域内是按照同一方向飞行,也就是说在该正方形区域内,闭塞区域是按直线行走的。运动初期,各圆的位置已知,飞行方向可调整(先不考虑调整角度的限制问题),而调整的目的是为了让这些圆在以相同速度且沿直线的行走过程中不存在重合部分。将模型简化,仅看两架飞机的运动形式,以一个圆的圆心出发,引发出两条相对于另一个圆的切线,而这两条切线在包含另一个圆的一侧所构成的夹域是禁飞域,夹角为禁飞角。在这种解法中不仅要考虑飞机之间的相对运动速度还要考虑飞机之间初始相对位置的具体关系。讨论初期,我们确定飞机之间相对运动的速度还是
引用原有的坐标系,飞机之间相对位置关系也在原有的坐标系中讨论,该模型的约束条件就是相对运动方向与相对位置方向之间的相对夹角只要满足不在禁飞夹角中即可。首先我们选择通过几何方法计算相对夹角,而在进一步计算过程中,我们发现相对夹角若以角度制表示要分不同种情况讨论,过程复杂难以计算。之后我们又将模型进一步改进,对于坐标系的选取,若以飞机之间相对位置所在向量确定相对坐标则使约束条件简化。在计算方面,为了避免角度的多情况影响,我们选取以复数形式表示角度,这样就将相对夹角的形式统一起来,计算也有了进一步突破。
在计算出具体结果后,就需要对结果的有效性进行检验,而我们选取的检验条件就是飞机以调整后的角度飞行,在离开飞行区域时,是否与飞行区域内的其他飞机的距离大于60公里,若都满足大于60公里,则认为结果可行性较强,若有出入,在对模型进行相应调整。
五、模型建立与求解
模型一
设飞机初始状态下的坐标为),(i i y x 易知任意时刻飞机的坐标为
)cos(i i i i vt x X θθ?++= )sin(i i i i vt y Y θθ?++= 则两架飞机的距离:
222)()(j i j i ij ij Y Y X X s d -+-==
讨论飞机在规定区域内的飞行状态,则飞机有飞行时间的限制:若飞机沿对角线飞行,则飞行距离最远,飞行时间最长,此时800
2
160max =
t ,即()2828.0,0∈t ,
在飞行时间内,如果任意两架飞机的距离都超过8公里,则认为不会发生碰撞事件,即:
()2828
.0,0∈?t )61,61(64 ==>j i d ij (1)
首先讨论如果不改变飞机偏转角是否有两架飞机相撞,即满足(1)式的j i 和
当小时公里/800
=v ,()i i y x ,和i θ已知的情况下,经计算: ()()
64
52sin 5
.220sin 80015552cos 5.220cos 8001502
2
2
36<=-++-+=?
?
?
?t t s 和
()()
6452sin 230
sin 80015052cos 230cos 8001302
2
2
56<=-++-+=?
?
?
?t t s
在()2828
.0,0∈t 内有实解,则如果按照初始方向飞行会发生碰撞事件,需要
做出相应调整。
为了确保飞行时不会发生碰撞事件,则任意两架飞机的飞行距离不得小于公里8,在此条件下寻求6架飞机能够安全通过飞行区域的最小调整方式。而对于“最小的调整方式”,本文采取下述两种理解方式。
1.偏转总和最小。因为每架飞机偏转角可正可负,所以,本模型取偏转角的平方和作为目标函数。结合上述分析与条件,可建立如下非线性规划模型
∑=?6
12)(min i i θ
..t s ()??
?
??∈===>2828.0,0)61(30
)61,61(64
t i j i d i ij θ 其中:
222)()(j i j i ij ij Y Y X X s d -+-== 2.每架飞机偏转都最小。 T =min
..t s ???
??
?
?∈==>)2828.0,0(30)
61,61(64t T j i d i i ij θθ
编程时,我们将t 离散化,分为200等分,
即00094.0300
2828
.0==
dt 误差km vdt ds 7541.0== (1)求得的最终结果为:
01=?θ;02=?θ;?=?540629
.23θ;04=?θ;05=?θ;?=?074526.16θ。 (2)求得的最终结果为:
?=?2906.21θ;?=?2500.02θ;?
=?5406
.23θ;04=?θ;?=?2500.05θ;0745.16=?θ
对比上述结果,结合题中误差在?01.0之内,可只取第一种方案,这样既满
足飞机偏转总和最小,又满足每架飞机偏转最小。
下面检验对于上述结果是否满足“进入该区域的飞机在到达区域边缘时,与区域内飞机的距离应在60公里以上;”
飞机的飞行路径方程为:
()()()i i i i x x y y θθ?+=--tan /;
求得每架飞机飞出飞行区域的坐标分别为
()0,6190.129;()0,6668.27;()1227.15,0; ()6603.105,0;()0,1351.4;()160,1314.120
飞机飞出飞行区域的时间为:
()()v
y y x x t i i i i i 2
121-+-=
求得每架飞机的飞行时间为:
0561.01=t ;1282.02=t ;2564.03=t ; 1941.04=t ;2448.05=t ;2501.06=t 。
飞机飞出飞行区域的先后顺序为:1,2,4,5,6,3。按照此飞行时间顺序
由少到多重新排序。
而飞机飞出区域要保证和飞行区域内的其他飞机距离都大于60公里; 所以:
()()()
()
()6~1,6~136002
12
1+==>-+-i j i t y y t x x i j i i j i
即:
()()3600sin cos 2
12
1>--+--j i j i j i j i vt y y vt x x
θθ成立,
根据上述思想编写C++程序(见附录)
经检验。只当第5架飞机飞出区域是与第3架飞机的距离小于60公里,再缩短时间步长使结果更加精确。 模型二
先将模型简化为两架飞机,考虑两架飞机的相对运动,若保证飞机在飞行区域内不发生碰撞则相对飞行角度有一定的范围限制。以第j 架飞机坐标为圆心,以8为半径作圆,此圆即为第j 架飞机的闭塞区域。第k 架飞机相对于第j 架飞机的速度为kj v 。由第k 架飞机坐标引出的相对于此圆的两条切线的夹角kj α为第
k 架飞机做相对运动的禁飞方向。
规定:kj β为第k 架飞机相对于第j 架飞机的夹角:
以k 为原点,j k →为极轴kj f ,将极轴kj f 转至kj v 的角度即为kj β,逆时针为正,顺时针为负。(按照这种规定方法,求出的值与MATLAB 语言中angle 命令产生的结果一致)
设A v =
则第k 架飞机的速度矢量:
()k k i k Ae Z θθ?+=
则第k 架飞机相对于第j 架飞机的相对速度矢量: )()
()
(j j k k i i j k kj
e
e
A Z Z v θθθθ?+?+-=-=
第k 架飞机与第j 架飞机相对位置矢量(极轴):
)(k j k j kj y y i x x f -+-=
则按照上述规定: )arg(
kj
kj kj f v =β
图4-2
由图1-2知:
kj
kj s 8
sin =α
以各个偏转角绝对值的和作为目标函数建立如下非线性规划:
∑=?6
1min i i θ
..t s ??
???=<
?==≥)
61(6)61,61( i j k i kj kj πθαβ
对于此类非线性规划,通常运用罚函数法。
罚函数法基本思想是通过构造罚函数把约束问题转化为一系列无约束最优化问题,进而用无约束最优化方法去求解.这类方法称为序列无约束最小化方法.简称为SUMT 法.
其一为SUMT 外点法,其二为SUMT 内点法.本文运用SUMT 内点法(障碍函数法):
考虑问题:()())1( ,...,2,1i 0 ..min ?
??=≥m X g t s X f i
设集合(){}00,,2,1,0| D m i X g X D i ,φ≠=>= 是可行域中所有严格内点的集合。
构造障碍函数
()()()()()
∑
∑==+=+=m
i i m i i X g r X f r X I X g r X f r X I r X I 1
1
1
)(),( ln ,,或: 其中称()()
∑
∑==m
i i m i i X g r X g r 1
1
1
ln 或为障碍项,为障碍因子,这样问题(1)就转化为求一系列极值问题:())(得k k
k D
X r X r X I ,min 0
∈ 内点法的迭代步骤:
S1:给定允许误差0>ε,取10,01<<>βr ; S2:求出约束集合D 的一个内点00D X ∈,令1=k ;
S3:以01
D X
k ∈-为初始点,求解()k D
X r X I ,min 0
∈,其中0
D X ∈的最优解,设为()0D r X X k k ∈=;
S4:检验是否满足()ε≤-∑=m
i k
i X
g r
1
ln 或()ε≤∑=m
i i
k X g r 1
1,
满足,停止迭代,k X X ≈*
,
否则取k k r r 1β=+,令1+=k k ,返回S3.
根据题意可取?=01.0ε,?
??
???=∈=)61(),6,0()(|0 i i X X D π的六维样本空间,根
据算法进行迭代。
考虑到实现此算法的复杂性以及当今计算机的发展,本文采取直接运用MATLAB 编程得到如下计算结果:
01=?θ;02=?θ;?
=?540629.23θ;04=?θ;05=?θ;074526.16=?θ
五、模型的评价与改进
模型一
优点:
(1)将约束条件简化,以条件检验结果,使模型简单化;
(2)该模型比较传统,便于对模型结构的理解和把握;
(3)将时间作为已知条件,选取特定时间点进行条件约束,计算量大大减少;缺点:
(1)对时间步长的选取要求较高;步长选取不当直接影响结果的有效性;(2)计算软件的局限性会导致结果误差较大。
模型二
优点:
(1)将距离之间的关系转化为角度之间的关系,使模型更加直观;
(2)相对运动和相对位置的选取有利于约束条件的进一步简化;
(3)运用复数形式表示矢量,避免讨论多种情况;
(4)闭塞区间的引用使模型更加形象,便于理解;
(5)内点法计算最优化问题使约束条件在结果中的选取也更具优越性;
缺点:
(1)对角度的处理要求较高。
对于该模型的讨论,我们只考虑了人为因素的影响,具有一定的局限性。我们假定人员的素质极高,能够按照指示迅速做出正确的反映,且飞机也能够迅速按照指示执行,但是现实生活中这种情况很难实现。飞机驾驶员素质参差不齐、飞机设备和技术不够先进等等因素都会导致飞机事故的发生。
六、模型的推广
该模型对于避免飞机事故的各种研究可以推广到其他交通领域,正如最近发生的7.23高铁追尾事件,如果信息系统能够正常工作,人员素质较高,合理的调配必然能够有效避免不必要事故的发生。高铁系统黑匣子的设计原理正好符合飞机安全飞行相对于其他飞机飞行的闭塞区域,轨道设计也正符合飞机短区域内直线飞行的特点,而该模型对飞机飞行方向的调整正是对高铁安全轨道设计有一定的启发。当然,不仅交通领域可以运用此模型设计合理的运行路线对运河改道、维修公路等等方面都能有一定的运用。
七、参考文献
[1]作者王文波《数学建模及其基础知识详解》武汉大学出版社2006年出版
[2]作者1 姜启源作者2 谢金星《数学建模案例选集》高等教育出版社2006年出版
[3]作者1 冯杰作者2 黄力伟《数学建模原理与案例》科学出版社2009年出版
[4]作者章绍辉《数学建模》科学出版社2010年出版
[5]作者卓金武《MATLAB在数学建模中的应用》北京航空航天大学出版社2011年出版
[6]作者1 董文永作者2 刘进作者3丁建立《最优化技术与数学建模》清华大学出版社
2010年出版
附录一lingo:
model:
sets:
ii/1..6/:x,y,q,qt,xx,yy;
time/1..200/:t;
pair(ii,ii):d;
endsets
data:
x=150 85 150 145 130 0;
y=40,85,155,50,150,0;
q=243,236,220.5,159,230,52;
enddata
min=@sum(ii:qt^2);
@for(time(k):t(k)=0.2*2^0.5*k/300);
@for(pair(i,j)|i#lt#j:@for(time(k):(((x(i)+800*t(k)*@cos(q(i)*3.14159 /180+qt(i)*3.14159/180))-(x(j)+800*t(k)*@cos(q(j)*3.14159/180+qt(j)*3 .14159/180)))^2+
((y(i)+800*t(k)*@sin(q(i)*3.14159/180+qt(i)*3.14159/180))-(y(j)+800*t (k)*@sin(q(j)*3.14159/180+qt(j)*3.14159/180)))^2)>=64));
@for(ii(i):@abs(qt(i))<30);
end
附录2:
sets:
ii/1..6/:x,y,q,qt,xx,yy;
time/1..200/:t;
pair(ii,ii):d;
endsets
data:
x=150 85 150 145 130 0;
y=40,85,155,50,150,0;
q=243,236,220.5,159,230,52;
enddata
!min=@sum(ii:qt^2);
@for(time(k):t(k)=0.2*2^0.5*k/300);
@for(pair(i,j)|i#lt#j:@for(time(k):(((x(i)+800*t(k)*@cos(q(i)*3.14159 /180+qt(i)*3.14159/180))-(x(j)+800*t(k)*@cos(q(j)*3.14159/180+qt(j)*3 .14159/180)))^2+
((y(i)+800*t(k)*@sin(q(i)*3.14159/180+qt(i)*3.14159/180))-(y(j)+800*t (k)*@sin(q(j)*3.14159/180+qt(j)*3.14159/180)))^2)>=64));
@for(ii(i):@abs(qt(i))<30);
@for(ii:qt End 附录3: clc; clear all x=[150 85 150 145 130 0]; y=[40,85,155,50,150,0]; q=[4.2412 4.1190 3.8485 2.7751 4.0143 0.9076]; for k=1:6 for j=1:6 if(k>=j) b(k,j)=0; d(k,j)=0; a(k,j)=0; else z1(k,j)=x(j)-x(k)+i*(y(j)-y(k)); z2(k,j)=exp(i*q(k))-exp(i*q(j)); b(k,j)=angle(z2(k,j)/z1(k,j)); d(k,j)=sqrt((x(k)-x(j))^2+(y(k)-y(j))^2); a(k,j)=asin(8/d(k,j)); end end end b=b*180/pi a= 附录4: function f=myfun(qt); f=0; for k=1:6 f=f+abs(qt(k)); end function [g,ceq]=mycon(qt) x=[150 85 150 145 130 0]; y=[40,85,155,50,150,0]; q=[4.2412 4.1190 3.8485 2.7751 4.0143 0.9076]; for k=1:6 for j=1:6 if(k>=j) d(k,j)=0; a(k,j)=0; else z1(k,j)=x(j)-x(k)+i*(y(j)-y(k)); z2(k,j)=exp(i*(q(k)+qt(k)))-exp(i*(q(j)+qt(j))); d(k,j)=sqrt((x(k)-x(j))^2+(y(k)-y(j))^2); a(k,j)=asin(8/d(k,j)); end end end g=[a(1,2)-abs(angle((z2(1,2)/z1(1,2))));a(1,3)-abs(angle((z2(1,3)/z1( 1,3))));a(1,4)-abs(angle((z2(1,4)/z1(1,4)))); a(1,5)-abs(angle((z2(1,5)/z1(1,5))));a(1,6)-abs(angle((z2(1,6)/z1(1,6 ))));a(2,3)-abs(angle((z2(2,3)/z1(2,3)))); a(2,4)-abs(angle((z2(2,4)/z1(2,4))));a(2,5)-abs(angle((z2(2,5)/z1(2,5 ))));a(2,6)-abs(angle((z2(2,6)/z1(2,6)))); a(3,4)-abs(angle((z2(3,4)/z1(3,4))));a(3,5)-abs(angle((z2(3,5)/z1(3,5 ))));a(3,6)-abs(angle((z2(3,6)/z1(3,6)))); a(4,5)-abs(angle((z2(4,5)/z1(4,5))));a(4,6)-abs(angle((z2(4,6)/z1(4,6 ))));a(5,6)-abs(angle((z2(5,6)/z1(5,6))))]; ceq=[]; 附录6: %主函数 clc clear all qt0=[0 0 0 0 0 0]; vlb=[0 0 0 0 0 0]; vub=[30 30 30 30 30 30]*pi/180; [qt,fval]=fmincon('myfun',qt0,[],[],[],[],vlb,vub,'mycon'); qt=qt*180/pi; fprintf('各架飞机的偏转角为:\n'); qt 附录7:c++判断飞过去 #include using namespace std; #include #define PI 3.1415926535897932384626 float juli(float x1,float x2,float y1,float y2) { if(x2==y2&&y2==0) { return 100000000; } else{ float c; c=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); return c; } } float istime(float x,float y,float angle) { float c; if(x==y&&angle==0&&y==angle) return 100000000; else