建模案例飞行管理问题VIP

数学建模竞赛中的优化问题 95A 飞行管理问题 一个飞行管理问题 1 1问题描述1995年全国大学生数学建模竞赛中的A题 一个飞行管理问题 在约10000米高空的某边长为160km的正方形区域内 经常有若干架飞机做水平飞行 区域内每架飞机的位置和速度向量均由计算机记录其数据 以便进行飞行管理 当一架欲进入该区域的飞机到达区域边缘时 记录其数据后 要立即计算并判断是否会与区域内的飞机发生碰撞 现假定条件如下 不碰撞的标准为任意两架飞机的距离大于8km 飞机飞行方向角调整幅度不应超过30 所有飞机飞行速度均为800km h 进入该区域的飞机在到达该区域边缘时 与区域内飞机的距离应在60km以上 最多需考虑6架飞机 不必考虑飞机离开此区域后的情况 请你对这个避免碰撞的飞行管理问题建立数学模型 列出计算步骤 对以下数据进行计算 方向角误差不超过0 01 要求飞机飞行方向角调整的幅度尽量小 该区域四个定点的坐标为 0 0 160 0 160 160 0 160 记录数据见表2 1 表2 1飞机位置和方向角记录数据 飞机编号横坐标纵坐标方向角 飞机编号横坐标x纵坐标y方向角 说明 方向角指飞行方向与x轴正向的夹角 试根据实际应用背景对你的模型进行评价和推广 对问题仔细阅读 首先抓住题目中的关键词 管理 进行联想 抓住诸如 碰撞 调整 避免碰撞 立即 判断 等等词语 联系解决问题的方案 不加约束继续联想 再将关键词搭配起来 1 问题的前期分析 160km 160km 飞行位置示意图 立即判断 碰撞 条件 实时 算法 避免碰撞 调整方向角 实时 幅度尽量小 相对 距离 优化问题 优化算法 优化调整方案 问题的初步理解和想法 飞行管理问题是优化问题 在调整方向角的幅度尽量小的同时 还必须注意调整方案及算法的实时性 2 问题探究 1 优化问题的目标函数为何 方向角调整的尽量小 方向角如何表示 方向角的概念是什么 方向角的平方和 目标函数 任意两架飞机的距离大于8公里 飞机飞行方向角调整的幅度不应超过30 2 优化问题的约束条件为何 两点间距离表示方法 判断避免碰撞的依据 把飞机视为点 总结 目标函数和约束条件中都含有方向角 3 分析6架飞机目前碰撞情况 描点作图 分析飞机在飞行区域的时间 判断5架飞机此刻的情况 判断新进入飞机与其他5架飞机碰撞的情况 4 求解方法 特殊到一般 先考虑2架 然考虑3架 优化问题为非线性规划问题 编程求解 总结 初等算法和高等算法都可 创新之处 算法 模型一及求解 模型建立这个问题显然是一个优化问题 设第i架飞机在调整时的方向角为 题目中已给出 调整后的方向为 题目中就是要求飞机飞行方向角调整的幅度尽量小 因此有化的目的函数可以是 1 为了建立这个问题的优化模型 只需要明确约束条件就可以了 一个简单的约束是飞机飞行方向角调整的幅度不应超过30 即 2 题目中要求进入该区域的飞机在到达该区域边缘时 与区域内的飞机的距离应在60km以上 这个条件是个初始条件 很容易验证目前所给的数据是满足的 因此本模型中可以不予考虑 剩下的关键是要满足题目中描述的任意两架位于该区域内的飞机的距离应该大于8km 但这个问题的难点在于飞机是动态的 这个约束不好直接描述 为此我们首先需要描述每架飞机的飞行轨迹 记飞机飞行速率为 800km h 以当前时刻为0时刻 设第架飞机在调整时的位置坐标为 已知条件 时刻的位置坐标为 则 3 如果要严格表示两架位于该区域内的飞机的距离应大于8km 则需要考虑每架飞机在区域内的飞行时间的长度 记Ti为第架飞机飞出区域的时间 即 4 记时刻第架飞机与第架飞机的距离为 并记 这时在区域内飞机不相撞的约束条件就变成了 5 其中 6 此外 经过计算可以得到 7 8 9 10 所以是一个关于t的二次函数 表示的是一条开口向上的抛物线 当即 记为 时 函数取最小值 注意到 初始时刻不相撞 如果 即 则此时约束条件 5 一定成立 所以如果且 只要在右端点的函数值非负即可 即 11 如果且 只需要求最小值即可 即 12 实际上 约束 11 表示的是在右端点的函数值非负 这个约束在 12 的条件下也是自然成立的 所以可以是对约束 11 不再附加且的条件 于是我们的模型就是 13 14 15 16 模型求解 上面这是一个非线性规划模型 虽然是严格满足题目要求的模型 但得到的模型逻辑关系比较复杂 约束 16 是在一定条件下才成立的约束 而且其中的计算式 4 也含有相当复杂的关系式 使用LINGO软件不太容易将模型很方便的输入 因为逻辑处理不是LINGO的优势所在 即使想办法把这个模型输入到LINGO 也不一定能求出好的解 笔者尝试过 但是LINGO运行时有时会出现系统内部错误 可能是系统有问题 无法继续求解 而且 在实时飞行调度中显然需要快速求解 所以下面我们想办法简化模型 这个模型麻烦之处就在于 要求严格表示两架飞机的飞行距离应大于8km 所以需要考虑每架飞机在区域内的飞行时间的长度 比较繁琐 注意到区域对角线的长度只有 任何一架飞机在所考虑的区域内停留的时间不会超过 因此这里我们简化一下问题 不再单独考虑每架飞机在区域内停留的时间 而是以最大时间 这是已经是一个常数 代替之 此时所有 这实际上强化了问题的要求 即考虑了有些飞机可能已经飞出区域 但仍不允许两架飞机的距离小于8km 程序 MODEL TITLE飞行管理问题的非线性规划模型 SETS Plane 1 6 x0 y0 cita0 cita1 d cita cita0表示初始角度 cita1为调整后的角度 d cita为调整的角度 link plane plane ENDSETSDATA x0y0cita0 max cita 30 T max 0 283 V 800 ENDDATAINIT d cita 000000 ENDINIT for plane cita1 cita0 d cita 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 避免碰撞的条件 右端点非负 for link i j Right 2 V T max sin cita1 i cita1 j 3 14159265 360 2 b i j 2 V T max sin cita1 i cita1 j 3 14159265 360 c i j 0 最小点非负 for link i j Minimum 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 max b i j 2 4 c i j 1 0 for link i j if b i j lt 0 b i j 2 4 c i j 1 0 for link free b 调整角度上下限 单位为角度 for plane bnd max cita d cita max cita obj MIN SUM plane d cita 2 obj MIN SUM plane abs d cita END 运行这个程序 结果得到的是一个局部极小点 调整角度较大 能找到更好的解吗 如果不用全局求解程序 通常很难得到稍大规模的非线性规划问题的全局最优解 所以我们启动LINGO全局求解程序求解这个模型 过程省略 可以得到全局最优解 可以看出 在0 01度的误差要求下 需要调整第3 4 6三架飞机的角度 分别调整2 06度 0 5度 1 57度 调整量的平方和为6 95 其实 使用全局变量求解程序 通常也不一定要等到全局最优解 而是观察求解状态窗口 看到一个较好的当前解 或当前最好解在较长时间内不发生变化 时 就可以中止程序 用当前最好的局部最优解作为最后结果 例如对于本例 LINGO求出最优解大约需要一分 而实际上5秒内LINGO得到了与全局最优解类似的解 此外 上面的模型还可以进一步简化 例如可以假设要求飞机永远不相撞 即认为为无穷大 这时显然约束条件 15 是多余的 而且约束 16 中只需要的条件就可以了 也就是说上面程序中的对应部分 约束和可以改写为更简单的形式 右端点非负 不再需要 最小点非负 需化为以下形式 实际计算结果显示 此时得到的结果与前面计算的结果几乎没有差别 备注 优化的目标函数除了外 也可以假定为或等 用LINGO求解的过程的是完全类似的 计算结果略有差异 这里就不再对这个目标函数进行计算了 甚至可以考虑让参与调整的飞机的数量尽量少 这种想法在实际中也不能说没有道理 但与题目的要求不符 而且解题难度并没有减小 意义似乎不大 在实际调度中 由于计算上面的调度方案 需要时间将调度信息告知飞机驾驶员并做出调整方向角的操作也需要时间 因此如果考虑一定的反应滞后时间应该是比较合理的 也就是说 如果反应时间是10秒 则计算中应采用飞机沿当前方向角飞行10秒以后的位置作为计算的基础 问题分析 根据题目的条件 可将飞机飞行的空域视为二维平面xoy中的一个正方形 顶点在 0 0 160 0 160 160 0 160 各架飞机的飞行方向角为飞行方向与x轴正向夹角 转角 根据两飞机不碰撞的标准为二者距离大于8km 可将每架飞机视为一个以飞机为圆心 以4为半径的圆状物体 每架飞机在空域中的状态由圆心的位置矢量和飞行速度矢量确定 这样两架飞机是否碰撞就化为两圆在运动中是否相交的问题 两圆是否相交只要讨论它们的相对运动即可 C 建模时补充假定条件 1 飞机在所定区域内作直线飞行 不偏离航向 2 飞机管理系统内不发生意外 如发动机失灵 或其他意外原因迫使飞机改变航向 4 飞机管理系统发出的指令应被飞机立即执行 即认为转向是瞬间完成的 忽略飞机转向的影响 即转弯半径和转弯时间的影响 3 飞机进入区域边缘时 立即作出计算 每架飞机按照计算后的指示立即作方向角改变 5 每架飞机在在整个过程中指点改变一次方向 6 新飞机进入区域时 已在区域内部的飞机的飞行方向已调整合适 不会碰撞 7 对每架飞机方向角的相同调整量的满意程度是一样的 建模方案二 模型的建立 1 圆状模型 采用相对速度作为研究对象 符号说明 i j 第i 第j架飞机的圆心 ij 第i 第j架飞机的碰撞角 ij ji vij 第i架飞机相对第j架飞机的相对飞行速度 lij 第i 第j架飞机的圆心距 i 第i架飞机的飞行方向与x轴正向夹角 逆时针为正 xi 第i架飞机的位置矢量 vi 第i架飞机的的速度矢量 ij 第i飞机对第j架飞机的相对速度与两架飞机圆心连线的夹角 逆时针为正 不碰撞 ij ij 2 由圆状模型导出的方程 讨论 ij的改变量与第i第j两架飞机飞行方向角改变量 i j的关系 由题目条件知 vi A 800 可用复数表示速度 设第i j飞机飞行方向改变前的速度分别为 改变后的速度分别为 改变前后相对速度分别为 两者之商的幅角就是 ij 定理 对第i 第j两架飞机 其相对速度方向 ij的改变量 ij等于两飞机飞行方向角改变量之和的一半 即 模型 目标函数 Min 其中 为各飞机方向角调整量的最大值 或为 约束条件 调整方向角时不能超过300 调整飞行方向后飞机不能碰撞 模型为 化为线性规划模型 由于 i可正可负 为使各变量均非负 引入新变量 模型化为 模型求解 ij的计算 model sets plane 1 6 x0 y0 link plane plane alpha sin2 endsets for link i j i ne j sin2 i j 64 x0 i x0 j 2 y0 i y0 j 2 for link i j i ne j sin alpha 3 14159265 180 0 2 sin2 data x0 150 85 150 145 130 0 y0 140 85 155 50 150 0 enddataend ALPHA 1 1 1 234568ALPHA 1 2 5 391190ALPHA 1 3 752 2310ALPHA 1 4 5 091816ALPHA 1 5 2000 963ALPHA 1 6 2 234507ALPHA 2 1 5 391190ALPHA 2 2 1 234568ALPHA 2 3 4 804024ALPHA 2 4 6 613460ALPHA 2 5 5 807866ALPHA 2 6 3 815925ALPHA 3 1 752 2310ALPHA 3 2 4 804024ALPHA 3 3 1 234568ALPHA 3 4 4 364672ALPHA 3 5 1102 834ALPHA 3 6 2 125539ALPHA 4 1 5 091816 ALPHA 4 2 6 613460ALPHA 4 3 4 364672ALPHA 4 4 1 234568ALPHA 4 5 4 537692ALPHA 4 6 2 989819ALPHA 5 1 2000 963ALPHA 5 2 5 807866ALPHA 5 3 1102 834ALPHA 5 4 4 537692ALPHA 5 5 1 234568ALPHA 5 6 2 309841ALPHA 6 1 2 234507ALPHA 6 2 3 815925ALPHA 6 3 2 125539ALPHA 6 4 2 989819ALPHA 6 5 2 309841ALPHA 6 6 1 234568 整理可得 ij的值 单位角度 也可以用MATLAB计算 ij的值 x 150 85 150 145 130 0 y 140 85 155 50 150 0 k length x alpha zeros k fori 1 kforj 1 kifi jalpha i j 0 elsealpha i j 180 3 14159265 asin 8 sqrt x i x j 2 y i y j 2 endendendalpha 计算 ij的值程序为 计算结果为 05 3911902372235 391190237223032 2309526723314 8040239337975 0918164485506 61346048987220 9633608931285 8078662434212 2345067369953 815924775399 32 2309526723315 0918164485504 8040239337986 61346048987204 3646718991114 364671899111022 8336542040094 5376924624022 1255388575512 989819139045 20 9633608931282 2345067369955 8078662434213 81592477539922 8336542040092 1255388575514 5376924624032 98981913904502 3098413654052 3098413654050 ij的计算 a 150 85 150 145 130 0 b 140 85 155 50 150 0 x a b i c 243 236 220 5 159 230 52 pi 180 v exp i c k length a fori 1 kforj 1 kbeita i j angle v i v j angle x j x i 180 pi endendbeita 用matlab程序编写 beita 0109 2636 128 250024 1798 186 934914 4749109 26360 88 8711 42 2436 92 30489 0000231 7500271 1289012 4763301 21380 310824 1798 42 243612 476305 9692 3 5256173 0651267 6952 58 78625 969201 914414 47499 00000 3108 3 52561 91440 运算结果 最优解的计算 用LINGO求解 程序如下 model sets plane 1 6 cita link plane plane alpha beta endsetsmin sum plane abs cita for plane i bnd 30 cita i 30 for link i j i ne j abs beta i j 0 5 cita i 0 5 cita j alpha i j data alpha 0 000000 5 391190 32 230953 5 091816 20 963361 2 234507 5 391190 0 000000 4 8040024 6 813460 5 807866 3 815925 32 230953 4 804024 0 000000 4 364672 22 833654 2 125539 5 091816 6 613460 4 363673 0 000000 4 537692 2 989819 20 963361 5 807866 22 833654 4 537692 0 000000 2 309841 2 234507 3 815925 2 125539 2 989819 2 309841 0 000000 beta 0 000000109 263642 128 25000024 179830173 06505113 474934109 2636420 000000 88 871096 42 243563 92 3048479 000000 128 250000 88 870960 00000012 476311 58 7862430 31080924 179830 42 24356312 4763110 0000005 969234 3 525606174 065051 92 304846 58 7862445 9692340 0000001 91438314 4749349 0000000 310809 3 5256061 9133830 000000 enddataend 用MATLAB计算编程如下 function f g plane x alph 0 000000 5 391190 32 230953 5 091816 20 963361 2 234507 5 391190 0 000000 4 8040024 6 813460 5 807866 3 815925 32 230953 4 804024 0 000000 4 364672 22 833654 2 125539 5 091816 6 613460 4 363673 0 000000 4 537692 2 989819 20 963361 5 807866 22 833654 4 537692 0 000000 2 309841 2 234507 3 815925 2 125539 2 989819 2 309841 0 000000 bet 0 000000109 263642 128 25000024 179830173 06505113 474934109 2636420 000000 88 871096 42 243563 92 3048479 000000 128 250000 88 870960 00000012 476311 58 7862430 31080924 179830 42 24356312 4763110 0000005 969234 3 525606174 065051 92 304846 58 7862445 9692340 0000001 91438314 4749349 0000000 310809 3 5256061 9133830 000000 f abs x 1 abs x 2 abs x 3 abs x 4 abs x 5 abs x 6 g 1 alpha 1 2 abs beta 1 2 0 5 x 1 0 5 x 2 g 2 alpha 1 3 abs beta 1 3 0 5 x 1 0 5 x 3 g 3 alpha 1 4 abs beta 1 4 0 5 x 1 0 5 x 4 g 4 alpha 1 5 abs beta 1 5 0 5 x 1 0 5 x 5 g 5 alpha 1 6 abs beta 1 6 0 5 x 1 0 5 x 6 g 6 alpha 2 3 abs beta 2 3 0 5 x 2 0 5 x 3 g 7 alpha 2 4 abs beta 2 4 0 5 x 2 0 5 x 4 g 8 alpha 2 5 abs beta 2 5 0 5 x 2 0 5 x 5 g 9 alpha 2 6 abs beta 2 6 0 5 x 2 0 5 x 6 g 10 alpha 3 4 abs beta 3 4 0 5 x 3 0 5 x 4 g 11 alpha 3 5 abs beta 3 5 0 5 x 3 0 5 x 5 g 12 alpha 3 6 abs beta 3 6 0 5 x 3 0 5 x 6 g 13 alpha 4 5 abs beta 4 5 0 5 x 4 0 5 x 5 g 14 alpha 4 6 abs beta 4 6 0 5 x 4 0 5 x 6 g 15 alpha 5 6 abs beta 5 6 0 5 x 5 0 5 x 6 执行程序 x0 0 0 0 0 0 0 v1 30 ones 1 6 v2 30 ones 1 6 opt x constr plane x0 opt v1 v2 结果 x 0 00000576637983 0 000005766379832 58794980234726 0 000012434879850 000036204730951 04151019765274 最优解 模型检验 各飞行方向按此方案调整后 系统各架飞机均满足 ij I j 2 ij 结果是正确的 模型的评价与推广 1 此模型采用圆状模型分析碰撞问题是合理的 同时采用相对速度作为判别标准 既体现了碰撞的本质 相对运动 又简化了模型的计算 2 建模中用了适当的简化 将一个复杂的非线性规划问题简化为线性规划问题 既求到合理的解 又提高了运算速度 这对解决高速飞行的飞机碰撞问题是十分重要的 此模型对题目所提供的例子计算得出的结果是令人满意的 3 由对称性知模型中的约束个数是 n是飞机数 所有约束条件数是 计算量增加不大 建模方案二 1 问题分析目的 不碰撞手段 调整飞行方向角要求 调整的幅度尽量小求解思路 1 找出不碰撞的条件 可用matlab画出各架飞机在t 0时刻的位置和飞行方向 2 求调