Benchmarking a fast, satisficing vehicle routing algorithm for public health emergency planning and response: “Good Enough for Jazz”

Related work

The unpredictable nature of disaster situations requires a modified approach to the capacitated vehicle routing problem (CVRP) to include a number of context specific constraints that are inherently variable and often only loosely defined. This can present some unique challenges within the vehicle routing problems (VRPs) within emergency planning, response, and MCM distribution. During our literature review, we were unable to find any algorithms that matched the same criteria or solved the exact same problem as the RSSD algorithm; however, we did discover algorithms in a similar application domain. See Table 1 for the strengths and weaknesses of some of these algorithms and RSSD algorithm. Table 1 contains check marks and × marks to indicate some of the unique aspects of the chosen algorithms. Brief descriptions of the algorithms and article contributions are included below.

Shen, Dessouky & Ordóñez (2009) and Shen, Ordóñez & Dessouky (2009), in two articles, propose a mix integer heuristic model with a focus of mitigating some uncertainties in the road network. Their assumptions include: (1) that the demands and time between PODs are unknown at the time of plan creation and, (2) when PODs are activated the demands at each POD location will not be met by the MCMs delivered. As a result of these assumptions, they derived a two-step algorithm approach for the planning stage and the operational phases. In the planning stage, their first algorithm optimally solves the routes (Shen, Dessouky & Ordóñez, 2009) and in the activation stage, they use a stochastic routing algorithm that optimizes the amount of demand that is not met (Shen, Ordóñez & Dessouky, 2009). Qin et al. (2017) propose an update to a Genetic Algorithm utilizing a single-depot VRP after an emergency, where supplies are limited. They attempt to minimize the distance and number of vehicle costs. Gharib et al. (2018) propose an artificial neural fuzzy inference system of clustering the demand points by crisis severity and other criteria. Once these clusters are formed, the points in the cluster are prioritized. They then used two commodity distribution models, MO Firefly and NSGA-II, for their heterogeneous fleet and multi-depot problem. Li & Chung (2018) purport that stochastic programming may not work in emergency situations and therefore deterministic models like CVRP and SDVRP are more robust in applications where the probability distribution is unknown. They propose a robust CVRP and SDVRP using robust optimization techniques to be used in disaster relief routing. Goli & Malmir (2020) use a covering tour approach with fuzzy demand to minimize the arrival time of the last vehicle to nodes. They utilize a mathematical model solved by credibility theory and harmony search algorithm to decide which points are not mandatory, using this they are able to send vehicles to nodes that have the most need.

The described algorithms create solutions that deliver products from a depot to nodes in emergency situations with time and/or capacity constraints. These algorithms, including RSSD, have formulated the VRP differently in each case. As a result, we determined that it would not be adequate to directly compare algorithms and their solutions due to the differences in the problems being solved. We chose to benchmark the solutions of RSSD against best known solutions for CVRP, treating these best known solutions as a form of ground truth.

Algorithm overview

The RSSD algorithm utilizes a greedy approach to identify time and capacity compliant delivery routes for the distribution of MCMs to PODs. The RSSD algorithm attempts to solve a version of the vehicle routing problem (VRP) that can be described as an open, capacitated, homogeneous multi-vehicular, MCM distribution, and multi-depot delivery VRP. The RSSD algorithm guarantees that the generated delivery routes are both time and capacity compliant without requiring the number of routes (or number of vehicles) as input. RSSD has been created to satisfice the constraints of MCM delivery during the planning and activation stages. Although, in this article, we focus on the use of RSSD after a deliberate or accidental release of biochemicals, the routing algorithm can be used in situations where: the number of vehicles required is unknown, a time and capacity constraint must be met, and a need exists for a quick convergence (i.e., solutions in seconds). A version of this algorithm was previously reported in Urbanovsky (2018). A flow chart of the algorithm is included below, Fig. 2. Appendix A includes further modes of explaining the algorithm: (1) RSSD in algorithm notation; and (2) text version of the flowchart, Fig. 2. The version of RSSD as described below and used for running the experiments has been implemented in a Python3 Jupyter Notebook and published at www.doi.org/10.17605/OSF.IO/5H3GS.

The RSSD algorithm performs under five assumptions: (1) all demand can be met; (2) capacity of the vehicles is homogeneous; (3) vehicles leave the depot at the same time; (4) routes have no minimum time; and (5) vehicles do not return to the depot. Due to these assumptions, RSSD can rely on creating solutions that will serve all the demand of the PODs, use the same capacity and time constraint for all routes, and create open routes.

