|
15 | 15 | # by using a different backend.
|
16 | 16 | if sys.platform == "darwin": # OS X
|
17 | 17 | import matplotlib as mpl
|
18 |
| - mpl.use('Agg') |
| 18 | + |
| 19 | + mpl.use("Agg") |
19 | 20 | del mpl
|
20 | 21 |
|
| 22 | +import matplotlib.pyplot as plt |
21 | 23 | from math import sqrt, log
|
22 | 24 | from itertools import product
|
23 | 25 | from mip import Model, xsum, minimize, OptimizationStatus
|
|
35 | 37 | C = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
|
36 | 38 |
|
37 | 39 | # position of clients
|
38 |
| -pc = {1: (94, 10), 2: (57, 26), 3: (74, 44), 4: (27, 51), 5: (78, 30), 6: (23, 30), |
39 |
| - 7: (20, 72), 8: (3, 27), 9: (5, 39), 10: (51, 1)} |
| 40 | +pc = { |
| 41 | + 1: (94, 10), |
| 42 | + 2: (57, 26), |
| 43 | + 3: (74, 44), |
| 44 | + 4: (27, 51), |
| 45 | + 5: (78, 30), |
| 46 | + 6: (23, 30), |
| 47 | + 7: (20, 72), |
| 48 | + 8: (3, 27), |
| 49 | + 9: (5, 39), |
| 50 | + 10: (51, 1), |
| 51 | +} |
40 | 52 |
|
41 | 53 | # demands
|
42 | 54 | d = {1: 302, 2: 273, 3: 275, 4: 266, 5: 287, 6: 296, 7: 297, 8: 310, 9: 302, 10: 309}
|
43 | 55 |
|
44 |
| -dist = {(f, c): round(sqrt((pf[f][0] - pc[c][0]) ** 2 + (pf[f][1] - pc[c][1]) ** 2), 1) |
45 |
| - for (f, c) in product(F, C) } |
| 56 | +# plotting possible plant locations |
| 57 | +for i, p in pf.items(): |
| 58 | + plt.scatter((p[0]), (p[1]), marker="^", color="purple", s=50) |
| 59 | + plt.text((p[0]), (p[1]), "$f_%d$" % i) |
| 60 | + |
| 61 | +# plotting location of clients |
| 62 | +for i, p in pc.items(): |
| 63 | + plt.scatter((p[0]), (p[1]), marker="o", color="black", s=15) |
| 64 | + plt.text((p[0]), (p[1]), "$c_{%d}$" % i) |
| 65 | + |
| 66 | +plt.text((20), (78), "Region 1") |
| 67 | +plt.text((70), (78), "Region 2") |
| 68 | +plt.plot((50, 50), (0, 80)) |
| 69 | + |
| 70 | +dist = { |
| 71 | + (f, c): round(sqrt((pf[f][0] - pc[c][0]) ** 2 + (pf[f][1] - pc[c][1]) ** 2), 1) |
| 72 | + for (f, c) in product(F, C) |
| 73 | +} |
46 | 74 |
|
47 | 75 | m = Model()
|
48 | 76 |
|
|
67 | 95 | D = 6 # nr. of discretization points, increase for more precision
|
68 | 96 | v = [c[f] * (v / (D - 1)) for v in range(D)] # points
|
69 | 97 | # non-linear function values for points in v
|
70 |
| - vn = [0 if k == 0 else 1520 * log(v[k]) for k in range(D)] |
| 98 | + vn = [0 if k == 0 else 1520 * log(v[k]) for k in range(D)] |
71 | 99 | # w variables
|
72 | 100 | w = [m.add_var() for v in range(D)]
|
73 | 101 | m += xsum(w) == 1 # convexification
|
|
83 | 111 |
|
84 | 112 | # objective function
|
85 | 113 | m.objective = minimize(
|
86 |
| - xsum(dist[i, j] * x[i, j] for (i, j) in product(F, C)) + xsum(y[i] for i in F) ) |
| 114 | + xsum(dist[i, j] * x[i, j] for (i, j) in product(F, C)) + xsum(y[i] for i in F) |
| 115 | +) |
87 | 116 |
|
88 | 117 | m.optimize()
|
89 | 118 |
|
| 119 | +plt.savefig("location.pdf") |
| 120 | + |
90 | 121 | if m.num_solutions:
|
91 | 122 | print("Solution with cost {} found.".format(m.objective_value))
|
92 | 123 | print("Facilities capacities: {} ".format([z[f].x for f in F]))
|
93 | 124 | print("Facilities cost: {}".format([y[f].x for f in F]))
|
94 |
| - |
| 125 | + |
| 126 | + # plotting allocations |
| 127 | + for i, j in [(i, j) for (i, j) in product(F, C) if x[(i, j)].x >= 1e-6]: |
| 128 | + plt.plot( |
| 129 | + (pf[i][0], pc[j][0]), (pf[i][1], pc[j][1]), linestyle="--", color="darkgray" |
| 130 | + ) |
| 131 | + |
| 132 | + plt.savefig("location-sol.pdf") |
| 133 | + |
95 | 134 | # sanity checks
|
96 | 135 | opt = 99733.94905406
|
97 | 136 | if m.status == OptimizationStatus.OPTIMAL:
|
|
0 commit comments