粒子群算法(PSO)求解标准车辆路径问题(VRP)的MATLAB实现

1. 标准VRP问题描述

标准VRP问题可以描述为：

• 一个仓库(Depot)和多个客户点(Customers)
• 车队有相同容量的车辆
• 每个客户有确定的货物需求量
• 目标：设计车辆路线，使得总行驶距离最小，同时满足：
  • 每条路线从仓库出发并返回仓库
  • 每个客户只被访问一次
  • 车辆载重不超过容量限制

2. 粒子群算法求解VRP的MATLAB实现

function pso_vrp_solution()
%% PSO求解标准VRP问题
clear; clc; close all;

%% 问题参数设置
rng(42); % 设置随机种子，确保结果可重现

% 客户数量和数据生成
num_customers = 50;
customer_demands = randi([1, 10], 1, num_customers); % 客户需求
vehicle_capacity = 100; % 车辆容量
num_vehicles = 5; % 车辆数量

% 生成客户坐标（仓库在中心）
depot = [50, 50];
customer_positions = [depot; 
                     rand(num_customers, 2) * 100]; % 客户位置

% 计算距离矩阵
distance_matrix = calculate_distance_matrix(customer_positions);

%% PSO参数设置
pso_params.num_particles = 50;    % 粒子数量
pso_params.max_iterations = 200;  % 最大迭代次数
pso_params.w = 0.7;              % 惯性权重
pso_params.c1 = 1.5;             % 个体学习因子
pso_params.c2 = 1.5;             % 社会学习因子
pso_params.w_damp = 0.99;        % 惯性权重衰减系数

%% 初始化粒子群
particles = initialize_particles(pso_params, num_customers, num_vehicles, ...
                                customer_demands, vehicle_capacity, distance_matrix);

% 初始化全局最优
global_best.position = [];
global_best.cost = inf;
global_best.routes = [];

% 更新全局最优
for i = 1:pso_params.num_particles
    if particles(i).cost < global_best.cost
        global_best = particles(i);
    end
end

% 存储收敛曲线
best_costs = zeros(pso_params.max_iterations, 1);

%% PSO主循环
fprintf('开始PSO优化...\n');
for iter = 1:pso_params.max_iterations
    
    % 更新每个粒子
    for i = 1:pso_params.num_particles
        % 更新速度
        particles(i).velocity = pso_params.w * particles(i).velocity + ...
                               pso_params.c1 * rand() * (particles(i).best.position - particles(i).position) + ...
                               pso_params.c2 * rand() * (global_best.position - particles(i).position);
        
        % 限制速度范围
        particles(i).velocity = max(min(particles(i).velocity, 1), -1);
        
        % 更新位置
        new_position = particles(i).position + particles(i).velocity;
        
        % 处理边界（位置映射到客户索引）
        new_position = mod(new_position, num_customers) + 1;
        
        % 评估新位置
        [new_cost, new_routes] = evaluate_solution(new_position, num_vehicles, ...
                                                  customer_demands, vehicle_capacity, distance_matrix);
        
        % 更新个体最优
        if new_cost < particles(i).cost
            particles(i).position = new_position;
            particles(i).cost = new_cost;
            particles(i).routes = new_routes;
            
            if new_cost < particles(i).best.cost
                particles(i).best.position = new_position;
                particles(i).best.cost = new_cost;
                particles(i).best.routes = new_routes;
            end
        end
        
        % 更新全局最优
        if particles(i).cost < global_best.cost
            global_best = particles(i);
        end
    end
    
    % 记录当前最优解
    best_costs(iter) = global_best.cost;
    
    % 惯性权重衰减
    pso_params.w = pso_params.w * pso_params.w_damp;
    
    % 显示进度
    if mod(iter, 20) == 0
        fprintf('迭代 %d/%d, 当前最优距离: %.2f\n', iter, pso_params.max_iterations, global_best.cost);
    end
end

%% 结果显示
fprintf('\n=== 优化结果 ===\n');
fprintf('最优总行驶距离: %.2f\n', global_best.cost);
fprintf('使用的车辆数: %d\n', length(global_best.routes));

% 显示每条路径详情
for v = 1:length(global_best.routes)
    route = global_best.routes{v};
    route_demand = sum(customer_demands(route(2:end-1)));
    fprintf('车辆 %d: 仓库->%s->仓库, 载重: %d/%d, 距离: %.2f\n', ...
            v, num2str(route(2:end-1)), route_demand, vehicle_capacity, ...
            calculate_route_distance(route, distance_matrix));
end

%% 可视化结果
visualize_solution(global_best, customer_positions, depot, best_costs);

