The Vehicle Routing Problem (VRP) is a superset of the traveling salesman problem, which attempts to find the shortest route between cities without revisiting any of them. The twist with VRP is that the route must also minimize cost. 

VRPs must also adhere to actual-world constraints such as driver specialties, delivery windows and vehicle capacities. VRPs must come up with a solution that honors these constraints, while minimizing an objective function composed of operating costs and user preferences.

VRPs typically follow a set formula:

1. Generate an origin-destination matrix of shortest-path costs between all order and depot locations along a given network. 
2. Construct an initial solution by inserting the orders one at a time onto the most appropriate route. 
3. Improve upon the initial solution by re-sequencing the orders on each route, as well as moving orders from one route to another, and exchanging orders between routes. 

Depending on the complexity of the problem, optimized solutions can take minutes to hours to be obtained.

Using ArcGIS API for Python to Minimize Vehicle Routing

The ArcGIS API for Python is designed to make it easy for developers to work with maps and geospatial data. The library provides a tool called solve_vehicle_routing_problem designed to solve (obviously) vehicle routing problems, but it also includes other relevant tools shown in the table below:

OPERATION	METHOD: Network.analysis	METHOD: Features.use_proximity
Route	Find_routes	Plan_routes
ServiceArea	Generate_service_areas	Create_drive_time_areas
ClosestFacility	Find_closest_facilities	Find_nearest
OD Cost Matrix	Generate_origin_destination_cost_matrix	Connect_origins_to_destinations
Location Allocation	Solve_location_allocation	Choose_best_facilities
Vehicle Routing Problem	Solve_vehicle_routing_problem	Plan_routes

The two methods are defined in different modules of the ArcGIS package, and will make distinct REST calls to the back end: 

• Network.analysis provides full solver capabilities and runs fast
• Features.use_proximity is workflow-driven and provides a service-to-service I/O approach

The data preparation, implementation, and visualization of output is all going to be done here.

Vehicle Routing Problem Statement 

The main objective is to find the best routes for a fleet of vehicles operated by a distribution company to deliver goods from a distribution center to a set of 25 grocery stores. 

But each store has a specific quantity of demand for the goods and each truck has a limited capacity for carrying the goods, so we’ll need to:

• Assign trucks in the fleet to a subset of the stores to service
• Sequence the deliveries in a way that minimizes the overall transportation costs

We’ll solve the problem by approaching it as a VRP: 

1. Prepare the data.
2. Determine the pickup and delivery points.
3. Visualize the results by plotting the pickup->delivery routes on a map.
4. Optimize the pickup and delivery routes by solving the VRP.

The first step, as always, is to prepare the required input data.

1 — Vehicle Routing Data Preparation 

Import required libraries and establish a connection to an organization which could be an ArcGIS Online organization or an ArcGIS Enterprise. The code is shown below:

from arcgis.gis import GIS
import arcgis.network as network
from arcgis.features import FeatureLayer, Feature, FeatureSet, FeatureCollection, analysis
import pandas as pd
import time
import datetime as dt

To solve the VRP we’ll need: 

• An orders layer with stop information
• A depots layer with the warehouse location information from where the routes start
• A routes table with constraints on routes like maximum total time the driver can work

To provide all of this information, we’ll need a set of inputs:

• An “Input Feature” service that contains information for orders (grocery stores) and depots (the distribution center)
• CSV files that contain information about self-defined routes
• A set  of JSON variables for hand-picked prohibited/restricted areas

Define The Input Feature Class

The existing Feature Service item contains the sublayer (id=0) for distribution centers and sublayer(id=1) for all 25 grocery stores. The following code will:

1. Search for the item
2. Create FeatureLayer object per sublayer
3. Create a FeatureSet class object using query ()

try:
 sf_item = my_gis.content.get("fa809b2ae20a4c18959403d87ffdc3a1")
 display(sf_item)
except RuntimeError as re:
 print ("You don't have access to the item.")

Orders Layer

First, we’ll define the orders feature class (in this case, grocery stores) and use this parameter to specify the orders the routes should visit. An order can represent:

