飞行管理_数学建模

1、b题：飞行管理问题摘要：飞行管理问题是一个既现实又重耍的课题，本文利用偏转角度尽可能的小建 立两个非线性规划模型。模型一：时间模型。考虑到各架飞机的偏转角有正有负，在此模型屮，对于 各架e机调整选取各个偏转角的绝对值的和作为目标函数，要求任意两架0机任 意时刻的距离大于8公里，则可以求出任意两架飞机的距离如。此吋，两架飞机距离是时间/与各个飞机偏转角a0的函数，编写程序时将/离散化，且/有 最大值0. 2828 5 （沿对角线飞过的时间）,这样可得到表1t的结果:表1t飞机序号123456总和偏转角（度）002.540629001.0745263.615155模型二：闭塞区域模型。在两架飞机中

2、，将其屮一架看成“静止”，另一架 相对于它而运动。而以“静止”飞机为圆心，8如7为半径的圆形区域构成该飞机 的闭塞区域，任意一架飞机的方向角均不能在此区域内，则为不相撞。为此，木 文用复变函数的知识表示各架飞机的速度，从而算出相对速度，再求出相对位移, 以相对速度与相对位移的夹角大丁每两架飞机的临界夹角來刻画不相撞。冃标函 数为每架飞机偏转角的平方和。利用计算机编程得到表1-2的结果：表1-2飞机序号123456总和偏转角（度）002.5566001.06843.6250对于上述两个非线性规划，在理论方面，木文利用sutm内点法（障碍函数法） 进行算法描述，在操作方面，分别利用lingo语言与

3、matlab语言宜接编写程序 进行计算关键词：非线性规划、复变函数、sumt内点法、闭塞区域、禁飞角一、问题重述1. 背景知识与其他交通工具相比，飞机以其速度快、安全舒适等特点在交通领域占据了 绝对地位。而近年来飞机事故的频繁发生也预示着飞机存在一定的安全隐患。经 调杳造成飞机相撞事故的原因主要是人、飞机（设备）、环境，而人的因素是事 故屮通常起主体作用的因素，直接影响事故的发生和结局。飞机事故的发生难以 预测口死亡率极高，所以航空安全机制的健全，航空人员素质的提高已变得刻不 容缓。2问题重述在约10000米的高空某边长为160公里的止方形区域内，经常有若干架飞机 作水平飞行。区域内每架飞机的

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

5、算（方向角谋差不超过001度），要求飞机飞行方向角调整的幅度 尽量小。设该区域内4个顶点的坐标为（0, 0）, （160, 0）, （160, 160）, （0, 160）。记录数据为：飞机编号横坐标x纵坐标y方向角（度）115040243285852363150155220.54145501595130150230新进入0052 1二、基本假设1飞机在规定区域内飞行，不会发生任何自然灾害；2. 飞行人员能够探知飞行区域内所有飞机的飞行轨迹，并能够及吋作出相应调 整；3. 飞机在规定区域内最多只调整一次，且在第六架飞机刚进入区域时调整；4飞机在调整过程屮是瞬间完成，飞行轨迹可看成直线；5. 不

6、碰撞的标准是任何两架飞机的距离大于8如2 ；6. 飞机飞行方向的调整幅度不应超过30° ；7. 所有飞机调整角度总和最小为最优方案；&暂时不考虑飞机离开此区域的情况；9. 0机在调整之后按固定速度飞行；10. 飞机在飞行过程屮可看为质点。三、参数说明1. （心,兀）表示笫i架飞机的初始位置；2. 0表示第i架飞机初始飞行角度；3. 表示第i架飞机的调整角度；4. 5.表示第，架飞机和第丿架飞机的距离；5. ©表示第i架飞机和第/架飞机距离的平方；5. 岭表示笫i架飞机的飞行速度；6. /表示飞机的飞行吋间；7. 厶表示第/架飞机飞出飞行区域的时间；10. 0表示第/

7、架飞机和对于第j架飞机的飞行角度；11. z&表示第r架飞机的速度矢量；12. &表示第k架飞机指向第/飞机的飞行矢量；13. 叫/表示第r架飞机对第/架飞机的禁飞角的大小。四、问题分析第？三/第二架_飞机四廟1 1 1 1 114012010080架60 新进入的02040603010012014016020图4-1：飞机初始状态从距离约束条件来说，在初始状态下，我们可以通过调整飞行方向，口的 使飞机在飞行区域内满足任意两架飞机的距离都大于8公里。调整角度我们可以 假设已知，将上述初始角替换成调整角，通过计算，我们知道存在调整方法能够 使飞机安全飞出飞行区域，月调整方式多样。

