In this notebook we examine the problem of last mile delivery routing with multiple depots and a heterogeneous vehicle fleet. In this problem, each customer has ordered a specific quantity of each product, and each depot has a certain quantity of each product type available. There is a single, heterogeneous fleet of trucks available for shipping, and each truck is based at exactly one depot (meaning its route always begins and ends at its base depot). Each truck may have a different volume capacity and cost per mile. We also note that since not all products may be available at the same depot, customer orders may need to be split and fulfilled by multiple depots. The objective is to simultaneously determine 1) the products to be shipped from each depot to each customer, 2) how to assign trucks to customers, and 3) how to route each truck to its customers, all in a way that achieves lowest total delivery cost possible.
This notebook will demonstrate how to model and solve the problem using integer programming (IP). We will use the PuLP library to build and solve the IP model, and VeRoViz to visualize the truck routes.
# Import the libraries needed.
from pulp import *
import pandas as pd
import veroviz as vrv
import time
In this scenario, we will assume a furniture company has two depots in the Fredericksburg, VA area with eight customer orders to be delivered. The scenario location and transportation data was prepared with the VeRoViz Sketch Tool. To save time, we have pre-computed the distance matrix and saved it to a pickle file. However, the code to build this matrix using VeRoViz is shown commented out.
# Units of each product available at each depot
depot_product_availabilities = {'D1':{'P1':3, 'P2': 5, 'P3': 3, 'P4':0},
'D2':{'P1':2, 'P2': 4, 'P3': 0, 'P4':4}}
# The volume (cubic feet) required per unit of each product
product_volumes = {'P1':12, 'P2':8, 'P3':9, 'P4':16}
# Truck data (capacity is cubic feet available).
trucks = {'T1':{'base':'D1', 'capacity':50, 'cost_per_mile':0.47},
'T2':{'base':'D1', 'capacity':40, 'cost_per_mile':0.41},
'T3':{'base':'D2', 'capacity':30, 'cost_per_mile':0.35},
'T4':{'base':'D2', 'capacity':90, 'cost_per_mile':0.55}}
# Number of units of each product ordered by each customer.
customer_orders = {'C1':{'P1':0, 'P2':1, 'P3':1, 'P4':0},
'C2':{'P1':0, 'P2':0, 'P3':0, 'P4':1},
'C3':{'P1':1, 'P2':0, 'P3':0, 'P4':0},
'C4':{'P1':0, 'P2':0, 'P3':0, 'P4':2},
'C5':{'P1':2, 'P2':0, 'P3':0, 'P4':0},
'C6':{'P1':0, 'P2':2, 'P3':1, 'P4':0},
'C7':{'P1':0, 'P2':1, 'P3':1, 'P4':0},
'C8':{'P1':1, 'P2':1, 'P3':0, 'P4':0}}
# Labels for trucks, depots, customers, and products.
truck_labels = sorted(trucks.keys())
depot_labels = sorted(depot_product_availabilities.keys())
customer_labels = sorted(customer_orders.keys())
product_labels = sorted(product_volumes.keys())
# These is the scenario node (i.e., locations of depots and customers) data, initially generated by the Veroviz sketch tool.
# The leafletIconType was then manually set to 'home' for depot nodes and 'star' for customer nodes.
nodesArray = [
{'id': 0, 'lat': 38.3064174, 'lon': -77.5188606, 'altMeters': 0.0, 'nodeName': 'D1', 'nodeType': 'Depot', 'popupText': 'D1', 'leafletIconPrefix': 'glyphicon', 'leafletIconType': 'home', 'leafletColor': 'blue', 'leafletIconText': '0', 'cesiumIconType': 'pin', 'cesiumColor': 'blue', 'cesiumIconText': '0', 'elevMeters': None},
{'id': 1, 'lat': 38.1876887, 'lon': -77.619984, 'altMeters': 0.0, 'nodeName': 'D2', 'nodeType': 'Depot', 'popupText': 'D2', 'leafletIconPrefix': 'glyphicon', 'leafletIconType': 'home', 'leafletColor': 'blue', 'leafletIconText': '1', 'cesiumIconType': 'pin', 'cesiumColor': 'blue', 'cesiumIconText': '1', 'elevMeters': None},
{'id': 2, 'lat': 38.2510008, 'lon': -77.5844244, 'altMeters': 0.0, 'nodeName': 'C1', 'nodeType': 'Customer', 'popupText': 'C1', 'leafletIconPrefix': 'glyphicon', 'leafletIconType': 'star', 'leafletColor': 'orange', 'leafletIconText': '2', 'cesiumIconType': 'pin', 'cesiumColor': 'orange', 'cesiumIconText': '2', 'elevMeters': None},
{'id': 3, 'lat': 38.3494599, 'lon': -77.6573983, 'altMeters': 0.0, 'nodeName': 'C2', 'nodeType': 'Customer', 'popupText': 'C2', 'leafletIconPrefix': 'glyphicon', 'leafletIconType': 'star', 'leafletColor': 'orange', 'leafletIconText': '3', 'cesiumIconType': 'pin', 'cesiumColor': 'orange', 'cesiumIconText': '3', 'elevMeters': None},
{'id': 4, 'lat': 38.3706877, 'lon': -77.5007117, 'altMeters': 0.0, 'nodeName': 'C3', 'nodeType': 'Customer', 'popupText': 'C3', 'leafletIconPrefix': 'glyphicon', 'leafletIconType': 'star', 'leafletColor': 'orange', 'leafletIconText': '4', 'cesiumIconType': 'pin', 'cesiumColor': 'orange', 'cesiumIconText': '4', 'elevMeters': None},
{'id': 5, 'lat': 38.3199284, 'lon': -77.4104575, 'altMeters': 0.0, 'nodeName': 'C4', 'nodeType': 'Customer', 'popupText': 'C4', 'leafletIconPrefix': 'glyphicon', 'leafletIconType': 'star', 'leafletColor': 'orange', 'leafletIconText': '5', 'cesiumIconType': 'pin', 'cesiumColor': 'orange', 'cesiumIconText': '5', 'elevMeters': None},
{'id': 6, 'lat': 38.1691186, 'lon': -77.8724912, 'altMeters': 0.0, 'nodeName': 'C5', 'nodeType': 'Customer', 'popupText': 'C5', 'leafletIconPrefix': 'glyphicon', 'leafletIconType': 'star', 'leafletColor': 'orange', 'leafletIconText': '6', 'cesiumIconType': 'pin', 'cesiumColor': 'orange', 'cesiumIconText': '6', 'elevMeters': None},
{'id': 7, 'lat': 38.5123587, 'lon': -78.0333479, 'altMeters': 0.0, 'nodeName': 'C6', 'nodeType': 'Customer', 'popupText': 'C6', 'leafletIconPrefix': 'glyphicon', 'leafletIconType': 'star', 'leafletColor': 'orange', 'leafletIconText': '7', 'cesiumIconType': 'pin', 'cesiumColor': 'orange', 'cesiumIconText': '7', 'elevMeters': None},
{'id': 8, 'lat': 38.4804832, 'lon': -77.4413181, 'altMeters': 0.0, 'nodeName': 'C7', 'nodeType': 'Customer', 'popupText': 'C7', 'leafletIconPrefix': 'glyphicon', 'leafletIconType': 'star', 'leafletColor': 'orange', 'leafletIconText': '8', 'cesiumIconType': 'pin', 'cesiumColor': 'orange', 'cesiumIconText': '8', 'elevMeters': None},
{'id': 9, 'lat': 38.5125157, 'lon': -77.7038413, 'altMeters': 0.0, 'nodeName': 'C8', 'nodeType': 'Customer', 'popupText': 'C8', 'leafletIconPrefix': 'glyphicon', 'leafletIconType': 'star', 'leafletColor': 'orange', 'leafletIconText': '9', 'cesiumIconType': 'pin', 'cesiumColor': 'orange', 'cesiumIconText': '9', 'elevMeters': None},
]
nodesDF = pd.DataFrame(nodesArray)
# Get the raw distance matrix from Veroviz. NOTE: This takes ~ 30 seconds - not instantaneous.
[_, dist] = vrv.getTimeDist2D(nodes = nodesDF,
matrixType = 'all2all',
routeType = 'fastest',
dataProvider = 'OSRM-online',
outputDistUnits = 'miles')
# This is the "Big M" in linear programming terminology - a large constant
# used to impose prohibitively high costs on infeasible choices.
BIG_M = 9999.9
# Set distances from each node to itself as BIG_M to ensure these paths are
# never taken.
dist = {(nodesArray[i]['nodeName'], nodesArray[j]['nodeName']):dist[i,j] if i != j else BIG_M
for i in range(len(nodesArray)) for j in range(len(nodesArray))}
# View the nodes on a map - blue/home markers are depots and orange/star markers are customers.
vrv.createLeaflet(nodes=nodesDF, mapBackground='Arcgis Topo')
Since we want this to be generic and able to be reused with different data, we will construct a function to create the IP model using PuLP.
# The inputs for this function are as follows:
# trucks: A list of distinct truck labels (e.g., ['T1', 'T2', ...])
# depots: A list of distinct depot labels (e.g., ['D1', 'D2', ...])
# customers: A list of distinct customer labels (e.g., ['C1', 'C2', ...])
# products: A list of distinct product labels (e.g., ['C1', 'C2', ...])
# prod_availability: A dict of form {<depot label>:{<product label>:<units available>, ...}, ...}
# prod_volume: A dict of form {<product label>:<volume>}
# truck_cap: A dict of form {<truck label>:<volume capacity>}
# demand: A dict of form {<customer label>}:{<product label>:<units ordered>, ...}, ...}
# cost_per_mile: A dict of form {<truck label>:<cost per mile>}
# dist_matrix: dict of form {(<location label>, <location label>):<distance>, ...}. Locations are depots + customers
# truck_base: A dict of form {<truck label>:<depot label>}
# big_m: An arbitrarily large constant.
def build_model(trucks, depots, customers, products, prod_availability,
prod_volume, truck_cap, demand, cost_per_mile,
dist_matrix, truck_base, big_m):
H = trucks
I = depots
J = customers
K = products
L = I + J # Set of all distinct locations
a = prod_availability
e = prod_volume
E = truck_cap
r = demand
d = dist_matrix
c = cost_per_mile
b = {h:{i:1 if truck_base[h] == i else 0 for i in I} for h in H}
# u[h][j][k] = Units of product k delivered to customer j on truck h
u = LpVariable.dicts("u", (H, J, K), 0, None, LpInteger)
# x[h][i][j] == 1 if truck h travels directly from location i to j, 0 otherwise
x = LpVariable.dicts("x", (H, L, L), 0, 1, LpInteger)
prob = LpProblem("MultiDepot_VRP", LpMinimize)
prob += (lpSum([c[h] * d[i,j] * x[h][i][j] for h in H for i in L for j in L]),
'Total_Cost') # Objective Function
for h in H:
prob += (lpSum([e[k] * u[h][j][k] for j in J for k in K]) <= E[h]) # Ensure no truck capacity exceeded
for i in L:
prob += (lpSum([x[h][j][i] for j in L]) == lpSum([x[h][i][j] for j in L])) # For each loc, truck in -> truck out
for i in I:
prob += (lpSum([x[h][i][j] for j in L]) <= b[h][i]) # Ensure no truck leaves a non-base depot
for W in allcombinations(J, len(J)):
prob += (lpSum([x[h][i][j] for i in W for j in W]) <= len(W) - 1) # Ensure no subset of customers contains a circuit.
for k in K:
for i in I:
prob += (lpSum([b[h][i] * u[h][j][k] for h in H for j in J]) <= a[i][k]) # No depot ships more of a product than it has
for j in J:
prob += (lpSum(u[h][j][k] for h in H) == r[j][k]) # Each customer gets the products they ordered
for h in H:
prob += (u[h][j][k] <= big_m * lpSum(x[h][i][j] for i in L)) # No truck carries products for a customer unless it visits the customer.
return prob
Below we first construct the IP model and solve it. If you are using the CBC solver that comes with PuLP, it may take 15-30 seconds to find a solution. If you have access to the CPLEX solver, solutions can be obtained in a couple of seconds - just comment out the model.solve() line and uncomment the two below it.
# Construct the model
model = build_model(truck_labels,
depot_labels,
customer_labels,
product_labels,
depot_product_availabilities,
product_volumes,
{t:trucks[t]['capacity'] for t in trucks.keys()},
customer_orders,
{t:trucks[t]['cost_per_mile'] for t in trucks.keys()},
dist,
{t:trucks[t]['base'] for t in trucks.keys()},
BIG_M)
# Set the initial time
start_time = time.time()
# Solve the model using the CBC solver that comes with PuLP
model.solve()
# If you have the faster CPLEX solver and want to use it, uncomment the two
# lines below and comment out the one above.
# solver = getSolver('CPLEX_CMD')
# model.solve(solver)
# Compute the elapsed time
elapsed_time = time.time() - start_time
# Output the total cost of transportation and computational time
print("Total Cost of Transportation = ", value(model.objective))
print("Computational Time = ", round(elapsed_time, 1), 'seconds')
Total Cost of Transportation = 105.69849876346827 Computational Time = 24.2 seconds
In this step we need to extract the truck itineraries from the decision variables. For each truck, we want to know its stops and which products to deliver at each stop. To get this information, we need to sift through the non-zero decision variables from the solved model.
# Populate truck path and customer order dicts by extracting variable values
# from the solved model.
path, order = {}, {}
for v in model.variables():
var, truck, i, j, val = v.name.split('_') + [v.varValue]
if var == 'x' and val == 1:
if not(truck in path):
path[truck] = {}
path[truck][i] = j
if i.startswith('D'):
path[truck]['depot'] = i
elif var == 'u' and val > 0:
if not(truck in order):
order[truck] = {}
if not(i in order[truck]):
order[truck][i] = []
order[truck][i].append((j, val))
# Utility function - given a truck label, build its route as a list (excluding the depot).
def build_truck_route(truck):
h = path[truck]
depot = h['depot']
curr_loc = depot
route = []
while h[curr_loc] != depot:
route.append(h[curr_loc])
curr_loc = h[curr_loc]
return route
# Build a dict for all truck routes including the depots. We will use this later.
truck_routes = {}
# Build Itinerary Table
itinerary_list = []
for t in truck_labels:
if t in path :
depot = path[t]['depot']
route = build_truck_route(t)
truck_routes[t] = [depot] + route + [depot]
stop_num = 1
for stop in route:
delivery = str(order[t][stop]).replace('[', '').replace(']', '').replace("'", '')
itinerary_list.append({'Truck':t, 'Stop Number':str(stop_num),
'Customer':stop, 'Delivery (Product, Quantity)':delivery})
stop_num += 1
# Create a data frame to hold the truck iteneraries and show it
itinerary_df = pd.DataFrame(itinerary_list)
itinerary_df
| Truck | Stop Number | Customer | Delivery (Product, Quantity) | |
|---|---|---|---|---|
| 0 | T1 | 1 | C6 | (P2, 2.0), (P3, 1.0) |
| 1 | T1 | 2 | C8 | (P1, 1.0), (P2, 1.0) |
| 2 | T2 | 1 | C1 | (P3, 1.0) |
| 3 | T2 | 2 | C3 | (P1, 1.0) |
| 4 | T2 | 3 | C7 | (P2, 1.0), (P3, 1.0) |
| 5 | T3 | 1 | C5 | (P1, 2.0) |
| 6 | T4 | 1 | C4 | (P4, 2.0) |
| 7 | T4 | 2 | C2 | (P4, 1.0) |
| 8 | T4 | 3 | C1 | (P2, 1.0) |
Below we use VeRoViz to build the route visualizations.
# Construct a mapping from node labels (depots + customers) to node objects
node_dict = {x['nodeName']:x for x in nodesArray}
# Specify the truck colors
truck_colors = {'T1':'blue', 'T2':'red', 'T3':'black', 'T4':'green'}
# Construct an empty data frame to hold the routes
assignmentsDF = vrv.initDataframe('assignments')
# For each truck, add each one of its stops in order to the assignmentsDF
for t in truck_routes:
for i in range(1, len(truck_routes[t])):
start_node = truck_routes[t][i - 1]
end_node = truck_routes[t][i]
shapepointsDF = vrv.getShapepoints2D(
objectID = t,
modelFile = 'veroviz/models/ub_truck.gltf',
startLoc = [node_dict[start_node]['lat'], node_dict[start_node]['lon']],
endLoc = [node_dict[end_node]['lat'], node_dict[end_node]['lon']],
routeType = 'fastest',
leafletColor = truck_colors[t],
dataProvider = 'OSRM-online')
assignmentsDF = pd.concat([assignmentsDF, shapepointsDF],
ignore_index=True, sort=False)
# Construct the leaflet map and view. NOTE: This also takes a little time - not instantaneous.
vrv.createLeaflet(nodes=nodesDF, arcs=assignmentsDF, mapBackground='Arcgis Topo')
Message: The destination point (lat: 38.3706877, lon: -77.5007117) is 12.8 meters away from the road. You might find a gap between destination point and the route. Message: The origin point (lat: 38.3706877, lon: -77.5007117) is 12.8 meters away from the road. You might find a gap between the origin point and the route.