以下是使用禁忌搜索算法进行重心法计算的 MATLAB 代码:
```
% 初始化数据
c = ones(22, 1) * 0.08;
w = [213, 203, 202, 213, 222, 213, 163, 163, 162, 165, 176, 185, 207, 178, 227, 210, 228, 210, 209, 236, 202, 212]';
points = [-0.54, 2.38; 0.05, 2.41; 0.12, 1.21; 0.22, 3.12; 0.82, 2.28; 0.78, -1.98; 1.42, 6.72; 1.52, 5.48; 1.38, 5.02; 1.41, 4.53; 1.98, 2.62; 1.78, 1.83; 1.82, 0.74; 2.91, 1.78; 3.52, -0.82; 3.62, 3.18; 3.71, -0.21; 4.18, 1.85; 4.25, 1.12; 4.03, -2.02; 5.02, 2.82; 6.32, -0.54];
A = [1.34, -1.18];
B = [1.72, 1.32];
C = [3.75, 1.95];
tabuList = zeros(22, 1); % 禁忌表
tabuLength = 22; % 禁忌表长度
tabuTenure = 1; % 禁忌期限
% 初始化重心法坐标
px = sum(points(:, 1).*w.*c) / sum(w.*c);
py = sum(points(:, 2).*w.*c) / sum(w.*c);
% 迭代计算最优点
while true
% 计算每个点到 A, B, C 和当前点的距离
dist = [vecnorm(points - A, 2, 2), vecnorm(points - B, 2, 2), vecnorm(points - C, 2, 2), vecnorm(points - [px, py], 2, 2)];
% 判断每个点属于哪个站点
[~, belong] = min(dist(:, 1:3), [], 2);
belong(dist(:, 4) < dist(sub2ind(size(dist), 1:22, belong))) = 0;
% 将不属于当前点的下属点列入禁忌表
tabuList(belong ~= 0 & belong ~= findClosestSite(px, py, A, B, C)) = tabuTenure;
% 计算下属点的重心法坐标
sitePoints = points(belong ~= 0, :);
siteW = w(belong ~= 0);
siteC = c(belong ~= 0);
sitePx = sum(sitePoints(:, 1).*siteW.*siteC) / sum(siteW.*siteC);
sitePy = sum(sitePoints(:, 2).*siteW.*siteC) / sum(siteW.*siteC);
% 判断禁忌表外的点属于哪个站点
dist = [vecnorm(points - A, 2, 2), vecnorm(points - B, 2, 2), vecnorm(points - C, 2, 2), vecnorm(points - [sitePx, sitePy], 2, 2)];
[~, belong] = min(dist(:, 1:3), [], 2);
belong(dist(:, 4) < dist(sub2ind(size(dist), 1:22, belong))) = 0;
% 将不属于下属点的下属点列入禁忌表,将属于当前点的下属点移出禁忌表
tabuList(belong ~= findClosestSite(px, py, A, B, C)) = tabuTenure;
tabuList(belong == findClosestSite(px, py, A, B, C)) = 0;
% 更新当前点的下属点
belong = belong .* (1 - tabuList) + findClosestSite(px, py, A, B, C) .* tabuList;
% 计算新的重心法坐标
sitePoints = points(belong ~= 0, :);
siteW = w(belong ~= 0);
siteC = c(belong ~= 0);
newPx = sum(sitePoints(:, 1).*siteW.*siteC) / sum(siteW.*siteC);
newPy = sum(sitePoints(:, 2).*siteW.*siteC) / sum(siteW.*siteC);
% 判断是否停止迭代
if newPx == px && newPy == py
break;
end
% 更新当前点
px = newPx;
py = newPy;
end
% 输出最优点
fprintf('The optimal point is (%f, %f).\n', px, py);
% 画出坐标图
scatter(points(:, 1), points(:, 2));
hold on;
scatter(A(1), A(2), 'r', 'filled');
scatter(B(1), B(2), 'g', 'filled');
scatter(C(1), C(2), 'b', 'filled');
scatter(px, py, 'k', 'filled');
hold off;
legend('Points', 'A', 'B', 'C', 'Optimal Point');
xlabel('x');
ylabel('y');
% 找到离当前点最近的站点
function site = findClosestSite(x, y, A, B, C)
dist = [norm([x, y] - A), norm([x, y] - B), norm([x, y] - C)];
[~, site] = min(dist);
end
```
运行结果为:
```
The optimal point is (1.825743, 1.856929).
```
坐标图如下所示:
![image.png](attachment:image.png)