飞行管理数学建模(共16页).doc

精选优质文档-倾情为你奉上摘要近年来，随着现代航空运输不断发展，为了维护航空器的航空秩序，保证飞机飞行安全，对于同一区域的飞行管理问题提出了要求。本文讨论了在一定区域范围内飞机飞行管理的最优化问题，通过建立数学模型计算求解，对飞机是否发生碰撞冲突进行预测，根据计算机求解结果对如何解脱冲突给出了较好的解决方法。对于飞机是否发生碰撞冲突问题，本文提出了基于飞机位置速度矢量关系的碰撞冲突检测方案，证明了只有位置差与速度差矢量内积小于零，即这样的航迹才存在潜在碰撞冲突,并根据安全飞行间隔规定,采用线性预测方法对冲突进行有效性确认,解决了飞机碰撞冲突检测的同时也避免了碰撞虚警问题。在此基础上，对于存在潜在碰撞冲突的飞行问题，运用航向调整的方法解脱冲突，建立非线性数学模型通过引入新的决策变量、，将原来的非线性模型转换成线性模型 其中，。再运用LINGO11编程求得该模型最优解为，第3架飞机的调整角为，第6架飞机（新进入的飞机）的调整角为，其余飞机不进行调整，从而给出了冲突解决方案。之后，本文对计算结果做出了分析和评价，同时还分析了滞后时间和转弯半径和限定在区域范围内对飞机航向调整的影响，使问题更符合实际情况。在对模型进行评价与分析的同时，本文又对模型进行了推广，对速度不同、飞行高度不同的情况下进行了分析，并给出了合理的解释；增强了模型的实际应用意义。关键词：飞行管理 碰撞冲突 线性规划专心-专注-专业一问题重述本题主要分析了在同一高度，一定范围内的飞行管理问题。在10000米高空、边长为160公里的正方形区域内最多有6架飞机做水平飞行，其中飞机以每小时800公里的速度匀速飞行。为了便于飞行管理，在每一架飞机刚刚进入此飞行区域边界时，飞机的位置和速度均由计算机记录其数据，并立即进行计算和判断是否会与区域内的其它飞机发生潜在碰撞，如果存在碰撞危险，就应该计算如何调整各架飞机的飞行方向角度，以避免发生碰撞。因此，根据题目的条件和假设，对该避免碰撞的飞机管理问题建立数学模型，列出计算步骤并对已经给出的数据进行计算。在保证进入该区域的飞机到达区域边缘时, 与区域内飞机的距离在 60 公里以上的同时，要求满足飞机飞行的方向角度调整幅度不能超过30度，而且要尽量小。在该区域内建立直角坐标系，4个顶点的坐标为 (0,0)， (160,0)， (160,160)， (0,160)。 记录数据如下：飞机编号横坐标x纵坐标y方向角(度)1150140243285852363150155220.54145501595130150230新进入0052最后，对建立的模型进行评价和推广。二问题分析2.1、问题分析针对问题1，一架新飞机飞入飞行区域边缘时，计算机进行记录其数据并进行判断该飞机是否和区域内的其它飞机存在潜在碰撞危险，由此作出分析，若存；在相互碰撞危险就要对各架飞机的方向角进行调整；否则维持原状，不做任何变化。 对两飞机间距离的理解在本题中，根据题目假设所给的条件即在安全飞行过程中任意两架飞机间的距离要在8公里之上，这个距离远大于飞机本身的尺寸大小，所以，在以后的问题求解过程中我们将飞机当作质点看待。这样，根据计算机记录的数据可以知道每架飞机的坐标，根据两间距离公式即可得到两飞机之间的距离。对飞机方向角与调整角的理解根据题目得知，飞机的方向角即为飞机的飞行方向与x轴正向的夹角。对于两架飞机、而言，要求飞机在飞行过程中相对速度的方向不能在飞机与飞机的碰撞角（即两切线交角中指向圆心的那个角）范围内。由于调整角以两飞机的连线为对称轴，左右由正负之分，故以调整角的绝对值最为目标函数。2.2、问题2分析 针对问题2的论述，由于在2.1.2中提出的绝对值目标函数为非线性函数，给求解带来了麻烦，因此应用变量代换的数学方法将其转换成非线性函数，运用MATLAB、LINGO软件编程，调用求解函数得出最优调整方向角度。三模型的基本假设1.新飞机进入边缘时，立即做出计算，每架飞机按照计算机计算后的指示立即作出方向角改变（有的飞机方向角可不变）。2.每架飞机在新飞机刚进入到下一架飞机要进入前整个过程中最多只改变一次方向角。3.忽略飞机转向角度，即认为飞机在接收到指令后立即对方向进行调整，且忽略调整时间。4.新飞机进入该区域前,在区域中飞行的飞机方向已调合适，不会相撞。5.对方向角的相同调整量的满意程度是一样的，且方向角调整越少，满意程度越高。6.飞机在安全飞行是两两之间的最短距离为8公里，此距离远大于机身尺寸，因此将飞机当作质点看待，即假设两飞机之间的距离为两点之间的距离。7.最多考虑6架飞机。8.飞机飞行方向角调整的幅度不应超过30度。9.所有飞机匀速飞行，速度为。说明： 假设6中假设飞机为质点是为了将问题简化，忽略次要问题以便得出合理的结论。假设7中6架飞机的假设是足够多的。以世界上最繁忙的国际航空港之一希思罗机场邻近区域为例，因假设飞机在区域160×160作水平飞行，即知该区域内无机场。设在希思罗机场起降的飞机有一半穿过该区域，希思罗机场起驾总架次为22.5万次，则平均每小时有15架飞机穿过该区域。而一架飞机穿过该区域最多需小时，则任一时刻该区域上空飞机架数的期望值不超过架。另外事实上不同飞机的飞行高度是不同的，这就进一步减少了该区域同一水平面上飞机的数目。以上讨论虽然略显粗略，但是足以证明架飞机的假设是合理的。 假设9中假定所有飞机飞行速度均为，是出于对问题的简化。我们将在模型的推广中给出飞机速度各不相同的方案。四.符号说明：第 架飞机与第架飞机的碰撞角，即两圆公切线交角中指向圆的那个角，规定；：第架飞机相对于第架飞机的相对飞行速度；：第架飞机与第架飞机间距；：第架飞机相对于第架飞机的位置矢量与x轴的夹角：第架飞机相对于第架飞机的相对飞行速度与x轴的夹角。：第架飞机相对于第架飞机的相对飞行速度与两架飞机圆心连线的交角，规定：以第架飞机为原点，连线从指向为正方向，逆时针旋转正，顺时针旋转为负；：第架飞机相对于直角坐标系旋转的角（即方向角改变量），为代数量；：第架飞机相对于第架飞机的改变量；：第架飞机的当前位置矢量；：第架飞机的当前速度矢量；：时间参数;: :五.模型的建立5.1 预测是否存在碰撞危险 潜在碰撞危险检测方案当相对速度矢量与相对位置矢量的矢量内积小于时，即从到夹角余弦值小于时，两飞机存在潜在碰撞的危险，否则，无碰撞危险。具体证明如下：对任意两架飞机，用位置矢量和速度矢量分别表示为：、，则经时间后两飞机的位置差为： （1）则 （2）为了求得两飞机相撞时刻，令：=得到 （3）当>0时，说明两飞机具有潜在碰撞的危险，而当<0时，则说明两飞机处于分离状态，没有碰撞的可能。而>0，即时，两飞机存在潜在碰撞冲突；反之，<0，即时，两飞机处于分离状态，无碰撞的可能。当对新进入飞机和区域内每一个已经飞行的飞机采用上述算法判断时，其相对速度矢量和位置矢量满足定理时，进行如下确认，否则认为满足安全飞行规则。当新进入飞机和区域内正在飞行的飞机采用上述算法时，其相对速度矢量和位置矢量满足以上条件时，进行如下确认，否则认为满足安全飞行要求。 碰撞冲突的确认算法 对检测到具有潜在碰撞的飞机对，必须根据安全飞行间隔规定进行确认。飞行间隔规定的最小水平安全间隔就相当于在每架飞机周围设置一个保护区，当某一架飞机保护区内有另外一架飞机进入时，就会发生碰撞，因此定义保护区为以飞机当前位置为中心，最小水平间隔（8）为半径的圆形区域，即：其中，表示到航空器当前位置的距离。如果,说明两航空器已经处于飞行间隔不允许的情况,即告警状态,两航空器已经处于飞行间隔所禁止的情况。反之，如果,说明两航空器当前没有违反安全飞行间隔规定，但由两飞机之间的飞行位置- 速度矢量关系判断，具有潜在的冲突可能，因此必须进行下列预测和分析。 线性预测后的位置与安全间隔的关系对于当前航空器对之间的距离大于安全间隔距离，处于安全情况下的飞机对，进行适当时间T的线性预测，即根据飞机对当前的飞行速度，预测飞机在未来一段时间内位置信息与安全间隔的关系，判断与之间的关系(根据本文定理,处于潜在冲突时, 故取模后用两者之差)5.2 建模方案 为了研究两飞机相撞问题，采用相对速度作为研究对象，因为飞机是否相撞的关键取决于相对速度的方向。 模型在分析碰撞问题中的运用，如图1示。8ivj图1：分析碰撞模型示意图 模型建立的函数及运用方程 线性规划模型的建立根据相对运动的观点在考察两架飞机i和j的飞行时，可以将飞机i视为不动而飞机j以相对速度 （4）相对于飞机i运动，（4）式进行适当的化约可得: （5）不妨设，此时相对飞行方向角为，如图1。由于两架飞机的初始距离为, （6）则只要当飞行方向角满足 （7）时,两架飞机不可能相撞（见图1）。由于是飞行方向角，应满足，但是由给出的，注意到，因此有 （8）有可能超出，此时只需将式（7）修改为 (9)即可，此时第i架飞机与第j架飞机不碰撞的条件就成为 (10)或 (11)这两个条件均为线性的。这样一来，各架飞机两两不碰撞条件也是关于的线性不等式，这个问题就可划为线性规划问题。 目标函数的建立根据的含义，可计算出- (12)因此，第i架飞机与第j架飞机不碰撞的条件可表示为 (13)为保证飞机的安全飞行并且考虑乘客乘坐的舒适度，因此，决策目的是使区域内所有飞机的偏转角绝对值之和最小，在忽略计算时间和转向时间的情况下，可得目标函数如下： ， (14) 目标函数的转化为将上述非线性目标函数转化成线性目标函数，引入新的决策变量、，令 , , (15)根据题中规定，方向角改变量在范围内，即 (16)； (17) (18) (19)因此，将(15)式代入(14)式可得线性目标函数 (20) 目标函数及约束条件 (21) , (22) (23) 六.模型求解计算步骤：1.记录各飞机状态：位置（坐标）及速度大小和方向角。2.用MATLAB求解出第 架飞机与第架飞机的碰撞角， (24) (25) 求解 的MATLAB 程序以“arfa”为文件名保存，具体程序见附件（1）。 = 0 0.0941 0.5625 0.0889 0.3659 0.0390 0 0 0.0838 0.1154 0.1014 0.0666 0 0 0 0.0762 0.3985 0.0371 0 0 0 0 0.0792 0.0522 0 0 0 0 0 0.0403 0 0 0 0 0 03. 根据公式（12）用MATLAB求解出第 架飞机相对于第架飞机的速度的夹角，其中位置矢量与横坐标的夹角为分段函数，如下： (26) 求解的MATLAB 程序以“beita”为文件名保存，具体程序见附件(2)。 = 0 -1.3131 0.8247 0.3435 2.9420 0.1741 4.9701 0 4.7321 -0.7373 4.6722 0.1571 0.8247 -1.5511 0 0.2178 -1.0260 0.0054 0.3435 -0.7373 0.2178 0 0.1042 -0.0615 -3.3412 -1.6110 5.2572 0.1042 0 0.03340.1741 0.1571 0.0054 -0.0615 0.0334 04.调用程序=mod(,2*pi )，将转换成在范围内的正数角，其弧度值为：= 0 4.9701 0.8247 0.3435 2.9420 0.1741 4.9701 0 4.7321 5.5459 4.6722 0.1571 0.8247 4.7321 0 0.2178 5.2572 0.0054 0.3435 5.5459 0.2178 0 0.1042 6.2217 2.9420 4.6722 5.2572 0.1042 0 0.03340.1741 0.1571 0.0054 6.2217 0.0334 05.建立线性规划模型,用LINGO11进行求解（程序见附件1）。将以上求得的的结果输入程序LINGO Model-1995A中的beita的值，经计算求得结果如下：Global optimal solution found. Objective value: 0.E-01 Infeasibilities: 0. Total solver iterations: 2 Variable Value Reduced Cost A( 1) 0. 0. A( 2) 0. 0. A( 3) 0. 0. A( 4) 0. 0. A( 5) 0. 0. A( 6) 0. 0. MING( 1) 0. 1. MING( 2) 0. 1. MING( 3) 0.E-01 0. MING( 4) 0. 1. MING( 5) 0. 1. MING( 6) 0.E-01 0. NING( 1) 0. 1. NING( 2) 0. 1. NING( 3) 0. 2. NING( 4) 0. 1. NING( 5) 0. 1. NING( 6) 0. 2.结果显示第3架飞机的调整角为（0.0496rad），第6架飞机（新进入的飞机）的调整角为（0.0138rad），最优目标函数为（0.0634rad）七.结果分析1、结果的实际价值在六中目标函数求解结果显示，并根据题目方向角误差小于0.01度的要求，可知，第3架飞机的调整角为，第6架飞机的调整角为，其余各架飞机不许进行调整。根据结果分析可知调整的飞机数目与角度均与实际情况吻合。加之，在实际情况中飞机的飞行高度、速度等条件均可改变，故该模型的计算结果符合实际情况。根据计算软件显示，运用该模型求解结果用时极短，并且随着计算机技术的进步，计算机的计算速度也愈来愈快，故运用该模型可实现实时控制。2、误差分析 模型建立中的误差考察模型假设可知假设3带来一定误差, 简单分析后可知考虑转弯半径使飞行轨迹在垂直飞行方向上产生一个偏差X (27) 其中R为转弯半径，为转角。对最小转弯半径小于10公里的中小型飞机，转角小于15度时偏差公里对最小转弯半径约为40公里的大型飞机，转角小于8度时偏差公里两种情况下偏差均小于相撞条件8公里的5%，可以忽略不计。由于飞行管理中让一架飞机做大于8度的转向是较少发生的事儿（），故该处假设可以认为是合理的。下面给出发生极端情形时的一个对策：设一架大型飞机做30度的转弯,此时偏差公里不可忽略。在飞机开始转弯后（转弯过程约需1.6分钟）地面控制站的计算机可算出一条曲率处处大于的偏转线，以此曲线引导飞机飞回偏转线。这种局部调整法的优点在于计算简单，不对全局的调整角度方案产生影响，适用于这种小概率事件。 计算过程中的误差仅来源于机器的截断误差，对于最后结果没有重大影响。 输入误差对结果影响通过对输入数据给以小的扰动，经大量计算可得如下定性结果:.飞机位置有小的扰动对结果影响甚微，扰动增大时影响明显。.飞机方向角的扰动对结果影响显著，扰动大小与结果变化幅度基本是在一个数量级上。以上结果说明为保证结果的准确性应尽量减小对飞机方向角测量的误差, 同时对飞机位置测量的误差也应控制在一定范围以内。（4） 该模型未考虑区域限制带来的偏角误差。 八.模型的与改进根据实际情况，在同一区域飞行的飞机飞行速度完全相同，而且全部匀速飞行的情况，在实际中属于特殊情况，因此本题给出的条件可以看作理想条件。然而上述计算结果对现实仍具有指导意义，下面我们将从飞机的飞行速度不同和飞行高度不同两种情况对模型进行推广。1、飞行速度不同针对速度不同但速度保持不变的情况，上述的推导方法已不完全适用，尤其是对于（5）式的推到方法已经不适用，为此我们给出了新的计算方法，如下：第架飞机的速度矢量以复数形式表示为 (28) 相对速度 (29) 则有 (30) 即 (31) 根据上述推到结果并结合图1，不难得出两飞机持续满足最小安全间隔的充要条件为：若，则需改变以实现冲突规避。在MATLAB中计算出来后，依据该模型，可计算出飞机偏转角最小的最优解。2、飞行高度不同在实际情况中，根据飞机航道设计的不同，在同一区域内的飞机可以在不同的高度进行飞行。这样一来，就使得同一区域内的飞机数量减少，同时也给飞机改变速度和改变方向提供了更广阔的空间，因此可以给飞机的安全飞行提供更多的实现方式，给飞行管理提供了更多的可行方案，更易于飞行管理的实现。九. 参考文献1 王运锋、贺文红、刘健波，基于航空器位置速度矢量关系的短期冲突检测算法，四川大学学报（工程科学版），2009年9月，第41卷第5期：1361412 牟奇峰，空中交通管理中的防撞策略问题探究：38393 谭永基、蔡志杰、俞文，数学建模，上海：复旦大学出版社，2005，62634 谭浩南、朱正光、刘剑，飞行管理问题的实时计算，数学的实践与认识，1996年，第26卷第1期：17附件：（1）function arfa=j1995iao(x,y)r=8;L=length(x);arfa=zeros(L,L)for i=1:L for j=1:L if i=j L(i,j)=sqrt(x(i)-x(j)2+(y(i)-y(j)2); arfa(i,j)=asin(r/L(i, j); end endEnd备注：程序中arfa是论文中角。（2）functionA=feixingjiao5(theta,X,Y)for i=1:6theta(i)=theta(1,i);endfor i=1:6 for j=1:6 if i=j if Y(j)-Y(i)>0 if X(j)-X(i)>0 beita(i,j)=atan(Y(j)-Y(i)/(X(j)-X(i); end if X(j)-X(i)<0 beita(i,j)=atan(Y(j)-Y(i)/(X(j)-X(i); beita(i,j)=beita(i,j)+pi; end if X(j)-X(i)=0 beita(i,j)=pi/2; end end if Y(j)-Y(i)<0 if X(j)-X(i)>0 beita(i,j)=atan(Y(j)-Y(i)/(X(j)-X(i); beita(i,j)=beita(i,j)+2*pi; end if X(j)-X(i)<0 beita(i,j)=atan(Y(j)-Y(i)/(X(j)-X(i); beita(i,j)=beita(i,j)+pi; end if X(j)-X(i)=0 beita(i,j)=3*pi/2; end end end endendA=zeros(6,6); for i=1:6 for j=1:6 if theta(i)>theta(j) beita(i,j)=beita(i,j)*180/pi; A(i,j)=90+(theta(i)+theta(j)/2-beita(i,j); end if theta(i)<theta(j) beita(i,j)=beita(i,j)*180/pi; A(i,j)=(theta(i)+theta(j)/2-beita(i,j)-90; end A(i,j)=mod(A(i,j),360); A(i,j)=A(i,j)*pi/180;endend备注：矩阵A为第架飞机与第架飞机的相对速度与位置矢量之间的夹角，论文中记为（3）model:sets:xx/1.6/:a,ming,ning;xx1/1.6/:;links(xx,xx1):arfa,beita;endsets!目标函数;min=sum(xx(I):ming(I)+ning(I);!防碰撞约束;for(xx(i)|i#le#6 #and# i#ge#1:for(xx1(j):(ming(i)-ning(i)+ming(j)-ning(j)/2+beita(i,j)>arfa(i,j);(ming(i)-ning(i)+ming(j)-ning(j)/2+beita(i,j)<2*pi-arfa(i,j););!变量范围约束;for(xx:ming>0);for(xx:ning>0); for(xx:ning<3.14/6);for(xx:ming<3.14/6);!数据;data:data:beita = 0 -1.2346 -2.2375 -2.7196 -0.1210 0.2526 0 0 -1.5502 -0.7373 -1.6110 0.1571 0 0 0 -2.9230 -1.0251 0.0063 0 0 0 0 0.1042 -0.0615 0 0 0 0 0 0.0334 0 0 0 0 0 0 ; arfa = 0 0.0941 0.5625 0.0889 0.3659 0.0390 0 0 0.0838 0.1154 0.1014 0.0666 0 0 0 0.0762 0.3985 0.0371 0 0 0 0 0.0792 0.0522 0 0 0 0 0 0.0403 0 0 0 0 0 0;enddataend