Transport Analytics Training Series - Last Revision: October 2022

Capacitated Facility Network 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   = 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 [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]:
array([568, 927, 893, 984, 934])

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)})')
sum_facil (4306) > sum_cust (2222) = True)

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')

image.png

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')
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/d70614f588c1491daeeba5c1ab29df92-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/kk/14m2mg8s0w9c_tln749srnrw0000gp/T/d70614f588c1491daeeba5c1ab29df92-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 52 COLUMNS
At line 683 RHS
At line 731 BOUNDS
At line 737 ENDATA
Problem MODEL has 47 rows, 205 columns and 415 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 306819 - 0.00 seconds
Cgl0004I processed model has 46 rows, 205 columns (5 integer (5 of which binary)) and 410 elements
Cbc0038I Initial state - 4 integers unsatisfied sum - 1.24033
Cbc0038I Pass   1: suminf.    0.33298 (1) obj. 357493 iterations 22
Cbc0038I Solution found of 419793
Cbc0038I Relaxing continuous gives 395296
Cbc0038I Before mini branch and bound, 1 integers at bound fixed and 175 continuous
Cbc0038I Mini branch and bound did not improve solution (0.02 seconds)
Cbc0038I Round again with cutoff of 386448
Cbc0038I Pass   2: suminf.    0.33298 (1) obj. 357493 iterations 0
Cbc0038I Pass   3: suminf.    0.09473 (1) obj. 386448 iterations 25
Cbc0038I Pass   4: suminf.    0.36687 (1) obj. 369079 iterations 17
Cbc0038I Pass   5: suminf.    0.09473 (1) obj. 386448 iterations 19
Cbc0038I Pass   6: suminf.    0.36687 (1) obj. 367604 iterations 18
Cbc0038I Pass   7: suminf.    0.29628 (1) obj. 386448 iterations 51
Cbc0038I Pass   8: suminf.    0.40854 (1) obj. 344413 iterations 24
Cbc0038I Pass   9: suminf.    0.00000 (0) obj. 386448 iterations 8
Cbc0038I Solution found of 386448
Cbc0038I Relaxing continuous gives 373806
Cbc0038I Before mini branch and bound, 1 integers at bound fixed and 114 continuous
Cbc0038I Full problem 46 rows 205 columns, reduced to 22 rows 75 columns
Cbc0038I Mini branch and bound did not improve solution (0.03 seconds)
Cbc0038I Freeing continuous variables gives a solution of 373806
Cbc0038I Round again with cutoff of 360409
Cbc0038I Pass  10: suminf.    0.33298 (1) obj. 357493 iterations 0
Cbc0038I Pass  11: suminf.    0.37352 (1) obj. 360409 iterations 26
Cbc0038I Pass  12: suminf.    0.56491 (2) obj. 360409 iterations 23
Cbc0038I Pass  13: suminf.    0.00000 (0) obj. 360409 iterations 10
Cbc0038I Solution found of 360409
Cbc0038I Relaxing continuous gives 359983
Cbc0038I Before mini branch and bound, 1 integers at bound fixed and 150 continuous
Cbc0038I Full problem 46 rows 205 columns, reduced to 16 rows 47 columns
Cbc0038I Mini branch and bound did not improve solution (0.03 seconds)
Cbc0038I Freeing continuous variables gives a solution of 359983
Cbc0038I Round again with cutoff of 344034
Cbc0038I Pass  14: suminf.    0.33298 (1) obj. 344034 iterations 8
Cbc0038I Pass  15: suminf.    0.45115 (1) obj. 344034 iterations 19
Cbc0038I Pass  16: suminf.    0.66834 (2) obj. 344034 iterations 16
Cbc0038I Pass  17: suminf.    0.45115 (1) obj. 344034 iterations 8
Cbc0038I Pass  18: suminf.    0.33298 (1) obj. 344034 iterations 3
Cbc0038I Pass  19: suminf.    0.40854 (1) obj. 344034 iterations 42
Cbc0038I Pass  20: suminf.    0.32175 (2) obj. 344034 iterations 20
Cbc0038I Pass  21: suminf.    0.33298 (1) obj. 344034 iterations 26
Cbc0038I Pass  22: suminf.    0.33298 (1) obj. 344034 iterations 0
Cbc0038I Pass  23: suminf.    0.66834 (2) obj. 344034 iterations 16
Cbc0038I Pass  24: suminf.    0.33298 (1) obj. 344034 iterations 4
Cbc0038I Pass  25: suminf.    0.05246 (1) obj. 344034 iterations 62
Cbc0038I Pass  26: suminf.    0.26118 (1) obj. 344034 iterations 60
Cbc0038I Pass  27: suminf.    0.16612 (1) obj. 344034 iterations 17
Cbc0038I Pass  28: suminf.    0.47192 (2) obj. 344034 iterations 18
Cbc0038I Pass  29: suminf.    0.16612 (1) obj. 344034 iterations 16
Cbc0038I Pass  30: suminf.    0.16612 (1) obj. 344034 iterations 0
Cbc0038I Pass  31: suminf.    0.37004 (2) obj. 344034 iterations 66
Cbc0038I Pass  32: suminf.    0.57337 (3) obj. 344034 iterations 27
Cbc0038I Pass  33: suminf.    0.32175 (2) obj. 344034 iterations 11
Cbc0038I Pass  34: suminf.    0.33298 (1) obj. 344034 iterations 29
Cbc0038I Pass  35: suminf.    0.33298 (1) obj. 344034 iterations 0
Cbc0038I Pass  36: suminf.    0.45115 (1) obj. 344034 iterations 16
Cbc0038I Pass  37: suminf.    0.80127 (4) obj. 344034 iterations 49
Cbc0038I Pass  38: suminf.    0.40854 (1) obj. 344034 iterations 18
Cbc0038I Pass  39: suminf.    0.32175 (2) obj. 344034 iterations 21
Cbc0038I Pass  40: suminf.    0.33298 (1) obj. 344034 iterations 28
Cbc0038I Pass  41: suminf.    0.33298 (1) obj. 344034 iterations 0
Cbc0038I Pass  42: suminf.    0.45115 (1) obj. 344034 iterations 20
Cbc0038I Pass  43: suminf.    0.95341 (3) obj. 344034 iterations 31
Cbc0038I Rounding solution of 348719 is better than previous of 359983