• A delivery, such as a furniture delivery
• A pickup, such as an airport shuttle bus picking up a passenger
• A service, such as a tree trimming job 
• An inspection, such as a building inspection, etc

When specifying the orders, we’ll specify additional properties for orders using attributes, such as their names, service times, time windows, pickup or delivery quantities, and more. 

Below is the orders layer code:

stores_fl = sf_item.layers[1]
try:
 stores_fset = stores_fl.query(where="1=1", as_df=False)
 display(stores_fset)
except RuntimeError as re:
 print ("Query failed.")
--------------------------------------------------------------------------
for f in stores_fset:
 tmp1 = f.get_value("TimeStart1")
 tmp2 = f.get_value("TimeEnd1")
 f.attributes.update({"TimeWindowStart1":tmp1, "TimeWindowEnd1":tmp2})

Depots Layer 

Depots in this case can be interpreted as the distribution center. We’ll use this parameter to specify the location that a vehicle departs from at the beginning of its workday and returns to at the end of the workday. Vehicles are loaded (for deliveries) or unloaded (for pickups) at depots at the start of the route. 

Here’s the code for the depots layer:

distribution_center_fl = sf_item.layers[0]

try:
 distribution_center_fset = distribution_center_fl.query(where="1=1", as_df=False)
 display(distribution_center_fset)
except RuntimeError as re:
 print ("Query failed.")

Routes Table

A route specifies vehicle and driver characteristics, including: 

• Start and end depot service times
• A fixed or flexible starting time
• Time-based operating costs
• Distance-based operating costs
• Constraints on a driver’s workday, and so on. 

When specifying the routes, we’ll set properties for each one by using attributes. We’ll generate all of this information in a CSV file that includes: 

• Name – The name of the route
• StartDepotName – The name of the starting depot for the route. This field is a foreign key to the Name field in Depots. 
• EndDepotName – The name of the ending depot for the route. This field is a foreign key to the Name field in the Depots class.
• EarliestStartTime – The earliest allowable starting time for the route.
• LatestStartTime – The latest allowable starting time for the route.
• Capacities – The maximum capacity of the vehicle.
• CostPerUnitTime – The monetary cost incurred per unit of work time, for the total route duration, including travel times as well as service times and wait times at orders, depots, and breaks.
• MaxOrderCount – The maximum allowable number of orders on the route.
• MaxTotalTime – The maximum allowable route duration.

To get a FeatureSet we’ll convert the CSV file into a pandas data frame using panda’s read_csv function. Note that in this CSV, EarliestStartTime and LatestStartTime values are represented as strings denoting time in the local time zone of the computer. So in the next section, we’ll parse these values as date-time values, which is accomplished by specifying the to_datetime function as the datetime parser.

2 — Vehicle Routing Pickup and Delivery Timing

Now that we have all our orders, routes and depot data, we need to be able to calculate the time taken to obtain an order from a depot and route it to a destination.

When calling the arcgis.network.analysis.solve_vehicle_routing_problem function, we’ll pass the datetime values in milliseconds since epoch. The routes_df dataframe stores these values as a datetime type. We’ll need to convert from the datetime type to an int64 data type, which stores the values in nanoseconds, and then convert those to milliseconds:

