非线性规划

飞行管理问题

在约 10，000m 高空的某边长 160km 的正方形区域内，经常有若干架飞机作水平
飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。
当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会
与区域内的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机
飞行的方向角，以避免碰撞。现假定条件如下：
1）不碰撞的标准为任意两架飞机的距离大于 8km;
2）飞机飞行方向角调整的幅度不应超过 30 度；
3）所有飞机飞行速度均为每小时 800km;
4）进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在 60km 以上；
5）最多需考虑 6 架飞机；
6）不必考虑飞机离开此区域后的状况。
请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进
行计算（方向角误差不超过 0.01 度），要求飞机飞行方向角调整的幅度尽量小。
设该区域 4 个顶点的座标为(0,0)，(160,0)，(160,160)，(0,160)。记录数据见表 1。

表1

飞机编号	横坐标 $\, x$	纵坐标 $\, y$	方向角（度）
1	150	140	243
2	85	85	223
3	150	155	220
4	145	50	159
5	130	150	230
新进入	0	0	52

注：方向角指飞行方向与 $x$ 轴正向的夹角

为了方便以后的讨论，引进如下记号：

$D$ 为飞行管理区域的边长

$\Omega$ 为飞行管理区域，取直角坐标系使其为 $[0,D]\times[0,D]$

$a$ 为飞机飞行速度， $a=800km/h$

$(x_i^0,y_i^0)$ 为第 $i$ 架飞机的初始位置

$(x_i(t),y_i(t))$ 为第 $i$ 架飞机在 $t$ 时刻的位置

$\theta_i^0$ 为第 $i$ 架飞机的的原飞行角度，即飞行方向与 $x$ 轴夹角， $0\le\theta_i^0<2\pi$

$\Delta\theta_i$ 为第 $i$ 架飞机的方向角调整， $-\frac\pi6\le\Delta\theta_i\le\frac\pi6$

$\theta_i=\theta_i^0+\Delta\theta_i$ 为第 $i$ 架飞机调整后的飞行方向角

根据相对运动观点在考察两架飞机 $i$ 和 $j$ 的飞行时，可以将飞机 $i$ 视为不动，而飞机 $j$ 以相对速度

v_{i j} = v_{j} - v_{i} = (a \cos θ_{j} = a \cos θ_{i}, a \sin θ_{j} - a \sin θ_{i}) = a (- 2 \sin \frac{θ_{j} + θ_{j}}{2} \sin \frac{θ_{j} - θ_{i}}{2}, 2 \cos \frac{θ_{j} + θ_{j}}{2} \sin \frac{θ_{j} - θ_{i}}{2}) = 2 a \sin \frac{θ_{j} - θ_{i}}{2} (\cos (\frac{θ_{j} + θ_{j}}{2} + \frac{π}{2}), \sin (\frac{θ_{j} + θ_{j}}{2} + \frac{π}{2}))

$v_{ij}=v_j-v_i=(a\cos\theta_j=a\cos\theta_i,a\sin\theta_j-a\sin\theta_i) \\=a(-2\sin\frac{\theta_j+\theta_j}2\sin\frac{\theta_j-\theta_i}2,2\cos\frac{\theta_j+\theta_j}2\sin\frac{\theta_j-\theta_i}2) \\=2a\sin\frac{\theta_j-\theta_i}2(\cos(\frac{\theta_j+\theta_j}2+\frac\pi2),\sin(\frac{\theta_j+\theta_j}2+\frac\pi2))$

不妨设 $\theta_j\ge\theta_i$ ，此时相对飞行方向角为 $\beta_{ij}=\frac\pi2+\frac{\theta_i+\theta_j}2$

由于两架飞机的初始距离为

图1

r_{i j} (0) = \sqrt{(x_{i}^{0} - x_{j}^{0})^{2} + (y_{i}^{0} - y_{j}^{0})^{2}} α_{i j}^{0} = \arcsin \frac{8}{r_{i j} (0)}

$r_{ij}(0)=\sqrt{(x_i^0-x_j^0)^2+(y_i^0-y_j^0)^2}\\ \alpha_{ij}^0=\arcsin\frac8{r_{ij}(0)}$

只要相对飞行方向角满足

α_{i j}^{0} < β_{i j} < 2 π - α_{i j}^{0}

$\alpha_{ij}^0<\beta_{ij}<2\pi-\alpha_{ij}^0$

时，两架飞机不能碰撞

记 $\beta_{ij}^0$ 为调整前第 $j$ 架飞机相对于第 $i$ 架飞机的相对速度与这两架飞机连线（从 $j$ 指向 $i$ 的矢量）的夹角（以连线矢量为基准，逆时针方向为正，顺时针方向为负），两架飞机不碰撞的条件为

| β_{i j}^{0} + \frac{1}{2} (Δ θ_{i} + Δ θ_{j}) | > α_{i j}^{0}

$|\beta_{ij}^0+\frac12(\Delta\theta_i+\Delta\theta_j)|>\alpha_{ij}^0$

其中

