Weighted Constrained K-Means

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 at minimum cost. While wiring cost is considered, the optimal assignment has to satisfy the capacities of varying 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 solve this problem.

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 wT_{i,h} ||x^i-C^h||^2 + \Sigma _{h=1}^k \Sigma _{g=1}^n \rho _gS_{h,g} \] \[s.t. \quad \Sigma _{h=1}^k T_{i,h}=1, \quad i=1,...,m\] \[\quad \quad \quad \quad \quad \quad \quad \quad \Sigma _{i=1}^m d _{i,t}T_{i,h}\leq \Sigma _{g=1}^n U _gS_{h,g}, \quad h=1,...,k \quad t=1,...,\tau\] \[\quad \quad \quad \quad \quad \quad \quad \quad \Sigma _{i=1}^m d _{i,t}T_{i,h} \geq \Sigma _{g=1}^n L _gS_{h,g}, \quad h=1,...,k \quad t=1,...,\tau\] \[\quad \quad \quad \quad \Sigma _{g=1}^n S_{h,g}=1, \quad h=1,...,k.\] \[\quad \quad \quad \quad \quad \quad T_{i,h} \geq 0, \quad i=1,...,m, h=1,...,k.\] I adopt the following notations: \(i \in \{1,2,3,...,m\}\): index of buildings

\(h \in \{1,2,3,...,k\}\): index of clusters

\(t \in \{1,2,3,...,\tau\}\): index of time stamp

\(g \in \{1,2,3,...,n\}\): index of generator types

\(C^h\): the centroid of cluster \(h\)

\(x^i\): the position of building \(i\)

\(w\): wiring cost per kilometer

\(\rho _g\): generator set-up cost of type \(g\)

\(L_g\) and \(U_g\): lower and upper limit of type \(g\) generator

\(d_{i,t}\): electricity consumption (demand) of building \(i\) at time \(t\) in Kw

\(S_{h,g}\): decision variable, equals one if cluster \(h\) adopts type \(g\) generator

\(T_{i,h}\): decision variable, equals one if data point \(x^i\) is closest to centroid \(C^h\) and zero otherwise

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)


d1 = c(283, 9,  30, 67, 56, 530,    24, 37, 98,   2,    6,  48) #demand for month 1
d2 = c(310, 9,  34, 60, 65, 491,    26, 32, 118,    5,  7,  43)
d3 = c(336, 8,  31, 50, 57, 465,    22, 34, 100,    5,  7,  41)
d4 = c(323, 8,  26, 49, 31, 472,    23, 33, 130,    4,  8,  45)
d5 = c(329, 7,  30, 28, 26, 385,    18, 40, 112,    6,  7,  40)
d6 = c(351, 10, 41, 16, 27, 436,    19, 67, 325,    12, 11, 47)
d7 = c(363, 12, 88, 17, 31, 422,    18, 132,    31, 15, 13, 54)
d8 = c(351, 12, 61, 21, 23, 371,    7,  31, 174,    14, 13, 52)
d9 = c(329, 12, 36, 26, 27, 404,    11, 121,    311,    14, 12, 57)
d10 = c(309,    10, 32, 20, 21, 297,    24, 62, 168,    4,  9,  47)
d11 = c(301,    9,  33, 32, 33, 352,    17, 7,  36, 12, 8,  44)
d12 = c(322,    8,  34, 46, 47, 347,    21, 11, 48, 2,  6,  40)

d = cbind(d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12)

time = dim(d)[2]
building = dim(d)[1]

cost_g = c(10000,10000,10000)
lower = c(50,100,300)
upper = c(100,300,1000)

cost_g = matrix(cost_g,nrow=1)
generator = length(upper)
upper = matrix(upper,nrow=1)
lower = matrix(lower,nrow=1)


demand <- c(363,12,88,67,65,530,26,132,325,15,13,57)  # this is peak demand calculated as max for each building during the 12-month period.

data = cbind(coord,demand)
data = cbind(data,d)
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

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

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

Run the weighted constrained K-means algorithm

# constraints on demand capacity using vector of decision variables x1 (mk by 1). Be aware that x2 will be included later to accommodate varying capacities.
km <- function(num_cluster){

W = matrix(t(d),nrow=time)
Au = kronecker(diag(num_cluster), W) 
Al = kronecker(diag(num_cluster), W) 
#bu = c(rep(bu[1],time),rep(bu[2],time),rep(bu[3],time))
#bl = c(rep(bl[1],time),rep(bl[2],time),rep(bl[3],time))

# set up constraints to make sure each data point is assigned to one cluster only using vector of decision variables x1 (mk by 1)
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)){
  print(x)
  zero = rep(0,length(demand))
  zero[x] = 1
  temp = kronecker(diag(num_cluster),matrix(zero,nrow=1))
  temp = one%*%temp
  B = rbind(B,temp)
}


