Notebook 5.4 - Capacitated Facility Network Location Problem¶
In [1]:
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)
In [1]:
import os
os.system('jupyter nbconvert --to html yourNotebook.ipynb')
Out[1]:
In [2]:
center_box = (100, 200)
n_customer = 40
n_facility = 5
facil_cap_low = 400
facil_cap_high = 1200
custm_dem_low = 30
custm_dem_high = 80
transp_unit_cost = 2
capacity_unit_cost = 100
For this problem we are going to be randomly generating quite a few values to complement our dataset of customers and facilities. These are defined in the cell above.
In [3]:
np.random.seed(2)
We are using a random number generator to populate the values for facility capacities and customer demands. If we want to ensure that we always get the same instance if we rerun the notebook, we can set a random seed (here: 2) using the command above.
Different seed values will result in differnt datasets.
In [4]:
facil_capacities = np.random.randint(low=facil_cap_low,
high=facil_cap_high,
size=n_facility)
facil_capacities
Out[4]:
The np.random.randint will create an array of random integers (here, of size n_facility), between the values defined by the low and high parameters.
In [5]:
custm_demands = np.random.randint(low=custm_dem_low,
high=custm_dem_high,
size=n_customer)
In [6]:
print(f'sum_facil ({sum(facil_capacities)}) > sum_cust ({sum(custm_demands)}) = {sum(facil_capacities) > sum(custm_demands)})')
As a sanity check, we are testing whether we have enough facility capacity to serve our customers.
In [7]:
cust_coord , _ = make_blobs(n_samples=n_customer,
centers=2,
cluster_std=20,
center_box=center_box,
random_state = 2)
facil_coord, _ = make_blobs(n_samples=n_facility,
centers=1,
cluster_std=30,
center_box=center_box,
random_state = 50)
In [8]:
custm_node_size = plot_size*10*((custm_demands/custm_dem_high)**2)
facil_node_size = plot_size*300*(((facil_capacities-facil_cap_low)/facil_cap_high)**2)
plt.scatter(cust_coord[:, 0],
cust_coord[:, 1],
s=custm_node_size,
cmap='viridis',
zorder=1500);
plt.scatter(facil_coord[:, 0],
facil_coord[:, 1],
s=facil_node_size,
cmap='viridis',
zorder=2000);
We are using nodes of varying sizes to illustrate the differences in capacities and demands. Adjust the parameters as necessary.
In [9]:
from scipy.spatial import distance
dist_matrix = distance.cdist(cust_coord, facil_coord, 'euclidean')
In [10]:
from pulp import *
prob = LpProblem('prob', LpMinimize)
In [11]:
S = [('S'+str(i+1)) for i in range(n_facility)]
D = [('C'+str(i+1)) for i in range(n_customer)]
P = {facility:(facil_capacities[j]) for j,facility in enumerate(S)}
In [12]:
transp_cost_matrix = dist_matrix * transp_unit_cost
t = {customer:{facility:(transp_cost_matrix[i][j]) for j,facility in enumerate(S)} for i,customer in enumerate(D)}
c = {facility:(facil_capacities[j]*capacity_unit_cost) for j,facility in enumerate(S)}
Q = {customer:(custm_demands[i]) for i,customer in enumerate(D)}
In [13]:
n_min = 1
n_max = n_facility
In [14]:
x = LpVariable.dicts('x', (D,S), lowBound = 0, cat = LpContinuous)
y = LpVariable.dicts('y', S, lowBound = 0, upBound = 1, cat = LpInteger)
In [15]:
prob += lpSum([lpSum([x[i][j]*t[i][j] for j in S]) for i in D]) + lpSum([c[j]*y[j] for j in S])
In [16]:
for i in D:
prob += lpSum([x[i][j] for j in S]) == Q[i]
for j in S:
prob += lpSum([x[i][j] for i in D]) <= P[j]*y[j]
prob += lpSum([y[j] for j in S]) >= n_min
prob += lpSum([y[j] for j in S]) <= n_max
In [17]:
status = prob.solve()
print(f'STATUS\n{LpStatus[status]}\n')
In [18]:
chosen_facil = []
for j,facility in enumerate(S):
if y[facility].varValue == 1:
chosen_facil.append(j)
print(f'We will be establishing a facility in {facility}')
chosen_facil
Out[18]:
Summary¶
from pulp import *
from scipy.spatial import distance
# Define sets of customers and facilities
S = [('S'+str(i+1)) for i in range(n_facility)]
D = [('C'+str(i+1)) for i in range(n_customer)]
# Populate cost, demand and capacity arrays
P = {facility:(facil_capacities[j]) for j,facility in enumerate(S)}
c = {facility:(facil_capacities[j]*capacity_unit_cost) for j,facility in enumerate(S)}
Q = {customer:(custm_demands[i]) for i,customer in enumerate(D)}
trc_matrix = dist_matrix * transp_unit_cost
t = {customer:{facility:(trc_matrix[i][j]) for j,facility in enumerate(S)} for i,customer in enumerate(D)}
# Define minimum and maximum number of open facilities
n_min = 1
n_max = n_facility
# Define decision variables
x = LpVariable.dicts('x', (D,S), lowBound = 0, cat = LpContinuous)
y = LpVariable.dicts('y', S, lowBound = 0, upBound = 1, cat = LpInteger)
# Objective function: minimize the sum of operating and transportation costs
prob += lpSum([lpSum([x[i][j]*t[i][j] for j in S]) for i in D]) + lpSum([c[j]*y[j] for j in S])
# Constraint 1: We must meet the demand of all customers
for i in D:
prob += lpSum([x[i][j] for j in S]) == Q[i]
# Constraint 2: The logical constraint - facilities supply material only if they are open
for j in S:
prob += lpSum([x[i][j] for i in D]) <= P[j]*y[j]
# Constraint 3: Define the minimum and maximum number of facilities
prob += lpSum([y[j] for j in S]) >= n_min
prob += lpSum([y[j] for j in S]) <= n_max
# Solve the model
status = prob.solve()
In [19]:
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 > 0:
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 [20]:
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');
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')
plt.scatter(facil_coord[chosen_facil, 0],
facil_coord[chosen_facil, 1],
s=plot_size*10,
cmap='viridis',
zorder=99999)
Out[20]:
