飞行器管理问题

飞机碰撞调整优化模型

摘要

本文研究的是在一定区域内水平飞行的n(n<=6)架飞机的飞行管理问题。本文通过建立一个非线性规划模型来求解出飞机飞行的调整角度。

首先，根据题目中“要立即计算并判断新进入飞机是否会与区域内的飞机发生碰撞“的要求，我们在编程给出判断是否碰撞的算法之后，需要求出进行判断碰撞所需要的时间。我们运用题目中所给假设中的条件作为约束条件，通过判断在区域内任两架飞机之间在区域内飞行时任意时刻的距离是否小于8公里作为碰撞标准，进行判断是否会发生碰撞。利用MATLAB中tic和toc函数进行计时。以题目中飞机的初始状态求解，实验得出判断结果的用时均在0.3s以内，且存在飞机会发生碰撞。

为了得出飞机的调整角度，我们首先运用蒙特卡洛法随机产生200组飞行调整角度。其中每组数据包含6个数据项，分别作为架飞机的调整角度。然后，我们通过碰撞条件筛选出这200组数据中能使飞机进行角度调整后不发生碰撞的数据。继而将上一步中所得数据中的调整角度绝对值之和最小的那组数据作为

进一步优化的对象。接着，我们以o

0.01为步长对上一步得到的优化对象按飞机编号顺序进行逐角优化，使每一个调整角度逐步减小。每优化一步进行一次碰撞判断，直到该飞机的调整角度绝对值为0或将导致飞机发生碰撞为止，然后再按同样的算法优化下一个飞机的调整角度。从而得出近似最优解。通过模型求解后

4.5之内。

结果显示各飞机调整角度绝对值之和均在o

最后，为了对我们所建的模型进行推广，说明我们所建的模型是具有普适性的，我们又自行随机设计了两组符合初始状态要求即进入该区域的飞机在到达该区域边缘时，与区域内飞机的距离应在60公里以上的数据来作为6架飞机的初始状态，并利用我们所建的模型得出了调整角度结果。实验结果比较理想，各飞机调整角度绝对值之和均在o

1.5之内。

关键词：飞行管理判断调整蒙特卡罗逐步优化

一、问题重述

飞行管理问题在现代生活中是一个很值得考虑和解决的问题。

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

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

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

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

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

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

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

设该区域4个顶点的坐标为：（0，0），（160，0），（160，160），（0，160）。

试根据实际应用背景对所建模型进行评价和推广。

二、模型假设与符号说明

2.1 模型假设

1）模型的目标是在满足各架飞机不碰撞的条件下，每架飞机的调整角度绝对值之和最小；

2）不碰撞的标准为任意两架飞机在区域内飞行的任何时刻，它们之间的距离大于8 公里；

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

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

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

6）考虑至多六架飞机的飞行管理；

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

8）在模型中先运用蒙特卡洛法进行初选时，我们令调整角β～()20，

N 分布。在“建模前准备”中我们给出了解释。最后模型结果分析中，也对此假设的进行

了验证并通过。在优化算法中，我们还会对调整角度进行优化，以求更接近“调整角度尽可能小”的建模目标； 9）忽略飞机转换角度的时间；

10)假设飞机在进行角度调整之后将沿直线水平飞行且速率保持800km/h 不变。 2.2 符号说明

三、问题分析

本文研究的是在一定区域内水平飞行的n(n<=6)架飞机的飞行管理问题。通过初步分析可知，我们可以通过建立一个非线性规划模型求解出飞机飞行的方向角应该进行如何合理的调整来满足在该区域内的飞行满足要求。

首先，我们需要考虑的是碰撞判断的算法的即时性。题目要求，当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，需要立即计算是否会与区域内的飞机发生碰撞，这就要求实现该判断算法的即时性。这样我们就需要编程给出判断是否碰撞的算法并记录其花费时间。判断是否碰撞的思路是求出第六架（即新进入）飞机与其他五架飞机中任意一架飞机在t 时刻之间的距离，若求出的距离都始终大于8公里，则不需要调整。否则，需要调整。

若飞机的飞行角度需要调整，我们尝试运用蒙特卡洛法随机产生一定数量组的飞机调整角度。假设2）和4）给出了飞机的碰撞条件。通过这些条件判断这些组数据中有哪些组数据可以满足飞机不发生碰撞，然后可以对得到调整角绝对值之和最小的可行解进行优化，将使飞机的调整角度和绝对值尽量小为目的进而

将可行解在其原来数据的基础上一步步减小，从而得到近似最优解。

四、建模前的准备

我们的目标是使飞机调整角度的绝对值之和尽可能小，因此我们具体采用的

蒙特卡罗方法是根据以μ=0为条件的正态分布概率密度来随机产生每一架飞机的调整角度。这样，随机产生的角度数据服从以0度为均值的正态分布。为兼顾模型求解效率，我们利用程序循环规定每一次循环通过蒙德卡洛法产生1组数，循环200次。每组数由6个数据项构成，分别作为每一架飞机的调整角度。我们先令σ=1，以样例中飞机的飞行状态作为初始条件，经过多次实验，发现每次产生的200组数据中，满足飞机不碰撞条件的调整角度方案很少，大多数情况下产生的200组数据中没有一组数据能满足使所有飞机不碰撞的条件，因此，我们尝试令σ=2，经过上百次试验，发现通过蒙特卡罗法得出的数据，平均每次有10组解满足使飞机不发生碰撞的要求。因此我们认为σ=2的条件下，用此方法得

到可行解的概率相当高，确定首先通过产生服从()20，

N 分布随机数的方法对样例数据进行求解。

为了提高算法的效率，尝试使每次求解产生的随机数减小为100组，仍以

()20，N 的规则产生随机数，经过多次实验发现，有超过一半的情况没有获得可

行解，因此我们决定在对样例数据求解时，确定运用蒙特卡罗具体实现算法中按

照()20，

N 的规则每次求解时产生200组随机数。 五、模型的建立和求解

5.1模型建立

假设()00,i i y x 表示第i 架飞机的初始位置，()i i y x ,表示第i 架飞机在规定

区域内飞行时任一时刻的位置， i α表示第i 架飞机的初始飞行角（单位°），

i β表示第i 架飞机的调整角度（单位°），t 表示飞行时刻（单位s ），v 表示飞机速度（单位km/s ）。则有：

()i i i i vt x x βα++=cos 0 63,2,1 =i (1) ()i i i i vt y y βα++=sin 0 63,2,1 =i (2)

那么第i 架飞机和第j 架飞机的距离为：()()2

2

j i

j i y y

x x -+-

由碰撞条件（即假设2））得:

()()642

2

<-+-j i j i

y y x x

j i j i ≠=且，63,2,1 (3)

将式（1）和式（2）代入式（3）得：

(()())()()()

64 )(sin - )(sin y )(cos - )(cos 2j j 0i i 02

j j 0i i 0<+++++++++βαβαβαβαvt y vt vt x vt x j i j i

化简得：

()()[

]()

[

]{}++++-+2

2

22sin cos cos j i j i j i t v βαβαβα

()()[]()()()[]

(){}vt y y x x j i j i j i j i j i i i 0000sin sin cos cos 2-+-++-+-+βαβαβαβα()()642

002

00<-+-+j i j i y y x x （4）

式（4）可以表示为：

64At 2<++C Bt (5)

式（5）在飞机处于规定区域内时成立。由此题目可以转化为在定区间[0,tmin]上对二次函数求极值的问题。其中min t 为第i 架飞机和第j 架飞机飞出规定区域时的时刻的较小值。 假设()C Bt At t g ++=2

则飞机在规定区域内的飞行约束条件即为：()64m in 0。因此，

当-B/(2A)<0时, ()()640m in <=g t g

当0<-B/(2A)< min t 时，()()()642/m in <-=A B g t g 当-B/(2A)> min t 时, ()()64m in m in <=t g t g

为了避免复杂运算，提高效率，我们在用以上模型的约束条件进行判断之前，先进行了一步碰撞的必要条件判断。考虑到两架飞机的飞行轨迹为两条射线，若两条射线所在直线交点在射线的反向延长线上时，因为初始位置两架飞机距离一定大于8公里（否则此题不可解），所以这两架飞机一定不会相撞。如图1所示：

图1 两架飞机的飞行轨迹相交于延长线

考虑，第i 架飞机和第j 架飞机：

第i 架飞机和第j 架飞机的轨迹曲线为：

解得两直线交点横坐标为

)

an(-)an()

an(-)an(-x j j i i j j 0i i 000c βαβαβαβα+++++=

t t t x t x y y j i j j

对上式分析得，若出现下列四种情况，第i 架飞机和第j 架飞机一定不会碰撞：

（1）当[][]????∈+360270900i i ，， βα且[][]????∈+360270900j j ，， βα时，

若

0j c i c x x x x <<且则一定不会碰撞；

（2）当[][]????∈+360270900i i ，， βα且[]??∈+27009j j ，βα时，

若

0j c i c x x x x ><且则一定不会碰撞；

（3）当[]??∈+27090i i ，βα且[][]????∈+360270900j j ，， βα时，

若

0j c i c x x x x <>且则一定不会碰撞；

（4）当[]??∈+27090i i ，βα且[]??∈+27009j j ，βα时，

若

0j c i c x x x x >>且则一定不会碰撞；

如果出现以上这四种情况，则不必用二次函数极值的约束条件来判断两架飞机是否会相撞了。

另，由假设5）可得：

()()22

6002

600

60<-+-y y x x

i i （6）

)(sin y )(cos x i

i 0i i 0βαβα++=++=vt y vt x i i {

综上所述

我们把数学模型建立成为一个非线性规划问题 目标函数：

()∑==

6

1

654321,,,,,f min i i

β

ββββββ

由式（1），（2），（3），（6）可得约束条件：

经过合并化简，等价的约束条件为：

[][]{}[

]

[

]

{}[]均成立

，对tmin 0t 64) y -(y )x -(x )y -(y )(sin -)(sin )x -(x )(cos - )(cos 2v )(sin -)(sin )(cos - )(cos 20020000j j i i 00j j i i 2

22

j

j i

i 2

j

j i i ∈<++++++++++++j i j i j i j i vt t

βαβαβαβαβαβαβαβα

5.2模型求解

考虑到直接搜索法花费时间代价会很大，所以我们考虑用蒙特卡洛方法求解一个满足约束条件的粗略的可行解，再对该可行解进行逐步调整，最终找到一个近似最优解。

在上百次试验中，通过蒙特卡洛法随机取得200组调整角度β矩阵，在β

矩阵的值服从()20，

N 分布的条件下，基本上能找到10组左右的可行解。这说明了产生200组数据一般来说是足够的。这一点我们已在建模前准备里详细进行

了说明。在优化模型中，为了提高算法通用性，我们还对调整角度搜索方法进行

的优化，即在随机数服从()20，

N 分布的条件下，在得到的200组调整矩阵中不能找到可行解时，就在随机数服从()30，

N 分布的条件下取200组调整矩阵。如果仍然找不到可行解，则再把σ增大1，直到σ>10时停止。考虑到3σ原理，

.

.t s []时间的最小值

架飞机飞出j 架飞机和第i 表示第tmin 中， 64 )y -(y )x -(x 比如域的最短时间，为式中出现飞机飞出区tmin 均成立，其中tmin ，0t 对22<+∈j i j i j

i and 1,2,3...6j i, 64 )y -(y )x -(x 1,2,3...5i 3600 )y -(y )x -(x 1,2,3...6i )(sin y 1,2,3...6i )(cos x 222

0602

060i i 0i i 0≠=<+=<+=++==++=j i j i i i i i i i vt y vt x βαβα{

我们就认为在调整角范围[-30°,30°]的约束下找不到可行解，即在所给定的初始状态条件下无解。

在用蒙特卡洛法找到的粗略的可行解中，取调整角度绝对值之和最小的一组作为局部调整的对象和调整的起点，对其进行调整角度优化，即依次使

654321,,,,,ββββββ向绝对值减小的方向调整。考虑到题目要求角度精度为

0.01°，调整步长设定为0.01°，改变调整角度直到满足o

i 005.0|0|<-β的

条件或者满足再调整一步后即发生相撞的条件为止。根据此规则按顺序对每一个角度进行调整。对以上优化结果四舍五入保留两位小数，即为本次求解得到的近似最优解。模型求解程序见附录。

我们把最后所得的近似最优解进行对时间t 的遍历验证来检验解的正确性。验证程序见附录。

我们以题目中所给的六架飞机的初始状态作为条件得出了每架飞机的调整角度。我们从所进行的试验中给出了5次试验的结果，如表1所示。

注：以下结果的测试条件：

测试用计算机——联想M490笔记本（CPU Intel i5 3320，默认主频2.6GHz ） 操作系统——windows 8 Matlab 版本——7.0

表1 飞机调整角度结果

为了对我们所建的模型进行推广，说明我们所建的模型是具有普适性的，我们又自行随机设计了两组符合初始状态要求（进入该区域的飞机在到达该区域边缘时，与区域内飞机的距离应在60公里以上）的数据来作为6架飞机的初始状态，并利用我们所建的模型得出了调整角度结果。

如图2所示，这是在我们自行随机设计的第一组符合初始状态要求的数据条件下，进行了三次实验所得的结果。图3表示的是在我们自行随机设计的第二组符合初始状态要求的数据条件下，进行了三次实验所得的结果

图2 自行设计参数条件下的模型求解情况（1）

图3 自行设计参数条件下的模型求解情况（2）

*注：请老师注意下，图2和图3中各有一个错别字：每个图右边部分的测量数据表格中，第10行中的“优化后方向调整角绝对值之和”改为“优化前方向调整角绝对值之和”。即把“前”改为“后”。

六、模型评价和模型改进

6.1模型评价

优点：

蒙德卡洛法求可行解的效率非常高，对找到的可行解的优化也非常有效。在参数为样例数据的条件下运行，最好的记录可以达到调整角度绝对值之和小于1.00度。

缺点：

通过实践发现对某些样例，模型算法平均用时情况与最坏用时情况的差距较大（最大曾出现两者比值约为1:5）。考虑到这个模型应用于保障飞行安全，故应考虑即使在最坏用时情况下也能保障飞机的安全。导致最坏用时情况的主要原因是通过蒙德卡洛法找出的某些可行解对模型的优化算法来说，可优化的幅度太大，故算法优化部分执行时间长。据此，在实际运用中，需要使用高性能计算机

求解或者限制优化执行次数以保证得出结果的及时性。

6.2模型改进：

1)在算法最后的模型优化步骤，我们只进行了一步，即对由蒙特卡罗法得到的一组角度绝对值之和最小的可行解，我们只进行了一次遍历：将这组数据中每个飞机的飞行角度依次以0.01的角度改变，遍历结束后使得最终得出的飞机调整角度和最小。

考虑到飞机上计算机性能我们用的笔记本电脑性能优越的多，在算法运算时间允许的条件下，其实还可以再次将上步得出的最后结果再次进行遍历，可以使飞机的调整角度和更小。

2)对于在某些飞机的初始状态条件下可能会出现的无解的情况，若能找出在约束条件下，造成方向角调整方案无解的必要条件，用于对初始状态的检查，则可避免盲目尝试求解造成的时间浪费。或者尽可能快得知在这样的初始状态下，某些飞机即使调整飞行角度也会发生碰撞，我们可以将模型转换为制定对应的应急决策，如求得尽可能减少损失的方向角调整方案。

七、参考文献

[1]毕开波，王晓东，刘智平，飞行器制导与控制及其MATLAB仿真技术技术：国防工业出版社，2009。

[2]郑昌文，严平，丁明跃，孙富春，飞行器航迹规划：国防工业出版者，2008。

八、附件清单

附件 1 具有初始状态碰撞判断和计算近似最优调整角的的主函数（文件名main.m）

附件 2 检验是否调整角是否合理的MATLAB子函数（文件名Confirmation.m）

附件 3 判断初始条件下是否存在飞机会发生碰撞的MATLAB子函数（文件名CrashDetection.m）

附件 4 自行设定初始状态时，检查设定是否满足题目要求的MATLAB函数（文件名OriginDisCheck.

附件1

%简介%————————————————————————————————————————

% 功能：95国赛A题求解主程序

% 作者：庞平&张航&宁培阳

% 日期：2014-4-5 %————————————————————————————————————————

%变量命名前缀约定（大致遵守“驼峰命名法”）%————————————————————————————————————————

% 前缀含义

% c 常量

% d 浮点性变量

% i 整形变量

% m 矩阵

% s 符号变量

% v 向量%————————————————————————————————————————

%变量说明（按字母升序排列）%————————————————————————————————————————

% 变量名含义或作用

% cMiu 指定正态分布随机数均值

% dSigmaMax 对应最大调整角度的随机数标准差

% dError 浮点数比较误差限度

% dSigma 指定正态分布随机数标准差

% dStep 近似最优解的每次优化步长

% i、j、k、m for循环指针

% iAdjustNum 对每个调整角度的优化次数

% iCount 记录可行解个数

% iFlightNum 指定飞机数量（题设为6）

% iMotCalNum 指定蒙德卡洛试验次数

% mBetaGood 存储通过蒙德卡洛法收集到的可行解

% t 时间符号变量

% vAdjustDirect 存储优化最优解时各角度的调整方向（正负号）

% vAlpha 存储尝试调整前各架飞机飞行的角度

% vIdealBeta 存储程序所求的近似最优解

% vX0 存储调整前各架飞机位置横坐标

% vY0 存储调整前各架飞机位置纵坐标%————————————————————————————————————————

function main(vAlpha,vX0,vY0)

%初始化数据

if nargin==0 %若无参数输入，缺省使用题目所给样例数据

vAlpha=[243 236 220.5 159 230 52]';

vX0=[150 85 150 145 130 0]';

vY0=[140 85 155 50 150 0]';

elseif nargin~=3 || length(vAlpha)~=length(vX0) || length(vAlpha)~=length(vX0) || length(vX0)~=length(vY0)

error('参数调用错误！'); %若输入参数不齐全，报错退出elseif OriginDisCheck(vX0,vY0)

error('使用数据不满足题设初态要求！');

%若输入参数不满足题设初态要求（零时刻，区域内任意两架飞机距离至少8km，且新进入飞机与区域内飞机距离至少60km）

end

tic %定时器计时开始

if(~CrashDetection(vAlpha,zeros(6,1),vX0,vY0))

%若CrashDetection返回值为0，在输入参数条件下，无飞机发生碰撞，否则存在发生碰撞的飞机

disp('在本初始状态下不会发生碰撞，即无需调整。'); %zeros(6,1)的作用是将全部调整角置0，以判断是否需要调整

return

else

disp('在本初始状态下会发生碰撞，必须调整！');

end

sprintf('判断是否会发生碰撞用时：%fs',toc)

tic %重启定时器

disp('调整角方案求解中')

mBetaGood=[];

iCount=0;

iMotCalNum=200;

iFlightNum=length(vAlpha);

cMiu=0;

iAdjustNum=1000;

dStep=0.01;

dError=0.005;

dSigmaMax=30/3;

%通过蒙德卡罗法求近似最优解。找到可行解前，逐步增加dSigma的值以扩大随机数取值范围，对一个dSigma值，试验次数为iMotCalNum

for dSigma=2:dSigmaMax

for i=1:iMotCalNum

vBeta=normrnd(cMiu,dSigma,[iFlightNum 1]);

%产生均值为dMiu，方差为dSigma的iFlightNum行1列正态分布随机数

if(~CrashDetection(vAlpha,vBeta,vX0,vY0))

%若CrashDetection返回值为0，说明按参数进行调整，无飞机发生碰撞，否则存在发生碰撞的飞机

mBetaGood=[mBetaGood;vBeta']; %存储可行解

iCount=iCount+1; %记录可行解个数 end

end

if(iCount~=0) %找到可行解，退出循环 break;

end

if(dSigma==10 && iCount==0)

%当dSigma=10时，3*dSigma=30，P(-3*dSigma<=angle<=3*dSigma)=95.3%,可认为达到调整角极限

error('未找到可行解！');

end

end

vIdealBeta=mBetaGood(sum(abs(mBetaGood),2)==min(sum(abs(mBetaGood),2) ),:); %在可行解中选出角度绝对值之和最小的解向量

%对所得近似最优解进行优化

vIdealBeta=vIdealBeta'; %下文计算需要转置

vAdjustDirect=ones(iFlightNum,1); %定义一个iFlightNum行1列的全1矩阵，作为调整方向的向量的基础

vAdjustDirect(find(vIdealBeta>0))=-1; %联系下文第12行，vAdjustDirect向量确定角度优化方向，使调整方向角的值向0靠近。

for

j=1:iFlightNum %对iFlightNum架飞机分别进行最多iAdjustNum次的优化

for k=1:iAdjustNum

if(CrashDetection(vAlpha,vIdealBeta,vX0,vY0))

if(k~=1)

vIdealBeta(j)=vIdealBeta(j)-vAdjustDirect(j)*dStep; %修正做多一次for循环的错误

end

break;

end

if((abs(vIdealBeta(j)-0)

break;

end

vIdealBeta(j)=vIdealBeta(j)+vAdjustDirect(j)*dStep; %按步长为dStep、方向为vAdjustDirect的策略优化调整方向角

end

end

%输出运算结果

disp('角度角度方案(单位：°)：');

for m=1:(iFlightNum-1)

sprintf('%d号飞机：%.2f',m,vIdealBeta(m))

end

sprintf('新进入飞机：%.2f',vIdealBeta(m+1))

sprintf('调整角度绝对值之和：%.2f',sum(roundn(abs(vIdealBeta),-2))) %通过abs()取各角度绝对值，再通过roundn()保留2位小数，最后通过sum()求调整角之和sprintf('调整角方案求解用时：%fs',toc)

附件2

function [iCrashTimeNum,mRecord]=Confirmation(vAlpha,vBeta,vX0,vY0) if(nargin==0)

vAlpha=[

243

236

220.5

159

230

52];

vX0=[

150

85

150

145

130

0];

vY0=[

140

85

155

50

150

0];

vBeta=zeros(6,1);

end

v0=800/3600;

syms t;

sX=vX0+v0*cos(pi*(vBeta+vAlpha)/180)*t;

sY=vY0+v0*sin(pi*(vBeta+vAlpha)/180)*t;

mRecord=zeros(16000,4);

k=1;

for i=1:6

for j=i+1:6

sDis=(sX(i)-sX(j))^2+(sY(i)-sY(j))^2;

iArrivEdge=0;

iT=0;

while iArrivEdge==0

CrashFlag=0;

dXi=subs(sX(i),iT);

dYi=subs(sY(i),iT);

dXj=subs(sX(j),iT);

dYj=subs(sY(j),iT);

dDisValue=subs(sDis,iT);

if dXi>=0 && dXi<=160 && dYi>=0 && dYi<=160 && dXj>=0 && dXj<=160 && dYj>=0 && dYj<=160

iT=iT+1;

if dDisValue<64

CrashFlag=1;

end

mRecord(k,1)=i;

mRecord(k,2)=j;

mRecord(k,3)=k;

mRecord(k,4)=CrashFlag;

k=k+1;

else

iArrivEdge=1;

end

end

end

end

iCrashTimeNum=length(find(mRecord(:,4)==1));

if(iCrashTimeNum==0)

disp('不存在发生碰撞飞机。');

else

disp('存在发生碰撞飞机!');

end

附件3

%简介%————————————————————————————————————————

% 功能：95国赛A题碰撞判定函数

% 作者：庞平&张航&宁培阳

% 日期：2014-4-5 %————————————————————————————————————————

%变量命名前缀约定（大致遵守“驼峰命名法”）%————————————————————————————————————————

% 前缀含义

% c 常量

% d 浮点性变量

% i 整形变量

% m 矩阵

% s 符号变量

% v 向量%————————————————————————————————————————

%变量说明（按字母升序排列）%————————————————————————————————————————

% 变量名含义或作用

% vAlpha （输入参数）存储尝试调整前各架飞机飞行的角度

% vBeta （输入参数）考察的调整角度向量

% vX0 （输入参数）存储调整前各架飞机位置横坐标

% vY0 （输入参数）存储调整前各架飞机位置纵坐标

% iCrashFlag （输出参数）碰撞存在标志符

% cBorder 区域边界（题设160km）

% cFlightNum 指定所需调整飞机数量

% cLimitDis 两架飞机发生碰撞的极限距离的平方（题设8^2=64）% cV0 记录各飞机匀速飞行的速率（题设800km/h）

% dA 两飞机距离平方的二次函数表达式的二次项系数

% dB 两飞机距离平方的二次函数表达式的一次项系数

% dDis1 两架飞机距离可能的最小值（当两架飞机中的一架在到达边缘时距离最小）

% dDis2 两架飞机距离可能的最小值（若两架飞机的距离最小值产生在区域内）

% dMinDis 两架飞机距离的最小值

% dXc 假设两条射线交点存在，交点横坐标的值

% sDis 两飞机距离随时间变化的符号变量表达式

% sX 横坐标随时间变化的符号变量表达式

% sY 纵坐标随时间变化的符号变量表达式

% t 时间符号变量

% vAngle 所考察的两架飞机试图调整到的角度值

% vEdgeT 两架飞机到达区域的横、纵边缘的时间

% vIdealBeta 存储程序所求的近似最优解

% vRadAngle 经单位换算后的角度数值

% vVx 两架飞机飞行速度的水平分量绝对值

% vVy 两架飞机飞行速度的垂直分量绝对值

% vXb 所考察的两架飞机的起始横坐标

% vYb 所考察的两架飞机的起始纵坐标%————————————————————————————————————————

function iCrashFlag=CrashDetection(vAlpha,vBeta,vX0,vY0)

%初始化数据

cBorder=160;

cV0=800/3600;

cFlightNum=length(vAlpha);

syms t;

sX=vX0+cV0*cos(pi*(vBeta+vAlpha)/180)*t;

sY=vY0+cV0*sin(pi*(vBeta+vAlpha)/180)*t;

iCrashFlag=0;

cLimitDis=64;

%对cFlightNum架飞机两两进行碰撞判定

for i=1:cFlightNum

for j=i+1:cFlightNum

vAngle=[vAlpha(i)+vBeta(i),vAlpha(j)+vBeta(j)];

%提取本次循环要考察的参数

vXb=[vX0(i),vX0(j)];

vYb=[vY0(i),vY0(j)];

for k=1:2 %若角度超出[0,360]的范围，需进行修正 if vAngle(k)<0

vAngle(k)=vAngle(k)+360;

elseif vAngle(k)>360

vAngle(k)=vAngle(k)-360;

end

end

vRadAngle=pi*vAngle/180;

%将角度值换算为弧度制，便于表达下文的计算式

%利用简单的射线相交判定方法，若考察的两架飞机不满足碰撞的必要条件，则无需进行复杂运算即可判定两者不会碰撞，结束本次循环

dXc=(vYb(2)-vYb(1)+tan(vRadAngle(1))*vXb(1)-tan(vRadAngle(2))*vXb(2)) /(tan(vRadAngle(1))-tan(vRadAngle(2))); %假设两射线交点存在，计算交点 if (vAngle(1)>0&&vAngle(1)<90) || (vAngle(1)>270 && vAngle(1)<360)

%分情况讨论，判断两者的情况是否属于不会发生碰撞的情形

if (vAngle(2)>0&&vAngle(2)<90) || (vAngle(2)>270 && vAngle(2)<360)

if dXc

continue

end

elseif (vAngle(2)>90 && vAngle(2)<270)

if dXcvXb(2)

continue

end

end

elseif (vAngle(1)>90 && vAngle(1)<270)

if (vAngle(2)>0&&vAngle(2)<90) || (vAngle(2)>270 && vAngle(2)<360)

if dXc>vXb(1) && dXc

continue

end

elseif (vAngle(2)>90 && vAngle(2)<270)

if dXc>vXb(1) && dXc>vXb(2)

continue

end

end

end

%若仍存在发生碰撞的可能，则通过表达式为抛物线一段的轨迹方程，求解两架飞机的距离最小值判定

%以下准备判定所需数据，按计算顺序为：距离表达式、抛物线二次项系数、一次项系数、两架飞机到达题设边界时间

sDis=(sX(i)-sX(j))^2+(sY(i)-sY(j))^2;

dA=(cV0*cos(vRadAngle(1))-cV0*cos(vRadAngle(2)))^2+(cV0*sin(vRadAngle (1))-cV0*sin(vRadAngle(2)))^2;

dB=2*(vXb(1)-vXb(2))*(cV0*cos(vRadAngle(1))-cV0*cos(vRadAngle(2)))+2*

(vYb(1)-vYb(2))*(cV0*sin(vRadAngle(1))-cV0*sin(vRadAngle(2)));

vEdgeT=zeros(1,2);

vVx=zeros(1,2);

vVy=zeros(1,2);

for n=1:2

%根据方向角进行分类讨论，分别得到两架到达区域边界的时间

vVx(n)=abs(cV0*cos(vRadAngle(n)));

vVy(n)=abs(cV0*sin(vRadAngle(n)));

if vAngle(n)<=90

vEdgeT(n)=min((cBorder-vXb(n))/vVx(n),(cBorder-vYb(n))/vVy(n));

elseif vAngle(n)<=180

vEdgeT(n)=min(vXb(n)/vVx(n),(cBorder-vYb(n))/vVy(n)); elseif vAngle(n)<=270

vEdgeT(n)=min(vXb(n)/vVx(n),vYb(n)/vVy(n));

else

vEdgeT(n)=min((cBorder-vXb(n))/vVx(n),vYb(n)/vVy(n)); end

end

dMinEdgeT=min(vEdgeT(1),vEdgeT(2));

%取两架飞机到达区域边界的时间较小值

dDis1=subs(sDis,dMinEdgeT);

%求在dMinEdgeT时刻，两架飞机距离的数值

if(-dB/(2*dA)<0)

%若抛物线轴线位置不合理（即时间变量 t=-dB/(2*dA) 为负值），将dDis设置为无穷大，否则代入sDis求数值

dDis2=inf;

%这种情况在t=0时取最短距离，在题目可解的情况下均成立

else

dDis2=subs(sDis,-dB/(2*dA));

end

dMinDis=min([dDis1,dDis2]);

%取dDis1,dDis2较小者，即为两架飞机同在区域内时的距离最小值

iCrashFlag=0;

if dMinDis

%判定距离最小值的平方是否小于发生碰撞的极限距离的平方，若是，碰撞存在标志符置为“真”

iCrashFlag=1;

return

%凡发现存在会发生碰撞的飞机，立即退出子函数

end

end

end

end

附件4

%简介%————————————————————————————————————————

% 功能：95国赛A题输入参数检查子函数

% 作者：庞平&张航&宁培阳

% 日期：2014-4-5 %————————————————————————————————————————

%变量命名前缀约定（大致遵守“驼峰命名法”）%————————————————————————————————————————

% 前缀含义

% c 常量

% d 浮点性变量

% i 整形变量

% m 矩阵

% s 符号变量

% v 向量%————————————————————————————————————————

%变量说明（按字母升序排列）%————————————————————————————————————————

% 变量名含义或作用

% cCheckPass 坐标参数检查通过（增强可读性）

% cCheckUnpass 坐标参数检查不通过（增强可读性）

% cFlightNum 考虑飞机数量

% cLimitDis 两架飞机发生碰撞的极限距离的平方（题设8^2=64）

% cNewLimitDis 区域内飞机与新进入飞机的极限距离的平方（题设60^2=3600）

% iCheckFlag 坐标参数检查结果标志位

% vX0 存储调整前各架飞机位置横坐标

% vY0 存储调整前各架飞机位置纵坐标%————————————————————————————————————————

function iCheckFlag=OriginDisCheck(vX0,vY0)

%数据初始化