Consider a graph G defined on a Euclidean plane $G=(V,E)$, composed of a vertex set of PODs $V=\{v_1,v_2,...v_n\}$ which includes a depot set, $p\subset V$, and a directed, symmetrical edge set from which a distance matrix, D, is made. The distance matrix, D contains all distances between each pair of nodes, $D_{ij}=t(v_i,v_j)$. In a non-Euclidean graph implementation, the edges may not be symmetrical because of one-way roads and/or inclusion of traffic conditions. For every $v_i$, there is a demand $c(v_i)$, except for the depot/s $p$, where the demand is 0, $c(p)=0$. When implemented within a biomedical emergency context, each node’s demand can be equal to the population of the catchment area that is being served by the POD. The time constraint for the length of the route, $\tau$, and the capacity constraint of the vehicles carrying the supplies, $\psi$, are also given. In the RE-PLAN implementation, the time constraint is chosen by the public health preparedness planners, and can be a range of time, where for every value within the range, solutions will be produced. The capacity constraint in the RE-PLAN implementation is the amount of individual MCMs an eighteen-wheeler can carry, which for an anthrax outbreak is around 211,200 doses. The output of the algorithm is a set of ordered lists made up of a number of $v\in V$, called routes $R=\{r_1,r_2,...r_m\}$ where $r_1=v_1,v_2,...v_k$. For each $r_m$, there is the distance that the route covers $t(r_m)$, and the total demand being served on the route is $c(r_m)$. This graph and its variables are input of the RSSD algorithm, which is depicted in Fig. 2.

The first step of the algorithm is to check if routes can be created with the given time constraint, $\mathrm{\forall }{v}_i\in \{V$\ $p\}$, $max(t(v_i,p))\le \tau$. This checks to see if the maximum distance from the vertices to their closest depot is less than the time constraint. If one of the distances from a vertex to the closest depot is greater than the time constraint, routes cannot be created and RSSD will return an empty set. After this initial check, another check on the capacity constraint is executed for every node. For every node that has a greater demand than the capacity constraint, $\mathrm{\forall }{v}_i\in \{V$\ $p\},c(v_i)\ge \psi$, a route will be created from the depot to that node $R\cup r_i:v_i$. The demand needed to be met at that node is then updated $c(v_i)-=\psi$.

The Time Compliant (yellow) section ensures created routes are within the given time constraint $\tau$ by surrounding a function, Spatial Partition (green), with a while loop that runs until all routes are time compliant $\mathrm{\forall }{r}_i\in R,t(r_i)\le \tau$. The Spatial Partition function, first identifies the nodes with the maximum distance from one another, $v_i,v_j=max(D$\ $p)$, not including the depot/s. These two nodes $v_i,v_j$ are selected as the first nodes of two routes, $v_i\in r_a$ and $v_j\in r_b$. This metaheuristic was implemented because it is unlikely that the two furthest apart nodes will be on the same route given a time constraint. Next, these two routes, each composed of one node, will add the closest node to itself. The route with the shortest distance $t(r_a)\le t(r_b)$ will then add the next closest node to the either end of the route $min(D_{r_c^0,V},D_{r_c^{-1},V})$. The routes will iteratively add the closest node to either end of the route, with the route that is currently the shortest being the one to add a node next. This continues until all of the nodes have been added to one of the routes, $|r_a\cup r_b|=n-|p|$. When all nodes are added to either $r_a$ or $r_b$, the closest depot will be added to the closest end of the routes. Then $t(r_a)$ and $t(r_b)$ will be checked against the time constraint, $\tau$. If either route is over $\tau$ that route will be sent through the Time Compliant section once again, but only with the nodes within the route as part of the distance matrix D and vertex set V is updated to only contain nodes from the route $r_i$. This Time Compliant section recursively repeats until every route meets the given time constraint.

The Capacity Compliant (purple) portion of the algorithm takes in the time compliant routes and ensures they meet the capacity constraint, $\psi$. The set of routes R that do not meet the capacity constraint, $c(r_i)\ge \psi$ have the closest nodes to the depot pruned $v_x=min(t(r_i,p))$ until $r_i$ meets the capacity constraint. The nodes that are pruned from the non-capacity compliant routes are added to a new route, $pr$. This metaheuristic is implemented because pruning nodes that are closest to the depot will more likely result in a route that is time compliant. If the new route of pruned nodes do not meet the time constraint, it is sent through the Time Compliant portion of the algorithm until every route made up of the pruned nodes meets the time constraint. Then the capacity constraint is again checked. If the routes composed of the pruned nodes do not meet the capacity constraint, these routes will go through the capacity compliant portion again. This process is repeated until all routes are time and capacity compliant. When all routes, R, are compliant the final solutions are returned.

