Weighted Constrained K-Means

In this demo, the electricity consumption of 12 buildings is collected. The project aims to setup a number of microgrids powered by solar panel generators and provide electricity to these buildings. While minimization of wiring cost is considered, the optimum assignment has to satisfy the capacity of solar panel generators. In other words, the research question is to group the buildings into a number of clusters such that the aggregate inter-cluster distance is minimized within each cluster while keeping the aggregate electricity consumption (demand) of each cluster at a moderate level defined by upper/lower bound. I develop a weighted constrained k-means approach to find the optimum solution.

Inspired by Bradley et al. (2000), I formulate the constrained k-means clustering as a linear programming problem: \[min_{C,T} \quad \Sigma _{i=1}^m \Sigma _{h=1}^k T_{i,h} ||x^i-C^h||^2\] \[s.t. \quad \Sigma _{h=1}^k T_{i,h}=1, \quad i=1,...,m\] \[\quad \quad \quad \quad \quad \quad L_h \leq \Sigma _{i=1}^m \sigma _iT_{i,h}\leq U_h, \quad h=1,...,k\] \[\quad \quad \quad \quad \quad \quad T_{i,h} \geq 0, \quad i=1,...,m, h=1,...,k.\] Note that \(T_{i,h}=1\) if data point \(x^i\) is closest to centroid \(C^h\) and zero otherwise.\(\sigma _i\) is the demand for data point \(x^i\). The second constraint requires that the aggregate demand of each cluster should be bounded by the upper bound \(U_h\) and lower bound \(L_h\).

It’s worth mentioning that centroid \(C^{h,t+1}\) at iteration \(t+1\) is calculated using centroid at iteration \(t\). The algorithm contains two steps: 1) initiate \(C^h\) and run the LP problem to determine the cluster assignment \(T_{ih}^t\), and 2) update \(C^{h,t+1}=\frac{\Sigma _{i=1}^m T_{i,h}^t x^i}{\Sigma _{i=1}^m T_{i,h}^t}\) until convergence.

An Example

I deploy the model on a real-world dataset which consists the demand and location of 12 buildings.

coord <- matrix(c(40.555181921935095, -74.27603516627407,    # Create longitude/latitude matrix
                  40.55620901014459, -74.28006920857614,
                  40.555915557978196,-74.27901778265698,
                  40.560415016521716, -74.27509102871402,
                  40.560480223842745, -74.2754128937913,
                  40.55449718770875, -74.26743063987445,
                  40.55607858711832, -74.28137812655713,
                  40.560741052491544, -74.2889312270376,
                  40.559208669626244, -74.27987608952975,
                  40.55720347736548, -74.28328785934906,
                  40.55609489001049, -74.27689347314684,
                  40.55444827785318, -74.28251538316356),
                  ncol = 2,byrow = TRUE)

demand <- c(363,12,88,67,65,530,26,132,325,15,13,57)

data = cbind(coord,demand)
colnames(data)[1]="x"
colnames(data)[2]="y"
colnames(data)[3]="demand"
data = as.data.frame(data)

A scatter plot is created to visualize the dataset.

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
graph <- ggplot(data, aes(x=x, y=y,label=demand)) +geom_point(size=5)
graph +  geom_text(hjust=-0.33, vjust=0)

Set the number of clusters \(k\) and upper/lower bound for aggregate demand in each cluster

num_cluster =3
bu = c(600,600,600)
bl = c(100,100,100)

Run the weighted constrained K-means algorithm

library(lpSolve)
## Warning: package 'lpSolve' was built under R version 4.2.2
library(geodist) 
## Warning: package 'geodist' was built under R version 4.2.3
W = matrix(demand,nrow=1)
Au = kronecker(diag(num_cluster), W) 
Al = kronecker(diag(num_cluster), W) 

