Guide to Network Analysis (Part 7 - Vehicle Routing Problem)¶
Table of Contents
- Introduction
- Data Preparation
- Solution 1: A Basic Scenario
- The basic scenario
- Optional Attributes
- Solve the VRP
- Tabularizing the response from solve_vehicle_routing_problem
- Visualizing the response from from solve_vehicle_routing_problem
- Animating the response from from solve_vehicle_routing_problem
- Saving the response from from solve_vehicle_routing_problem to online
- A second solver in arcgis.features module
- Solution 2: A modified scenario
- Solution 3: Delineates work territories
- Conclusions
- References
Introduction¶
Now we have learned about Network Datasets and Network Analysis services in Part 1, how to find routes from one point to another, and among multiple points in Part 2, how to generate service area in Part 3, how to find closest facility in Part 4, how to create an Origin Destination Cost Matrix in Part 5, how to solve location allocation in Part 6, let's move onto the seventh topic - how to perform Vehicle Routing Problem
service. Please refer to the road map below if you want to revisit the previous topics or jump to the next topic -
- Network Dataset and Network Analysis services (Part 1)
- Find Routes (Part 2)
- Generate Service Area (Part 3)
- Find Closest Facility (Part 4)
- Generate Origin Destination Cost Matrix (Part 5)
- Solve Location Allocation (Part 6)
- Vehicle Routing Problem Service (You are here!)
What is a Vehicle Routing Problem?¶
The vehicle routing problem
(VRP) is a superset of the traveling salesman problem (TSP). In a TSP, one set of stops is sequenced in an optimal fashion. In a VRP, a set of orders needs to be assigned to a set of routes or vehicles such that the overall path cost is minimized. It also needs to honor real-world constraints including vehicle capacities, delivery time windows, and driver specialties. The VRP produces a solution that honors these constraints while minimizing an objective function composed of operating costs and user preferences, such as the importance of meeting time windows [1].
The VRP solver starts by generating an origin-destination matrix of shortest-path costs between all order and depot locations along the network. Using this cost matrix, it constructs an initial solution by inserting the orders one at a time onto the most appropriate route. The initial solution is then improved upon by re-sequencing the orders on each route, as well as moving orders from one route to another, and exchanging orders between routes. The heuristics used in this process are based on a tabu search metaheuristic and are proprietary, but these have been under continual research and development in-house at Esri for many years and quickly yield good results [1].
When is the VRP service applicable?¶
Various organizations service orders with a fleet of vehicles. For example, a large furniture store might use several trucks to deliver furniture to homes. A specialized grease recycling company might route trucks from a facility to pick up used grease from restaurants. A health department might schedule daily inspection visits for each of its health inspectors. The problem that is common to these examples is the vehicle routing problem (VRP) [2].
Each organization needs to determine which orders (homes, restaurants, or inspection sites) should be serviced by each route (truck or inspector) and in what sequence the orders should be visited. The primary goal is to best service the orders and minimize the overall operating cost for the fleet of vehicles. The VRP service can be used to determine solutions for such complex fleet management tasks. In addition, the service can solve more specific problems because numerous options are available, such as matching vehicle capacities with order quantities, providing a high level of customer service by honoring any time windows on orders, giving breaks to drivers, and pairing orders so they are serviced by the same route [2].
About the Async execution mode¶
The maximum time an application can use the vehicle routing problem service when using the asynchronous execution mode is 4 hours (14,400 seconds). If your request does not complete within the time limit, it will time out and return a failure. When using the synchronous execution mode, the request must complete within 60 seconds. If your request takes longer, the web server handling the request will time out and return the appropriate HTTP error code in the response [2].
Work with ArcGIS API for Python¶
The ArcGIS API for Python provides a tool called solve_vehicle_routing_problem
to solve the vehicle routing problems, which is shown in the table below, along with other tools we have learned so far from previous chapters. Or user can still use plan_routes
for VRP analysis.
Operation | network.analysis | 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 |
These two methods are defined in different modules of the arcgis package, and will make distinct REST calls in the back end. A key separation from network.analysis
to features.use_proximity
is that the former provides full capabilities of solvers and runs faster, and the latter is workflow-driven and provides service-to-service I/O approach.
Defined in the network.analysis
module, solve_vehicle_routing_problem
supports full capabilities of operations; while plan_routes
provides a workflow approach that user can input a feature service and get returned a feature service. We will walk through the data preparation, implementation, and visualization of output here. Remember that if you run the solve_vehicle_routing_problem
with ArcGIS Online, 2 credits will be consumed per usage.
Problem statement¶
The goal of part 7 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. Each store has a specific quantity of demand for the goods, and each truck has a limited capacity for carrying the goods. The main objective is to assign trucks in the fleet a subset of the stores to service and to sequence the deliveries in a way that minimizes the overall transportation costs.
This can be achieved by solving a vehicle routing problem (VRP). Once the delivery sequence is determined, you will generate the turn-by-turn directions for the resulting routes, which can be electronically distributed or printed and given to the drivers to make the deliveries [4].
Three examples will be demonstrated in the following sections, covering three commonly seen scenarios, and they are namely:
- Basic scenario, given the stores to visit, the distribution center to load supplies, and the vehicle(s) to deliver goods;
- Modified scenario, when one of the truck drivers go on vacation, and overtime is required;
- With work territories delineated, assuming that certain areas cannot be visited on the route (or under certain penalties if visited).
Before diving into the implementation, let's first prepare the required input data.
Data Preparation¶
As a first step, let's import required libraries and establish a connection to your organization which could be an ArcGIS Online organization or an ArcGIS Enterprise.
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
If you have already set up a profile to connect to your ArcGIS Online organization, execute the cell below to load the profile and create the GIS class object. If not, use a traditional username/password log-in e.g. my_gis = GIS('https://www.arcgis.com', 'username', 'password', verify_cert=False, set_active=True)
my_gis = GIS('home')
To solve the Vehicle Routing Problem, we need orders layer
with stop information, depots layer
with the warehouse location information from where the routes start and routes table
with constraints on routes like maximum total time the driver can work etc. To provide this information to the service, different types of inputs are supported as listed below:
- An existing feature service that contains information for
orders
(grocery stores) anddepots
(the distribution center) - CSV files for self defined routes
- JSON variables for hand-picked prohibited/restricted areas
Let's see how to extract the feature classes from the existing service:
Define Input Feature Class¶
The existing Feature Service item contains the sublayer (id=0) for distribution center, and sublayer(id=1) for all 25 grocery stores. We will search for the item, create FeatureLayer
object per sublayer, and then create a FeatureSet
class object using query()
.
try:
sf_item = my_gis.content.get("fa809b2ae20a4c18959403d87ffdc3a1")
display(sf_item)
except RuntimeError as re:
print("You dont have access to the item.")
orders layer¶
First, we need to get the orders
feature class (in this case, grocery stores) - Use this parameter to specify the orders the routes should visit. An order can represent a delivery (for example, furniture delivery), a pickup (such as an airport shuttle bus picking up a passenger), or some type of service or inspection (a tree trimming job or building inspection, for instance). When specifying the orders, you can specify additional properties for orders using attributes, such as their names, service times, time windows, pickup or delivery quantities etc.
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. Use this parameter to specify a 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.
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¶
Next, we will create routes feature class with csv file. A route specifies vehicle and driver characteristics. A route can have start and end depot service times, a fixed or flexible starting time, time-based operating costs, distance-based operating costs, multiple capacities, various constraints on a driver’s workday, and so on. When specifying the routes, you can set properties for each one by using attributes. Attributes in the csv are explained below.
Name
- The name of the routeStartDepotName
- 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
from dataframe, we convert the CSV to a pandas data frame using read_csv
function. Note that in our CSV, EarliestStartTime
and LatestStartTime
values are represented as strings denoting time in the local time zone of the computer. So we need to parse these values as date-time values which we accomplish by specifying to_datetime
function as the datetime parser.
When calling arcgis.network.analysis.solve_vehicle_routing_problem
function we need to pass the datetime values in milliseconds since epoch. The routes_df
dataframe stores these values as datetime type. We convert from datetime to int64 datatype which stores the values in nano seconds. We then convert those to milliseconds [4].
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)
Visualize the problem set¶
Before moving onto the solution, let's take a look at the visualization of what input data 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
#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)
# 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.
map0 = my_gis.map('San Francisco, CA')
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)
Once you have all the inputs as featuresets, you can 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 featureset geocoded from address, orders and routes could be read from csv files to convert to featureset.
Now, we are ready to explore the implementations with three practical examples:
Solution 1: A Basic Scenario¶
The basic scenario¶
Assuming that the requirements for the basic scenario is solving the problem of how to dispatch the three trucks in San Francisco (working from 8AM to 5PM) in delivering goods to 25 different stores. In the basic scenario, the distributor is given three required input parameters:
orders
You will add the grocery store locations to theOrders
feature class. You can think of orders as orders to be filled, 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 attributes of Stores contain information about the total weight of goods (in pounds) required at each store, the time window during which the delivery has to be made, and the 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 theDistributionCenter
feature class. The distribution center operates between 8:00 a.m. and 5:00 p.m.routes
The distribution center has three trucks, each with a maximum capacity to carry 15,000 pounds of goods. You will add three routes (one for each vehicle) and set the properties for the routes based on the center's operational procedures.
Optional Attributes¶
Other optional attributes include:
- If we need driving directions for navigation, populate_directions must be set to true.
Time Attribute = TravelTime (Minutes)
The VRP solver will use this attribute to calculate time-based costs between orders and the depot. 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. Use the default here.Default Date
is set to be the day of today (i.e. Monday)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 is 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 in minutes. Distance Field Units
is set to Miles. This specifies that all distance-based attributes, such as MaxTotalDistance and CostPerUnitDistance for Routes, are 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
orTrue Shape option
for theOutput Shape Type
. Note that this option only affects the display of the routes, not the results determined by the VRP solver. - Using
Use Hierarchy
as default here (a.k.a. True).
You can set save_route_data
to True if you want to save the route data from result to local disk, which would then be used to upload to online to share with drivers eventually and share the routes in ArcGIS online or Enterprise. Individual routes are saved as route layers which could then be opened in navigator with directions(if you solve with populate_directions
=True) [4].
Solve the VRP¶
The following operations can help you sort out the basic scenario - how to dispatch the three trucks in San Francisco (working from 8AM to 5PM) in delivering goods to 25 different stores. The output will also include the driving directions in Spanish.
Also note that you can set the if_async
variable to True, when you need to execute multiple solvers in parallel.
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="es",
future=if_async)
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.
Only when the job is finished and shown as succeeded can we proceed to explore the results. Otherwise, skip the rest of this section and check out the solution 2 instead.
if if_async:
if result1.done():
result1 = result1.result()
print("Async job done!")
else:
print("Async job not done yet!")
print('Analysis succeeded? {}'.format(result1.solve_succeeded))
Here result1
is a arcgis.geoprocessing._support.ToolOutput
Class object, and contains multiple objects - out_routes (FeatureSet), out_stops(FeatureSet), etc. Since that we have enabled save_route_data
, out_route_data
will appear in the resulting tool output as a dictionary object that is the url pointing to the zipped file of the route data (saved on the GIS object).
result1
Tabularizing the response from solve_vehicle_routing_problem¶
Now, let's explore the tabularized output from solve_vehicle_routing_problem
. What will be useful for distributor and the drivers will be the summarized route information, and sequences of stops per route.
# Display the analysis results in a pandas dataframe.
out_routes_df = result1.out_routes.sdf
out_routes_df[['Name','OrderCount','StartTime','EndTime',
'TotalCost','TotalDistance','TotalTime','TotalTravelTime','StartTimeUTC','EndTimeUTC']]
Based on the dataframe display of the out_routes object, we can tell the optimal routing option provided by solve_vehicle_routing_problem
is for Truck_1 to visit 8 stops, Truck_2 to visit 6 stops, and Truck_3 to visit 11 stops. Upon this selection, the total cost will be 162.13 + 72.36 + 186.84 = 421.33, the total distance is 55.13 + 13.22 + 63.70 = 132.05, and the total travel time will be 149.15 + 55.65 + 145.42 = 350.22.
Scenario | Total Cost | Total Distance | Total Travel Time | Scheduled Stops |
---|---|---|---|---|
#1 | 421.33 | 132.05 | 350.22 | [8,6,11] |
out_stops_df = result1.out_stops.sdf
out_stops_df[['Name','RouteName','Sequence','ArriveTime','DepartTime']].sort_values(by=['RouteName',
'Sequence'])
Visualizing the response from from solve_vehicle_routing_problem¶
In order to improve the re-usability of codes, we will define a method called visualize_vehicle_routing_problem_results
to render the map, and visualize the orders
, depots
and the routing results calculated by the VRP solver. This method will be reused in scenarios 2 and 3 as well.
# 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, CA')
map1.basemap = 'dark-gray'
map1.layout.height = '650px'
map1
# Call custom function defined earlier in this notebook to
# display the analysis results in the map.
visualize_vehicle_routing_problem_results(map1, result1,
orders_fset=stores_fset, depots_fset=distribution_center_fset, zoom_level=8)
Judging from what's displayed in map1
, Truck_1 (blue) tends to take care of the stores located at the east side of San Francisco, while the Truck_2 (red) and Truck_3 (green) are responsible for delivering goods to stores located at the west. Also, the difference between Truck_2 and Truck_3 is that the former handles the downtown area, and the latter focuses on the outer rim.
Animating the response from from solve_vehicle_routing_problem¶
In order to show a stronger sequential relationship between origin, stops and destination of each solved route, we can also use animate_vehicle_routing_problem_results
function to be defined below, to animate each stop along the route sequentially:
# Display the analysis results in a map.
# Create a map of SF, California.
map1a = my_gis.map('San Francisco, CA')
map1a.basemap = 'dark-gray'
map1a.layout.height = '650px'
map1a