Transport Analytics Training Series - Last Revision: October 2022

Multiple Facility Location Problems¶

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   = 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)
In [2]:
num_customer = 40
num_facility = 5
In [3]:
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)
In [4]:
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');
In [5]:
from pulp import *

prob = LpProblem('prob', LpMinimize)
In [6]:
from scipy.spatial import distance

dist_matrix = distance.cdist(cust_coord, facil_coord, 'euclidean')
In [7]:
S = [('S'+str(i+1)) for i in range(num_facility)]
D = [('C'+str(i+1)) for i in range(num_customer)]
In [8]:
n = 2

We are using the above parameter variable to set the number of facilities to be opened.

In [9]:
d = {customer:{facility:dist_matrix[i][j] for j,facility in enumerate(S)} for i,customer in enumerate(D)}
In [10]:
x = LpVariable.dicts('x', (D,S), lowBound = 0, upBound = 1, cat = LpInteger)
y = LpVariable.dicts('y', S, lowBound = 0, upBound = 1, cat = LpInteger)
In [11]:
prob += lpSum([lpSum([x[i][j]*d[i][j] for j in S]) for i in D])
In [12]:
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.

In [13]:
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
In [14]:
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
Out[14]:
[1, 3]

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)

# 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()
In [15]:
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)
Out[15]:
<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?

In [16]:
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.

In [17]:
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.