zero = rep(0,length(demand))
zero[1] = 1
B = kronecker(diag(num_cluster),matrix(zero,nrow=1))
one = matrix(rep(1,num_cluster),nrow=1)
B = one%*%B
for (x in 2:length(demand)){
  zero = rep(0,length(demand))
  zero[x] = 1
  temp = kronecker(diag(num_cluster),matrix(zero,nrow=1))
  temp = one%*%temp
  B = rbind(B,temp)
}

b= rep(1,length(demand))

left = rbind(Au,Al,B)
right =c(bu,bl,b)
direction = c(rep("<=",num_cluster),rep(">=",num_cluster),rep("=",length(demand)))

set.seed(1)
cx <- runif(n = num_cluster, min = min(coord[,1]), max = max(coord[,1]))
cy <- runif(n = num_cluster, min = min(coord[,2]), max = max(coord[,2]))
cluster = cbind(cx,cy)


adj = as.matrix(geodist(rbind(cluster[1,],coord)))
coef = adj[1,2:(length(demand)+1)]
for (x in 2:num_cluster){
  adj = as.matrix(geodist(rbind(cluster[x,],coord)))
  coef = c(coef,adj[1,2:(length(demand)+1)])
}

optimum <-  lp(direction="min",
               objective.in = coef,
               const.mat = left,
               const.dir = direction,
               const.rhs = right,
               all.bin = TRUE)

result = optimum$solution

for (i in 1:80){
  #print(i)
  
  xnew = coord[,1]*result
  ynew = coord[,2]*result
  cnewx = c()
  cnewy = c()
  for (x in 1:num_cluster){
    tmpx = sum(xnew[((x-1)*length(demand)+1):((x)*length(demand))])/sum(xnew[((x-1)*length(demand)+1):((x)*length(demand))] != 0)
    tmpy = sum(ynew[((x-1)*length(demand)+1):((x)*length(demand))])/sum(ynew[((x-1)*length(demand)+1):((x)*length(demand))] != 0)
    cnewx = c(cnewx,tmpx)
    cnewy = c(cnewy,tmpy)
  }
  clusternew = cbind(cnewx,cnewy)
  
  adj = as.matrix(geodist(rbind(clusternew[1,],coord)))
  coef = unname(adj[1,2:(length(demand)+1)])
  
  for (j in 2:num_cluster){
    adj = as.matrix(geodist(rbind(clusternew[j,],coord)))
    coef = c(coef,unname(adj[1,2:(length(demand)+1)]))
  }
  
  optimum <-  lp(direction="min",
                 objective.in = coef,
                 const.mat = left,
                 const.dir = direction,
                 const.rhs = right,
                 all.bin = TRUE)
  result = optimum$solution
  
}

group =rep(0,length(demand))
for (i in 1:num_cluster){
  group = group + result[((i-1)*length(demand)+1):((i)*length(demand))]*i
}

data$k <- kmeans(data[,1:2],3)$cluster
data$group <- group

The buildings have been grouped into the following three clusters. It can be shown that the capacity constraints are satisfied for each of the cluster,i.e., \(demand_{\ group1}=530\), \(demand_{\ group2}=567\), and \(demand_{\ group3}=596\). The distribution of electricity supply has achieved a balance across mocrogrids.

graph <- ggplot(data, aes(x=x, y=y,color=group,label=demand)) +geom_point(size=5)
graph +  geom_text(hjust=-0.33, vjust=0)

We can also compare the result with original K-means clustering. It can be shown that the capacity constraint is not satisfied for cluster 2 (dark blue) where aggregate demand exceeds 600.

graph <- ggplot(data, aes(x=x, y=y,color=k,label=demand)) +geom_point(size=5)
graph +  geom_text(hjust=-0.33, vjust=0)

Extensions

There are several directions to extend the model:

  1. Introduce heterogeneous solar panel generators with different capacity and setup cost.

  2. Consider the electricity gap/surplus between microgrids and state grids and model the electricity trading.

  3. Introduce temporal features such as seasonality to model the dynamics of electricity consumption.

References

Bradley, P. S., Bennett, K. P., & Demiriz, A. (2000). Constrained k-means clustering. Microsoft Research, Redmond, 20(0), 0.