遗传算法求解物流路径优化(CVRP)的MATLAB实现

1. 物流路径优化与CVRP问题概述

最近在物流配送领域，容量受限的车辆路径问题（CVRP）一直是个既经典又实用的课题。简单来说，就是如何用有限数量的运输车辆（每辆车都有载重限制）来服务一组地理位置分散的客户点，同时满足以下两个核心目标：第一，使用的车辆数尽可能少；第二，所有车辆行驶的总距离最短。这在实际物流运营中直接关系到运输成本的控制。

举个例子，假设你管理着一个区域配送中心，每天需要向50家便利店补货。每辆货车的最大载重是5吨，而每家店的订单重量从100公斤到1吨不等。如何安排车辆路线才能既不多派车（减少固定成本），又让司机们不走冤枉路（降低油耗和工时）？这就是典型的CVRP问题。

这类问题的复杂性在于：

• 组合爆炸：随着客户点数量增加，可能的路径组合呈阶乘级增长
• 双重约束：既要考虑单车的载重限制，又要优化全局路径效率
• 现实扰动：实际场景还需考虑时间窗、路况等更多变量

2. 遗传算法求解CVRP的核心设计

2.1 染色体编码方案

遗传算法的第一步是如何将路径方案编码为"染色体"。在我们的MATLAB实现中，采用了一种直观的编码方式：

• 用数字1-N表示N个客户点
• 用0作为车辆路径的分隔符
• 例如染色体[0,3,1,0,2,4]表示：
  • 第一辆车：仓库 → 客户3 → 客户1 → 仓库
  • 第二辆车：仓库 → 客户2 → 客户4 → 仓库

这种表示法的优势在于：

1. 直观反映车辆调度顺序
2. 分隔符自动包含仓库往返路径
3. 便于实施遗传操作（交叉、变异）

关键细节：初始化种群时，分隔符的位置需要随机生成，但要确保每段路径的累计载重不超过车辆最大容量。这是解的有效性保证。

2.2 适应度函数设计

适应度函数是遗传算法的导航仪，我们的设计需要同时考虑：

• 总行驶距离（越小越好）
• 使用车辆数（越少越好）

经过多次测试，最终采用的适应度计算公式为：

code复制fitness = 1 / (total_distance + vehicle_count * distance_weight)

其中distance_weight是个调节参数，经过实验设置为100效果较好。这意味着：

• 每多用一辆车，相当于增加100公里的距离惩罚
• 算法会优先减少车辆使用，其次优化路径长度

这种加权方式在实践中很有效，因为：

1. 车辆固定成本（司机工资、折旧等）通常高于可变成本（燃油费）
2. 减少车队规模能显著降低管理复杂度
3. 符合多数物流企业的实际成本结构

2.3 遗传操作实现

2.3.1 顺序交叉(OX)实现

交叉操作是遗传算法的核心，我们采用顺序交叉(OX)来生成子代：

1. 从父代1随机选取一段连续基因段
2. 在父代2中删除这些基因对应的客户点
3. 将父代1的基因段插入到父代2剩余基因的对应位置

matlab复制function child = crossover(parent1, parent2)
    % 去除分隔符纯客户序列
    p1 = parent1(parent1~=0); 
    p2 = parent2(parent2~=0);
    
    % 随机选取交叉段
    startPos = randi(length(p1));
    endPos = randi([startPos, length(p1)]);
    segment = p1(startPos:endPos);
    
    % 构造子代
    remaining = p2(~ismember(p2, segment));
    child = [remaining(1:startPos-1), segment, remaining(startPos:end)];
    
    % 重新插入智能分隔符
    child = insertSplitPoints(child, maxLoad);
end

2.3.2 三重变异策略

为提高种群多样性，我们实现了三种变异方式随机触发：

1. 交换变异：随机选择两个位置交换客户
2. 逆转变异：随机选择一段路径进行反转
3. 分隔符变异：重新计算载重分布并插入分隔符

matlab复制function mutated = mutate(chromosome)
    if rand() < 0.3 % 交换变异
        pos = randperm(length(chromosome),2);
        mutated = chromosome;
        mutated(pos) = mutated(fliplr(pos));
    elseif rand() < 0.3 % 逆转变异
        start = randi(length(chromosome));
        finish = randi([start, length(chromosome)]);
        mutated = chromosome;
        mutated(start:finish) = fliplr(mutated(start:finish));
    else % 分隔符变异
        mutated = insertSplitPoints(chromosome, maxLoad);
    end
end

3. 关键实现细节与避坑指南

3.1 分隔符插入算法

insertSplitPoints函数是整个算法中最容易出错的环节，其核心逻辑是：

