Matlab【路径规划】—— 无人机药品配送路线最优化_无人机配送数学建模

问题描述

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

请制订无人机的飞行路线，使尽可能多的病人尽早得到救治。

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

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

具体思路

由题目给出的条件和附件数据（见文末），根据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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

求解出的距离矩阵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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

%%寻找无人机配送路线规划算法
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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

权重系数$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值越小越符合最优解所需的权重）：

%%飞行路线总距离不包含返航路线的最优解
x_2 = data_all(:,4);  %病人数量
%计算决策指标
obj_trk = zeros(625,1);  %决策指标初始化
dist_index = zeros(25,1);
pat_index = [data_all(1,4);zeros(24,1)];
pat_indexs = zeros(25,1);  %累计救治病人数量
dist_dxs = zeros(25,1);  %累计巡航距离
best_index = zeros(100,1); 
for index = 1:100 %自定义权重范围 
    obj_weight = index;
    disp("当前权重系数为：" + num2str(obj_weight))
    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
    %根据决策指标从起点开始规划路线
    
    %%寻找无人机配送路线规划算法
    obj_index = way_back(obj_trk,25,25);  %输入参数依次为决策指标、节点个数、起点标号
    
    for i=1:24
        pat_index(i+1,1) = data_all(obj_index(i),4);
        dist_index(i,1) = dist_A(25*(obj_index(i)-1)+obj_index(i+1),7);
        dist_index(25,1) = dist_A(25*(obj_index(i)-1)+25,7);  %返航距离
    end

    for i = 1:24
        pat_indexs(i,1) = sum(pat_index(1:i+1,1));
        dist_dxs(i,1) = sum(dist_index(1:i+1,1));
    end
    pat_indexs(25,1) = sum(pat_index) + pat_index(25,1);

    %拉格朗日中值定理确定变化率
    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(24,1)/600;  %不考虑返航路线
end
[f_x,~] = find(best_index == min(best_index), 1);
best_weight = f_x;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

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

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

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

支撑材料&工程附件

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

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