优化算法 | 混合K-Means蚁群算法求解CVRP问题（附Matlab ...

redhat9i

文案：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 [clusterRes] = 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)';
    tempCustomerSort(:, 2) = customer.demand;
    tempCustomerSort = sortrows(tempCustomerSort, -2);
    sortedCustomerIndex = tempCustomerSort(:, 1)';  % 排序后的客户索引
   
    % 遍历并根据距离和容量限制分配客户至各个簇中
    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);
        
        % 分配到对应距离"相对最小" 且 容量限制满足 的聚类中心中
        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,'ks',...
    'LineWidth',2,...
    'MarkerSize',14,...
    'DisplayName','物流中心');
hold on

% 遍历绘制各个聚类中心及其邻域顾客
shape = ['.', '*', '+', 'o', 'd', 'p', 'h', 'x', 's', '^'];
for j=1:K
    tempCustomerIndex = clusterRes.neighborCollection{j};   
    plot(customerVertex(tempCustomerIndex, 1), customerVertex(tempCustomerIndex, 2), shape(j), ...
                'color', 'k', 'MarkerSize',5, 'DisplayName',['cluster', num2str(j)]);
    hold on
end

xlabel('横坐标','FontSize',14)
ylabel('纵坐标','FontSize',14)
title('客户聚类结果','FontSize',14)
legend()
axis equal
end
ACO_TSP函数代码如下所示：
function [route] = 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 = [customer.allVertex(1, :); customerVertex];                     % 配送中心 + 卡车服务客户的坐标
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 = [randPos, randperm(customerCount2)];
    end
    Tabu(:, 1) = (randPos(1:m))';
   
    % m只蚂蚁按概率选择下一个客户，完成各自的周游
    for j=2:customerCount2
        for i=1:m
            visited = Tabu(i, 1:(j-1));                     % 已访问的客户
            unVisited = setdiff((1:customerCount2), visited, 'stable'); % 待访问客户
            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
   
    % 记录并更新最优解
    [minValue, minIndex] = 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(['第' num2str(iter) '代, 最优适应度:' num2str(bestRouteLength)]);
    iter = iter + 1;
end

% 将卡车相对路径还原为实际序号
[carRoute] = indexConversion(bestRoute, customerIndex);

% 邻域搜索算法优化最优路线
initRoute = carRoute(2:end);
initPath = [0, initRoute, 0] + 1;
initPathLen = calPathLength(initPath, customer.Dist2);
route = initRoute;
minRoute = initRoute;
minPathLen = initPathLen;
for iter2=1:1000
    [newRoute] = ACO_VNS_neighbor(route);
    newPath = [0, newRoute, 0] + 1;
    [newPathLen] = calPathLength(newPath, customer.Dist2);
    % 更新最短距离
    if newPathLen <= minPathLen
        route = newRoute;
        minRoute = newRoute;
        minPathLen = newPathLen;
    end
end

carRoute = [0, minRoute];

% 保存结果
route = carRoute;


%% 绘图验证
% 绘制网点中心
x0=customer.allVertex(1,1);
y0=customer.allVertex(1,2);
plot(x0,y0,'ks',...
    'LineWidth',2,...
    'MarkerSize',14);
hold on
shape = ['*', '+', 'o', 'd', 'p', 'h', 'x', 's', '^'];
for j=1:clusterRes.K
    tempCustomerIndex = clusterRes.neighborCollection{j};   
    plot(customer.vertex(tempCustomerIndex, 1), customer.vertex(tempCustomerIndex, 2), shape(j), ...
                'color', 'k', 'MarkerSize',7, 'DisplayName',['cluster', num2str(j)]);
    hold on
end

% 绘制卡车行驶路线
tempCarRoute = carRoute(2:end);
tempCarRoute = [0 tempCarRoute 0];
tempCarRoute = tempCarRoute + 1;
axis_x = customer.allVertex(tempCarRoute, 1)';        % 获取坐标
axis_y = customer.allVertex(tempCarRoute, 2)';
plot(axis_x, axis_y, '-', 'color', 'k', 'LineWidth', 1.2);
hold on

xlabel('横坐标','FontSize',14)
ylabel('纵坐标','FontSize',14)
title('配送路径结果','FontSize',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