#set up constraints for generator-cluster decision variables using vector of decision variables x2 (nk by 1)
#a cluster can only choose a generator of one type
one = matrix(rep(1,generator),nrow=1)
S = kronecker(diag(num_cluster),one)
s = rep(1,num_cluster)

################################################
select = matrix(0, ncol = num_cluster*(generator), nrow = time*num_cluster) #upper 
left11 = cbind(Au,select)
one = matrix(1,nrow=time,ncol=1)
select = matrix(0, ncol = num_cluster*(building), nrow = time*num_cluster)
left12 = cbind(select,kronecker(diag(num_cluster),kronecker(one,upper)))
left1 = left11-left12
right1 = rep(0,(time*num_cluster))
direction1 = rep("<=",(time*num_cluster))


select = matrix(0, ncol = num_cluster*(generator), nrow = time*num_cluster) #lower
left21 = cbind(Au,select)
one = matrix(1,nrow=time,ncol=1)
select = matrix(0, ncol = num_cluster*(building), nrow = time*num_cluster)
left22 = cbind(select,kronecker(diag(num_cluster),kronecker(one,lower)))
left2 = left21-left22
right2 = rep(0,(time*num_cluster))
direction2 = rep(">=",(time*num_cluster))

select = matrix(0, ncol = num_cluster*(generator), nrow = building)
left3 = cbind(B,select)
right3 = rep(1,length(demand))
direction3 = rep("=",length(demand))

select = matrix(0, ncol = num_cluster*(building), nrow = num_cluster)
left4 = cbind(select,S)
right4 = s
direction4 = rep("=",num_cluster)



left = rbind(left1,left2,left3,left4)
right =c(right1,right2,right3,right4)
direction = c(direction1,direction2,direction3,direction4)

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_init = cbind(cx,cy)


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

one = matrix(1,nrow=1,ncol=num_cluster)
coef = c(coef,kronecker(one,cost_g))

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)
  position = result[1:(building*num_cluster)]
  xnew = coord[,1]*position
  ynew = coord[,2]*position
  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)
  }
  cluster = cbind(cnewx,cnewy)
  
  adj = as.matrix(geodist(rbind(cluster[1,],coord)))
  coef = unname(adj[1,2:(length(demand)+1)])
  
  for (j in 2:num_cluster){
    adj = as.matrix(geodist(rbind(cluster[j,],coord)))
    coef = c(coef,unname(adj[1,2:(length(demand)+1)]))
  }
  one = matrix(1,nrow=1,ncol=num_cluster)
  coef = c(coef,kronecker(one,cost_g))
  
  optimum <-  lp(direction="min",
                 objective.in = coef,
                 const.mat = left,
                 const.dir = direction,
                 const.rhs = right,
                 all.bin = TRUE)
  result = optimum$solution
  
  position = result[1:(building*num_cluster)]
  xnew = coord[,1]*position
  ynew = coord[,2]*position
  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)
  
  if (norm(cluster-clusternew)<0.00001){
    break
  }
  
  #print(i)
  
}
  
return (result%*%coef)}

Run the algorithm with k values ranging from 2 to 4. Please note that the LP problem is not feasible if \(k \geq 5\) because the lower bound of capacity constraints can not be satisfied if there are too many clusters.

obj = c()
for (i in 2:4){
  obj = c(obj,km(i))
}

## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12

k_cost = cbind(obj,c(2,3,4))
colnames(k_cost)[2]="k"
k_cost = as.data.frame(k_cost)
ggplot(k_cost, aes(x=k)) + 
  geom_line(aes(y = obj), color = "green") + scale_x_continuous(breaks =seq(1,3,1)) + ylab("total cost")

The result indicates 2 clusters in the optimal assignment.

num_cluster = 2
W = matrix(t(d),nrow=time)
Au = kronecker(diag(num_cluster), W) 
Al = kronecker(diag(num_cluster), W) 
#bu = c(rep(bu[1],time),rep(bu[2],time),rep(bu[3],time))
#bl = c(rep(bl[1],time),rep(bl[2],time),rep(bl[3],time))

# set up constraints to make sure each data point is assigned to one cluster only using vector of decision variables x1 (mk by 1)
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)){
  print(x)
  zero = rep(0,length(demand))
  zero[x] = 1
  temp = kronecker(diag(num_cluster),matrix(zero,nrow=1))
  temp = one%*%temp
  B = rbind(B,temp)
}

## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12

