Notebook 5.2 - Single Facility Network Location Problem¶
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 = 14
plot_width = 5
plot_height = 5
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
We set here the number of customers and facilities that will be used by the data generators and later in the model.
center_box = (100, 200)
cust_coord, _ = make_blobs(n_samples=num_customer,
centers=2,
cluster_std=20,
center_box=center_box,
random_state = 2)
We are again using make_blobs to generate a random set of customers.
You might notice that we are saving predicted cluster indices (which we don't need) to a variable named _. This is a real variable name, which is commonly used whenever we want to disregard a certain value.
facil_coord, _ = make_blobs(n_samples=num_facility,
centers=1,
cluster_std=30,
center_box=center_box,
random_state = 50)
facil_coord
To solve the facility location problem, we also need to provide a set of candidate locations. We use make_blobs again for this.
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)
We are back again in using PuLP.
from scipy.spatial import distance
dist_matrix = distance.cdist(facil_coord, facil_coord, 'euclidean').round(2)
dist_matrix
We are going to use the distance.cdist() command to determine the euclidean distances between our sets of customers and facility coordinates.
In a real application, we will instead by using a pathfinder, or a travel time database to obtain realistic values that reflect the structure of the network.
The cell above calculates a distance matrix between all facilitiy locations. You might also notice that we are also using the round() function, as the original values will have a very large number of decimals, and would not render well when printed out.
dist_matrix = distance.cdist(cust_coord, facil_coord, 'euclidean')
dist_matrix.shape
With the above command we have now calculated the distance matrix between the set of customers and facilities.
The shape property is used check the size of the matrix. This appears to be 40 x 5, as expected.
S = [('S'+str(i+1)) for i in range(num_facility)]
D = [('C'+str(i+1)) for i in range(num_customer)]
With the above in-line for loop, we are generating the indices for each customer and facility location element, to be saved in D and S, respectively.
You can check out their contents int he cell below:
print(f'The S array contains the following: \n{S}\n')
print(f'While D, contains: \n{D}')
d = {customer:{facility: round(dist_matrix[i][j],2) for j,facility in enumerate(S)} for i,customer in enumerate(D)}
With another inline for loop (a nested one, in this case), we have generated a dictionary of distances between customers i and facilities j. For visualisation purposes, we have rounded again the values to two decimals.
d
x = LpVariable.dicts('x', (D,S), lowBound = 0, upBound = 1, cat = LpInteger)
y = LpVariable.dicts('y', S, lowBound = 0, upBound = 1, cat = LpInteger)
With the above we have created our sets of decision variables, both defined as Booleans (integers with values between 0 and 1).
prob += lpSum([lpSum([x[i][j]*d[i][j] for j in S]) for i in D])
The above expression defines our objective, which minimises the product of distances with the decision variables that indicate whether a facility and customer are connected.
for i in D:
prob += lpSum([x[i][j] for j in S]) == 1
Contraint 1: Ensure that only one facility serves each customer.
prob += lpSum([y[j] for j in S]) == 1
Contraint 2: Ensure that only one facility is open.
for i in D:
for j in S:
prob += x[i][j] <= y[j]
Contraint 3: A facility will not serve a customer if it has not been opened.
status = prob.solve()
print(f'STATUS\n{LpStatus[status]}\n')
for j,facility in enumerate(S):
if y[facility].varValue == 1:
chosen_facil = j
print(f'We will be opening a facility in {facility}')
With the above loop, we are iterating through all decision variables in y to establish whether a facility has been opened.
tran_cost = 0;
for j,facility in enumerate(S):
for i,customer in enumerate(D):
if x[customer][facility].varValue == 1:
tran_cost += dist_matrix[i][j]
print(f'The objective value of the problem (total transportation cost) is {value(prob.objective)}');
Summary¶
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)
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');
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');
In green color, the above graph demonstrates the location of the chosen facility.