非线性规划

飞行管理问题

在约 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$轴正向的夹角

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

$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

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

表3

上述飞行管理数学规划模型可如下输入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°