#set up constraints for generator-cluster decision variables using vector of decision variables x2 (nk by 1)
#a cluster can only choose a generator of one type
one = matrix(rep(1,generator),nrow=1)
S = kronecker(diag(num_cluster),one)
s = rep(1,num_cluster)

################################################
select = matrix(0, ncol = num_cluster*(generator), nrow = time*num_cluster) #upper 
left11 = cbind(Au,select)
one = matrix(1,nrow=time,ncol=1)
select = matrix(0, ncol = num_cluster*(building), nrow = time*num_cluster)
left12 = cbind(select,kronecker(diag(num_cluster),kronecker(one,upper)))
left1 = left11-left12
right1 = rep(0,(time*num_cluster))
direction1 = rep("<=",(time*num_cluster))


select = matrix(0, ncol = num_cluster*(generator), nrow = time*num_cluster) #lower
left21 = cbind(Au,select)
one = matrix(1,nrow=time,ncol=1)
select = matrix(0, ncol = num_cluster*(building), nrow = time*num_cluster)
left22 = cbind(select,kronecker(diag(num_cluster),kronecker(one,lower)))
left2 = left21-left22
right2 = rep(0,(time*num_cluster))
direction2 = rep(">=",(time*num_cluster))

select = matrix(0, ncol = num_cluster*(generator), nrow = building)
left3 = cbind(B,select)
right3 = rep(1,length(demand))
direction3 = rep("=",length(demand))

select = matrix(0, ncol = num_cluster*(building), nrow = num_cluster)
left4 = cbind(select,S)
right4 = s
direction4 = rep("=",num_cluster)



left = rbind(left1,left2,left3,left4)
right =c(right1,right2,right3,right4)
direction = c(direction1,direction2,direction3,direction4)

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_init = cbind(cx,cy)


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

one = matrix(1,nrow=1,ncol=num_cluster)
coef = c(coef,kronecker(one,cost_g))

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)
  position = result[1:(building*num_cluster)]
  xnew = coord[,1]*position
  ynew = coord[,2]*position
  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)
  }
  cluster = cbind(cnewx,cnewy)
  
  adj = as.matrix(geodist(rbind(cluster[1,],coord)))
  coef = unname(adj[1,2:(length(demand)+1)])
  
  for (j in 2:num_cluster){
    adj = as.matrix(geodist(rbind(cluster[j,],coord)))
    coef = c(coef,unname(adj[1,2:(length(demand)+1)]))
  }
  one = matrix(1,nrow=1,ncol=num_cluster)
  coef = c(coef,kronecker(one,cost_g))
  
  optimum <-  lp(direction="min",
                 objective.in = coef,
                 const.mat = left,
                 const.dir = direction,
                 const.rhs = right,
                 all.bin = TRUE)
  result = optimum$solution
  
  position = result[1:(building*num_cluster)]
  xnew = coord[,1]*position
  ynew = coord[,2]*position
  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)
  
  if (norm(cluster-clusternew)<0.00001){
    break
  }}

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],num_cluster)$cluster
data$group <- group

generator_assign = result[(1+num_cluster*building):(num_cluster*(building+generator))]*rep(c(1,2,3),num_cluster)
generator_assign = kronecker(diag(num_cluster),matrix(1,nrow=1,ncol=generator))%*%generator_assign
print(generator_assign)

##      [,1]
## [1,]    3
## [2,]    3

In the optimal assignment, both cluster 1 and cluster 2 adopt type-3 generator.

We create the time series of electricity consumption for each cluster as follows.It can be shown that the variation of electricity consumption in each cluster is constrained by the pre-defined upper and lower bound during the 12-month period. Both cluster 1 and cluster 2 are bounded by the capacity of type-3 generator from 300 to 1000 Kw.

library(doBy)

## Warning: package 'doBy' was built under R version 4.2.3

ts = t(summaryBy(d ~ group, data=data,FUN=sum)[,2:(time+1)])
ts = cbind(ts,as.integer(seq(1,time,1)))
colnames(ts)[num_cluster+1] = "month"
colnames(ts)[1]="g1"
colnames(ts)[2]="g2"
ts = as.data.frame(ts)

ggplot(ts, aes(x=month)) + 
  geom_line(aes(y = g1), color = "green") + 
  geom_line(aes(y = g2), color="black", linetype="twodash") + scale_x_continuous(breaks =seq(1,time,1)) + ylab("demand")

We also visualize the clustering result using the electricity consumption of each building at \(t=1\).

library(gridExtra)

## Warning: package 'gridExtra' was built under R version 4.2.3

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

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

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

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

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

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

#grid.arrange(plot1, plot2, plot3,plot4,plot5,plot6, ncol=3, nrow = 2)

We can compare the result with original K-means clustering.

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

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.

References

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