function [optimal_k, L1_star, L2_star, min_cost, primaries] = pipeline_optimization()
% 主函数 - 管道最优铺设解决方案
% 加载数据
[A, points, labels] = load_data();
% STEP-1: 估算k范围
[k_min, k_max, L_total] = estimate_k_range(points);
fprintf('MST总长: %.2fkm, k_min=%d, k_max=%d\n', L_total/1000, k_min, k_max);
% 初始化结果存储
costs = inf(k_max, 1);
L1_values = zeros(k_max, 1);
L2_values = zeros(k_max, 1);
all_primaries = cell(k_max, 1);
all_clusters = cell(k_max, 1);
feasible_flags = false(k_max, 1);
% 遍历k值
for k = k_min:k_max
try
% STEP-2: 分簇与选站 (使用MST约束聚类)
[clusters, primaries] = mst_constrained_clustering(points, A, k, labels);
% STEP-3: 费用计算 (考虑一级站间连接)
[L1, L2, cost] = calculate_cost_with_level1_connections(A, primaries, clusters);
% 存储结果
costs(k) = cost;
L1_values(k) = L1;
L2_values(k) = L2;
all_primaries{k} = primaries;
all_clusters{k} = clusters;
feasible_flags(k) = true;
fprintf('k=%2d: L1=%.2fkm, L2=%.2fkm, Cost=%.2f万元\n', ...
k, L1/1000, L2/1000, cost/10000);
catch ME
fprintf('k=%d 计算失败: %s\n', k, ME.message);
end
end
% 选择最优k值 (只考虑可行解)
feasible_k = find(feasible_flags);
if isempty(feasible_k)
error('没有找到可行解');
end
[min_cost_val, idx] = min(costs(feasible_k));
optimal_k = feasible_k(idx);
L1_star = L1_values(optimal_k);
L2_star = L2_values(optimal_k);
min_cost = min_cost_val;
primaries = all_primaries{optimal_k};
clusters = all_clusters{optimal_k};
% 输出最终结果
fprintf('\n=== 最优解 ===\n');
fprintf('最优k*: %d\n', optimal_k);
fprintf('I型管总长L1*: %.2fkm\n', L1_star/1000);
fprintf('II型管总长L2*: %.2fkm\n', L2_star/1000);
fprintf('最小总费用C*: %.2f万元\n', min_cost/10000);
% 显示一级站列表
fprintf('\n一级站坐标列表 (k*=%d):\n', optimal_k);
disp('序号 X坐标 Y坐标 标签');
for i = 1:size(primaries, 1)
fprintf('%2d %5.1f %5.1f %s\n', ...
i, primaries{i, 1}(1), primaries{i, 1}(2), primaries{i, 3});
end
% 可视化
visualize_solution(A, points, primaries, clusters, optimal_k, min_cost);
end
function [clusters, primaries] = mst_constrained_clustering(points, A, k, labels)
% 构建全局MST
dist_matrix = pdist2(points, points);
G = graph(dist_matrix, 'upper');
T = minspantree(G);
% 初始化簇
clusters = {1:size(points, 1)}; % 初始为一个大簇
% 迭代切割直到满足约束且达到k个簇
while numel(clusters) < k
new_clusters = {};
for i = 1:numel(clusters)
cluster_indices = clusters{i};
cluster_points = points(cluster_indices, :);
% 计算簇的MST长度
if numel(cluster_indices) > 1
cluster_dist = pdist2(cluster_points, cluster_points);
G_cluster = graph(cluster_dist, 'upper');
T_cluster = minspantree(G_cluster);
cluster_length = sum(T_cluster.Edges.Weight);
else
cluster_length = 0;
end
% 若长度超限则切割
if cluster_length > 30000
[subclusters, success] = split_cluster_by_mst(cluster_points, cluster_indices);
if success
new_clusters = [new_clusters, subclusters];
else
new_clusters = [new_clusters, {cluster_indices}];
end
else
new_clusters = [new_clusters, {cluster_indices}];
end
end
% 检查簇数量
if numel(new_clusters) > k
% 合并小簇
[new_clusters, merged] = merge_small_clusters(new_clusters, points, k);
if ~merged
error('无法调整簇数量');
end
end
clusters = new_clusters;
end
% 选择一级站 (考虑中心站距离和簇内中心性)
primaries = cell(numel(clusters), 3); % {坐标, 原始索引, 标签}
for i = 1:numel(clusters)
cluster_indices = clusters{i};
cluster_points = points(cluster_indices, :);
% 计算加权得分
dist_to_A = vecnorm(cluster_points - A, 2, 2);
centrality = zeros(size(cluster_points, 1), 1);
for j = 1:size(cluster_points, 1)
centrality(j) = mean(vecnorm(cluster_points - cluster_points(j, :), 2, 2));
end
scores = 0.7 * dist_to_A + 0.3 * centrality;
% 选择得分最小的点
[~, min_idx] = min(scores);
selected_point = cluster_points(min_idx, :);
point_idx = cluster_indices(min_idx);
primaries(i, :) = {selected_point, point_idx, labels{point_idx}};
end
end
function [subclusters, success] = split_cluster_by_mst(points, indices)
% 构建簇的MST
dist_matrix = pdist2(points, points);
G = graph(dist_matrix, 'upper');
T = minspantree(G);
% 寻找最长边
[~, max_idx] = max(T.Edges.Weight);
edge = T.Edges.EndNodes(max_idx, :);
% 移除最长边
T_temp = rmedge(T, edge(1), edge(2));
comp = conncomp(T_temp);
% 创建子簇
subclusters = {};
for i = 1:max(comp)
sub_indices = indices(comp == i);
subclusters{end+1} = sub_indices;
end
success = true;
end
function [clusters, merged] = merge_small_clusters(clusters, points, target_k)
merged = false;
while numel(clusters) > target_k
% 计算簇大小和簇间距离
cluster_sizes = cellfun(@numel, clusters);
[min_size, min_idx] = min(cluster_sizes);
% 找到最近邻簇
min_dist = inf;
merge_idx = 0;
for j = 1:numel(clusters)
if j == min_idx, continue; end
dist = min(pdist2(points(clusters{min_idx}, :), points(clusters{j}, :), [], 'all');
if dist < min_dist
min_dist = dist;
merge_idx = j;
end
end
% 尝试合并簇
merged_cluster = [clusters{min_idx}, clusters{merge_idx}];
merged_points = points(merged_cluster, :);
% 检查约束
if numel(merged_cluster) > 1
merged_dist = pdist2(merged_points, merged_points);
G_merged = graph(merged_dist, 'upper');
T_merged = minspantree(G_merged);
merged_length = sum(T_merged.Edges.Weight);
else
merged_length = 0;
end
if merged_length <= 30000
% 合并簇
clusters{min_idx} = merged_cluster;
clusters(merge_idx) = [];
merged = true;
else
% 无法合并
break;
end
end
end
function [L1, L2, total_cost] = calculate_cost_with_level1_connections(A, primaries, clusters)
% 提取一级站坐标
primary_coords = cell2mat(cellfun(@(x) x, primaries(:, 1), 'UniformOutput', false));
% I型管道: 中心站和所有一级站构成的MST
all_nodes = [A; primary_coords];
dist_I = pdist2(all_nodes, all_nodes);
G_I = graph(dist_I, 'upper');
T_I = minspantree(G_I);
L1 = sum(T_I.Edges.Weight); % 米单位
% II型管道: 每个簇的MST长度
L2 = 0;
for i = 1:size(primaries, 1)
% 获取簇内所有点 (包括一级站)
cluster_indices = clusters{i};
cluster_points = points(cluster_indices, :);
% 计算簇的MST
if size(cluster_points, 1) > 1
cluster_dist = pdist2(cluster_points, cluster_points);
G_cluster = graph(cluster_dist, 'upper');
T_cluster = minspantree(G_cluster);
cluster_length = sum(T_cluster.Edges.Weight);
else
cluster_length = 0;
end
% 检查约束
if cluster_length > 30000
error('簇%d违反长度约束: %.2fkm > 30km', i, cluster_length/1000);
end
L2 = L2 + cluster_length;
end
% 总成本 (米单位)
total_cost = 1291 * L1 + 445 * L2;
end
function visualize_solution(A, points, primaries, clusters, k, cost)
% ... (与之前相同的可视化代码)
% 添加一级站间连接的绘制
primary_coords = cell2mat(cellfun(@(x) x, primaries(:, 1), 'UniformOutput', false));
all_nodes = [A; primary_coords];
dist_I = pdist2(all_nodes, all_nodes);
G_I = graph(dist_I, 'upper');
T_I = minspantree(G_I);
% 绘制一级站间连接
for j = 1:size(T_I.Edges, 1)
edge = T_I.Edges.EndNodes(j, :);
if edge(1) == 1 || edge(2) == 1 % 跳过中心站的连接
continue;
end
idx1 = edge(1)-1; % 调整索引 (第一个节点是中心站)
idx2 = edge(2)-1;
p1 = primaries{idx1, 1};
p2 = primaries{idx2, 1};
plot([p1(1), p2(1)], [p1(2), p2(2)], ...
'b-', 'LineWidth', 2.0, 'Color', [0, 0.5, 0]);
end
% 更新图例
h6 = plot(nan, nan, 'b-', 'Color', [0, 0.5, 0], 'LineWidth', 2.0);
legend_items = [legend_items, {'一级站间连接'}];
legend([h1, h2, h3, h4, h5, h6], legend_items, 'Location', 'Best');
end
% 其他辅助函数 (load_data, estimate_k_range) 保持不变,帮我优化这个代码,同时帮我看一下,如何能够让这个K在合理的范围内最小
最新发布