routes_csv = "data/vrp/routes.csv"
# Read the csv file
routes_df = pd.read_csv(routes_csv, parse_dates=["EarliestStartTime", LatestStartTime"],
                   date_parser=pd.to_datetime)
routes_df["EarliestStartTime"] = routes_df["EarliestStartTime"].astype("int64") / 10 ** 6
routes_df["LatestStartTime"] = routes_df["LatestStartTime"].astype("int64") / 10 ** 6
routes_df
routes_fset = FeatureSet.from_dataframe(routes_df)
display(routes_fset)

3 — Vehicle Routing Visualization

Let’s visualize the routes we currently have, namely the depots and the orders:

#Define a function to display the problem domain in a map
def visualize_vehicle_routing_problem_domain(map_widget, orders_fset,     depots_fset, zoom_level, route_zones_fset = None):
#The map widget
map_view_outputs = map_widget
 #Visualize the inputs with different symbols
 map_view_outputs.draw(orders_fset, symbol={"type": "esriSMS",
 "style": "esriSMSCircle",
 "color": [76,115,0,255],"size": 8})
 map_view_outputs.draw(depots_fset, symbol={"type": "esriSMS",
 "style": "esriSMSSquare",
 "color": [255,115,0,255], "size": 10})
 if route_zones_fset is not None:
    route_zones_sym = {
	      "type": "esriSFS",
      "style": "esriSFSSolid",
      "color": [255,165,0,0],
      "outline": {
	      "type": "esriSLS",
      "style": "esriSLSSolid",
      "color": [255,0,0,255],
      "width": 4}
    }
 map_view_outputs.draw(route_zones_fset, symbol=route_zones_sym)
#Zoom out to display all of the allocated census points.
 map_view_outputs.zoom = zoom_level
#Display the analysis results in a map.
#Create a map of Nairobi, Kenya.
map0 = my_gis.map ('Nairobi, KEN')
map0.basemap = 'dark-gray'
map0.layout.height = '650px'
map0
# Call custom function defined earlier in this notebook to
# display the analysis results in the map.
visualize_vehicle_routing_problem_domain (map0, orders_fset=stores_fset,
 depots_fset=distribution_center_fset, zoom_level=8)

After setting all the inputs as feature sets, it is possible to pass inputs converted from different formats. The preparation step shown above is not the only way to do it. For example, “depot” could be a feature set geo-coded from address, while orders and routes could be read from a CSV file which is then converted to a feature set.

Now that we have all the data and general routes, we can explore how those routes might change when subject to various practical scenarios. 

4 — How to Optimize Deliveries by Solving VRP 

Now let’s create a practical problem involving three trucks in San Francisco (working from 8AM to 5PM) that need to deliver goods to 25 different grocery stores. 

In this scenario, the distributor is given three required input parameters: 

• orders – here we’ll add the grocery store locations to the Orders feature class. Think of orders as orders to be fulfilled, since each grocery store has requested goods to be delivered to it from the distribution center. Members of the Orders class will eventually become stops along the vehicles’ routes.
  The attribute of “Stores” contains information about: 
  • Total weight of goods (in pounds) required at each store
  • Time window during which the delivery has to be made
  • Service time (in minutes) incurred while visiting a particular store. The service time is the time required to unload the goods.
• depots – the goods are delivered from a single distribution center whose location is shown in the DistributionCentre feature class. The distribution center operates between 8:00 a.m. and 5:00 p.m.
• routes – the distribution center has three trucks, each of which can carry a maximum of 15,000 pounds of goods. We’ll add three routes (one for each vehicle), and set the properties for the routes based on the center’s operational procedures.

Optional Attributes:

• populate_directions must be set to true if driving directions are needed for navigation.
• Time Attribute = TravelTime (Minutes) – the VRP solver will use this attribute to calculate time-based costs between orders and the depot. We’ll use the default here.
• Distance Attribute = Meters – this attribute is used to determine travel distances between orders and the depot for constraint purposes and creating directions. However, the VRP solver’s objective is to minimize time costs so we’ll use the default here.
• Default Date is set to be the day of today (i.e. Tuesday).
• Capacity Count is set to 1. This setting indicates that the goods being delivered have only one measurement. In this case, that measurement is weight (pounds). If the capacities were specified in terms of two measurements, such as weight and volume, then the capacity count would be set to 2.
• Minutes are selected for Time Field Units. This specifies that all time-based attributes, such as ServiceTime and MaxViolationTime1 for Orders and MaxTotalTime, MaxTotalTravelTime, and CostPerUnitTime for Route are all designated in minutes.
• Distance Field Units are set to miles. This specifies that all distance-based attributes, such as MaxTotalDistance and CostPerUnitDistance for Routes, are denoted in miles.
• Since it is difficult for these delivery trucks to make U-turns, set U-Turns at Junctions to Not Allowed.
• Select between Straight Line, True Shape with Measures or True Shape option for the Output Shape Type. Note that this option only affects the display of the routes, not the results determined by the VRP solver.
• Using Hierarchy as default here (a.k.a. True).

Individual routes are saved as route layers, which could then be opened in the navigator with directions (solve with populate_directions=True).

We can now solve the VRP:

if_async = False
%% time
current_date = dt.datetime.now().date()
result1 = network.analysis.solve_vehicle_routing_problem(orders=stores_fset,                  depots=distribution_center_fset,                                                            default_date=current_date,                                                            routes=routes_fset, populate_route_lines=True,                                                save_route_data=True,                                                            populate_directions=True,                                                            directions_language="en",

The VRP solver calculates the three routes required to service the orders, and draws lines connecting the orders. Each route begins and ends at the distribution center and serves a set of orders along the way.

The results can also be visualized on a map. In order to improve reusability, we’ll define a method called visualize_vehicle_routing_problem_results to render the map on which to visualize the orders, depots and the routing results calculated by the VRP solver:

# Define the route symbols as blue, red and green
route_symbols = [{"type": "esriSLS", "style": "esriSLSSolid",
                  "color": [0, 100, 240, 255], "size":10},
                 {"type": "esriSLS", "style": "esriSLSSolid",
                  "color": [255, 0, 0, 255], "size":10},
                 {"type": "esriSLS", "style": "esriSLSSolid",
                  "color": [100, 240, 0, 255], "size":10}
                 ]
# Define a function to display the output analysis results in a map
def visualize_vehicle_routing_problem_results(map_widget,                                            solve_vehicle_routing_problem_result,                                            orders_fset, depots_fset, zoom_level,                                               route_zones_fset=None):
    # The map widget
    map_view_outputs = map_widget
    
     # The solve_vehicle_routing_problem analysis result
    results = solve_vehicle_routing_problem_result
    # Visusalize the inputs with different symbols
    map_view_outputs.draw(orders_fset, symbol={"type": "esriSMS",
                                               "style": "esriSMSCircle",
                                               "color": [76, 115, 0, 255], "size": 8})
    map_view_outputs.draw(depots_fset, symbol={"type": "esriSMS",
                                               "style": "esriSMSSquare",
                                               "color": [255, 115, 0, 255], "size": 10})
    if route_zones_fset is not None:
         route_zones_sym = {
	            "type": "esriSFS",
            "style": "esriSFSSolid",
            "color": [255, 165, 0, 0],
            "outline": {
	                "type": "esriSLS",
                "style": "esriSLSSolid",
                "color": [255, 0, 0, 255],
                "width": 4}
        }
    map_view_outputs.draw(route_zones_fset, symbol=route_zones_sym)
   #Visualize each route
    for i in range(len(results.out_routes.features)):
        out_routes_flist = []
        out_routes_flist.append(results.out_routes.features[i])
        out_routes_fset = []
        out_routes_fset = FeatureSet(out_routes_flist)
        map_view_outputs.draw(out_routes_fset,
                              symbol=route_symbols[i % 3])
        # Zoom out to display all of the allocated census points.
        map_view_outputs.zoom = zoom_level
# Display the analysis results in a map.
# Create a map of SF, California.
map1 = my_gis.map('San Francisco, California')
map1.basemap = 'dark-gray'
map1.layout.height = '650px'
map1

The above results can also be animated on the map in order to show a stronger sequential relationship between origin, stops and destination of each solved route.

Conclusions – Why VRP is Critical

With the pandemic having normalized online ordering, delivery trucks have become more and more of a common sight on all of our streets. Optimizing delivery routes with a technique like VRP has a number of positive benefits, including:

• Minimizing gas consumption, which is key to profitability for shipping companies given the high price of gas these days.
• Minimizing miles driven, which cuts down on truck wear and tear in order to minimize repair costs.
• Speeding delivery, which can improve customer satisfaction.

While the ArcGIS API for Python is quite complex, once the data and basic routing points are input, it becomes an extremely flexible tool that allows you to solve realistic VRP issues that business face every day, from:

• Suddenly sick drivers, which decreases the number of trucks that can make deliveries on any given day. 
• The shifting need for overtime and work breaks on specific days.
• The addition of satellite distribution centers, which means trucks can pick up orders elsewhere than the main distribution center.

And so on. In each case, solving VRPs is the best way to minimize total cost, decrease travel time, and minimize distance traveled. 

Read Similar Stories