end

%% 计算距离矩阵
function dist_matrix = calculate_distance_matrix(positions)
    n = size(positions, 1);
    dist_matrix = zeros(n);
    for i = 1:n
        for j = 1:n
            dist_matrix(i,j) = norm(positions(i,:) - positions(j,:));
        end
    end
end

%% 初始化粒子群
function particles = initialize_particles(params, num_customers, num_vehicles, ...
                                         demands, capacity, dist_matrix)
    particles = [];
    
    for i = 1:params.num_particles
        % 随机生成客户序列
        customer_sequence = randperm(num_customers);
        
        % 评估初始解
        [cost, routes] = evaluate_solution(customer_sequence, num_vehicles, ...
                                          demands, capacity, dist_matrix);
        
        % 创建粒子
        particle.position = customer_sequence;
        particle.velocity = zeros(1, num_customers);
        particle.cost = cost;
        particle.routes = routes;
        
        % 个体最优
        particle.best.position = particle.position;
        particle.best.cost = particle.cost;
        particle.best.routes = particle.routes;
        
        particles = [particles, particle];
    end
end

%% 评估解的质量
function [total_cost, routes] = evaluate_solution(customer_sequence, num_vehicles, ...
                                                 demands, capacity, dist_matrix)
    num_customers = length(customer_sequence);
    routes = cell(1, num_vehicles);
    current_vehicle = 1;
    current_load = 0;
    current_route = [1]; % 从仓库开始（索引1）
    
    total_cost = 0;
    penalty = 10000; % 违反约束的惩罚项
    
    for i = 1:num_customers
        customer = customer_sequence(i) + 1; % 转换为位置索引（仓库是1）
        
        % 检查是否超过车辆容量
        if current_load + demands(customer-1) <= capacity
            % 可以加入当前路径
            current_route = [current_route, customer];
            current_load = current_load + demands(customer-1);
        else
            % 结束当前路径并开始新路径
            if current_vehicle < num_vehicles
                current_route = [current_route, 1]; % 返回仓库
                routes{current_vehicle} = current_route;
                total_cost = total_cost + calculate_route_distance(current_route, dist_matrix);
                
                current_vehicle = current_vehicle + 1;
                current_route = [1, customer];
                current_load = demands(customer-1);
            else
                % 没有可用车辆，施加惩罚
                total_cost = total_cost + penalty;
                break;
            end
        end
    end
    
    % 处理最后一条路径
    if ~isempty(current_route) && current_vehicle <= num_vehicles
        current_route = [current_route, 1]; % 返回仓库
        routes{current_vehicle} = current_route;
        total_cost = total_cost + calculate_route_distance(current_route, dist_matrix);
    end
    
    % 移除空路径
    routes = routes(1:current_vehicle);
end

%% 计算单条路径距离
function distance = calculate_route_distance(route, dist_matrix)
    distance = 0;
    for i = 1:length(route)-1
        distance = distance + dist_matrix(route(i), route(i+1));
    end
end

%% 可视化结果
function visualize_solution(best_solution, positions, depot, convergence)
    figure('Position', [100, 100, 1200, 500]);
    
    % 子图1: 路径可视化
    subplot(1,2,1);
    hold on;
    grid on;
    
    % 绘制客户点
    scatter(positions(2:end,1), positions(2:end,2), 50, 'b', 'filled');
    
    % 绘制仓库
    scatter(depot(1), depot(2), 100, 'r', 's', 'filled');
    text(depot(1), depot(2), ' 仓库', 'FontSize', 12, 'FontWeight', 'bold');
    
    % 绘制路径
    colors = lines(length(best_solution.routes));
    for v = 1:length(best_solution.routes)
        route = best_solution.routes{v};
        route_coords = positions(route, :);
        
        plot(route_coords(:,1), route_coords(:,2), 'o-', ...
             'Color', colors(v,:), 'LineWidth', 2, 'MarkerSize', 6);
        
        % 标注路径编号
        text(mean(route_coords(2:end-1,1)), mean(route_coords(2:end-1,2)), ...
             sprintf('V%d', v), 'FontSize', 10, 'BackgroundColor', 'white');
    end
    
    xlabel('X坐标');
    ylabel('Y坐标');
    title(sprintf('VRP路径规划 (总距离: %.2f)', best_solution.cost));
    legend('客户点', '仓库', 'Location', 'best');
    
    % 子图2: 收敛曲线
    subplot(1,2,2);
    plot(convergence, 'b-', 'LineWidth', 2);
    grid on;
    xlabel('迭代次数');
    ylabel('总行驶距离');
    title('PSO收敛曲线');
    
    % 添加统计信息
    annotation('textbox', [0.02, 0.02, 0.3, 0.1], 'String', ...
               sprintf('客户数量: %d\n车辆数量: %d\n最优距离: %.2f', ...
                       length(best_solution.position), length(best_solution.routes), best_solution.cost), ...
               'BackgroundColor', 'white', 'EdgeColor', 'black');
