"""
kcenter_binary_search.py: use bisection for solving the k-center problem
bisects the interval [0, max facility-customer distance] until finding a
distance such that all customers are covered, but decreasing that distance
by a small amount delta would leave some uncovered.
Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
"""
from gurobipy import *
def kcover(I,J,c,k):
"""kcover -- minimize the number of uncovered customers from k facilities.
Parameters:
- I: set of customers
- J: set of potential facilities
- c[i,j]: cost of servicing customer i from facility j
- k: number of facilities to be used
Returns a model, ready to be solved.
"""
model = Model("k-center")
z,y,x = {},{},{}
for i in I:
z[i] = model.addVar(vtype="B", name="z(%s)"%i)
for j in J:
y[j] = model.addVar(vtype="B", name="y(%s)"%j)
for i in I:
x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))
model.update()
for i in I:
model.addConstr(quicksum(x[i,j] for j in J) + z[i] == 1, "Assign(%s)"%i)
for j in J:
model.addConstr(x[i,j] <= y[j], "Strong(%s,%s)"%(i,j))
model.addConstr(quicksum(y[j] for j in J) == k, "k_center")
model.setObjective(quicksum(z[i] for i in I), GRB.MINIMIZE)
model.update()
model.__data = x,y,z
return model
def solve_kcenter(I,J,c,k,delta):
"""solve_kcenter -- locate k facilities minimizing distance of most distant customer.
Parameters:
I - set of customers
J - set of potential facilities
c[i,j] - cost of servicing customer i from facility j
k - number of facilities to be used
delta - tolerance for terminating bisection
Returns:
- list of facilities to be used
- edges linking them to customers
"""
model = kcover(I,J,c,k)
x,y,z = model.__data
facilities,edges = [],[]
LB = 0
UB = max(c[i,j] for (i,j) in c)
while UB-LB > delta:
theta = (UB+LB) / 2.
# print "\n\ncurrent theta:", theta
for j in J:
for i in I:
if c[i,j]>theta:
x[i,j].UB = 0
else:
x[i,j].UB = 1.0
model.update()
# model.Params.OutputFlag = 0 # silent mode
model.Params.Cutoff = .1
model.optimize()
if model.status == GRB.Status.OPTIMAL:
# infeasibility = sum([z[i].X for i in I])
# print "infeasibility=",infeasibility
UB = theta
facilities = [j for j in y if y[j].X > .5]
edges = [(i,j) for (i,j) in x if x[i,j].X > .5]
# print "updated solution:"
# print "facilities",facilities
# print "edges",edges
else: # infeasibility > 0:
LB = theta
return facilities,edges
import math
import random
def distance(x1,y1,x2,y2):
return math.sqrt((x2-x1)**2 + (y2-y1)**2)
def make_data(n,m,same=True):
if same == True:
I = range(n)
J = range(m)
x = [random.random() for i in range(max(m,n))] # positions of the points in the plane
y = [random.random() for i in range(max(m,n))]
else:
I = range(n)
J = range(n,n+m)
x = [random.random() for i in range(n+m)] # positions of the points in the plane
y = [random.random() for i in range(n+m)]
c = {}
for i in I:
for j in J:
c[i,j] = distance(x[i],y[i],x[j],y[j])
return I,J,c,x,y
if __name__ == "__main__":
random.seed(67)
n = 200
m = n
I,J,c,x_pos,y_pos = make_data(n,m,same=True)
k = 20
delta = 1.e-4
facilities,edges = solve_kcenter(I,J,c,k,delta)
print "Selected facilities:",facilities
print "Edges:",edges
print "Max distance from a facility to a customer: ", max([c[i,j] for (i,j) in edges])
try: # plot the result using networkx and matplotlib
import networkx as NX
import matplotlib.pyplot as P
P.clf()
G = NX.Graph()
facilities = set(facilities)
unused = set(j for j in J if j not in facilities)
client = set(i for i in I if i not in facilities and i not in unused)
G.add_nodes_from(facilities)
G.add_nodes_from(client)
G.add_nodes_from(unused)
for (i,j) in edges:
G.add_edge(i,j)
position = {}
for i in range(len(x_pos)):
position[i] = (x_pos[i],y_pos[i])
NX.draw(G,position,with_labels=False,node_color="w",nodelist=facilities)
NX.draw(G,position,with_labels=False,node_color="c",nodelist=unused,node_size=50)
NX.draw(G,position,with_labels=False,node_color="g",nodelist=client,node_size=50)
P.show()
except ImportError:
print "install 'networkx' and 'matplotlib' for plotting"