#!/usr/bin/env python # coding: utf-8 # # Multi-Depot Last Mile Delivery Routing # # 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](https://medium.com/hackernoon/mixed-integer-programming-a-straight-forward-tutorial-41cc50fb9c23) (IP). We will use the [PuLP](https://coin-or.github.io/pulp/) library to build and solve the IP model, and [VeRoViz](https://veroviz.org) to visualize the truck routes. # In[1]: # Import the libraries needed. from pulp import * import pandas as pd import veroviz as vrv import time # ## Step 1: Prepare the Scenario Data # # 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](https://veroviz.org/sketch.html). 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. # In[2]: # 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') # ## Step 2: Create the Model # # 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. # In[3]: # 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 {:{:, ...}, ...} # prod_volume: A dict of form {:} # truck_cap: A dict of form {:} # demand: A dict of form {}:{:, ...}, ...} # cost_per_mile: A dict of form {:} # dist_matrix: dict of form {(, ):, ...}. Locations are depots + customers # truck_base: A dict of form {:} # 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 # ## Step 3: Solve the Model # # 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. # In[4]: # 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') # ## Step 4: Extract and View the Truck Itineraries # # 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. # In[5]: # 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 # ## Step 5: Build and Display the Truck Route Visualizations # # Below we use VeRoViz to build the route visualizations. # In[6]: # 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') # In[ ]: