Matlab【路径规划】—— 无人机药品配送路线最优化

问题描述

某市引进一架专业大型无人机用于紧急状态下的药品投递。已知该市设有25处可用于在紧急状态接纳病人的医疗机构。其地理位置坐标（单位为公理）如下图所示。具体数据及可容纳病人数量见附件1。现要求通过数学建模，提供药品紧急配送策略，具体问题如下：

已知该市唯一的药品仓库兼设在地理位置x,y坐标分别为（82，55）的医疗机构内部，请制订无人机的飞行路线，使尽可能多的病人尽早得到救治。

上图为医疗机构及病人数量分布

题目假设：建模过程不考虑其他运载工具，也不考虑无人机的续航能力、巡航时间及承载容量限制。

具体思路

由题目给出的条件和附件数据（见文末），根据demo1计算出各节点之间的直线距离（假设直升机最短路线就是走直线）；并根据demo2进行一个简单的数学建模，实现无人机配送的路径规划——通过调整目标函数中的权重系数obj_weight，获得不同的路线方案，并以累计路径最短、优先救治病人最多作为优化目标，通过demo3的优化算法给出最优解。

demo1:求解所有节点间的距离

上图为添加了各医疗机构病人数量的重制地图，黄色节点为兼设药品仓库的医疗机构。

数据准备：
需要根据已有的数据，根据demo1将横纵坐标，全连接节点标号，计算出的距离值绘制成数据集dist_A

clc,clear,close all;
load data_all

%%求解所有节点间的距离
x_1 = [data_all(:,1),data_all(:,2),data_all(:,3)];  %节点标号，x坐标、y坐标
dist_A = zeros(25*25,7);  %用于存放计算出的距离矩阵

count = 1;
for i = 1:25
    for j = 1:25
        dist_A(j+25*(i-1),1) = count;  %第一列
        dist_A(j+25*(i-1),2) = data_all(i,2);  %起点x坐标
        dist_A(j+25*(i-1),3) = data_all(i,3);  %起点y坐标
        dist_A(j+25*(i-1),4) = j;  %第四列
        dist_A(j+25*(i-1),5) = data_all(j,2);  %终点x坐标
        dist_A(j+25*(i-1),6) = data_all(j,3);  %终点y坐标
        dist_A(j+25*(i-1),7) = sqrt((dist_A(j+25*(i-1),2)- ...
            dist_A(j+25*(i-1),5))^2+(dist_A(j+25*(i-1),3)-dist_A(j+25*(i-1),6))^2);
    end
    count = count+1;
end
save dist_A dist_A

求解出的距离矩阵dist_A——第1列为起点标号，第2、3列为起点x,y坐标，第4列为起点标号，第5、6列为起点x,y坐标，第7列为各节点之间的距离值（单位为公理）：

有了具体的距离值dist_A(:,7)和每个站点的病人数量x_2（题目附件数据已给），就可建立一个简单的目标函数，求出每一步路径规划依赖的指标——决策参数${\mathcal{L}}_{obj-trk}$：

${\mathcal{L}}_{obj-trk(i)}=\frac{1}{ob{j}_{-}weight}\cdot x_2(i)+max(dist{A}_{(25\ast (i-1)+1:25\ast (i-1)+25,7)})-dist{A}_{(25\ast (i-1)+j,7)}$

优先将决策参数${\mathcal{L}}_{obj-trk}$值较大的所对应的节点标号作为下一个预测配送点；

权重系数$ob{j}_{-}weight$值越大，则说明最短距离对路径规划的影响越大，反之病人数量的影响越大；

迭代计算过程中若遇到自身节点，则将${\mathcal{L}}_{obj-trk}$置为0——表示一轮计算中（如节点1对应节点1-25为第一轮，节点2对应节点1-25为第二轮…），无人机不重复配送同一站点；实际上，程序中还需加以条件判断是否在前往上一轮已经配送过的站点，并限制重复配送（见demo2），否则会出现两个站点间“反复横跳”的奇葩现象;

依据${\mathcal{L}}_{obj-trk}$每一轮预测出一个配送目标站点，下一轮中再将配送过的站点排除，以此类推，得出配送站点的规划顺序，其大致过程如下图所示：

demo2:直升机配送路线规划算法

%计算决策指标
obj_trk = zeros(625,1);  %决策指标初始化
obj_weight = 30; %设置影响决策的权重——值越大，最短距离影响越大，反之病人数量影响越大
for i = 1:25
    for j = 1:25
        if dist_A(25*(i-1)+j,7) ~= 0
            obj_trk(25*(i-1)+j,1) = x_2(j,1)/obj_weight + ...
                (max(dist_A(25*(i-1)+1:25*(i-1)+25,7)-dist_A(25*(i-1)+j,7)));
        else
            obj_trk(25*(i-1)+j,1) = 0;
        end
    end
