在众多实际场景中,诸如消防设施选址、基站布局规划以及充电桩站点部署等,都面临着如何利用最少的资源,实现对所有目标对象全面覆盖的难题。为有效解决这类问题,本文提出一种全新的组合算法模型 —— 基于聚类与引力斥力优化的选址算法。
算法简介
1. K - means 聚类初始化
-
目的:快速生成初始中心点,初步划分数据分布。
-
方法:通过 K - means 算法将二维点聚类为 k 组,每组中心作为初始候选点。结合肘部法则或轮廓系数确定最优 k,并使用 K - means++ 优化初始中心选择。
2. 贪心算法补充覆盖
-
目的:确保所有数据点被覆盖,补充初始聚类未覆盖的区域。
-
方法:设置中心点覆盖半径 r,标记未被覆盖的点。每次从未覆盖点中选择能覆盖最多未覆盖点的点作为新中心点,直至所有点被覆盖。
3. 引力搜索算法(GSA)优化
-
目的:动态调整中心点位置,剔除冗余点,实现最少中心点全覆盖。
-
方法:
-
引力与斥力:中心点受未覆盖点引力(向数据移动)和其他中心点斥力(避免密集)。
-
合力驱动:根据引力和斥力的矢量和移动中心点,验证覆盖后更新位置。
-
剪枝操作:定期尝试删除冗余中心点,确保覆盖性的同时最小化中心点数量。
4. 组合算法优势
-
高效性:K - means 快速聚类,贪心算法快速覆盖,GSA 精细调优。
-
鲁棒性:结合几何覆盖与物理启发式优化,适应复杂分布。
-
可解释性:每一步逻辑清晰,便于根据场景调整参数或替换算法。
算法步骤与公式
第一步:K - means 聚类
1. 初始聚类中心
使用算法初始化聚类中心,选择距离已有中心最远的点作为新中心:
2. 聚类分配与中心更新
计算每个点到中心的欧氏距离并分配:
重新计算中心坐标:
其中Sj是第j个聚类的数据点集合。
第二步:贪心算法覆盖未被覆盖点
1. 覆盖条件
点x被中心cj覆盖当且仅当:
2. 选择最优新中心
每次选择覆盖最多未覆盖点的点:
其中I(·)是指示函数。
第三步:引力搜索算法(GSA)优化中心点
1. 质量定义
中心点ci的质量mi由覆盖点数决定:
2. 引力计算
数据点x对中心ci的引力:
其中是时变引力常数。
3. 斥力计算
其他中心cj对ci的斥力:
4. 合力与位置更新
总合力:
位置更新:
其中 λ是方向系数,step_size是步长。
5. 中心点剔除条件
若剔除中心ci后所有点仍被覆盖:
代码
python代码
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans# 第一步:K - means 聚类
# 随机生成一些示例数据
points = np.random.rand(100, 2)
# 设置聚类中心数
k = 10
# 执行 K - means 聚类
kmeans = KMeans(n_clusters=k, random_state=42)
idx = kmeans.fit_predict(points)
centers = kmeans.cluster_centers_# 第二步:贪心算法覆盖未被覆盖点
# 设置覆盖半径
r = 0.2
# 标记未被覆盖的点
covered = np.zeros(points.shape[0], dtype=bool)
for center in centers:distances = np.linalg.norm(points - center, axis=1)covered[distances <= r] = True
uncovered_points = points[~covered]while len(uncovered_points) > 0:max_coverage = 0best_point = Nonefor point in uncovered_points:distances = np.linalg.norm(points - point, axis=1)current_coverage = np.sum((distances <= r) & ~covered)if current_coverage > max_coverage:max_coverage = current_coveragebest_point = pointcenters = np.vstack((centers, best_point))distances = np.linalg.norm(points - best_point, axis=1)covered[distances <= r] = Trueuncovered_points = points[~covered]# 第三步:引力搜索算法(GSA)优化中心点
# 设置参数
T = 100 # 迭代次数
G0 = 1 # 引力常数
alpha = 2 # 质量衰减系数# 初始化质量
masses = np.zeros(centers.shape[0])
for i in range(centers.shape[0]):distances = np.linalg.norm(points - centers[i], axis=1)masses[i] = np.sum(distances <= r)for t in range(T):G = G0 * np.exp(-alpha * t / T)for i in range(centers.shape[0]):force = np.zeros(2)# 计算斥力for j in range(centers.shape[0]):if i != j:distance = np.linalg.norm(centers[i] - centers[j])if distance > 0:force += G * masses[j] * (centers[i] - centers[j]) / distance# 计算引力for j in range(points.shape[0]):distance = np.linalg.norm(centers[i] - points[j])if 0 < distance <= 2 * r:force -= G * (centers[i] - points[j]) / distance# 移动中心点(先假设移动)new_center = centers[i] + force# 检查移动后的中心点是否能覆盖所有数据点temp_centers = centers.copy()temp_centers[i] = new_centertemp_covered = np.zeros(points.shape[0], dtype=bool)for center in temp_centers:distances = np.linalg.norm(points - center, axis=1)temp_covered[distances <= r] = Trueif np.all(temp_covered):centers[i] = new_center# 尝试剔除冗余的中心点for i in range(centers.shape[0]):temp_centers = np.delete(centers, i, axis=0)temp_covered = np.zeros(points.shape[0], dtype=bool)for center in temp_centers:distances = np.linalg.norm(points - center, axis=1)temp_covered[distances <= r] = Trueif np.all(temp_covered):centers = temp_centersmasses = np.delete(masses, i)break# 输出最终的中心点
print('最终的中心点坐标:')
print(centers)# 可视化结果
plt.figure()
plt.scatter(points[:, 0], points[:, 1], s=20, c='b', marker='s', edgecolors='b')
plt.scatter(centers[:, 0], centers[:, 1], s=20, c='r', marker='^', edgecolors='r')
# 为每个中心点绘制覆盖半径的透明圆
theta = np.linspace(0, 2 * np.pi, 100)
for center in centers:x = center[0] + r * np.cos(theta)y = center[1] + r * np.sin(theta)plt.fill(x, y, 'r', alpha=0.1, edgecolor='r')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
结果图
matlab代码
clear
clc
close all% 第一步:K - means 聚类
% 假设我们有一组二维点数据 points
% 随机生成一些示例数据
points = rand(100, 2);
% 设置聚类中心数
k = 10;
% 执行 K - means 聚类
[idx, centers] = kmeans(points, k);% 第二步:贪心算法覆盖未被覆盖点
% 设置覆盖半径
r = 0.2;
% 标记未被覆盖的点
covered = false(size(points, 1), 1);
for i = 1:size(centers, 1)distances = sqrt(sum((points - repmat(centers(i, :), size(points, 1), 1)).^2, 2));covered(distances <= r) = true;
end
uncovered_points = points(~covered, :);while ~isempty(uncovered_points)max_coverage = 0;best_point = [];for i = 1:size(uncovered_points, 1)distances = sqrt(sum((points - repmat(uncovered_points(i, :), size(points, 1), 1)).^2, 2));current_coverage = sum(distances <= r & ~covered);if current_coverage > max_coveragemax_coverage = current_coverage;best_point = uncovered_points(i, :);endendcenters = [centers; best_point];distances = sqrt(sum((points - repmat(best_point, size(points, 1), 1)).^2, 2));covered(distances <= r) = true;uncovered_points = points(~covered, :);
end% 第三步:引力搜索算法(GSA)优化中心点
% 设置参数
T = 100; % 迭代次数
G0 = 1; % 引力常数
alpha = 2; % 质量衰减系数% 初始化质量
masses = zeros(size(centers, 1), 1);
for i = 1:size(centers, 1)distances = sqrt(sum((points - repmat(centers(i, :), size(points, 1), 1)).^2, 2));masses(i) = sum(distances <= r);
endfor t = 1:TG = G0 * exp(-alpha * t / T);for i = 1:size(centers, 1)force = zeros(1, 2);% 计算斥力for j = 1:size(centers, 1)if i ~= jdistance = sqrt(sum((centers(i, :) - centers(j, :)).^2));if distance > 0force = force + G * masses(j) * (centers(i, :) - centers(j, :)) / distance;endendend% 计算引力for j = 1:size(points, 1)distance = sqrt(sum((centers(i, :) - points(j, :)).^2));if distance <= 2*r && distance > 0force = force - G * (centers(i, :) - points(j, :)) / distance;endend% 移动中心点(先假设移动)new_center = centers(i, :) + force;% 检查移动后的中心点是否能覆盖所有数据点temp_centers = centers;temp_centers(i, :) = new_center;temp_covered = false(size(points, 1), 1);for j = 1:size(temp_centers, 1)distances = sqrt(sum((points - repmat(temp_centers(j, :), size(points, 1), 1)).^2, 2));temp_covered(distances <= r) = true;endif all(temp_covered)centers(i, :) = new_center;endend% 尝试剔除冗余的中心点for i = 1:size(centers, 1)temp_centers = centers;temp_centers(i, :) = [];temp_covered = false(size(points, 1), 1);for j = 1:size(temp_centers, 1)distances = sqrt(sum((points - repmat(temp_centers(j, :), size(points, 1), 1)).^2, 2));temp_covered(distances <= r) = true;endif all(temp_covered)centers(i, :) = [];masses(i) = [];break;endend
end% 输出最终的中心点
disp('最终的中心点坐标:');
disp(centers);% 可视化结果
figure;
scatter(points(:, 1), points(:, 2), 20, 'bs', 'filled');
hold on;
scatter(centers(:, 1), centers(:, 2), 20,'r^', 'filled');
% 为每个中心点绘制覆盖半径的透明圆
theta = linspace(0, 2*pi, 100);
for i = 1:size(centers, 1)x = centers(i, 1) + r * cos(theta);y = centers(i, 2) + r * sin(theta);h = fill(x, y, 'r');set(h, 'FaceAlpha', 0.1); % 设置透明度为 0.1set(h, 'EdgeColor', 'r');
end
xlabel('X')
ylabel('Y')
title('基于聚类算法的选址结果')
legend('数据点', '中心点', '覆盖区域', 'Location', 'northeast')