8、若我们假设调整角度对于每架飞机 来说，其满意度是一样的，那在这一系列调整方式中，必定存在调整角度之和达 到最小的调整方式。在飞行过程屮，初始位置和飞行方向对飞机的飞行路径有决 定性作用，而在速度已知的情况下，时间的确定就可以直接决定飞机的具体位置。 由此，我们想到，如果取适当的一段时间为步长，确定检测时间点，若这些时间 点各架飞机都能满足不相撞的条件，那我们可以近似认为，飞机在这些飞行方向 上是符合安全飞行条件的，之后只需找出调整最小的方式即可。从角度约束条件来说，将飞机看做质点，在飞行区域内，以飞机为圆心，4 公里为半径的圆形区域内是不能存在其他飞机的。即每架飞机都冇其自身的闭寒 区域，而且

9、是动态的闭塞区域。在边长为160公里的正方形二维平面内，这6 个闭塞区域在运动过程屮是不允许存在重合部分的。飞机在飞行区域内是按照同 一方向飞行，也就是说在该正方形区域内，闭塞区域是按直线行走的。运动初期, 各圆的位置已知，飞行方向可调整（先不考虑调整角度的限制问题），而调整的 廿的是为了让这些圆在以相同速度且沿直线的行走过程中不存在重合部分。将模 型简化，仅看两架飞机的运动形式，以一个i员i的i员i心出发，引发出两条相对于另 一个圆的切线，而这两条切线在包含另一个圆的一侧所构成的夹域是禁飞域，夹 角为禁飞角o在这种解法中不仅要考虑飞机之间的相对运动速度还要考虑飞机之 间初始相对位置的具体关系

10、。讨论初期，我们确定飞机之间相对运动的速度还是引用原冇的坐标系，飞机z间相对位置关系也在原冇的坐标系屮讨论，该模型的 约束条件就是相对运动方向与相对位置方向之间的相对夹角只要满足不在禁飞 夹角屮即可。首先我们选择通过几何方法计算相对夹角，而在进一步计算过程屮, 我们发现相对夹角若以角度制表示要分不同种情况讨论，过程复杂难以计算。之 后我们又将模型进一步改进，对于坐标系的选取，若以飞机z间相对位置所在向 量确定相对坐标则使约束条件简化。在计算方而，为了避免角度的多情况彩响， 我们选取以复数形式表示角度，这样就将相对夹角的形式统一起来，计算也有了 进一步突破。在计算出具体结果后，就需要对结果的有效

