[学习笔记]非线性规划实例——飞行管理问题

一、题目阐述

二、题目分析
（1）假设条件：
1、忽略飞机的体积，转角时间的问题。
2、新飞机在进入时，所有飞机完成飞行路线调整。
3、一次就调整到位，直到下一个飞机进入区域。
（2）模型选择：飞机的位置由初始方向、调整方向、飞行速度，时间确定。判断是否相撞，可以通过判断任意两架飞机距离是否大于8公里。可以由此确定该问题是非线性规划问题。
（3）构建模型：为了保证调整幅度最小，目标函数应设置为：
Min z = ∑|△Θi| （i=1、2、3、4、5、6）
为了保证飞机在题设条件下不相撞，约束条件设置为：
Min dij(t)>8 距离大于8公里
△Θi<[-Π/6，Π/6] 角度不超过30°
构建Min dij(t)模型（纯数学问题）：设某架飞机的初始坐标是(x,y)、飞机的初始角度是theta、角度变化量是delta、飞机速度已知为v、记时间为t。飞机在t时间的位置是：(x+vtcos(theta+delta),yi+vtsin(theta+delta))
设置两架飞机的角标为i、j，两架飞机的向量是：
{[xi-xj+vt[cos(theta i+delta i )-cos(theta j+delta j )],[yi-yj+vt[sin(theta i+delta i )-sin(theta j+delta j )]}
并且很容易求得其欧式距离。
令：
A = [xi-xj+vt[cos(theta i+delta i )-cos(theta j+delta j )]
B = [yi-yj+vt[sin(theta i+delta i )-sin(theta j+delta j )]
A1 = cos(theta i+delta i )-cos(theta j+delta j )
A2 = sin(theta i+delta i )-sin(theta j+delta j )
L = A2+B2(平方不好打…)
在其他条件确定、只有t为变量的情况下对两架飞机之间的距离求导：
dL/dt = AvA1+BvB1;
Ⅰ、dL/dt > 0 :在t = 0 时，取得距离最小值，但是由于初始状态之下每一架飞机间隔距离都>8公里，所以不会相撞。
Ⅱ、dL/dt = 0 :在t使得dL/dt = 0时，取得距离最小值。
为了方便计算，移项化成“t”的不等式。如下：
fenzi = ((y(j)-y(i))(si(j)-si(i)) +(x(j)-x(i))(co(j)-co(i))) ; % 分子
fenmu = v * ((co(j)-co(i))^2 + (si(j)-si(i))^2); % 分母
t(i,j) =- fenzi / fenmu;
三、matlab代码实现
（1）绘制机场图：

%% 画出六架飞机
clear;clc
figure(1) %生成一个图像
%存储有关于飞机的数据坐标、角度
data = [150 140 243
    85  85  236
    150 155 220.5
    145 50  159 
    130 150 230
    0   0   52];
plot(data(:,1),data(:,2),'.r')
axis([0 160,0,160]);% 设置坐标轴刻度范围
hold on;
for i = 1:6
    txt = ['飞机',num2str(i)];
    text(data(i,1)+2,data(i,2)+2,txt,'FontSize',8)
end

（2）设置目标函数：

function f=fun1(delta)
    f = sum(abs(delta))*180/pi;
end

（3）设置约束条件

function [c,ceq] = nonlfun1(delta) %决策变量是六个飞机的变化角度
     x = [150 85 150 145 130  0]; % 飞机初始位置的横坐标
     y = [140 85 155  50 150  0]; % 飞机初始位置的纵坐标
     theta = [243 236 220.5 159 230 52] * pi / 180; % 飞机初始的飞行方向角
     v=800;
     co = cos(theta + delta);  % 包含6个元素的向量
     si = sin(theta + delta);  % 包含6个元素的向量
     d = zeros(6);  % 初始化飞机两两之间的最短距离矩阵
     for i = 2: 6
         for j = 1: i-1
             % 套用我们推导出来的公式计算飞机i和飞机j相距最近的时间
             fenzi = ((y(j)-y(i))*(si(j)-si(i)) +(x(j)-x(i))*(co(j)-co(i))) ;  % 分子
             fenmu =  v * ((co(j)-co(i))^2 + (si(j)-si(i))^2);  % 分母
             t(i,j) =- fenzi / fenmu;
             if t(i,j) <0  
                 d(i, j) = 1000; 
             else
                 d(i, j) = sqrt((x(j)-x(i)+v*t(i,j)*(co(j)-co(i)))^2+(y(j)-y(i)+v*t(i,j)*(si(j)-si(i)))^2); 
             end 
         end
     end
     %非线性不等式约束条件
     c = ones(15,1)*8.000001-[d(2,1);d(3,1:2)';d(4,1:3)';d(5,1:4)';d(6,1:5)'];
     % 15表示15种飞机匹配方式，取值为8.000001表示大于题设条件8
     %无非线性等式约束
     ceq=[];
end

几个处理细节：
1.当t>0时表示两架飞机越来越远，设置距离为1000表示无论如何两架飞机都不可能相撞。
2.设置距离8.000001表示距离即使是最小值也略大于8千米。
3.注意sin和cos函数的参数是弧度制，要注意角度制的转换。
4.名称尽量设计的短！尤其是遇到长的公式！！！
（3）主函数设置

%% 求解非线性规划问题
x0 = [0 0 0 0 0 0];  % 初始值
lb = -pi/6*ones(6,1);
ub = pi/6*ones(6,1);
[x,fval] = fmincon(@fun1,x0,[],[],[],[],lb,ub,@nonlfun1)
x = x * 180 / pi    % 将弧度转换为度数