The RSSD algorithm is composed of two main sections that separately ensure time and capacity compliance. The time compliance portion relies on the metaheuristic that nodes furthest apart in a graph, are less likely to meet a given time constraint when in the same route. This section of the algorithm can also be thought of as clustering nodes together to identify nodes for routes and could potentially be used in other domains when the number of routes is not provided. The capacity portion of the algorithm relies on the metaheuristic that the nodes closest to the depot when pruned from non-capacity compliant routes will more likely create a time compliant route. RSSD is a simple, quick, and greedy algorithm that was designed to fit the needs of public health preparedness planners by providing the number of routes and creating solutions that are time and capacity compliant.

Benchmarking methods

To benchmark the quality of RSSD solutions, we utilize best known solutions from the CVRPLIB (Uchoa et al., 2017) as if they are the ground truth. We were unable to find any algorithms that solve the exact same problem as RSSD, therefore, benchmarking the RSSD algorithm to such algorithms would not yield a fair comparison. Instead of comparing the RSSD algorithm to another algorithm, we endeavor to benchmark the solutions from RSSD to a best known set of solutions within the CVRP domain so to evaluate the overall performance of the RSSD’s route solutions.

As discussed in Related Works, the CVRPLIB repository contains regularly used CVRP benchmarking datasets. Many of the solutions have a designated best known solution within the problem space of CVRP. The algorithms that produce the best known solutions create closed, primarily capacity-optimized, and secondly distance-optimized routes. Although RSSD produces open routes and does not optimize for capacity nor time, the VRP that RSSD is ultimately solving is a CVRP, and thus using route solutions within this problem space is relevant to compare the solution results.

We focus our analysis on one dataset from CVRPLIB, which is listed as ‘Uchoa et al. (2014)’ and also designated with an ‘X’ in the online repository. This dataset was published along with the CVRPLIB repository in Uchoa et al. (2017). Every instance in dataset ‘X’ is on a [0, 1,000] $\times$ [0, 1,000] grid. The locations of the depot are varied from central (C), at the origin (E), or randomly placed in the grid (R). The positions of the nodes are either random (R), clustered (C), or partially random and clustered (RC). The demand distributions for the nodes are uniform (U); small values between 1–10 or 5–10; large values between 1–100 or 50–100; dependent on quadrant (Q); or many small values and a few large values (SL). The combinations of these attributes in the datasets are split close to uniformly across the 100 instances in the dataset. Each instance in the dataset includes a capacity constraint, the locations of all nodes (including the depot), and the demands at each node.

The first portion of the benchmarking section reports CPU time between algorithms that engage in the NP-Hard task of identifying an optimal solution and RSSD which instead satisfices for the above described instances. Although, we find that comparing algorithms not helpful in our benchmarking process, we felt that comparing the CPU times of the algorithms reported in Uchoa et al. (2017) and RSSD as evidence for choosing a satisficing approach over engaging with the NP-Hard nature of the of the problem. In Uchoa et al. (2017), they report three algorithms’ CPU time in minutes for each of the instances in dataset ‘X.’ These three algorithms are iterated local search-based matheuristic (ILS-SP), the unified hybrid genetic search (UHGS), and branch-cut-and-price (BRS). We will combine their table with the RSSD process run-time in the ‘Results’ section. The ILS-SP and UHGS tests were conducted on Xeon CPU with 3.07 gigahertz and 16 gigabytes of RAM on an Oracle Linux Server 6.4 as reported in Uchoa et al. (2017). UHGS was given a high termination criterion of up to 50,000 consecutive iterations without improvement (Uchoa et al., 2017). The BCP tests were run on a single core of an Intel i7-3960X 3.30 gigahertz processor with 64 gigabytes RAM as reported in Uchoa et al. (2017). The RSSD tests were conducted on an Intel i5-10500 3.10 gigahertz processor, with 32 gigabytes of RAM. The implementation of RSSD for this article was not parallelized, therefore, all processes were conducted on one core, see notebook published at at www.doi.org/10.17605/OSF.IO/5H3GS.

The goal of our research is to benchmark the solutions of RSSD to the published ‘X’ best known solutions from CVRPLIB (Uchoa et al., 2017). Two aspects of the best known solutions need to be mitigated: the closed routes in the ‘X’ solution set and the lack of a given time constraint. To fix the closed routes attribute of the solutions, we calculated the final metrics by removing the longest edge that joins the route to the depot, thus making the routes open. To resolve the absence of a time constraint for the instances, we use the best known solution’s longest route distance (again excluding the longest segment to the depot) as the time constraint, $\tau$.

