基于聚类与引力斥力优化的选址算法

在众多实际场景中，诸如消防设施选址、基站布局规划以及充电桩站点部署等，都面临着如何利用最少的资源，实现对所有目标对象全面覆盖的难题。为有效解决这类问题，本文提出一种全新的组合算法模型 —— 基于聚类与引力斥力优化的选址算法。

算法简介

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')

结果图