end

3. 改进的离散粒子群算法

对于VRP这类组合优化问题，标准的连续PSO需要改进：

function improved_pso_vrp()
%% 改进的离散PSO求解VRP
% 这里展示关键改进部分

%% 离散位置更新策略
function new_sequence = discrete_position_update(old_sequence, personal_best, global_best)
    % 基于交换操作的离散更新
    new_sequence = old_sequence;
    
    % 以一定概率学习个体最优
    if rand() < 0.3
        % 选择个人最优中的一段路径插入
        segment_length = randi([2, 5]);
        start_pos = randi(length(personal_best) - segment_length + 1);
        segment = personal_best(start_pos:start_pos+segment_length-1);
        
        % 插入到当前序列
        new_sequence = insert_segment(new_sequence, segment);
    end
    
    % 以一定概率学习全局最优
    if rand() < 0.3
        segment_length = randi([2, 5]);
        start_pos = randi(length(global_best) - segment_length + 1);
        segment = global_best(start_pos:start_pos+segment_length-1);
        
        new_sequence = insert_segment(new_sequence, segment);
    end
    
    % 确保解的有效性（所有客户都被访问）
    new_sequence = repair_solution(new_sequence);
end

%% 局部搜索优化
function improved_solution = local_search(solution, demands, capacity, dist_matrix)
    improved_solution = solution;
    
    % 2-opt局部搜索
    for i = 1:length(solution.routes)
        route = solution.routes{i};
        if length(route) > 4 % 需要足够长的路径
            improved_route = two_opt(route, dist_matrix);
            solution.routes{i} = improved_route;
        end
    end
    
    % 重新计算总成本
    total_cost = 0;
    for i = 1:length(solution.routes)
        total_cost = total_cost + calculate_route_distance(solution.routes{i}, dist_matrix);
    end
    improved_solution.cost = total_cost;
end

%% 2-opt算法
function new_route = two_opt(route, dist_matrix)
    new_route = route;
    improved = true;
    
    while improved
        improved = false;
        for i = 2:length(route)-2
            for j = i+1:length(route)-1
                % 尝试交换边
                if dist_matrix(route(i), route(j)) + dist_matrix(route(i+1), route(j+1)) < ...
                   dist_matrix(route(i), route(i+1)) + dist_matrix(route(j), route(j+1))
                    
                    new_route(i+1:j) = flip(new_route(i+1:j));
                    improved = true;
                end
            end
        end
    end
end

4. 算法性能分析

%% 性能对比分析
function performance_analysis()
    % 测试不同规模的VRP问题
    problem_sizes = [20, 50, 100];
    algorithms = {'标准PSO', '改进PSO', '遗传算法'};
    
    results = zeros(length(problem_sizes), length(algorithms));
    
    for i = 1:length(problem_sizes)
        n = problem_sizes(i);
        fprintf('测试 %d 个客户的VRP问题...\n', n);
        
        % 生成测试实例
        % [这里生成测试数据]
        
        % 运行不同算法并记录结果
        for j = 1:length(algorithms)
            % [运行相应算法]
            % results(i,j) = 最终距离;
        end
    end
    
    % 绘制性能对比图
    figure;
    bar(results);
    set(gca, 'XTickLabel', arrayfun(@(x) sprintf('%d客户', x), problem_sizes, 'UniformOutput', false));
    ylabel('总行驶距离');
    title('不同算法在VRP问题上的性能对比');
    legend(algorithms, 'Location', 'northeast');
    grid on;
end

参考代码 粒子群算法求解标准VRP问题 www.3dddown.com/cna/82797.html

5. 关键改进点总结

1. 离散编码策略：将连续PSO适应于组合优化问题
2. 路径构建算法：将客户序列转换为可行路径
3. 约束处理：通过罚函数处理容量约束
4. 局部搜索：结合2-opt等局部优化算法
5. 自适应参数：动态调整惯性权重和学习因子

6. 使用建议

• 小规模问题(<50客户)：标准PSO + 局部搜索
• 中等规模(50-200客户)：改进离散PSO + 多种局部搜索算子
• 大规模问题(>200客户)：分层PSO + 聚类预处理