Cbc0038I Before mini branch and bound, 0 integers at bound fixed and 76 continuous
Cbc0038I Full problem 46 rows 205 columns, reduced to 35 rows 118 columns
Cbc0038I Mini branch and bound did not improve solution (0.04 seconds)
Cbc0038I Round again with cutoff of 332869
Cbc0038I Pass  43: suminf.    0.33304 (2) obj. 332869 iterations 11
Cbc0038I Pass  44: suminf.    0.69839 (3) obj. 332869 iterations 18
Cbc0038I Pass  45: suminf.    0.44677 (2) obj. 332869 iterations 11
Cbc0038I Pass  46: suminf.    0.33304 (2) obj. 332869 iterations 48
Cbc0038I Pass  47: suminf.    0.87923 (3) obj. 332869 iterations 12
Cbc0038I Pass  48: suminf.    0.96168 (4) obj. 332869 iterations 28
Cbc0038I Pass  49: suminf.    0.97509 (2) obj. 332869 iterations 52
Cbc0038I Pass  50: suminf.    0.30764 (2) obj. 332869 iterations 30
Cbc0038I Pass  51: suminf.    0.17749 (1) obj. 332869 iterations 50
Cbc0038I Pass  52: suminf.    0.31387 (2) obj. 332869 iterations 73
Cbc0038I Pass  53: suminf.    0.33304 (2) obj. 332869 iterations 31
Cbc0038I Pass  54: suminf.    0.33304 (2) obj. 332869 iterations 0
Cbc0038I Pass  55: suminf.    0.69839 (3) obj. 332869 iterations 18
Cbc0038I Pass  56: suminf.    0.44677 (2) obj. 332869 iterations 11
Cbc0038I Pass  57: suminf.    0.33304 (2) obj. 332869 iterations 48
Cbc0038I Pass  58: suminf.    0.58928 (2) obj. 332869 iterations 21
Cbc0038I Pass  59: suminf.    0.30764 (2) obj. 332869 iterations 15
Cbc0038I Pass  60: suminf.    0.17749 (1) obj. 332869 iterations 50
Cbc0038I Pass  61: suminf.    0.31387 (2) obj. 332869 iterations 73
Cbc0038I Pass  62: suminf.    0.97509 (2) obj. 332869 iterations 27
Cbc0038I Pass  63: suminf.    0.52534 (2) obj. 332869 iterations 45
Cbc0038I Pass  64: suminf.    0.52534 (2) obj. 332869 iterations 0
Cbc0038I Pass  65: suminf.    0.52534 (2) obj. 332869 iterations 0
Cbc0038I Pass  66: suminf.    0.96168 (4) obj. 332869 iterations 21
Cbc0038I Pass  67: suminf.    0.97509 (2) obj. 332869 iterations 45
Cbc0038I Pass  68: suminf.    0.33304 (2) obj. 332869 iterations 46
Cbc0038I Pass  69: suminf.    0.33304 (2) obj. 332869 iterations 0
Cbc0038I Pass  70: suminf.    0.69839 (3) obj. 332869 iterations 18
Cbc0038I Pass  71: suminf.    0.44677 (2) obj. 332869 iterations 7
Cbc0038I Pass  72: suminf.    0.33304 (2) obj. 332869 iterations 47
Cbc0038I Before mini branch and bound, 0 integers at bound fixed and 109 continuous
Cbc0038I Full problem 46 rows 205 columns, reduced to 22 rows 79 columns
Cbc0038I Mini branch and bound did not improve solution (0.05 seconds)
Cbc0038I Round again with cutoff of 322449
Cbc0038I Pass  72: suminf.    0.33858 (2) obj. 322449 iterations 3
Cbc0038I Pass  73: suminf.    0.81508 (3) obj. 322449 iterations 17
Cbc0038I Pass  74: suminf.    0.56345 (2) obj. 322449 iterations 11
Cbc0038I Pass  75: suminf.    0.33858 (2) obj. 322449 iterations 29
Cbc0038I Pass  76: suminf.    0.33858 (2) obj. 322449 iterations 2
Cbc0038I Pass  77: suminf.    0.56345 (2) obj. 322449 iterations 24
Cbc0038I Pass  78: suminf.    0.33STATUS
Optimal

