4 minimize the total (weighted) travel cost for servicing 5 a set of customers from k facilities. 7 Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012 11 from pyscipopt
import Model, quicksum, multidict
14 """kmedian -- minimize total cost of servicing customers from k facilities 17 - J: set of potential facilities 18 - c[i,j]: cost of servicing customer i from facility j 19 - k: number of facilities to be used 20 Returns a model, ready to be solved. 23 model = Model(
"k-median")
27 y[j] = model.addVar(vtype=
"B", name=
"y(%s)"%j)
29 x[i,j] = model.addVar(vtype=
"B", name=
"x(%s,%s)"%(i,j))
32 model.addCons(quicksum(x[i,j]
for j
in J) == 1,
"Assign(%s)"%i)
34 model.addCons(x[i,j] <= y[j],
"Strong(%s,%s)"%(i,j))
35 model.addCons(quicksum(y[j]
for j
in J) == k,
"Facilities")
37 model.setObjective(quicksum(c[i,j]*x[i,j]
for i
in I
for j
in J),
"minimize")
43 def distance(x1,y1,x2,y2):
44 """return distance of two points""" 45 return math.sqrt((x2-x1)**2 + (y2-y1)**2)
48 def make_data(n,m,same=True):
49 """creates example data set""" 53 x = [random.random()
for i
in range(max(m,n))]
54 y = [random.random()
for i
in range(max(m,n))]
58 x = [random.random()
for i
in range(n+m)]
59 y = [random.random()
for i
in range(n+m)]
63 c[i,j] = distance(x[i],y[i],x[j],y[j])
68 if __name__ ==
"__main__":
73 I,J,c,x_pos,y_pos = make_data(n,m,same=
True)
80 edges = [(i,j)
for (i,j)
in x
if model.getVal(x[i,j]) > EPS]
81 facilities = [j
for j
in y
if model.getVal(y[j]) > EPS]
83 print(
"Optimal value:",model.getObjVal())
84 print(
"Selected facilities:", facilities)
85 print(
"Edges:", edges)
86 print(
"max c:", max([c[i,j]
for (i,j)
in edges]))
90 import matplotlib.pyplot
as P
94 facilities = set(j
for j
in J
if model.getVal(y[j]) > EPS)
95 other = set(j
for j
in J
if j
not in facilities)
96 client = set(i
for i
in I
if i
not in facilities
and i
not in other)
97 G.add_nodes_from(facilities)
98 G.add_nodes_from(client)
99 G.add_nodes_from(other)
104 for i
in range(len(x_pos)):
105 position[i] = (x_pos[i],y_pos[i])
107 NX.draw(G,position,with_labels=
False,node_color=
"w",nodelist=facilities)
108 NX.draw(G,position,with_labels=
False,node_color=
"c",nodelist=other,node_size=50)
109 NX.draw(G,position,with_labels=
False,node_color=
"g",nodelist=client,node_size=50)
112 print(
"install 'networkx' and 'matplotlib' for plotting")