To measure the quality of the best known and RSSD solutions, we utilize three metrics that quantify consequential components in the domain of MCM dispersal at the time of plan activation. When evaluating the RSSD solution set, it was particularly imperative that we analyzed whether the utilization of resources (e.g., vehicles) was excessive as a result of RSSD’s lack of optimization. We assumed that the utilization of resources would be higher using the non-optimized, satisficing RSSD algorithm to produce the routes but we endeavored to find out how much. Using our three metrics, we quantify this trade-off for not engaging with and/or sidestepping this NP-Hard problem.

The UnCap metric is a measure of the percentage of unused capacity in relation to the sum of the demand of all nodes in the instance. It is an indication of potentially excessive routes/vehicles in the solution than necessary.

(1) $$UnCap=\frac{(|R|\ast \psi )-c(V)}{c(V)}$$where R is the set of routes in a solution, $\psi$ is the capacity constraint, V is the set of nodes, and $c(V)$ is the demand of all nodes. RSSD does not optimize for capacity, but it does perform certain metaheuristics in order to make better-informed greedy decisions. The closer the metric is to $0$ indicates the less unused capacity in relation to the sum of all demands. Even with these meta-heuristics, we expected to see a significant difference between the best known and RSSD solutions in terms of unused capacity.

The RngMax metric measures the range between the solution’s route distances and evaluates the variance in route sizes. All PODs in a region must open at the same time in order to reduce the possibility of over-crowding at the PODs that open first. This amount measures potential excess resources of guards at POD locations. RngMax is a ratio between the maximum route distance and minimum route distance over the maximum route distance.

(2) $$RngMax=\frac{max(t(R))-min(t(R))}{max(t(R))}$$where R is the set of routes in a solution and $t(R)$ is route distance. A high value of the RngMax metric indicates a large variance in route lengths. If some routes are significantly longer than others, guards may be required for longer for the PODs that have already received deliveries on shorter routes. When routes have similar distances, PODs along each route will be served around the same time. The closer the RngMax score is to 0, the less distance in between the route ranges.

The NoR metric is a ratio of the solution’s number of routes over the total number of nodes. This is similar to unused capacity in that it measures excessive use of resources, but instead with focusing on the number of nodes in relation to number of routes.

(3) $$NoR=\frac{|R|}{n}$$where R is the set of routes and $n$ is the total number of nodes. This ratio is a measure of how much of a route per node is in the instance, and also allows for uniform comparison across all instances in a dataset. Dividing by the total number of nodes in an instance normalizes the metric and allows for comparison across all instances. The lower the number of routes the fewer resources (vehicles, etc.) that need to be reserved in the case of plan activation. A smaller value of this metric for an instance indicates that fewer routes are needed to serve the number of nodes.

For further evaluation, we will compare metrics across a set of characteristic attributes of the instance. These attributes of an instance include: depot location, node positioning, demand distribution, number of nodes in an instance, time constraint length, capacity constraint size, and the sum of all node demands. This analysis is in order to quantify whether there are significant differences of quality of the solutions produced when the graph has a specific attribute. These attributes are features of each dataset (location of nodes/depot, demand distribution, etc.) that can be binned across a continuum and compared across all instances.

The difference between the mean metric (UnCap, RngMax, NoR) for the best known and RSSD solutions is taken to indicate the compared quality of the solutions for that metric. When the difference is negative, the best known solutions have a lower value for that metric, and when the difference is positive the RSSD metric is lower. A difference of $0$ indicates the RSSD and best known solutions perform the same.

For the mean difference values, a trendline is computed to allow for quantitative analysis of potential patterns across the attributes continuum (increase in the attribute quantity or randomness). This trendline is a simple least squares fitting line. If the coefficient of determination value or $R^2$ is high, this indicates that the data fits the linear regression model. If $R^2$ value is high and the slope of the trendline is 0, this indicates that the RSSD and the best known set of solutions do not behave differently based on the continuum of the attribute.

To further compare the means, we perform statistical tests that can indicate if the means of the dependent variables (the attribute categories) and the independent variables (the mean difference for the three metrics) perform significantly differently. The One-Way ANOVA p-values were used when Levene’s test results were not significant; the Welch ANOVA is used when the result of Levene’s test were significant. If there is significance between the variables based on the p-values from the ANOVA tests, the Games-Howell posthoc test was performed. The Welch ANOVA and the Games-Howell posthoc test were chosen for their robustness against datasets that violate the homogeneity of varianc and equal sample sizes. With only 100 instances in the ‘X’ dataset and fewer than 20 instances in each attribute category, these results may be affected by the law of small numbers. These statistical tests were performed on SPSS Version 28 and the queries utilized and every instance’s results are published at www.doi.org/10.17605/OSF.IO/5H3GS.