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

文章目录

问题描述
- 具体思路
- - demo1:求解所有节点间的距离
  - demo2:直升机配送路线规划算法
  - demo3:寻找最优解的优化算法
支撑材料&工程附件

问题描述

某市引进一架专业大型无人机用于紧急状态下的药品投递。已知该市设有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:25for j = 1:25dist_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);endcount = 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（题目附件数据已给），就可建立一个简单的目标函数，求出每一步路径规划依赖的指标——决策参数Lobj−trk\mathcal{L}_{obj-trk}Lobj−trk：

Lobj−trk(i)=1obj−weight⋅x2(i)+max(distA(25∗(i−1)+1:25∗(i−1)+25,7))−distA(25∗(i−1)+j,7)\mathcal{L}_{obj-trk{(i)}}=\frac{1}{obj_{-}weight}·x_2(i)+max(distA_{(25*(i-1)+1:25*(i-1)+25,7)})-distA_{(25*(i-1)+j,7)}Lobj−trk(i)=obj−weight1⋅x2(i)+max(distA(25∗(i−1)+1:25∗(i−1)+25,7))−distA(25∗(i−1)+j,7)

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

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

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

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

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

%计算决策指标
obj_trk = zeros(625,1);  %决策指标初始化
obj_weight = 30; %设置影响决策的权重——值越大，最短距离影响越大，反之病人数量影响越大
for i = 1:25for j = 1:25if dist_A(25*(i-1)+j,7) ~= 0obj_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)));elseobj_trk(25*(i-1)+j,1) = 0;endend
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:ntemp = 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;elsetk = 0;endif  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);elsesort_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:nobj_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;elsetk2 = 0;endif tk2 && isempty(temp2)  %是否为已配送过的站点breakelsecontinueendendendend
end
obj_index(n+1) = n;  %配送完毕返回起点
end

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

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

权重系数obj−weight=30obj_{-}weight=30obj−weight=30时，规划出的路线obj−indexobj_{-}indexobj−index:

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

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

obj−weight=30obj_{-}weight=30obj−weight=30时，计算出的飞行路线评价指标：

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

由于累计救治病人数量的初值和终值保持不变，故可通过拉格朗日中值定理确定一个导数f′(ξ)f'(\xi)f′(ξ)，我们希望ξ\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;

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

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

支撑材料&工程附件

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

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