%% 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
根据上述代码，详细解释，提供第一口吻的解释代码的稿件
本文介绍了一个使用批处理命令实现从文本文件中按一定概率筛选并输出行到另一个文件的方法。通过设置随机数，可以控制筛选的比例。
扫一扫
1826
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
