%% 0. 环境清理及随机种子设置
clear;
clc;
close all;
% ========== 随机性控制开关 ==========
set_seed = true; % true: 固定结果;false: 每次随机
if set_seed
seed = 42;
rng(seed);
fprintf('【随机控制】已启用固定种子,seed = %d\n', seed);
else
rng('shuffle');
fprintf('【随机控制】已禁用种子,启用完全随机模式');
end
%% 1. 定义问题规模和生成合成数据
num_candidate_centers = 50;
num_demand_points = 200;
fprintf('正在生成一个带容量约束的 %d 个候选配送中心, %d 个需求点的大规模问题...\n', ...
num_candidate_centers, num_demand_points);
% 生成带容量限制的问题数据
problem = generate_capacitated_problem_data(num_candidate_centers, num_demand_points);
fprintf('解空间大小: 2^%d ≈ %.2e\n\n', num_candidate_centers, 2^num_candidate_centers);
% 将详细的随机数据保存到 CSV 文件
save_problem_data_to_csv(problem);
% 保存每个中心 - 需求点对的完整中心 - 需求点成本矩阵
% 生成的 CSV 文件将包含以下列:CenterID(配送中心 ID)、DemandID(需求点 ID)、Distance_Center_Demand(配送中心 - 需求点距离)
% Transport_Cost_Internal_A(内部运输成本_A)和 Transport_Cost_External_B(外部运输成本_B)
save_center_demand_cost_matrix(problem);
%% 2. 算法参数定义
nPop = 50; % 蝙蝠群规模
MaxIt = 3000; % 最大迭代次数
fmin = 0; % 最小频率
fmax = 2; % 最大频率
A_init = 1.5; % 初始响度
r_init = 0.5; % 初始脉冲率
alpha = 0.95; % 响度衰减系数
gamma = 0.1; % 脉冲率增长系数
D = num_candidate_centers; % 解向量维度等于候选中心数量
%% 3. 数据生成函数
function problem_data = generate_synthetic_problem_data(num_centers, num_demands)
% 随机生成配送中心和需求点坐标
city_size = 100;
center_coords = rand(num_centers, 2) * city_size;
demand_coords = rand(num_demands, 2) * city_size;
main_hub_coord = [city_size/2, city_size/2];
% 计算配送中心到需求点的距离矩阵
dist_matrix_cd = pdist2(center_coords, demand_coords);
% 计算配送中心到物流中心的距离
dist_to_hub = pdist2(center_coords, main_hub_coord);
% 内部需求运输成本 (成本 = 固定 + 单位距离成本 * 距离 * 随机扰动)
base_cost_a = 20;
km_cost_a = 1.5;
problem_data.cost_transport_internal_a = base_cost_a + km_cost_a * dist_matrix_cd ...
.* (1 + 0.2*(rand(num_centers, num_demands) - 0.5));
% 外部需求运输成本
base_cost_b = 30;
km_cost_b = 2.0;
problem_data.cost_transport_external_b = base_cost_b + km_cost_b * dist_matrix_cd ...
.* (1 + 0.2*(rand(num_centers, num_demands) - 0.5));
% 各需求点的内部和外部需求量
problem_data.demand_internal_x = rand(1, num_demands) * 0.2 + 0.05;
problem_data.demand_external_w = rand(1, num_demands) * 0.8 + 0.2;
% 每个中心的单位仓储成本
problem_data.unit_storage_cost_C = rand(num_centers, 1) * 0.5 + 1.5;
% 固定建设成本 (随机整数)
problem_data.annual_fixed_cost_F = randi([30, 60], num_centers, 1);
% 轨道建设成本 (按距离 hub 并摊销到年)
track_cost_per_km = 30000;
num_years = 40;
problem_data.annual_track_cost_P = (dist_to_hub * track_cost_per_km * 2) / num_years;
% 保存坐标信息
problem_data.center_coords = center_coords;
problem_data.demand_coords = demand_coords;
problem_data.main_hub_coord = main_hub_coord;
% 保存维度信息
problem_data.num_dist_centers = num_centers;
problem_data.num_demand_points = num_demands;
end
CostFunction = @(x) calculate_capacitated_cost(x, problem);
function problem_data = generate_capacitated_problem_data(num_centers, num_demands)
% 调用 generate_synthetic_problem_data 生成基本数据
problem_data = generate_synthetic_problem_data(num_centers, num_demands);
% 计算所有需求点的总需求量
total_demand = sum(problem_data.demand_internal_x) + sum(problem_data.demand_external_w);
% 定义容量区间
min_capacity_per_center = total_demand / num_centers * 1.5; % 下界
max_capacity_per_center = total_demand / num_centers * 3.0; % 上界
% 随机生成各中心的最大容量
problem_data.max_capacity = rand(num_centers, 1) * ...
(max_capacity_per_center - min_capacity_per_center) + min_capacity_per_center;
% 确保总体容量充足:若总容量不足总需求 1.2 倍,则扩大容量
while sum(problem_data.max_capacity) < total_demand * 1.2
problem_data.max_capacity = problem_data.max_capacity * 1.1;
end
end
%% 4. 运行两种算法
disp('----------------------------------------------------');
disp('正在运行【优化前】的基础蝙蝠算法 (BA)...');
[BestSol_BA, BestCost_BA, CostCurve_BA] = bat_algorithm(nPop, MaxIt, fmin, fmax, A_init, r_init, alpha, gamma, D, CostFunction);
fprintf('BA 找到的最优成本: %.2f (万元)\n', BestCost_BA);
disp('----------------------------------------------------');
disp('正在运行【文献】的改进蝙蝠算法 (IBA)...');
[BestSol_IBA, BestCost_IBA, CostCurve_IBA] = iba_algorithm(nPop, MaxIt, fmin, fmax, A_init, r_init, alpha, gamma, D, CostFunction);
fprintf('IBA 找到的最优成本: %.2f (万元)\n', BestCost_IBA);
%% 5. 绘制收敛曲线
% 统一所有图的 Y 轴范围
all_costs = [CostCurve_BA(CostCurve_BA < inf); CostCurve_IBA(CostCurve_IBA < inf)];
y_min_limit = min(all_costs) * 0.98;
y_max_limit = max(all_costs) * 1.02;
figure('Name', '算法收敛曲线对比');
hold on;
plot(CostCurve_BA, 'r-', 'LineWidth', 3);
plot(CostCurve_IBA, 'g-', 'LineWidth', 3);
hold off;
legend({'BA','IBA'}, 'Location', 'best');
title('两种算法收敛曲线比较');
xlabel('迭代次数');
ylabel('总成本 (万元)');
grid on;
ax = gca;
ax.YAxis.Exponent = 0;
ytickformat('%,.0f');
ylim([y_min_limit, y_max_limit]);
drawnow;
%% 6. 最终结果输出
disp('====================== 最终结果输出 ======================');
[BestCost, best_alg_idx] = min([BestCost_BA, BestCost_IBA]);
alg_names = {'【优化前】的BA算法', '【文献】的IBA算法'};
fprintf('在两个算法中,%s 找到了最优的解。\n', alg_names{best_alg_idx});
if best_alg_idx == 1
BestSol = BestSol_BA;
else
BestSol = BestSol_IBA;
end
if BestCost == inf
disp('警告:所有算法均未能找到任何可行解!');
else
fprintf('最优总成本: %.2f (万元)\n', BestCost);
fprintf('选择的配送中心数量: %d\n', sum(BestSol));
fprintf('选择的配送中心编号: %s\n\n', mat2str(find(BestSol)));
try
generate_allocation_summary(BestSol, problem, CostFunction);
catch ME
if strcmp(ME.identifier, 'MATLAB:TooManyOutputs')
disp(' ');
disp('提示:generate_allocation_summary 需要成本函数返回额外信息。');
disp('请确保 calculate_capacitated_cost 函数的返回形式为:');
disp('[total_cost, allocation_map, goods_per_active_center] = ...');
else
rethrow(ME);
end
end
% 导出最佳解决方案中的需求中心分配表
try
[~, allocation_map] = calculate_capacitated_cost(BestSol, problem);
demand_ids = (1:problem.num_demand_points)';
mapping_table = table(demand_ids, allocation_map', ...
'VariableNames', {'需求点编号','配送中心编号'});
writetable(mapping_table, 'best_allocation_map.csv');
fprintf('已将最优方案的需求点与配送中心对应关系保存至 best_allocation_map.csv\n');
catch ME
disp('无法导出最优方案的需求–中心对应表:');
disp(ME.message);
end
end
%% 7. 辅助函数
function [total_cost, allocation_map, goods_per_active_center] = calculate_capacitated_cost(solution, problem)
% 计算二进制解向量的总成本,并返回分配信息
allocation_map = [];
goods_per_active_center = [];
active_centers_idx = find(solution == 1);
% 若没有选择任何中心,则返回无穷成本以表示不可行
if isempty(active_centers_idx)
total_cost = inf;
return;
end
num_active_centers = length(active_centers_idx);
capacity_remaining = problem.max_capacity(active_centers_idx);
% 每个需求点总需求量(内部 + 外部)
total_goods_per_demand = problem.demand_internal_x + problem.demand_external_w;
% 1. 计算每个中心向每个需求点的单位运输成本
unit_transport_cost = (problem.cost_transport_internal_a(active_centers_idx, :) ...
.* problem.demand_internal_x + ...
problem.cost_transport_external_b(active_centers_idx, :) ...
.* problem.demand_external_w) ./ total_goods_per_demand;
% 2. 创建并排序所有可能的分配 (配送中心编号, 需求编号, 单位成本)
[center_indices, demand_indices] = meshgrid(1:num_active_centers, 1:problem.num_demand_points);
allocation_list = [center_indices(:), demand_indices(:)];
costs_flat = unit_transport_cost';
allocation_list(:, 3) = costs_flat(:);
% 按成本升序排列
[~, sorted_idx] = sort(allocation_list(:, 3));
sorted_allocations = allocation_list(sorted_idx, :);
% 3. 逐一进行分配,检查容量约束
demand_assigned = false(1, problem.num_demand_points);
goods_per_active_center = zeros(num_active_centers, 1);
total_transport_cost = 0;
allocation_map = zeros(1, problem.num_demand_points);
for k = 1:size(sorted_allocations, 1)
local_center_idx = sorted_allocations(k, 1);
demand_idx = sorted_allocations(k, 2);
if ~demand_assigned(demand_idx)
demand_volume = total_goods_per_demand(demand_idx);
if capacity_remaining(local_center_idx) >= demand_volume
% 执行分配
demand_assigned(demand_idx) = true;
capacity_remaining(local_center_idx) = capacity_remaining(local_center_idx) - demand_volume;
goods_per_active_center(local_center_idx) = goods_per_active_center(local_center_idx) + demand_volume;
allocation_map(demand_idx) = active_centers_idx(local_center_idx);
% 累加该次分配的运输成本
transport_cost_this_alloc = ...
(problem.cost_transport_internal_a(active_centers_idx(local_center_idx), demand_idx) ...
* problem.demand_internal_x(demand_idx) + ...
problem.cost_transport_external_b(active_centers_idx(local_center_idx), demand_idx) ...
* problem.demand_external_w(demand_idx));
total_transport_cost = total_transport_cost + transport_cost_this_alloc;
end
end
end
% 4. 检查可行性:若仍有未分配的需求,则解不可行
if ~all(demand_assigned)
total_cost = inf;
allocation_map = [];
goods_per_active_center = [];
return;
end
% 5. 计算仓储成本与建设成本
active_storage_costs = problem.unit_storage_cost_C(active_centers_idx);
total_storage_cost = sum(goods_per_active_center .* active_storage_costs);
active_fixed_costs = problem.annual_fixed_cost_F(active_centers_idx);
active_track_costs = problem.annual_track_cost_P(active_centers_idx);
total_construction_cost = sum(active_fixed_costs) + sum(active_track_costs);
% 6. 总成本(万元)
total_cost = total_transport_cost + total_storage_cost + total_construction_cost;
end
function generate_allocation_summary(solution, problem, CostFunction)
% 根据成本函数计算分配结果,输出每个选中配送中心的服务情况
[~, allocation_map, goods_per_active_center] = CostFunction(solution);
active_centers_idx = find(solution == 1);
if isempty(active_centers_idx)
disp('没有选择任何配送中心。');
return;
end
disp('各选定配送中心的分配情况 (带容量利用率):');
disp(repmat('-', 1, 80));
fprintf('%-10s | %-15s | %-15s | %-15s | %-15s\n', '配送中心ID', '服务的需求点数', '总处理货量', '最大容量', '容量利用率');
disp(repmat('-', 1, 80));
for i = 1:length(active_centers_idx)
center_id = active_centers_idx(i);
num_served = sum(allocation_map == center_id);
total_volume = goods_per_active_center(i);
max_cap = problem.max_capacity(center_id);
utilization = total_volume / max_cap * 100;
fprintf('%-10d | %-15d | %-15.2f | %-15.2f | %-15.2f%%\n', center_id, num_served, total_volume, max_cap, utilization);
end
disp(repmat('-', 1, 80));
end
%% 8. 基础的蝙蝠算法(BA)
function [BestSol, BestCost, CostCurve] = bat_algorithm(nPop, MaxIt, fmin, fmax, A_init, r_init, alpha, gamma, D, CostFunction)
% 基础蝙蝠算法(BA)的实现
vel = zeros(nPop, D);
pos = rand(nPop, D) > 0.5;
cost = zeros(nPop, 1);
for i = 1:nPop
cost(i) = CostFunction(pos(i, :));
end
[BestCost, BestIdx] = min(cost);
BestSol = pos(BestIdx, :);
A = ones(nPop, 1) * A_init;
r = ones(nPop, 1) * r_init;
CostCurve = zeros(MaxIt, 1);
for it = 1:MaxIt
for i = 1:nPop
beta = rand;
freq = fmin + (fmax - fmin) * beta;
vel(i, :) = vel(i, :) + (pos(i, :) - BestSol) * freq;
new_pos_continuous = pos(i, :) + vel(i, :);
% 采用 Sigmoid 变换将连续值转换为二进制值
new_pos = (1 ./ (1 + exp(-new_pos_continuous))) > rand(1, D);
% 在全局最优值周围进行局部随机移动
if rand < r(i)
k = randperm(D, max(1, round(0.1*D)));
new_pos_local = BestSol;
new_pos_local(k) = 1 - new_pos_local(k);
new_pos = new_pos_local;
end
new_cost = CostFunction(new_pos);
if new_cost < cost(i) && rand < A(i)
pos(i, :) = new_pos;
cost(i) = new_cost;
A(i) = alpha * A(i);
r(i) = r_init * (1 - exp(-gamma * it));
end
if cost(i) < BestCost
BestCost = cost(i);
BestSol = pos(i, :);
end
end
CostCurve(it) = BestCost;
if BestCost ~= inf
fprintf('BA Iteration %d/%d, Best Cost = %.2f\n', it, MaxIt, BestCost);
end
end
end
%% 9. 莱维飞行函数
function s = levy(d)
% 生成莱维飞行步长
beta = 1.5;
sigma_u = (gamma(1 + beta) * sin(pi * beta / 2) ...
/ (gamma((1 + beta) / 2) * beta * 2^((beta - 1) / 2)))^(1 / beta);
sigma_v = 1;
u = random('Normal', 0, sigma_u, 1, d);
v = random('Normal', 0, sigma_v, 1, d);
s = u ./ (abs(v).^(1 / beta));
end
%% 10. 基于莱维飞行改进的蝙蝠算法(IBA)
function [BestSol, BestCost, CostCurve] = iba_algorithm(nPop, MaxIt, fmin, fmax, A_init, r_init, alpha, gamma, D, CostFunction)
% 改进蝙蝠算法(IBA),引入莱维飞行策略
pos = rand(nPop, D) > 0.5;
cost = zeros(nPop, 1);
for i = 1:nPop
cost(i) = CostFunction(pos(i, :));
end
[BestCost, BestIdx] = min(cost);
BestSol = pos(BestIdx, :);
A = ones(nPop, 1) * A_init;
r = ones(nPop, 1) * r_init;
CostCurve = zeros(MaxIt, 1);
scale_factor = 0.01;
for it = 1:MaxIt
for i = 1:nPop
levy_step = levy(D);
new_pos_continuous = pos(i,:) + scale_factor * levy_step .* (BestSol - pos(i,:));
new_pos = (1 ./ (1 + exp(-new_pos_continuous))) > rand(1, D);
if rand < r(i)
k = randperm(D, max(1, round(0.1*D)));
new_pos_local = BestSol;
new_pos_local(k) = 1 - new_pos_local(k);
new_pos = new_pos_local;
end
new_cost = CostFunction(new_pos);
if new_cost < cost(i) && rand < A(i)
pos(i, :) = new_pos;
cost(i) = new_cost;
A(i) = alpha * A(i);
r(i) = r_init * (1 - exp(-gamma * it));
end
if cost(i) < BestCost
BestCost = cost(i);
BestSol = pos(i, :);
end
end
CostCurve(it) = BestCost;
if BestCost ~= inf
fprintf('IBA Iteration %d/%d, Best Cost = %.2f\n', it, MaxIt, BestCost);
end
end
end
%% 11. 将合成数据导出为csv
function save_problem_data_to_csv(problem)
% 配送中心信息
num_centers = problem.num_dist_centers;
center_id = (1:num_centers)';
center_x = problem.center_coords(:,1);
center_y = problem.center_coords(:,2);
max_capacity = problem.max_capacity;
unit_storage_cost = problem.unit_storage_cost_C;
annual_fixed_cost = problem.annual_fixed_cost_F;
annual_track_cost = problem.annual_track_cost_P;
centers_table = table(center_id, center_x, center_y, max_capacity, unit_storage_cost, annual_fixed_cost, annual_track_cost, ...
'VariableNames', {'配送中心编号','X轴','Y轴','最大容量','单位储存费用','年度固定成本','年度运营成本'});
writetable(centers_table, 'centers_data.csv');
% 需求点信息
num_demands = problem.num_demand_points;
demand_id = (1:num_demands)';
demand_x = problem.demand_coords(:,1);
demand_y = problem.demand_coords(:,2);
internal_demand = problem.demand_internal_x';
external_demand = problem.demand_external_w';
total_demand = internal_demand + external_demand;
demand_table = table(demand_id, demand_x, demand_y, internal_demand, external_demand, total_demand, ...
'VariableNames', {'需求点编号','X轴','Y轴','内部需求量','外部需求量','总需求量'});
writetable(demand_table, 'demands_data.csv');
fprintf('已将随机生成的配送中心和需求点数据分别保存至 centers_data.csv 和 demands_data.csv\n');
end
%% 12. 导出配送中心 - 需求成本矩阵
function save_center_demand_cost_matrix(problem)
num_centers = problem.num_dist_centers;
num_demands = problem.num_demand_points;
rows = num_centers * num_demands;
data = zeros(rows, 5);
idx = 1;
for c = 1:num_centers
for d = 1:num_demands
dist_cd = pdist2(problem.center_coords(c, :), problem.demand_coords(d, :));
internal_cost = problem.cost_transport_internal_a(c, d);
external_cost = problem.cost_transport_external_b(c, d);
data(idx, :) = [c, d, dist_cd, internal_cost, external_cost];
idx = idx + 1;
end
end
% Convert to table with descriptive headers
tbl = array2table(data, 'VariableNames', ...
{'配送中心编号','需求点编号','配送中心与需求点之间的距离', ...
'内部运输成本','外部运输成本'});
writetable(tbl, 'center_demand_cost_data.csv');
fprintf('已将配送中心–需求点成本矩阵保存为 center_demand_cost_data.csv\n');
end
根据上述代码,详细解释,提供第一口吻的解释代码的稿件