import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
from sklearn.datasets import make_blobs
%matplotlib inline
plot_size = 30
plot_width = 19
plot_height = 17
params = {'legend.fontsize': 'large',
'figure.figsize': (plot_width,plot_height),
'axes.labelsize': plot_size,
'axes.titlesize': plot_size,
'xtick.labelsize': plot_size*0.75,
'ytick.labelsize': plot_size*0.75,
'axes.titlepad': 25}
plt.rcParams.update(params)
num_customer = 40
num_facility = 5
center_box = (100, 200)
cust_coord, _ = make_blobs(n_samples=num_customer,
centers=2,
cluster_std=20,
center_box=center_box,
random_state = 2)
facil_coord, _ = make_blobs(n_samples=num_facility,
centers=1,
cluster_std=30,
center_box=center_box,
random_state = 50)
plt.scatter(cust_coord[:, 0],
cust_coord[:, 1],
s=plot_size*2,
cmap='viridis');
plt.scatter(facil_coord[:, 0],
facil_coord[:, 1],
s=plot_size*5,
cmap='viridis');
from pulp import *
prob = LpProblem('prob', LpMinimize)
from scipy.spatial import distance
dist_matrix = distance.cdist(cust_coord, facil_coord, 'euclidean')
S = [('S'+str(i+1)) for i in range(num_facility)]
D = [('C'+str(i+1)) for i in range(num_customer)]
n = 2
We are using the above parameter variable to set the number of facilities to be opened.
d = {customer:{facility:dist_matrix[i][j] for j,facility in enumerate(S)} for i,customer in enumerate(D)}
x = LpVariable.dicts('x', (D,S), lowBound = 0, upBound = 1, cat = LpInteger)
y = LpVariable.dicts('y', S, lowBound = 0, upBound = 1, cat = LpInteger)
prob += lpSum([lpSum([x[i][j]*d[i][j] for j in S]) for i in D])
for i in D:
prob += lpSum([x[i][j] for j in S]) == 1
prob += lpSum([y[j] for j in S]) == n
for i in D:
for j in S:
prob += x[i][j] <= y[j]
Note that in the 2nd constrained we use n to define the number of facilities to be opened.
status = prob.solve()
print(f'STATUS\n{LpStatus[status]}\n')
Welcome to the CBC MILP Solver Version: 2.10.3 Build Date: Dec 15 2019 command line - /Users/pa01/opt/anaconda3/lib/python3.9/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/kk/14m2mg8s0w9c_tln749srnrw0000gp/T/4108bd386fb04e63a1def741dac2ae51-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/kk/14m2mg8s0w9c_tln749srnrw0000gp/T/4108bd386fb04e63a1def741dac2ae51-pulp.sol (default strategy 1) At line 2 NAME MODEL At line 3 ROWS At line 246 COLUMNS At line 1462 RHS At line 1704 BOUNDS At line 1910 ENDATA Problem MODEL has 241 rows, 205 columns and 605 elements Coin0008I MODEL read with 0 errors Option for timeMode changed from cpu to elapsed Continuous objective value is 916.066 - 0.00 seconds Cgl0004I processed model has 241 rows, 205 columns (205 integer (205 of which binary)) and 605 elements Cbc0038I Initial state - 0 integers unsatisfied sum - 0 Cbc0038I Solution found of 916.066 Cbc0038I Before mini branch and bound, 205 integers at bound fixed and 0 continuous Cbc0038I Mini branch and bound did not improve solution (0.01 seconds) Cbc0038I After 0.01 seconds - Feasibility pump exiting with objective of 916.066 - took 0.00 seconds Cbc0012I Integer solution of 916.06615 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds) Cbc0001I Search completed - best objective 916.0661493823824, took 0 iterations and 0 nodes (0.01 seconds) Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost Cuts at root node changed objective from 916.066 to 916.066 Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Result - Optimal solution found Objective value: 916.06614938 Enumerated nodes: 0 Total iterations: 0 Time (CPU seconds): 0.00 Time (Wallclock seconds): 0.01 Option for printingOptions changed from normal to all Total time (CPU seconds): 0.01 (Wallclock seconds): 0.01 STATUS Optimal
chosen_facil = []
for j,facility in enumerate(S):
if y[facility].varValue == 1:
chosen_facil.append(j)
print(f'We will be opening facility {facility}')
chosen_facil
We will be opening facility S2 We will be opening facility S4
[1, 3]
from pulp import *
from scipy.spatial import distance
# Calculate distance matrix
dist_matrix = distance.cdist(cust_coord, facil_coord, 'euclidean')
# Initialise problem
prob = LpProblem('prob', LpMinimize)
# Define sets of customers and facilities
n_customers = len(cust_coord)
n_facil = len(facil_coord)
S = [('S'+str(i+1)) for i in range(num_facility)]
D = [('C'+str(i+1)) for i in range(num_customer)]
# Prepare distance matrix parameters
d = {customer:{facility:dist_matrix[i][j] for j,facility in enumerate(S)} for i,customer in enumerate(D)}
# Define decision variables
x = LpVariable.dicts('x', (D,S), lowBound = 0, upBound = 1, cat = LpInteger)
y = LpVariable.dicts('y', S, lowBound = 0, upBound = 1, cat = LpInteger)
# Objective function
prob += lpSum([lpSum([x[i][j]*d[i][j] for j in S]) for i in D])
# Constraint 1: Each customer will be served by 1 facility
prob += lpSum([y[j] for j in S]) == 1
# Constraint 2: Only 1 facility will be established
prob += lpSum([y[j] for j in S]) == 1
# Constraint 3: Only established facilities can serve customers
for i in D:
for j in S:
prob += x[i][j] <= y[j]
# Solve the model
status = prob.solve()
plt.scatter(cust_coord[:, 0],
cust_coord[:, 1],
s=plot_size*2,
cmap='viridis',
zorder = 10000);
plt.scatter(facil_coord[:, 0],
facil_coord[:, 1],
s=plot_size*5,
cmap='viridis');
plt.scatter(facil_coord[chosen_facil, 0],
facil_coord[chosen_facil, 1],
s=plot_size*10,
cmap='viridis',
zorder=99999)
<matplotlib.collections.PathCollection at 0x13e46fd90>
The above graph plots correctly the facilities to be opened, but how do we know which one serves any customer?
vect_orig = []
vect_dest = []
u = v = np.zeros((len(vect_orig),len(vect_orig)))
for j,facility in enumerate(S):
for i,customer in enumerate(D):
if x[customer][facility].varValue == 1:
vect_orig.append(facil_coord[j])
vect_dest.append(cust_coord[i])
vect_orig = np.array(vect_orig)
vect_dest = np.array(vect_dest)
In the above code we are iterating through our sets of facilities and customers.
Wherever we find a pair that is connected, we will append the origin and destination coordinates of a link between them to vect_orig
and vect_dest
, respectively.
plt.scatter(cust_coord[:, 0],
cust_coord[:, 1],
s=plot_size*2,
cmap='viridis',
zorder = 10000);
plt.scatter(facil_coord[:, 0],
facil_coord[:, 1],
s=plot_size*5,
cmap='viridis');
plt.scatter(facil_coord[chosen_facil, 0],
facil_coord[chosen_facil, 1],
s=plot_size*10,
cmap='viridis',
zorder=99999)
for i in range(len(vect_orig)):
plt.arrow(vect_orig[i][0],
vect_orig[i][1],
vect_dest[i][0]-vect_orig[i][0],
vect_dest[i][1]-vect_orig[i][1],
color='black')
We have now replotted the graph, but with the addition of arrows between facilities and the customers that they serve.
The eagled-eyed among you might have notices that we are now also using the zorder
parameter, which allows us to define the order with which the elements of the graph are served.