β_{m n}^{0} = 相 对 速 度 v_{m n} 的 辐 角 - 从 n 指 向 m 的 连 线 矢 量 的 辐 角 = \arg \frac{e^{i θ_{n}} - e^{i θ_{m}}}{(x_{m} + i y_{m}) - (x_{n} + i y_{n})}

$\beta_{mn}^0=相对速度v_{mn}的辐角-从n指向m的连线矢量的辐角\\ =\arg\frac{e^{i\theta_n}-e^{i\theta_m}}{(x_m+iy_m)-(x_n+iy_n)}$

$i$ 是虚数单位，这里利用复数的辐角，来计算角度 $\beta_{mn}^0\,(m,n=1,2,\cdots,6)$

本问题中的优化目标可以有不同的形式：如是所有飞机的最大调节量最小；所有飞机的调整量绝对值之和最小等。这里以所有飞机的调整量绝对值之和为目标函数，可以得到如下的数学规划模型：

min \sum_{i = 1}^{6} | Δ θ_{i} | s . t | β_{i j}^{0} + \frac{1}{2} (Δ θ_{i} + Δ θ_{j}) | > α_{i j}^{0} i, j = 1, 2 \dots, 6, i \neq j | Δ θ_{i} | \leq 30 ° ， i = 1, 2, \dots, 6

$\min\sum_{i=1}^6|\Delta\theta_i|\\ s.t\quad|\beta_{ij}^0+\frac12(\Delta\theta_i+\Delta\theta_j)|>\alpha_{ij}^0\,i,j=1,2\cdots,6,i\neq j\\ |\Delta\theta_i|\le30°，\qquad i=1,2,\cdots,6$

利用如下的程序，使用python：

import numpy as np

x0 = np.array([150, 85, 150, 145, 130, 0]) # 初始横坐标
y0 = np.array([140, 85, 155, 50, 150, 0]) # 初始纵坐标
theta_i = np.array([243, 236, 220.5, 159, 230, 52])
xy0 = np.concatenate([[x0], [y0]], axis=0)
d0 = np.zeros((6, 6))
for i in range(xy0.shape[1]):  # 计算矩阵各个列向量之间的距离
    for j in range(xy0.shape[1]):
        d0[i, j] = np.linalg.norm(xy0[:, i] - xy0[:, j])
d0[d0 == 0] = np.inf  # 设置点到自身的距离为无穷大
a0 = np.arcsin(8 / d0) / np.pi * 180  # arcsin()返回弧度制，换算成角度，获得不碰撞的最小角度
xy1 = np.array([])
for i, j in zip(x0, y0):  # 计算x+yi
    xy1 = np.append(xy1, np.complex(i, j))
xy2 = np.array([])
for i in theta_i:# 计算e_iθ
    xy2 = np.append(xy2, np.exp(np.complex(0, i / 180 * np.pi)))
b0 = np.zeros((6, 6))
for i in range(6):
    for j in range(6):
        if i != j: # 不计算对角线
            b0[i, j] = np.angle((xy2[j] - xy2[i]) / (xy1[j] - xy1[i]), deg=True)  # 范围角度值，范围在-180°到180°
            # 获得两点之间的相对角度
np.savetxt('./a0.csv', a0, delimiter=',', fmt='%f')  # 写入数据
np.savetxt('./b0.csv', b0, delimiter=',', fmt='%f')

求得 $\alpha_{ij}^0$ 的值如表2所示

表2

0	5.39119	32.230953	5.091816	20.963361	2.234507
5.39119	0	4.804024	6.61346	5.807866	3.815925
32.230953	4.804024	0	4.364672	22.833654	2.125539
5.091816	6.61346	4.364672	0	4.537692	2.989819
20.963361	5.807866	22.833654	4.537692	0	2.309841
2.234507	3.815925	2.125539	2.989819	2.309841	0

求得 $\beta_{ij}^0$ 如表3所示

表3

0	-70.736358	51.75	-155.82017	-6.934949	-165.525066
-70.736358	0	91.128904	137.756437	87.695154	-171
51.75	91.128904	0	-167.523689	121.213757	-179.689191
-155.82017	137.756437	-167.523689	0	-174.030766	176.474394
-6.934949	87.695154	121.213757	-174.030766	0	-178.085617
-165.525066	-171	-179.689191	176.474394	-178.085617	0

上述飞行管理数学规划模型可如下输入LINGO求解：

model:
sets: ! 变量设置;
plane/1..6/:delta;
link(plane,plane):alpha,beta;
endsets
data:
alpha=@file('C:\Users\reion\OneDrive\Documents\机器学习实战\a0.csv'); ! 导入alpha_ij的数据;
beta=@file('C:\Users\reion\OneDrive\Documents\机器学习实战\b0.csv');! 导入beta_ij的数据;
enddata
min=@sum(plane:@abs(delta)); ! 最小化角度的绝对值之和;
@for(plane:@bnd(-30,delta,30)); ! 角度的调整范围在-30°到30°之间;
@for(link(i,j)|i#ne#j:@abs(beta(i,j)+0.5*delta(i)+0.5*delta(j))>alpha(i,j)); ! 调整后的角度的绝对值要小于a_ij;
end