11、性进行检验，而我们选取的检验 条件就是飞机以调整后的角度飞行，在离开飞行区域时，是否与飞行区域内的其 他飞机的距离大于60公里，若都满足大于60公里，则认为结呆可行性较强，若 有岀入，在对模型进行相应调整。五、模型建立与求解模型一设飞机初始状态下的坐标为(“)易知任意吋刻飞机的坐标为xi = xi + vt cos(- + ag)yi = yi + wsin(q +&,)则两架飞机的距离：dj = 5%- =( x 厂 x'2+(y厂岭y讨论飞机在规定区威内的飞行状态，则飞机冇飞行时间的限制：若飞机沿对 角线飞行，则飞行距离最远，飞行时间最长，此时/唤二型返，即t e(0,0.

12、2828), max 800在飞行时间内，如杲任意两架飞机的距离都超过8公里，则认为不会发生碰撞事 件，即：vr g(0,0.2828) di >64(i = l.6j = l6)(1)首先讨论如果不改变飞机偏转角是否有两架飞机相撞，即满足(1)式的当v = 800公里/小时，(兀,兀)和&心知的情况下，经计算：5362 =(150+ 800/cos220.5“ cos52丁 + (155 + 800fsin220.5“ sin52<=645562 = (130+ 800fcos230°cos52°f+(150+ 800/sin230°sin5

13、2°$ <= 64在t e(0,0.2828)内有实解，贝lj如果按照初始方向飞行会发生碰撞事件，需要做出相应调整。为了确保飞行时不会发生碰撞事件，则任意叨架飞机的飞行距离不得小于 8公里，在此条件下寻求6架飞机能够安全通过飞行区域的最小调整方式。而对 于“最小的调整方式”，本文采取下述两种理解方式。1 偏转总和最小。因为每架飞机偏转角可止可负，所以，本模型取偏转角的 平方和作为目标函数。结合上述分析与条件，可建立如下非线性规划模型6min工(aq尸1=1dq >64(z=b-6,j = l-6)皿|旳<30 (心16)虫(0,0.282*其中:2.毎架飞机偏转都最

14、小。min = tdj> 64t g (0,0.2828)编程时，我们将/离散化，分为200等分,hn .0.2828八八八八八彳即 dt = 0.00094300误差 ds = vdt = 0.7541如？（1）求得的最终结果为：q=0； &2=0； 03 =2.540629°； a4 = 0； 比=0； 06 = 1.074526° o(2) 求得的最终结果为：= 2.2906° ； 込=0.2500° ；= 2.5406 ； a% = 0 ； kd, = 0.2500° ；佻=10745对比上述结果，结合题屮误差在0.01。

15、之内，可只取第一种方案，这样既满 足飞机偏转总和最小，又满足每架飞机偏转最小。下面检验对于上述结杲是否满足“进入该区域的飞机在到达区域边缘时，与 区域内飞机的距离应在60公里以上；”飞机的飞行路径方程为:（y _ x）/（兀xi）= tan（q + aq）；求得每架飞机飞出飞行区域的坐标分别为（129.6190,0）； （27.666&0）； （0,15.1227）；（0,105.6603）；（4351,0）； （120314,160）飞机飞岀飞行区域的时间为：_（兀一兀j+（”一儿尸v求得每架飞机的飞行时间为：斤=0.0561 ； t2 = 0.1282 ；= 0.2564 ；t4

16、=0.1941； t5 = 0.2448 ； t6 = 0.2501。飞机飞出飞行区域的先后顺序为：1, 2, 4, 5, 6, 3o按照此飞行时间顺序 由少到多重新排序。而飞机飞出区域要保证和飞行区域内的其他飞机距离都人于60公里；所以：u1 -认）2 + （儿-儿s）2）> 3600（心 1 6, j = z +1 6）b|j：（兀 x> vti cos0j + （）% yj - vti sin 0> 3600成立，根据上述思想编写c+程序（见附录）经检验。只当第5架飞机0出区域是与第3架0机的距离小于60公里，再 缩短时间步长使结果更加精确。模型二先将模型简化为两架飞机

17、，考虑两架飞机的相对运动，若保证飞机在飞行区 域内不发生碰撞则相对飞行如度有一定的范围限制。以第丿架飞机坐标为圆心， 以8为半径作圆，此圆即为第丿架飞机的闭塞区域。第鸟架飞机相对于第丿架飞 机的速度为%o由第k架飞机坐标引出的相对于此圆的两条切线的夹角叫为第 k架飞机做相对运动的禁飞方向。规定：0灯为第£架飞机相对于第/架飞机的夹介：以r为原点，k»为极轴/旷 将极轴几转至的角度即为0“逆 时针为正，顺时针为负。（按照这种规定方法，求出的值与matlab语 言中angle命令产生的结果一致）设 |v| = a则第r架飞机的速度矢量：乙=人/偸+现）则第r架飞机相对于第j架飞

18、机的相对速度矢量:vkj = zk -zj =严第£架飞机与第丿架飞机相对位置矢量（极轴）:fkj =xj xk +i（儿-儿）则按照上述规加：v,.a； = arg(4j kj由图1-2知:8sin如=kj s、kj以各个偏转角绝对值的和作为目标函数建立如下非线性规划:6min£|a.|z=1pk akj 仇=1 6j = 16)tt岡驚（日6）对于此类非线性规划，通常运用罚函数法。罚函数法基木思想是通过构造罚函数把约朿问题转化为一系列无约朿最优化 问题，进而用无约朿最优化方法去求解这类方法称为序列无约束最小化方法简 称为sumt法.考虑问题:其一为sumt外点法，其二为

19、sumt内点法.本文运用sumt内点法（障碍函 数法）：p.l g£x)0 i = 1,2,.,/n设集合d°=xi&(x)>(u = 1,2,mh0是可行域中所有严格内点的集合。构造障碍函数加"i 1/(xd /(x,r)=/(x)+rlng,.(x)或 z(x,r) =/(x) + rj； 篙/=! g(x)m仍1其屮称专ngg或客硕为障碍项，为障碍因子，这样问题（1）就转化为求一系列极值问题:min内点法的迭代步骤:s1：给定允许误差£0,取作0,0<0<1；s2：求岀约束集合d的一个内点x°gd°,令

20、k = ；s3：以xa-'ed°为初始点，求解其中x e d0的最优解，设为s4：检验是杏满足"或"，满足，停止迭代,否则取心1=07，令k=k + l,返回s3.根据题意可取£ = 0.01。，ix（i）w（0,彳）,（心16）的六维样木空间，根据算法进行迭代。考虑到实现此算法的复杂性以及当今计算机的发展，本文采取直接运用 matlab编程得到如下计算结果：q=0； 2=0； 03=2.540629°； a4 = 0； a5 = 0； a6 =1.074526五、模型的评价与改进模型一优点：（1）将约束条件简化，以条件检验结果，使模型

21、简单化；（2）该模型比较传统，便于对模型结构的理解和把握；（3）将时间作为己知条件，选取特圧时间点进行条件约束，计算量大大减少； 缺点：（1）对时间步长的选取要求较高；步长选取不当直接影响结果的有效性；（2）计算软件的局限性会导致结果误差较大。模型二优点：（1）将距离z间的关系转化为角度z间的关系，使模型更加直观；（2）相对运动和相对位置的选取有利于约朿条件的进一步简化；（3）运用复数形式表示矢量，避免讨论多种情况；（4）闭塞区间的引用使模型更加形象，便于理解；（5）内点法计算最优化问题使约束条件在结果中的选取也更具优越性； 缺点：（1）对介度的处理要求较高。对于该模型的讨论，我们只考虑了人为

22、因素的影响，具有一定的局限性。我 们假定人员的素质极高，能够按照指示迅速做出正确的反映，旦飞机也能够迅速 按照指示执行，但是现实生活中这种情况很难实现。飞机驾驶员素质参井不齐、 飞机设备和技术不够先进等等因素都会导致飞机事故的发牛。六、模型的推广该模型对于避免飞机事故的各种研究可以推广到英他交通领域，正如最近发 生的7. 23高铁追尾事件，如果信息系统能够止常工作，人员素质较高，合理的 调配必然能够有效避免不必要事故的发牛。高铁系统黑匣子的设计原理正好符合 飞机安全飞行相对于其他飞机飞行的闭塞区域，轨道设计也正符合飞机短区域内 直线飞行的特点，而该模型对飞机飞行方向的调整正是对高铁安全轨道设计

23、有一 建的启发。当然，不仅交通领域可以运用此模型设计合理的运行路线对运河改道、 维修公路等等方面都能有一定的运用。七、参考文献1作者王文波数学建模及其基础知识详解武汉大学出版社2006年出版2作者1姜启源作者2谢金星数学建模案例选集高等教育出版社200611!版 作者1冯杰作者2黄力伟数学建模原理与案例科学出版社2009年出版4作者 章绍辉数学建模科学出版社2010年出版作者 卓金武matlab在数学建模中的应用北京航空航天大学出版社2011年出版6作者1董文永作者2刘进作者3 丁建立授优化技术与数学建模清华大学出版社 2010年出版附录一 i i ngo：model:sets :ii/1.6

24、/:x,y,q,qtzxx,yy;time/i.200/:t;pair (iir ii) :d;endsetsdata :x=150 85 150 145 130 0;y=40, 85, 155, 50, 150, 0;q=243z 236z 220.5, 159, 230, 52;enddatamin=sum(ii:qta2);for(time(k) :t (k)=0.2*2a0.5*k/300);for (pair (i,j) |i#lt#j:0for (time(k) : ( ( (x(i)+800*t(k)*0cos(q(i)*3.14159 /180 + qt(i)*3.1415 9

25、/180)-(x(j)+800*t(k) *cos(q(j) *3.14159/180 + qt(j) *3 .14159/180)a2+(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)a2)>=64);for (ii(i) :abs(qt(i)<30);end附录2：sets :ii/1.6/:x,y,q,qt,xx,yy; time/1200/:t;pair (ii,ii) :d;endsetsd

26、ata :x=150 85 150 145 130 0;y=40, 85, 155, 50, 150, 0;q=243z 2 36, 22 0-5, 159, 2 30, 52;enddata!min=sum(ii:qta2);min=z;for(time(k):t(k)=0.2*2a0.5*k/300);for(pair(izj)代ime(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)a2+(y (i)+800*

27、t(k)*sin(q(i)*3.14159/180 + qt(i)*3.1415 9/180)-(y(j)+800*t (k)*sin(q(j)*3.14159/180+qt(j)*3.14159/180)a2)>=64);for (ii(i) :abs(qt(i)<30);for (ii:qt<z);end附录3：clc;clear allx=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=l:6for j=l:6if (k>=

28、j)b (k,j)=0;d(kz j)=0;a(k, j)=0;elsezl(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(kz j)/zl (k,j);d(k, j)=sqrt ( (x(k)-x(j)a2+(y(k)-y(j)a2);a(k, j)=asin(8/d (kz j);endendendb=b*180/pia=附录4：function f=myfun(qt);f=0;for k=l:6f=f+abs(qt (k);end附录5：function g,ceq=mycon(qt)

29、x=150 85 150 145 130 0;厂40, 85, 155, 50, 150, 0;q=4.2412 4.1190 3.8485 2.7751 4.0143 0.9076;for k=l:6for j=l:6if (k>=j)d(kz j)=0;a(k, j)=0;elsezl(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) )a2+(y(k)-y(j)a2);a(kz j)=asin(8/d(kz j);endendendg=a

30、(1,2)-abs(angle(z2(1,2)/zl(1,2);a(1,3)-abs(angle ( (z2 (1,3)/zl ( 1,3) );a(1,4)-abs(angle(z2 (1,4)/zl (1,4);a(1,5)-abs(mngle(z2(1,5)/zl(1,5);a(1,6)-abs(mngle(z2(1,6)/zl(lz 6 );a(2,3)-abs(angle( (z2 (2z3) /zl (2z3);a(2z 4)-abs(angle(z2(2,4)/zl(2,4);a(2,5)-abs(angle(z2 (2,5)/zl (2,5);a (2z 6) -abs (an

31、gle (z2(2,6)/zl(2,6);a(3,4)-abs(angle(z2(3,4)/zl(3,4);a(3z 5)-abs(angle(z2(3z 5)/zl(3z 5 )；a (3,6)-abs(angle(z2(3,6)/zl(3,6);a(4,5)-abs(angle(z2(4z 5)/zl(4z 5);a(4,6)-abs(angle(z2 (4,6)/zl (4,6)；a(5,6)-abs(angle(z2(5,6)/zl(5, 6);ceq=;附录6：労主函数clcclear allqt0=0 00000;vlb=0 00000;vub=30 30 30 30 30 30*

32、pi/180;qt,fval=fmincon(fmyfun1z qto, , , , ,vlb,vub, 1mycon *);qt=qt*180/pi;fprintf (»各架飞机的偏转角为:n»);qt附录7： c+判断飞过去include <iostream>using namespace std;#include<cmath>#define pi 3. 1415926535897932384626float juli (float xl, float x2, float yl, float y2)if (x2=y2&&y2=0)

33、return 100000000;else float c;c=sqrt (xl-x2)* (xlx2) + (yl-y2)*(yl-y2);return c;float istime(float x,float v,float angle)float c;if (x=y&&angle二二0&&y二二angle)return 100000000;elseif 仃80/pt* (2*piatan(y/(160-x)<angle&&angle<3601|0<angle&&angle<180 /pt*atan(y

34、/(160-x)c= (160-x)/800/fabs(cos(pi*angle/180);return c;elseif (180/3. 14*atan (y/(160-x) <angle&&angle< (180t80/pi*atan (160-y) /x) )c=(160-y)/800/fabs (sin (angle*pi/180);return c;elseif (180-180/pi*atan(160-y)/x)<angle&&angle<(180+180/pi*atan(y/x) c=x/800/fabs(cos(angle*pi/180);return c;elseif (180+180/pi*atan (y/x)<angle&&angle<180/3. 14*(2*pi-atan(y/(160-x) )c=y/800/fabs(sin(anglc*pi/180);return c;elsereturn 160/800/sin (angle*pt/180);bool ischujie()rcturn 1;void main()floata6 = 150