带约束的K-means 聚类算法

1. 前言

上一期学习了K-means聚类算法，聚类是不受到限制的，单纯的无监督学习，但是当存在一些约束时，比如对每一簇的聚类样本点数量有限制，或者每个样本点带需求，每一簇存在满足约束的最大能力，这时候应该怎么添加约束，并且不破坏聚类的效果。

2. 场景构建

假设场景：有一定数量的区域需要分配邮递员配送信件，怎么将区域划分给邮递员？每个区域有一定的需求量，一个邮递员一天的工作量只能满足一定的需求，那么应该怎么将区域划分给邮递员？考虑到成本问题，这里设置一个目标——在满足需求的情况下最小化邮递员数量。

这样的场景在快递行业应该很常见，不太清楚实际业务中是怎么分配的，但我觉得带约束的聚类算法应该能求解一个可行解，为了给K-means加约束构想的场景，场景简化了很多，可能不太符合实际，欢迎提建议。

3. 带约束的K-means 聚类设计

K-means聚类算法流程

①算法说明：由于存在需求约束，那么就存在最小的K（总需求/簇的能力），可以从最小值K开始遍历。而约束就加在（3）将对象分配到最相似的簇中，具体修改为：在满足约束的情况下，将对象分配到最相似的簇中，如果最相似的簇无法容纳，就找第二相似的簇，依次类推。

②遍历改进：对遍历对象进行一个小的修改，通常K-means是随机/按列表顺序遍历对象，将对象分配到最相似的簇，这里改进为：取所有对象中距离某个簇最近的点，优先分配（这算是随机/按列表顺序遍历对象的一个特殊情况）。

③初始解敏感：为了减小K-means对初始解敏感的问题，这里采用多次求解的方式保留最优的方式降低对初始解的依赖。

④效果评价：聚类效果使用轮廓系数进行评价。

3. python实现

# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
import random
import math
from sklearn.metrics import silhouette_score
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 添加这条可以让图形显示中文def calDistance(centers,nodes):''' 计算节点间距离输入：centers-中心，nodes-节点；输出：距离矩阵-dis_matrix'''dis_matrix = pd.DataFrame(data=None,columns=range(len(centers)),index=range(len(nodes)))for i in range(len(nodes)):xi,yi = nodes[i][0],nodes[i][1]for j in range(len(centers)):xj,yj = centers[j][0],centers[j][1]dis_matrix.iloc[i,j] = round(math.sqrt((xi-xj)**2+(yi-yj)**2),2)return dis_matrixdef scatter_diagram(clusters,nodes):'''#画路径图输入：nodes-节点坐标；输出：散点图'''for cluster in clusters:x,y= [],[]for Coordinate in cluster:x.append(Coordinate[0])y.append(Coordinate[1])plt.scatter(x, y,alpha=0.8,)plt.xlabel('x')plt.ylabel('y')plt.show()def distribute(center,nodes,K,demand,d_limit):'''将节节点分给最近的中心输入：center-中心,nodes-节点,K-类数量,demand-需求,d_limit-一个簇的满足需求的能力输出：新的簇-clusters，簇的需求-clusters_d'''clusters = [[] for i in range(K)]#簇label = [None for i in range(len(nodes))]clusters_d = [0 for i in range(K)]#簇需求dis_matrix = calDistance(center,nodes).astype('float64')for i in range(len(dis_matrix)):row, col = dis_matrix.stack().idxmin()node = nodes[row]j = 1while clusters_d[col] + demand[row] > d_limit:#检验约束if j < K: j += 1dis_matrix.loc[row,col] = math.pow(10,10)col = dis_matrix.loc[row,:].idxmin()else:#所有类都不满足需求量约束print("K较小")scatter_diagram(clusters,nodes)return None#满足约束正常分配clusters[col].append(node)label[row] = colclusters_d[col] += demand[row]dis_matrix.loc[row,:] = math.pow(10,10)return clusters,clusters_d,labeldef cal_center(clusters):'''计算簇的中心输入：clusters-类坐标；输出：簇中心-new_center'''new_center = []for cluster in clusters:x,y= [],[]for Coordinate in cluster:x.append(Coordinate[0])y.append(Coordinate[1])x_mean = np.mean(x)y_mean = np.mean(y)new_center.append((round(x_mean,2),round(y_mean,2)))return new_centerif __name__ == "__main__":#输入数据nodes = [(11, 36),(13, 35),(14, 12),(25, 22),(15, 11),(16, 24),(14, 30),(16, 32),(26, 10),(27, 31),(22, 29),(1, 29),(20, 29),(17, 15),(5, 38),(3, 15),(22, 20),(4, 13),(4, 22),(12, 18),(23, 12),(18, 15),(5, 13),(9, 27),(17, 36),(2, 14),(1, 16),(7, 15),(25, 15),(19, 26),(25, 40),(26, 34),(25, 35),(2, 36),(29, 24),(17, 17),(8, 26),(4, 14),(5, 25),(6, 37),(1, 14),(6, 39),(11, 13),(10, 20),(21, 11),(5, 19),(5, 35),(1, 34),(16, 39),(19, 24),(39, 31),(49, 31),(41, 50),(31, 33),(32, 40),(35, 30),(31, 39),(34, 48),(42, 32),(32, 35),(35, 33),(35, 34),(43, 41),(35, 47),(49, 36),(37, 41),(43, 46),(41, 41),(45, 50),(41, 35),(45, 44),(41, 30),(43, 33),(31, 45),(48, 32),(39, 49),(38, 42),(33, 39),(49, 33),(43, 44),(32, 30),(40, 47),(36, 46),(47, 47),(37, 33),(35, 31),(42, 38),(43, 47),(30, 47),(30, 30),(37, 34),(41, 45),(27, 33),(42, 39),(43, 43),(50, 43),(28, 40),(35, 41),(32, 41),(31, 30)]demand = [3,3,5,9,6,1,4,8,8,3,6,5,6,4,1,4,7,1,3,6,2,5,7,2,2,1,9,3,7,8,4,3,1,2,5,1,3,4,8,5,2,9,3,10,6,1,9,3,5,3,3,3,9,9,6,2,5,2,4,10,4,10,10,6,3,7,9,4,2,9,7,4,5,3,3,8,2,3,2,9,1,3,3,3,9,5,7,8,8,5,2,8,3,4,2,4,6,3,7,4]scatter_diagram([list(nodes)],nodes)#初始位置分布图#参数d_limit = 150 # 一个区的约束K = math.ceil(sum(demand)/d_limit)  #簇的最小边界#历史最优参数best_s = -1best_center = 0best_clusters = 0best_labels = 0best_clusters_d = 0for n in range(20):#多次遍历，kmeans对初始解敏感#初始化随机生成簇中心index = random.sample(list(range(len(nodes))),K)new_center = [nodes[i] for i in index]i = 1center_list = []while True:#节点——> 簇center = new_center.copy()center_list.append(center)#保留中心，避免出现A生成B，B生成C，C有生成A的情况clusters,cluster_d,label = distribute(center,nodes,K,demand,d_limit)new_center = cal_center(clusters)if (center == new_center) | (new_center in center_list):breaki += 1s = silhouette_score(nodes,label, metric='euclidean') # 计算轮廓系数，聚类的评价指标# scatter_diagram(clusters,nodes)#当前最优图if best_s < s:best_s = sbest_clusters_d = cluster_dbest_center = center.copy()best_clusters = clusters.copy()best_labels = label.copy()scatter_diagram(best_clusters,nodes)#历史最优图# print("各类需求量：",best_clusters_d)

【结果展示】

图1 需求点位置

图2 聚类效果

4 谈论

从结果上看，带约束的聚类效果还是可以的，但是求解过程中会出现图3的聚类效果，蓝色的点距离黄色的点反而更近一点。而且还存在另外一个问题，假设需求量是500，而每簇满足需求的能力是100，理论上最少需要5个簇就可以，但是由于约束的存在和需求不可拆分，可能会导致无法刚好满足五个簇的容量，遇到这种情况我写的算法中是直接报错，虽然这种情况下大概率是出现如图3蓝色点的情况存在，但是也会减小求解最优解的可能性，程序的健壮性较差。

图3 聚类效果