优化算法 | 混合K-Means蚁群算法求解CVRP问题(附Matlab ...
文案:Timelomo排版:随心390
今天为各位讲解混合K-Means蚁群算法求解CVRP,之前我们在聚类算法 | K-means聚类与DBSCAN原理及代码实现这篇推文讲解了K-Means算法的基本思想,在车辆路径问题(VRP)合集这篇推文中整理了公众号讲解CVRP和VRPTW问题的若干篇推文,需要复习的小伙伴可以选择性进行复习。
目录
1.算法设计思路
2.算法代码
3.算法实例验证
<hr/>1.算法设计步骤
算法分为两个阶段,具体流程如下:
阶段1:改进K-Means聚类
步骤 1:根据需求量总和totalDemand与车辆载重的carCap比值确定聚类数量K=\lceil\ totalDemand/carCap \rceil;
步骤2:随机选择K个需求点坐标,作为各聚类中心的初始值,并设置聚类簇容量为车辆载重;
步骤3:将需求点按照需求量由大至小的顺序排序,依次将需求点分配至相应的簇中。如需求点i的具体分配流程如下:计算需求点i与各聚类簇的距离,同时判断将需求点i分配至距离最小的簇k后,簇k的剩余容量是否满足,若满足,将需求点i分配给簇k,否则需求点i分配至**距离次优的簇k^{\prime}**,并判断是否满足簇k^{\prime}的剩余容量,重复以上步骤,直至完成需求点的i分配;
步骤4:完成所有需求点的分配后,重新计算K个簇的重心,更新各聚类中心坐标;
步骤5:判断新聚类中心坐标与原聚类中心坐标的差值是否大于设定的阈值,若大于,则转至步骤2;否则,保存聚类结果,改进K-Means聚类算法结束。
阶段2:配送路径规划
经过阶段1的聚类后,每个聚类簇内的需求点需求量总和均小于车辆载重,可以分别安排一辆车配送,即通过改进K-Means算法将CVRP转为MTSP。再使用蚁群算法(或其他经典启发式算法) 对每一个聚类簇分别优化配送路径。
<hr/>2.算法代码
整个算法共包含6个文件,在这里我们只展示其中部分代码,需要完整代码的小伙伴可以在优化算法 | 混合K-Means蚁群算法求解CVRP问题(附Matlab代码)提取代码。
kMeansCluster函数代码如下所示:
function = kMeansCluster(customer)
%% 改进kMeans聚类函数
% 输入
% customer class:struct mean:客户相关信息,详见main.m
% 输出
% clusterRes class:struct mean:聚类结果
% 思路
% 将客户根据距离+簇容量限制进行聚类
%% 数据初始化
K = ceil(customer.totalDemand/customer.carCapacity); % 需要的卡车数量,也为聚类K的值
customerVertex = customer.vertex; % 待聚类顾客坐标
tempIndex = randperm(customer.count, K); % 随机初始化聚类中心
clusterCenter = customerVertex(tempIndex, :); % 保存聚类中心坐标
%% 改进kMeans聚类算法
while 1
distanceArray = zeros(K,2); % 记录某个客户距离各聚类中心距离(第一列存储簇类别索引号,便于后续排序匹配)
clusterIndex = zeros(1, customer.count);% 标记当前客户所分配到的簇类别索引号(索引对应客户序号)
clusterLoad = zeros(1, K); % 记录各个类的容量(考虑容量限制)
clusterNeighborCount = zeros(1, K); % 记录各个簇中分别包含客户数量
newClusterCenter = zeros(K, 2); % 保存新得到的聚类中心
% 优先分配大需求的客户(客户根据需求量排序)
tempCustomerSort(:, 1) = (1:customer.count)&#39;;
tempCustomerSort(:, 2) = customer.demand;
tempCustomerSort = sortrows(tempCustomerSort, -2);
sortedCustomerIndex = tempCustomerSort(:, 1)&#39;;% 排序后的客户索引
% 遍历并根据距离和容量限制分配客户至各个簇中
for i=sortedCustomerIndex
% 计算当前客户i到各个聚类中心的距离
for j=1:K
distanceArray(j, 1) = j; % 第一列存储类索引号,为后续可能按列排序做准备
distanceArray(j, 2) = norm(clusterCenter(j, :) - customerVertex(i, :));
end
% 对distanceArray按照第二列的距离由小到大排序
distanceArray = sortrows(distanceArray, 2);
distanceArray = ceil(distanceArray);
% 分配到对应距离&#34;相对最小&#34; 且 容量限制满足 的聚类中心中
for j=1:K
tempClusterIndex = distanceArray(j, 1);
tempLoad = clusterLoad(tempClusterIndex) + customer.demand(i);
if tempLoad <= customer.carCapacity
% 容量限制满足 则加入到当前聚类中
clusterIndex(i) = tempClusterIndex; % 标记类索引
clusterLoad(tempClusterIndex) =clusterLoad(tempClusterIndex) + customer.demand(i);
clusterNeighborCount(tempClusterIndex) = clusterNeighborCount(tempClusterIndex) + 1;
break;
end
end
end
% 更新聚类中心坐标
acceptCount = 0;
for i=1:K
tempVertex = customerVertex(clusterIndex==i,:);
newClusterCenter(i,:) = sum(tempVertex)/clusterNeighborCount(i);
if norm(newClusterCenter(i,:) - clusterCenter(i,:)) < 0.1
% 是否接受当前新产生的聚类中心点
acceptCount = acceptCount + 1;
end
end
% 是否完成循环
if acceptCount == K
% 保存返回结果
clusterRes.K = K;
clusterRes.clusterIndex = clusterIndex;
clusterRes.clusterLoad = clusterLoad;
clusterRes.clusterNeighborCount = clusterNeighborCount;
clusterRes.clusterCenter = newClusterCenter;
% 跳出循环
break;
else
clusterCenter = newClusterCenter;
end
end
% 保存每个簇内具体客户序号
for i=1:K
neighborCollection{i} = find(clusterRes.clusterIndex==i); % 各聚类邻域集合
end
clusterRes.neighborCollection = neighborCollection;
%% 绘图验证
% 绘制网点中心
x0=customer.allVertex(1,1);
y0=customer.allVertex(1,2);
plot(x0,y0,&#39;ks&#39;,...
&#39;LineWidth&#39;,2,...
&#39;MarkerSize&#39;,14,...
&#39;DisplayName&#39;,&#39;物流中心&#39;);
hold on
% 遍历绘制各个聚类中心及其邻域顾客
shape = [&#39;.&#39;, &#39;*&#39;, &#39;+&#39;, &#39;o&#39;, &#39;d&#39;, &#39;p&#39;, &#39;h&#39;, &#39;x&#39;, &#39;s&#39;, &#39;^&#39;];
for j=1:K
tempCustomerIndex = clusterRes.neighborCollection{j};
plot(customerVertex(tempCustomerIndex, 1), customerVertex(tempCustomerIndex, 2), shape(j), ...
&#39;color&#39;, &#39;k&#39;, &#39;MarkerSize&#39;,5, &#39;DisplayName&#39;,[&#39;cluster&#39;, num2str(j)]);
hold on
end
xlabel(&#39;横坐标&#39;,&#39;FontSize&#39;,14)
ylabel(&#39;纵坐标&#39;,&#39;FontSize&#39;,14)
title(&#39;客户聚类结果&#39;,&#39;FontSize&#39;,14)
legend()
axis equal
end
ACO_TSP函数代码如下所示:
function = ACO_TSP(clusterIndex, clusterRes, customer, ACO)
%% 蚁群算法
% 输入
% clusterIndex class:double mean:当前聚类簇序号
% clusterRes class:struct mean:客户聚类结果
% customer class:struct mean:客户相关信息,详见main.m
% ACO class:struct mean:蚁群算法参数
% 输出
% route class:ceil mean:路径规划结果
% 思路
% 根据clusterIndex取出当前簇内客户,使用相对序号重新对客户进行编码
%% 数据初始化
customerIndex = clusterRes.neighborCollection{clusterIndex}; % 当前簇内客户序号
customerCount = length(customerIndex); % 卡车服务客户的数量
customerCount2 = customerCount + 1; % 加上配送中心
customerVertex = customer.vertex(customerIndex, :); % 卡车服务客户的坐标
allVertex = ; % 配送中心 + 卡车服务客户的坐标
cityDist = squareform(pdist(allVertex)); % 配送中心 + 卡车服务客户的直线
% 蚁群算法初始化
m = ACO.m;
Alpha = ACO.Alpha;
Beta = ACO.Beta;
Rho = ACO.Rho;
Q = ACO.Q;
G = ACO.G;
Eta = 1 ./ cityDist; % 初始化启发因子(距离倒数),包含了配送中心
Eta(find(Eta==Inf)) = eps;
Tau = ones(customerCount2, customerCount2); % 初始化信息素矩阵
Tabu = zeros(m, customerCount2); % 存储并记录路径的生成
% 记录变量
bestRoute = zeros(1, customerCount2); % 记录最佳路线
bestRouteLength = Inf; % 记录最短路径长度
BestRouteLength = Inf * ones(1, G); % 记录各代最短路径长度
iter = 1; % 迭代变量
% 注意: 此时配送中心的序号为1,客户的序号顺推一位
% 迭代寻优
while iter <= G
% 初始化m只蚂蚁的位置
randPos = [];
for i=1:(ceil(m/(customerCount2)))
randPos = ;
end
Tabu(:, 1) = (randPos(1:m))&#39;;
% m只蚂蚁按概率选择下一个客户,完成各自的周游
for j=2:customerCount2
for i=1:m
visited = Tabu(i, 1:(j-1)); % 已访问的客户
unVisited = setdiff((1:customerCount2), visited, &#39;stable&#39;); % 待访问客户
p = zeros(1, length(unVisited)); % 待访问客户的选择概率分布
% 计算待选择客户的概率分布
for k=1:length(unVisited)
p(k) = (Tau(visited(end), unVisited(k))^Alpha)*(Eta(visited(end), unVisited(k))^Beta);
end
% 按照概率原则选取下一个客户
p = p / sum(p);
p = cumsum(p);
selectIndex = find(p >= rand);
pickCustomer = unVisited(selectIndex(1));
Tabu(i,j) = pickCustomer;
end
end
% 计算各蚂蚁的路径长度(作为适应度)
L = zeros(1, m);
for i=1:m
tempPath = Tabu(i, :);
L(i) = calPathLength(tempPath, cityDist);
end
% 记录并更新最优解
= min(L);
BestRouteLength(iter) = minValue;
if minValue <= bestRouteLength
bestRouteLength = minValue;
bestRoute(1:end) = Tabu(minIndex, :);
end
% 更新信息素(修改为精英蚂蚁系统)
deltaTau = zeros(customerCount2, customerCount2);
for i=1:m
for j=1:(customerCount2-1)
% L(i)表示路径长度
deltaTau(Tabu(i,j), Tabu(i,j+1)) = deltaTau(Tabu(i,j), Tabu(i,j+1)) + Q/L(i);
% 精英种群
e = 1;
if abs(find(bestRoute == Tabu(i,j)) - find(bestRoute == Tabu(i,j+1)))
deltaTau(Tabu(i,j), Tabu(i,j+1)) = deltaTau(Tabu(i,j), Tabu(i,j+1)) + e*(1/bestRouteLength);
else
deltaTau(Tabu(i,j), Tabu(i,j+1)) = deltaTau(Tabu(i,j), Tabu(i,j+1)) + 0;
end
end
deltaTau(Tabu(i,customerCount2), Tabu(i,1)) = deltaTau(Tabu(i,customerCount2), Tabu(i,1)) + Q/L(i);
% 精英种群
e = 1;
if abs(find(bestRoute == Tabu(i,customerCount2)) - find(bestRoute == Tabu(i,1)))
deltaTau(Tabu(i,customerCount2), Tabu(i,1)) = deltaTau(Tabu(i,customerCount2), Tabu(i,1)) + e*(1/bestRouteLength);
else
deltaTau(Tabu(i,customerCount2), Tabu(i,1)) = deltaTau(Tabu(i,customerCount2), Tabu(i,1)) + 0;
end
end
Tau = (1-Rho).*Tau + deltaTau;
% 禁忌表重置
Tabu = zeros(m, customerCount2);
% 更新迭代
disp([&#39;第&#39; num2str(iter) &#39;代, 最优适应度:&#39; num2str(bestRouteLength)]);
iter = iter + 1;
end
% 将卡车相对路径还原为实际序号
= indexConversion(bestRoute, customerIndex);
% 邻域搜索算法优化最优路线
initRoute = carRoute(2:end);
initPath = + 1;
initPathLen = calPathLength(initPath, customer.Dist2);
route = initRoute;
minRoute = initRoute;
minPathLen = initPathLen;
for iter2=1:1000
= ACO_VNS_neighbor(route);
newPath = + 1;
= calPathLength(newPath, customer.Dist2);
% 更新最短距离
if newPathLen <= minPathLen
route = newRoute;
minRoute = newRoute;
minPathLen = newPathLen;
end
end
carRoute = ;
% 保存结果
route = carRoute;
%% 绘图验证
% 绘制网点中心
x0=customer.allVertex(1,1);
y0=customer.allVertex(1,2);
plot(x0,y0,&#39;ks&#39;,...
&#39;LineWidth&#39;,2,...
&#39;MarkerSize&#39;,14);
hold on
shape = [&#39;*&#39;, &#39;+&#39;, &#39;o&#39;, &#39;d&#39;, &#39;p&#39;, &#39;h&#39;, &#39;x&#39;, &#39;s&#39;, &#39;^&#39;];
for j=1:clusterRes.K
tempCustomerIndex = clusterRes.neighborCollection{j};
plot(customer.vertex(tempCustomerIndex, 1), customer.vertex(tempCustomerIndex, 2), shape(j), ...
&#39;color&#39;, &#39;k&#39;, &#39;MarkerSize&#39;,7, &#39;DisplayName&#39;,[&#39;cluster&#39;, num2str(j)]);
hold on
end
% 绘制卡车行驶路线
tempCarRoute = carRoute(2:end);
tempCarRoute = ;
tempCarRoute = tempCarRoute + 1;
axis_x = customer.allVertex(tempCarRoute, 1)&#39;; % 获取坐标
axis_y = customer.allVertex(tempCarRoute, 2)&#39;;
plot(axis_x, axis_y, &#39;-&#39;, &#39;color&#39;, &#39;k&#39;, &#39;LineWidth&#39;, 1.2);
hold on
xlabel(&#39;横坐标&#39;,&#39;FontSize&#39;,14)
ylabel(&#39;纵坐标&#39;,&#39;FontSize&#39;,14)
title(&#39;配送路径结果&#39;,&#39;FontSize&#39;,14)
axis equal
end
3.算法实例验证
以下算例运行时间均在10s左右。
(1)Solomon C201
(2)Solomon R201
(3)Solomon RC201
<hr/>咱们下期再见
近期你可能错过了的好文章
新书上架 | 《MATLAB智能优化算法:从写代码到算法思想》
优化算法 | 灰狼优化算法(文末有福利)
优化算法 | 鲸鱼优化算法
遗传算法(GA)求解带时间窗的车辆路径(VRPTW)问题MATLAB代码
粒子群优化算法(PSO)求解带时间窗的车辆路径问题(VRPTW)MATLAB代码
知乎 | bilibili | CSDN:随心390
页:
[1]