!pip install pulp
Collecting pulp Downloading PuLP-2.7.0-py3-none-any.whl (14.3 MB) |████████████████████████████████| 14.3 MB 2.5 MB/s eta 0:00:01 Installing collected packages: pulp Successfully installed pulp-2.7.0
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
array([[110.99138668, 82.98394341], [153.25029374, 148.67412191], [126.04608814, 154.9163425 ], [130.83231184, 78.87089589], [191.81854821, 108.50634614]])
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
array([[ 0. , 78.11, 73.49, 20.26, 84.76], [78.11, 0. , 27.91, 73.31, 55.69], [73.49, 27.91, 0. , 76.2 , 80.5 ], [20.26, 73.31, 76.2 , 0. , 67.81], [84.76, 55.69, 80.5 , 67.81, 0. ]])
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
(40, 5)
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}')
The S array contains the following: ['S1', 'S2', 'S3', 'S4', 'S5'] While D, contains: ['C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7', 'C8', 'C9', 'C10', 'C11', 'C12', 'C13', 'C14', 'C15', 'C16', 'C17', 'C18', 'C19', 'C20', 'C21', 'C22', 'C23', 'C24', 'C25', 'C26', 'C27', 'C28', 'C29', 'C30', 'C31', 'C32', 'C33', 'C34', 'C35', 'C36', 'C37', 'C38', 'C39', 'C40']
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
{'C1': {'S1': 104.73, 'S2': 26.92, 'S3': 47.4, 'S4': 98.4, 'S5': 63.35}, 'C2': {'S1': 83.04, 'S2': 25.2, 'S3': 52.49, 'S4': 72.12, 'S5': 32.17}, 'C3': {'S1': 53.66, 'S2': 42.74, 'S3': 58.85, 'S4': 39.57, 'S5': 32.6}, 'C4': {'S1': 43.83, 'S2': 34.55, 'S3': 38.59, 'S4': 39.68, 'S5': 55.31}, 'C5': {'S1': 69.77, 'S2': 15.1, 'S3': 14.4, 'S4': 68.75, 'S5': 66.12}, 'C6': {'S1': 78.67, 'S2': 1.4, 'S3': 29.31, 'S4': 73.52, 'S5': 54.57}, 'C7': {'S1': 53.82, 'S2': 24.71, 'S3': 33.55, 'S4': 48.82, 'S5': 51.69}, 'C8': {'S1': 33.55, 'S2': 69.01, 'S3': 76.91, 'S4': 13.67, 'S5': 55.18}, 'C9': {'S1': 90.97, 'S2': 33.25, 'S3': 61.02, 'S4': 78.86, 'S5': 29.29}, 'C10': {'S1': 4.28, 'S2': 77.66, 'S3': 71.54, 'S4': 24.11, 'S5': 87.11}, 'C11': {'S1': 11.54, 'S2': 71.27, 'S3': 70.6, 'S4': 10.06, 'S5': 73.44}, 'C12': {'S1': 44.62, 'S2': 36.53, 'S3': 45.04, 'S4': 36.85, 'S5': 48.47}, 'C13': {'S1': 61.2, 'S2': 48.4, 'S3': 67.4, 'S4': 44.97, 'S5': 23.63}, 'C14': {'S1': 83.34, 'S2': 16.31, 'S3': 16.77, 'S4': 82.26, 'S5': 71.99}, 'C15': {'S1': 29.81, 'S2': 53.39, 'S3': 58.06, 'S4': 19.92, 'S5': 56.01}, 'C16': {'S1': 17.06, 'S2': 62.17, 'S3': 56.45, 'S4': 25.92, 'S5': 77.17}, 'C17': {'S1': 37.74, 'S2': 66.01, 'S3': 75.58, 'S4': 18.33, 'S5': 50.16}, 'C18': {'S1': 94.83, 'S2': 33.39, 'S3': 21.46, 'S4': 96.38, 'S5': 88.95}, 'C19': {'S1': 26.52, 'S2': 52.08, 'S3': 52.21, 'S4': 24.0, 'S5': 63.44}, 'C20': {'S1': 44.06, 'S2': 58.02, 'S3': 70.22, 'S4': 26.34, 'S5': 41.47}, 'C21': {'S1': 100.93, 'S2': 22.91, 'S3': 40.18, 'S4': 96.01, 'S5': 66.83}, 'C22': {'S1': 71.39, 'S2': 7.0, 'S3': 28.26, 'S4': 66.31, 'S5': 52.66}, 'C23': {'S1': 77.33, 'S2': 12.07, 'S3': 16.27, 'S4': 75.7, 'S5': 67.06}, 'C24': {'S1': 61.58, 'S2': 19.02, 'S3': 35.42, 'S4': 54.77, 'S5': 46.02}, 'C25': {'S1': 69.41, 'S2': 8.81, 'S3': 24.98, 'S4': 65.32, 'S5': 55.61}, 'C26': {'S1': 68.06, 'S2': 15.25, 'S3': 15.65, 'S4': 66.84, 'S5': 64.87}, 'C27': {'S1': 20.28, 'S2': 73.34, 'S3': 76.23, 'S4': 0.03, 'S5': 67.81}, 'C28': {'S1': 22.32, 'S2': 56.23, 'S3': 55.45, 'S4': 21.14, 'S5': 66.4}, 'C29': {'S1': 42.99, 'S2': 70.99, 'S3': 82.02, 'S4': 22.86, 'S5': 49.05}, 'C30': {'S1': 59.44, 'S2': 51.94, 'S3': 27.63, 'S4': 69.22, 'S5': 95.97}, 'C31': {'S1': 67.96, 'S2': 10.25, 'S3': 24.87, 'S4': 63.98, 'S5': 55.63}, 'C32': {'S1': 79.77, 'S2': 1.66, 'S3': 28.55, 'S4': 74.91, 'S5': 56.06}, 'C33': {'S1': 54.44, 'S2': 48.39, 'S3': 26.56, 'S4': 63.32, 'S5': 90.27}, 'C34': {'S1': 71.76, 'S2': 7.04, 'S3': 23.54, 'S4': 68.02, 'S5': 57.45}, 'C35': {'S1': 83.32, 'S2': 14.56, 'S3': 42.45, 'S4': 75.02, 'S5': 43.6}, 'C36': {'S1': 42.73, 'S2': 62.97, 'S3': 74.5, 'S4': 23.81, 'S5': 44.47}, 'C37': {'S1': 92.89, 'S2': 15.5, 'S3': 39.55, 'S4': 86.59, 'S5': 56.29}, 'C38': {'S1': 43.98, 'S2': 35.94, 'S3': 43.01, 'S4': 37.4, 'S5': 50.79}, 'C39': {'S1': 61.56, 'S2': 39.32, 'S3': 15.83, 'S4': 67.67, 'S5': 85.58}, 'C40': {'S1': 27.14, 'S2': 72.97, 'S3': 58.96, 'S4': 44.34, 'S5': 96.98}}
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')
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/6e769006fd6547e694d78c7a4a868341-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/kk/14m2mg8s0w9c_tln749srnrw0000gp/T/6e769006fd6547e694d78c7a4a868341-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 1494.27 - 0.00 seconds Cgl0005I 41 SOS with 205 members Cgl0004I processed model has 241 rows, 205 columns (205 integer (205 of which binary)) and 605 elements Cutoff increment increased from 1e-05 to 0.00999 Cbc0038I Initial state - 0 integers unsatisfied sum - 0 Cbc0038I Solution found of 1494.27 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 1494.27 - took 0.00 seconds Cbc0012I Integer solution of 1494.27 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds) Cbc0001I Search completed - best objective 1494.27, 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 1494.27 to 1494.27 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: 1494.27000000 Enumerated nodes: 0 Total iterations: 0 Time (CPU seconds): 0.01 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
for j,facility in enumerate(S):
if y[facility].varValue == 1:
chosen_facil = j
print(f'We will be opening a facility in {facility}')
We will be opening a facility in S2
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)}');
The objective value of the problem (total transportation cost) is 1494.27
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.