Transport Analytics Training Series - Last Revision: October 2022
Transport Analytics Training Series - Last Revision: October 2022
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)
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.
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.
facil_capacities = np.random.randint(low=facil_cap_low,
high=facil_cap_high,
size=n_facility)
facil_capacities
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.
custm_demands = np.random.randint(low=custm_dem_low,
high=custm_dem_high,
size=n_customer)
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.
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)
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.
from scipy.spatial import distance
dist_matrix = distance.cdist(cust_coord, facil_coord, 'euclidean')
from pulp import *
prob = LpProblem('prob', LpMinimize)
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)}
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)}
n_min = 1
n_max = n_facility
x = LpVariable.dicts('x', (D,S), lowBound = 0, cat = LpContinuous)
y = LpVariable.dicts('y', S, lowBound = 0, upBound = 1, cat = LpInteger)
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])
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
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
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
[0, 1, 2]
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()
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)
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)
<matplotlib.collections.PathCollection at 0x16f734ac0>
