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")
44 """return distance of two points"""
45 return math.sqrt((x2-x1)**2 + (y2-y1)**2)
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__":
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")