1. 遍历染色体中的每个客户点
2. 累计当前路径段的载重量
3. 当超过maxLoad时，插入分隔符0并重置累计量

matlab复制function chromo = insertSplitPoints(raw, maxLoad)
    chromo = [];
    currentLoad = 0;
    
    for i = 1:length(raw)
        if raw(i) == 0 % 已有分隔符
            chromo = [chromo, 0];
            currentLoad = 0;
        else
            demand = getDemand(raw(i)); % 获取该客户需求量
            if currentLoad + demand > maxLoad
                chromo = [chromo, 0, raw(i)];
                currentLoad = demand;
            else
                chromo = [chromo, raw(i)];
                currentLoad = currentLoad + demand;
            end
        end
    end
end

血泪教训：这里必须严格检查边界条件！特别是当连续多个客户的累计需求刚好等于maxLoad时，容易漏插分隔符，导致后续超载。

3.2 路径距离计算优化

计算路径距离时，有几点性能优化技巧：

1. 预计算所有客户点间的距离矩阵
2. 仓库到各点的距离单独存储
3. 使用向量化运算替代循环

matlab复制% 预计算距离矩阵
distMatrix = zeros(numCustomers+1);
for i = 0:numCustomers
    for j = 0:numCustomers
        distMatrix(i+1,j+1) = norm(coords(i+1,:)-coords(j+1,:));
    end
end

% 快速计算单条路径距离
function dist = calcRouteDist(route)
    dist = distMatrix(1,route(1)+1); % 仓库到第一个客户
    for k = 1:length(route)-1
        dist = dist + distMatrix(route(k)+1, route(k+1)+1);
    end
    dist = dist + distMatrix(route(end)+1,1); % 最后客户回仓库
end

3.3 参数调优经验

经过大量实验，推荐以下参数组合作为起点：

典型收敛曲线特征：

• 前20代：适应度快速提升
• 20-50代：进入平台期
• 50代后：缓慢优化
• 100代后：基本稳定

4. 进阶优化方向

4.1 局部搜索混合策略

在每代最优个体上应用2-opt局部优化，可以显著提升解的质量：

1. 随机选择路径中的两个边
2. 尝试交换这两边的连接方式
3. 如果改进则保留新路径

matlab复制function improved = twoOpt(route)
    improved = route;
    for i = 1:length(route)-2
        for j = i+2:length(route)-1
            % 计算交换前后的距离差
            oldDist = distMatrix(route(i)+1,route(i+1)+1) + ...
                     distMatrix(route(j)+1,route(j+1)+1);
            newDist = distMatrix(route(i)+1,route(j)+1) + ...
                     distMatrix(route(i+1)+1,route(j+1)+1);
            if newDist < oldDist
                improved(i+1:j) = fliplr(improved(i+1:j));
            end
        end
    end
end

4.2 动态参数调整

智能调整遗传参数可以平衡探索与开发：

• 当种群多样性下降时（适应度方差变小），增加变异率
• 当连续多代无改进时，临时增大交叉概率
• 对精英个体采用更保守的变异策略

4.3 可视化监控

实现算法运行过程的可视化有助于调试：

1. 实时绘制最优路径图
2. 显示适应度变化曲线
3. 输出关键指标统计

matlab复制function plotRoutes(routes)
    figure; hold on;
    % 绘制仓库
    plot(coords(1,1), coords(1,2), 'rp', 'MarkerSize',15); 
    % 绘制客户点
    plot(coords(2:end,1), coords(2:end,2), 'bo');
    % 绘制路径
    colors = lines(length(routes));
    for v = 1:length(routes)
        path = [1, routes{v}+1, 1]; % 添加仓库
        plot(coords(path,1), coords(path,2), 'Color',colors(v,:));
    end
    title(sprintf('总距离: %.1f | 车辆数: %d', totalDist, length(routes)));
end

5. 实际应用中的挑战与对策

5.1 非对称距离处理

现实路网往往是非对称的（A→B ≠ B→A），需要：

1. 修改距离矩阵为不对称结构
2. 路径计算时区分方向
3. 考虑单行道等约束

5.2 时间窗约束扩展

加入服务时间窗要求后：

1. 染色体需要编码服务时间
2. 适应度函数增加时间惩罚项
3. 变异操作需维护时间可行性

5.3 大规模问题分解

当客户点超过200个时：

1. 先进行地理聚类分区
2. 各分区独立求解
3. 最后协调跨区运输

经过多次实际项目验证，这套方法在客户点50-150个的中等规模问题上，通常能在5分钟内找到与人工经验相当甚至更优的解决方案。对于需要更高精度的场景，建议结合模拟退火或禁忌搜索进行后优化。