""" center.py: Using Gurobi for solving the k-center problem: bisection on the minimum max distance. Problem: minimize the maximum travel distance from n customers to k facilities. Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2010 """ from gurobipy import * def kcenter(m, n, c, k, max_c): model = Model("k-center") z, y, x = {}, {}, {} for i in range(n): z[i] = model.addVar(obj=1, vtype="B", name="z[%s]"%i) for j in range(m): y[j] = model.addVar(obj=0, vtype="B", name="y[%s]"%j) for i in range(n): x[i,j] = model.addVar(obj=0, vtype="B", name="x[%s,%s]"%(i,j)) model.update() for i in range(n): coef = [1 for j in range(m+1)] var = [x[i,j] for j in range(m)] var.append(z[i]) model.addConstr(LinExpr(coef,var), "=", 1, name="Assign[%s]"%i) for j in range(m): for i in range(n): model.addConstr(x[i,j], "<", y[j], name="Strong[%s,%s]"%(i,j)) coef = [1 for j in range(m)] var = [y[j] for j in range(m)] model.addConstr(LinExpr(coef,var), "=", rhs=k, name="k_center") model.update() model.__data = x,y,z return model def solve_kcenter(m, n, c, k, max_c, delta): model = kcenter(m, n, c, k, max_c) x,y,z = model.__data nodes = [] edges = [] LB = 0 UB = max_c while UB-LB > delta: theta = (UB+LB) / 2. # print "\n\ncurrent theta:", theta for j in range(m): for i in range(n): 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.optimize() infeasibility = sum([z[i].X for i in range(m)]) # print "infeasibility=",infeasibility if infeasibility > 0: LB = theta else: UB = theta nodes = [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 "nodes", nodes # print "edges", edges return nodes, edges import math import random def distance(x1, y1, x2, y2): return math.sqrt((x2-x1)**2 + (y2-y1)**2) def make_data(n): # positions of the points in the plane x = [random.random() for i in range(n)] y = [random.random() for i in range(n)] c = {} max_c = 0. for i in range(n): for j in range(n): c[i,j] = distance(x[i],y[i],x[j],y[j]) max_c = max(c[i,j],max_c) return c, x, y, max_c if __name__ == "__main__": random.seed(67) n = 200 c, x_pos, y_pos, max_c = make_data(n) m = n k = 20 delta = 1.e-4 nodes, edges = solve_kcenter(m, n, c, k, max_c, delta) print "Selected nodes:", nodes print "Edges:", edges print [((i,j),c[i,j]) for (i,j) in edges] print "max c:", 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.ion() # interactive mode on G = NX.Graph() other = [j for j in range(m) if j not in nodes] G.add_nodes_from(nodes) G.add_nodes_from(other) for (i,j) in edges: G.add_edge(i,j) position = {} for i in range(n): position[i]=(x_pos[i],y_pos[i]) NX.draw(G, position, node_color='y', nodelist=nodes) NX.draw(G, position, node_color='g', nodelist=other) raw_input("press [enter] to continue") except ImportError: print "install 'networkx' and 'matplotlib' for plotting"