end
save obj_trk obj_trk
%根据决策指标从起点25开始规划路线
%%寻找无人机配送路线规划算法（具体代码见函数way_back）
obj_index = way_back(obj_weight,obj_trk,25,25);  %输入参数依次为权重系数、决策指标、节点个数
save obj_index obj_index

%%寻找无人机配送路线规划算法
function obj_index = way_back(obj_trk,n,num1)
%输入参数为决策指标、节点个数、起点标号
obj_index = [num1;zeros(n-1,1)];
step = n;
for i = 2:n
    temp = find(obj_index(1:i-1)==step, 1);  %是否为配送过的站点
    if isempty(temp) || i == 2 || i == 3
        [value,obj_index(i)] = max(obj_trk(n*(step-1)+1:n*(step-1)+n,1));
        step = obj_index(i);
    else
        %以前配送过就再寻找下一个更合适的决策点
        temp = find(obj_index(1:i-1)==step);
        if size(temp) == [0,1]
            tk = 1;
        elseif size(temp) == [0,0]
            tk = 1;
        else
            tk = 0;
        end
        if  tk && isempty(temp)
            sort_obj = sort(obj_trk(n*(obj_index(i-1)-1)+1:n*(obj_index(i-1)-1)+n,1), 'descend'); %降序排列
            disp("已规划第"+num2str(i)+"个配送点")
            obj_temp = find(obj_trk(n*(obj_index(i-1)-1)+1:n*(obj_index(i-1)-1)+n,1)==sort_obj(2));
            obj_index(i) = obj_temp(1);
            step = obj_index(i);
        else
            sort_obj = sort(obj_trk(n*(obj_index(i-1)-1)+1:n*(obj_index(i-1)-1)+n,1), 'descend');
            disp("已规划第"+num2str(i)+"个配送点")
            for k = 2:n
                obj_temp = find(obj_trk(n*(obj_index(i-1)-1)+1:n*(obj_index(i-1)-1)+n,1)==sort_obj(k));
                obj_index(i) = obj_temp(1);
                step = obj_index(i);
                temp2 = find(obj_index(1:i-1)==step);
                if size(temp2) == [0,1]  
                    tk2 = 1;
                elseif size(temp2) == [0,0]
                    tk2 = 1;
                else
                    tk2 = 0;
                end
                if tk2 && isempty(temp2)  %是否为已配送过的站点
                    break
                else
                    continue
                end
            end
        end
    end
end
obj_index(n+1) = n;  %配送完毕返回起点
end

权重系数$ob{j}_{-}weight=30$时，计算出的决策参数${\mathcal{L}}_{obj-trk}$：

同样是25个数为一轮（以dist_A中的节点对应关系共计25*25=625个元素），每一轮中自身节点对应的参数值都为0，其余为按目标函数的算法计算；

权重系数$ob{j}_{-}weight=30$时，规划出的路线$ob{j}_{-}index$:

按照不同权重系数规划出的路线：




其中，$ob{j}_{-}weight=30$时还是比较符合现实生活中直升机、无人机飞行的巡航盘旋机制，其余按最短距离与优先救治病人数量折衷决策规划：

$ob{j}_{-}weight=30$时，计算出的飞行路线评价指标：

demo3:寻找最优解的优化算法

由于累计救治病人数量的初值和终值保持不变，故可通过拉格朗日中值定理确定一个导数$f^{\prime }(\xi )$，我们希望$\xi$出现得越早越好（越早说明优先救治病人数量越多），再结合优先选择总飞行距离最短的，组合成一个新的最优参数向量best_index，并返回其最小值best_weight（best_index值越小越符合最优解所需的权重）：

%%代码片展示：
obj_weight = 1:100;  定义权重范围
for...
    %拉格朗日中值定理确定变化率
    a = 1; b = length(pat_indexs);
    fa = pat_indexs(1,1);
    fb = pat_indexs(b,1);
    fc = (fb-fa)/(b-a);
    %对累计优先救治病人数量进行插值，使其连续化
     step = 0.01;
    x1 = a:step:b;
    pat_indexs_1 = interp1(a:b,pat_indexs,x1,'spline');
    ff = diff(pat_indexs_1)/step;  %求一阶导
    [~,p_y] = min(abs(ff-fc));  %找到最接近中值fc的数值横坐标
    best_index(index,1) = floor(p_y/length(x1)) + dist_dxs(25,1)/600;  %归一化后构建新指标
end
[f_x,~] = find(best_index == min(best_index), 1);
best_weight = f_x;

分以下两种情况求得最优解：

1. 飞行路线总距离包含返航路线在内
   best_weight = 23，最优解飞行路线及观测参数如下：
   
   

2. 飞行路线总距离不包含返航路线：
   best_weight = 24，最优解飞行路线如下：
   
   

支撑材料&工程附件

附件脚本可直接运行：
Matlab【路径规划】—— 无人机药品配送路线最优化

【注意】：本程序全部在Matlab 2022a中编写，编码格式为"UTF-8"，为保障流畅体验，建议使用2019a及以上版本打开！