非线性规划
飞行管理问题
在约 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_{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
只要相对飞行方向角满足
\[\alpha_{ij}^0<\beta_{ij}<2\pi-\alpha_{ij}^0 \]时，两架飞机不能碰撞
记\(\beta_{ij}^0\)为调整前第\(j\)架飞机相对于第\(i\)架飞机的相对速度与这两架飞机连线（从\(j\)指向\(i\)的矢量）的夹角（以连线矢量为基准，逆时针方向为正，顺时针方向为负），两架飞机不碰撞的条件为
\[|\beta_{ij}^0+\frac12(\Delta\theta_i+\Delta\theta_j)|>\alpha_{ij}^0 \]其中
\[\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|\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
求得的最优解为\(\Delta\theta_5=-28.05682 °\)，其余角度为0°