858 (2) obj. 322449 iterations 34
Cbc0038I Pass  79: suminf.    0.33858 (2) obj. 322449 iterations 0
Cbc0038I Pass  80: suminf.    0.81508 (3) obj. 322449 iterations 18
Cbc0038I Pass  81: suminf.    0.56345 (2) obj. 322449 iterations 11
Cbc0038I Pass  82: suminf.    0.87670 (2) obj. 322449 iterations 16
Cbc0038I Pass  83: suminf.    0.74107 (3) obj. 322449 iterations 30
Cbc0038I Pass  84: suminf.    0.57075 (3) obj. 322449 iterations 10
Cbc0038I Pass  85: suminf.    0.72543 (2) obj. 322449 iterations 39
Cbc0038I Pass  86: suminf.    0.41594 (2) obj. 322449 iterations 24
Cbc0038I Pass  87: suminf.    0.78529 (4) obj. 322449 iterations 17
Cbc0038I Pass  88: suminf.    0.74107 (3) obj. 322449 iterations 36
Cbc0038I Pass  89: suminf.    0.41594 (2) obj. 322449 iterations 35
Cbc0038I Pass  90: suminf.    0.41594 (2) obj. 322449 iterations 0
Cbc0038I Pass  91: suminf.    0.92320 (4) obj. 322449 iterations 20
Cbc0038I Pass  92: suminf.    0.56345 (2) obj. 322449 iterations 29
Cbc0038I Pass  93: suminf.    0.33858 (2) obj. 322449 iterations 30
Cbc0038I Pass  94: suminf.    0.33858 (2) obj. 322449 iterations 0
Cbc0038I Pass  95: suminf.    0.81508 (3) obj. 322449 iterations 16
Cbc0038I Pass  96: suminf.    0.56345 (2) obj. 322449 iterations 11
Cbc0038I Pass  97: suminf.    0.81508 (3) obj. 322449 iterations 5
Cbc0038I Pass  98: suminf.    0.56345 (2) obj. 322449 iterations 7
Cbc0038I Pass  99: suminf.    0.33858 (2) obj. 322449 iterations 32
Cbc0038I Pass 100: suminf.    0.33858 (2) obj. 322449 iterations 0
Cbc0038I Pass 101: suminf.    0.81508 (3) obj. 322449 iterations 16
Cbc0038I No solution found this major pass
Cbc0038I Before mini branch and bound, 0 integers at bound fixed and 116 continuous
Cbc0038I Mini branch and bound did not improve solution (0.06 seconds)
Cbc0038I After 0.06 seconds - Feasibility pump exiting with objective of 348719 - took 0.04 seconds
Cbc0012I Integer solution of 348718.67 found by feasibility pump after 0 iterations and 0 nodes (0.06 seconds)
Cbc0031I 45 added rows had average density of 17.844444
Cbc0013I At root node, 45 cuts changed objective from 306819.04 to 336754.02 in 13 passes
Cbc0014I Cut generator 0 (Probing) - 69 row cuts average 2.0 elements, 0 column cuts (0 active)  in 0.001 seconds - new frequency is 1
Cbc0014I Cut generator 1 (Gomory) - 15 row cuts average 115.3 elements, 0 column cuts (0 active)  in 0.001 seconds - new frequency is 1
Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.001 seconds - new frequency is -100
Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.000 seconds - new frequency is -100
Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 40 row cuts average 32.6 elements, 0 column cuts (0 active)  in 0.001 seconds - new frequency is 1
Cbc0014I Cut generator 5 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.002 seconds - new frequency is -100
Cbc0014I Cut generator 6 (TwoMirCuts) - 46 row cuts average 47.6 elements, 0 column cuts (0 active)  in 0.001 seconds - new frequency is 1
Cbc0010I After 0 nodes, 1 on tree, 348718.67 best solution, best possible 336754.02 (0.08 seconds)
Cbc0001I Search completed - best objective 348718.6698739209, took 706 iterations and 6 nodes (0.10 seconds)
Cbc0032I Strong branching done 24 times (579 iterations), fathomed 2 nodes and fixed 0 variables
Cbc0035I Maximum depth 2, 1 variables fixed on reduced cost
Cuts at root node changed objective from 306819 to 336754
Probing was tried 36 times and created 117 cuts of which 0 were active after adding rounds of cuts (0.002 seconds)
Gomory was tried 36 times and created 15 cuts of which 0 were active after adding rounds of cuts (0.002 seconds)
Knapsack was tried 13 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.001 seconds)
Clique was tried 13 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 36 times and created 139 cuts of which 0 were active after adding rounds of cuts (0.003 seconds)
FlowCover was tried 13 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.002 seconds)
TwoMirCuts was tried 36 times and created 114 cuts of which 0 were active after adding rounds of cuts (0.003 seconds)
ZeroHalf was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                348718.66987392
Enumerated nodes:               6
Total iterations:               706
Time (CPU seconds):             0.08
Time (Wallclock seconds):       0.10

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.08   (Wallclock seconds):       0.11

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
We will be establishing a facility in S1
We will be establishing a facility in S2
We will be establishing a facility in S3
Out[18]:
[0, 1, 2]

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]:
<matplotlib.collections.PathCollection at 0x16f734ac0>