1 Introduction

Healthcare facilities and professionals were, for the most part, caught unprepared for the magnitude of the COVID-19 pandemic. On March 23–27, 2020, the Office of Inspector General (U.S. Department of Health and Human Services) conduced brief telephone interviews (“pulse surveys”) with a random sample of administrators from 323 hospitals across 46 States, the District of Columbia, and Puerto Rico.Footnote 1 The survey key takeaway was:

Hospitals reported that their most significant challenges centered on testing and caring for patients with known or suspected COVID-19 and keeping staff safe. Hospitals also reported substantial challenges maintaining or expanding their facilities’ capacity to treat patients with COVID-19. Hospitals described specific challenges, mitigation strategies, and needs for assistance related to personal protective equipment (PPE), testing, staffing, supplies and durable equipment; maintaining or expanding facility capacity; and financial concerns.

One of the significant findings in this report was the limited supplies of personal protective equipment:

Hospitals reported that widespread shortages of PPE put staff and patients at risk. Hospitals reported that heavier use of PPE than normal was contributing to the shortage and that the lack of a robust supply chain was delaying or preventing them from restocking PPE needed to protect staff. Hospitals also expressed uncertainty about availability of PPE from Federal and State sources and noted sharp increases in prices for PPE from some vendors.

While this survey was conducted in the United States, hospitals in Europe had faced the same situation a few weeks earlier. The hospitals and other healthcare facilities in the city and province of Burgos (Spain) were not an exception as they were scrambling to deal with the eruption of the fast-moving pandemic. The shortage of PPE was quickly identified as critical in the COVID-19 fight. In response, various groups in the community, including the “Unidad de Cultura Científica” (UCC), a unit in charge of supporting and promoting scientific research at the University of Burgos, collaborated to set up the production of face shields. This resulted in a “call for help” to which many volunteers (the face shield “makers”) quickly responded. At the same time, agencies such as “Protección Civil” (similar to the Federal Emergency Management Agency in the United Sates), the national and local police, and the Red Cross volunteered vehicles and drivers to deliver materials to the makers and the face shields to hospitals, health centers, clinics, residences for the elderly, pharmacies, and small and medium-size stores.

After an initial attempt by the local government of Burgos to deal with the logistics associated with the production and delivery of face shields, UCC contacted us to develop a decision support system. This is the chronology of this development:

  • March 23, 2020: Initial contact and description of the problem

  • March 24, 2020: Start of heuristic development based on the original problem description

  • March 28, 2020: Data are made available for testing and integration

  • March 29, 2020: Meeting to discuss adequacy of the solutions generated by the system

  • March 30, 2020: Delivery and installation of the final routing system

This chronology highlights one of the main benefits of heuristic technology, namely, the ability to quickly produce a solution to a variety of problems, including those in logistics. Many academic publications in optimization focus on the solution speed. However, the speed that is reported and compared is the run time of the procedure that is being evaluated. What these publications lack is to mention the time that it took to formulate, develop, debug, test, and fine tune the procedures that found the solutions that are reported. While the overall time may not be important in many applications, it was a critical in the emergency that we describe here.

While vehicle routing has a long tradition in the operations research literature, most of the solution methods and applications addressed in the literature relate to commercial problems for which efficiency is associated with minimizing costs. However, there exists a healthy body of vehicle routing publications for disaster relief, emergency, and humanitarian logistics. We focus on some relevant examples in that subset of the vehicle routing literature.

Campbell et al. (2008) focus on the logistics associated with relief efforts in the aftermath of a large disaster. Specifically, their work is concerned with the routing of vehicles carrying critical supplies. Such routing can greatly impact the arrival times to those in need. They believe that it is critical that the deliveries are both fast and fair to those being served, and question whether the classic cost-minimizing routing problems properly reflect the relevant priorities in disaster relief. They focus on two alternative objective functions for the traveling salesman problem (TSP) and the vehicle routing problem (VRP): one that minimizes the maximum arrival time and one that minimizes the average arrival time. They examine the potential increase in routing costs that results from using these objectives. They present solution approaches for these two variants of the TSP and VRP, which are based on well-known insertion and local search techniques.

Vehicle routing in emergency situations such as disaster relief typically involve pickup and delivery points (Özdamar and Wei 2008; Wohlgemuth et al. 2012). It is also typical for supplies to be limited and for the demand for these supplies to be uncertain (Qin et al. 2017). Travel time is another potential source of uncertainty in disaster relief vehicle routing (Li and Chung 2019). Imposing travel time limits may be necessary in situations associated with healthcare activities. In the sample collection problem, vehicles visit a collection points to pick up thermal containers to be delivered to a laboratory. The length of the route (including travel and collection time) is restricted and the objective is to minimize total travel time (Grasas et al. 2014).

As we describe in detail below, the problem that we faced includes a combination of the features addressed in vehicle routing for emergency situations and disaster relief. In particular, the problem deals with multiple pickup and delivery points. Like (Campbell et al. 2008), routes must be fair to both volunteer drivers and those in need of supplies. The minmax objective function that we proposed to achieve this goal has been employed in formulations dealing with pickup and delivery problems arising in emergency situations (Xu et al. 2010; Yakici 2017). Narasimha et al. (2013) state that “this problem [minimization of the maximum distance traveled by any vehicle] is often of interest when minimization of time taken to visit all points is more important than the total distance traveled. The application includes emergency management situations where the objective is to use all available vehicles to minimize the time taken to attend to all points needing emergency resources.”

2 Problem description

The problem consists of routing vehicles to perform the following three activities:

  1. 1.

    Deliver material to face shield makers

  2. 2.

    Pick up face shields from makers

  3. 3.

    Deliver face shields to demand points

At the beginning of the day, vehicles are located in various origins. Routes could be closed or open, in the sense that vehicles may or may not return to their origin. In this context, the concept of a depot does not apply. That is, there is no central location with supplies to be delivered to demand points. Instead, routes must be configured in such a way that vehicles visit pickup points first and delivery points second.

As mentioned earlier, vehicle drivers are volunteers, who coordinate their efforts to meet the demand for face shields. As such, they expect an even distribution of work and time commitment. This leads to the need for configuring routes whose length is not too different from one to another.

The VRP literature includes objective functions designed to produce solutions with balanced route lengths. In response to complaints by parents whose kids spent too much time in their school buses, Pacheco et al. (2013) developed a solution procedure for a bi-objective optimization model that trades off service (i.e., the minimization of the longest route) and operational cost (i.e., the minimization of the total distance traveled). The minimization of the longest route is one way of searching for solutions with balance routes. Another possible objective consists of minimizing the difference between the longest and the shortest route (Jozefowiez et al. 2009). This objective function is similar to the minimization of the longest route in the sense that both seek length balancing. However, the objectives are not the same, as illustrated by the following example (Pacheco et al. 2013). Consider the routes depicted in Fig. 1. The circles represent pickup locations (A, B, and C) and the square is a demand point (D). The numbers between parentheses associated with the arcs represent distances.

Fig. 1
figure 1

Two alternative routes for three pickup points (A, B, and C) and one demand point (D)

The total distance of the longest route in both solutions in Fig. 1 is three. Solution 2 is better than solution 1 in terms of the difference between the longest and the shortest route. This difference has a value of two for solution 1 and a value of zero for solution 2. However, in terms of total distance traveled, solution 1, with a total distance of four, is better than solution 2, with a total distance of six.

In our context. We considered that minimizing the longest route would result in an adequate workload fairness, and it would also afford the possibility of cost savings from the reduction in the total distance traveled.

In order to formalize the description of the problem and solution approach, we define the following terms:

  • \( m \): Number of vehicles

  • \( P \): Set of \( n_{1} \) pickup points, such that \( P = \left\{ {1,2, \ldots ,n_{1} } \right\} \)

  • \( Q \): Set of \( n_{2} \) delivery points, such that \( Q = \left\{ {n_{1} + 1,n_{1} + 2, \ldots ,n_{t} } \right\} \) and \( n_{t} = n_{1} + n_{2} \)

  • \( O \): Set of origins, such that \( O = \left\{ {n_{t} + 1,n_{t} + 2, \ldots ,n_{t} + m} \right\} \)

  • \( E \): Set of end-of-route points, such that \( E = \left\{ {n_{t} + m + 1,n_{t} + m + 2, \ldots ,n_{t} + 2m} \right\} \)

  • \( t_{ij} \): Travel time between \( i \) and \( j \)

  • dij: Distance between \( i \) and \( j \)

  • \( p_{i} \): Quantity picked up at point \( i \in P \)

  • \( q_{i} \): Quantity delivered at point \( i \in Q \)

A solution is described as follows:

  • \( R \): Solution (i.e., set of feasible routes)

  • \( \pi_{r} \left( i \right) \) Index of the point in position \( i \) in route \( r \)

  • \( n_{r} \): Number of points in route \( r \)

  • \( T\left( r \right) \): Time to complete route \( r \), i.e., \( T\left( r \right) = \sum\nolimits_{i = 1}^{{i = n_{r} - 1}} {t_{{\pi_{r} \left( i \right),\pi_{r} \left( {i + 1} \right)}} } \)

  • \( D\left( r \right) \): Total distance of route \( r \), i.e., \( D\left( r \right) = \sum\nolimits_{i = 1}^{{i = n_{r} - 1}} {d_{{\pi_{r} \left( i \right),\pi_{r} \left( {i + 1} \right)}} } \)

For a route \( r \in R \), \( \pi_{r} \left( 1 \right) = n_{t} + r \), \( \pi_{r} \left( {n_{r} } \right) = n_{t} + m + r \).

The problem can be stated as follows:

$$ {\text{Minimize}}\quad \quad T^{max} \left( R \right) = \mathop {\hbox{max} }\limits_{r \in R} \left\{ {T\left( r \right)} \right\} $$
(1)
$$ {\text{Minimize}}\quad \quad D^{sum} \left( R \right) = \mathop \sum \limits_{r \in R} D\left( r \right) $$
(2)

Subject to

$$ \mathop \sum \limits_{i = 3}^{{i = n_{r} - 1}} q_{{\pi_{r} \left( i \right)}} \le \mathop \sum \limits_{i = 2}^{{i = n_{r} - 1}} p_{{\pi_{r} \left( i \right)}} $$
(3)
$$ i < j\forall \pi_{r} \left( i \right) \in P\,\,{\text{and}}\,\,\pi_{r} \left( j \right) \in Q $$
(4)

The objective functions are treated hierarchically. The first objective minimizes the time of the longest route. The second objective minimizes the total distance traveled. As we describe below, the second objective is only used to break ties when searching for the optimal solution for the first objective. Therefore, this formulation should not be viewed as a biobjective problem in which the goal is to find a frontier of efficient solutions. Constraints (3) ensure that the amount of supplies (i.e., materials for the makers or face shields for the demand points) that are picked up in the route are at least as large as those that are delivered. The constraint assumes that the first position in route \( r \) is the origin of vehicle \( r \)\( \left( {{\text{i}}.{\text{e}}.,\, n_{t} + r} \right) \) and that the second position is a pickup point \( ({\text{i}}.{\text{e}}.,\,\pi_{r} \left( 2 \right) \in P) \). Constraints (4) ensure that pickup points are visited before delivery points. The capacity of the vehicles is not relevant because it can be assumed that vehicles can handle any pickup or delivery quantity.

This is a summary of the characteristics of the problem that we tackled:

  1. 1.

    Vehicles and drivers belong to several organizations and therefore their starting and ending points may be different.

  2. 2.

    Routes may be closed or open, that is, the ending point may be different from the starting point.

  3. 3.

    Routes include pickup and delivery of face shields. As opposed to some commercial applications, the items that are picked up do not go back to a depot, instead, they are delivered to demand points. Therefore, quantities to be picked up must be at least as large as those that are demanded by the delivery points.

  4. 4.

    Routes include the delivery of raw materials (to make the face shields) in addition to the delivery of the face shields that are being picked up in the route.

The solution approach that we developed takes into considerations these characteristics in order to produce feasible routes that are acceptable to the volunteer organizations.

3 Solution approach

Our initial approach was to adapt a commercial vehicle routing software. The adaptation consisted of a redefinition of the distance matrix to force the procedure to visit pickup locations before delivery locations. In the redefined matrix, all the distances from delivery points to pickup points were made a large number. These large numbers guaranteed that no optimal solution would include a link from a delivery point to a pickup point. All routes would have to visit pickup points first before transitioning to delivery points. This simple adaptation produced routes that, in principal, were feasible solutions to the problem. However, those in charge of coordinating the delivery efforts, quickly pointed out the issues of implementing the highly unbalanced routes that the system was producing. Also, some solutions violated the constraint of picking up enough face shields to cover the demand at the delivery points. We realized that we needed a search procedure that could produce feasible routes for which the primary objective was minimizing the longest route and with a secondary objective of minimizing the total distance traveled. The system that we had for commercial deliveries was not equipped to deal with a problem with such constraints and objectives. Therefore, we took elements from solution procedures for routing problems that we have tackled in the past and blended them in a multistart framework that constructs and improves solutions, as shown in Algorithm 1.

The multistart procedure runs until the number of iterations without improvement reaches \( maxiter \). The best solution starts as an empty set of routes. The objective function associated with this solution is a large number for both the longest route and the total distance traveled. The do-loop starts with the construction of a solution by a pseudo-greedy insertion procedure (line 2). A short-term memory tabu search is then executed with the goal of finding an improved solution (line 3). The solution found after the tabu search is compared with the best solution found so far (line 4). This comparison is such that \( R \) is better than \( R^{best} \) if:

$$ T^{max} \left( R \right) < T^{max} \left( {R^{best} } \right)\quad {\text{or}}\quad T^{max} \left( R \right) = T^{max} \left( {R^{best} } \right)\quad {\text{and}}\quad D^{sum} \left( R \right) < D^{sum} \left( {R^{best} } \right) $$

If \( R \) is better than \( R^{best} \), then the best solution is updated in line 5. The procedure continues to iterate until the termination criterion is satisfied.

figure a

The solution construction consists of an insertion procedure that balances quality and diversity. The value of a parameter \( \alpha \) controls the degree of randomization that allows the construction process to generate diverse solutions. It is well known that solutions that are generated totally at random achieve a high degree of diversity at the expense of poor quality. It is also known that improving methods require of a large amount of computational effort to find high-quality solutions when starting from random solutions of extremely poor quality (Campos et al. 2001). Pseudo-greedy construction of the type that we chose to use have been shown effective, particularly in the context of greedy randomized search procedures (Feo and Resende 1995). These constructions produce diverse solution of reasonable quality.

Algorithm 2 shows the pseudo-code for our solution constructions based on insertions. Due to the nature of our routing problem, we know that a solution must use all available vehicles and that their routes should start at the specified set of origins \( O \) and end at the specified end-of-route points \( E \). Therefore, all routes are initialized as follows:

$$ \pi_{r} \left( 1 \right) = n_{t} + r\quad {\text{and}}\quad \pi_{r} \left( 2 \right) = n_{t} + m + r\quad {\text{for}}\,{\text{all}}\quad r = 1, \ldots ,m $$

The initialization of the routes is the first step of the procedure (line 1 in Algorithm 1). The length of each route is the time to travel from the origin to the endpoint of the route (line 2). The construction procedure uses the set \( U \) to keep track of the points that must be visited but that so far have not been inserted in any route. At the beginning, all pickup and delivery points are in \( U \) (line 3).

In line 4, the procedure builds a list \( L \) of all feasible insertions. An insertion is feasible if the resulting route does not violate the problem constraints. The candidate list \( CL \) is the subset of insertions in \( L \) that do not increase the length of the longest route (line 5). If all insertions increase the length of the longest route, then the candidate list \( L \) is the set of all feasible insertions (line 6). A restricted candidate list \( RCL \) is built with the top \( \alpha \)% insertions. If \( \alpha \) is zero, then only the best insertion is considered, making the procedure totally deterministic. If \( \alpha \) is one, then all insertions in \( CL \) are considered, and thus, allowing for maximum randomization. The best insertion is the one that causes the smallest increase in the length of the route where the insertion occurs. The insertion to be executed next is randomly selected from \( RCL \) (line 8). Once the selected insertion is executed (line 9), the solution is updated (\( R \) and \( T^{max} \)), and the point inserted in the solution is removed from \( U \) (line 10). The procedure finishes when all points have been inserted into the solution (i.e., when \( U \ne \emptyset \)).

figure b

Solutions constructed with the insertion procedure are subjected to an improving process consisting of a short-term memory tabu search (Glover 1989, 1990; Glover and Laguna 1997). The search starts from the solution generated by the insertion procedure. The best move is executed at each iteration, where move quality is determined in the same way as in Algorithm 1. That is, moves are first judged by the amount in which they are able to reduce \( T^{max} \).The secondary objective function \( D^{sum} \) is used as the tiebreaker when there is more than one move that improves upon the best solution found during the current tabu search. However, if two moves result in the same \( T^{max} \) but the corresponding solutions are not better than the current best, then we use a different criterion to break ties. This criterion favors moves that result in the least increase of the maximum length of the routes involved in the move.

Our procedure uses four types of moves: Or Exchange, Generalized Or Exchange, CROSS, and Insertion. We chose these moves because we have verified their merit when dealing with a bi-objective problem involving \( T^{max} \) and \( D^{sum} \) (Pacheco et al. 2013). Figure 2 shows a graphical representation of the moves that we embedded in the tabu search. The Or Exchanges happen in a single route. The CROSS and Insertion moves involve two routes. In all cases, the moves consist of deleting and adding arcs. For example, in the OR Exchange in Fig. 2, the arcs (A,B), (C,D), and (F,G) are deleted (dashed lines), while (A,D), (F,B), and (C,G) are added. This has the effect of exchanging two consecutive segments, that is, the one from B to C with the one from D to F. The Generalized Or Exchange deletes and adds four arcs and is used to exchange two non-consecutive segments. In Fig. 2, the Generalized Or exchanges segment B to C with segment F to G. CROSS exchanges segments between two routes while Insertion moves a segment from one route to another. In all cases, we only perform feasible moves. That is, the neighborhood \( N\left( R \right) \) consists of all the feasible solutions that can be reached from solution \( R \) by applying the four moves.

Fig. 2
figure 2

Moves within the tabu search

As customary in tabu search, the procedure can make non-improving moves. In tabu search, short-term memory is used to prevent cycling caused by move reversals. Our short-term memory structure is such that it records the iteration number at which an arc was deleted from the solution. Deleted arcs are declared tabu-active for several iterations (given by the value of \( tenure \)). A move is tabu if it contains a tabu-active arc. The tabu status of a move is overridden if the move results in a solution that is better than the best found during the current call to the tabu search procedure. This is what is known in the tabu search methodology as the aspiration criterion. The solutions that can be reached with non-tabu moves or that meet the aspiration criterion form the reduced neighborhood \( N^{*} \left( R \right) \subset N\left( R \right) \). Algorithm 3 shows the pseudo-code of our tabu search.

figure c

The search starts with the initialization of the tabu structure (line 1). This initialization consists of assigning a value of zero to all the entries of a matrix that records the iteration at which an arc becomes tabu. Each entry of the matrix corresponds to an arc. The best solution (\( R^{*} \)) at the beginning of the search is the initial solution \( R \) constructed with the insertion procedure (line 2). Note that \( R^{*} \) denotes the best solution found during the current call to \( TabuSearch \) with the starting solution \( R \). This is different from \( R^{best} \), which represents the best overall solution (as shown in Algorithm 1). The next solution to be visited is the best solution in the reduced neighborhood of the current solution (line 3). The best neighbor is chosen with the criteria described above. After the updating of the memory structure (line 4), the current solution is compared to the best solution found in the current tabu search (line 5). The current best is updated in line 6. The search stops when \( maxts \) iterations are performed without improving \( R^{*} \).

The performance of the multistart procedure depends on four parameters:

  • \( maxiter \): Maximum number of starts without improvement of \( R^{best} \)

  • \( maxts \): Maximum number of tabu search iterations without improvement of \( R^{*} \)

  • \( \alpha \): Degree of randomization

  • \( tenure \): Short-term memory tabu tenure

Given a fixed amount of time to find a solution to a problem, it is important to determine how to distribute the computational budget. Unfortunately, there is no general rule that prescribes the proportion of time that should be spent constructing solutions versus the time spent trying to improve them. Typically, searching a neighborhood to find the best move to make requires more effort than construing a solution. Therefore, in most applications the number of moves (that is, iterations within tabu search) is much larger than the number of steps in the construction phase.

4 Delivery of face shields

As we mentioned in the introduction, we engaged in this project because the city and province of Burgos (Spain) needed a system to organize the delivery of face shields produced by volunteers. After the call for help, 236 individuals and small companies volunteered to make face shields. From March 21, 2020 to June 3, 2020, these volunteers made a total of 37,311 face shields. Table 1 shows the distribution of delivery points throughout the province of Burgos.

Table 1 Delivery points in the province of Burgos

The routing problem was solved daily for 51 days with an average number of pickup points of 30 and an average number of delivery points of 65. There was a total of eleven vehicles available. On average, eight vehicles were used each day. Figure 3 shows the level of delivery activity from the third week of March to the first week of June.

Fig. 3
figure 3

Number of face shields delivered in the province of Burgos

We used cartographic data for the province of Burgos to calculate point to point travel times and distances. We also adapted a geographical information system to show input and output data, such as the routes proposed by our procedure (see Fig. 4).

Fig. 4
figure 4

Screen capture of our GIS with one of the routes on April 2, 2020

One of the reasons for developing a customized heuristic for this routing problem was that our adaptation of a commercial code that minimizes \( D^{sum} \) resulted in unbalanced routes. Consider the data associated with the delivery of face shields on April 2, 2020 in the city of Burgos. On this day, there were 34 points served by 4 vehicles. Table 2 shows the comparison between the solution found when the objective is to minimize the total distance traveled (i.e., \( D^{sum} \)) and the approach of minimizing the time of the longest route (i.e., \( T^{max} \)) as the primary objective and the total distance traveled as the secondary objective.

Table 2 Solution comparison for the delivery on April 2, 2020

The routes found with the distance minimization approach are unbalanced. The difference between the longest route and the shortest route is 252.4–16.4 = 236 min (i.e., almost 4 h). The minmax approach generates balanced routes with a 27% increase in total distance. There was a strong preference for the minmax solution in terms of improved coordination and workload fairness.

During the first 2 weeks of operation, the routes were configured manually. During this time, for the entire province, the number of points varied from 80 to 90 and the number of vehicles from 9 to 11. The manual routing consisted of dividing the points by zones. These manual solutions resulted in routes that often could not be completed because they did not consider speed limits and street directions. In addition, the manual routing required time from the volunteer drivers, who did not have the knowledge or the tools to be effective route planners. We used our optimization system to generate routes for the 2 weeks of manual operation. Comparing our routes with the manual routes, we estimate that the optimization system was able to save an average of 1 h per route. This is significant, considering that a total of 420 routes were completed during the time that the system was used. In addition, all routes were feasible, and drivers were able to complete them within the estimated time.

5 Conclusions

The real-world application of heuristic optimization that we described here illustrates the practical usefulness of this area of operations research. In a remarkable short period of time we were able to capture the key features of a problem and deliver implementable solutions with a heuristic approach. The optimization system that we developed has the following advantages over the manual approach that was initially used:

  • Significant time savings in the planning of routes

  • Accurate estimated time to complete routes

  • Detailed maps of the balanced routes

The routing system resulted in a total deliver of face shields in the province of Burgos that was an order of magnitude larger than in the rest of Spain. During the period of late March to early June of 2020, the number of face shields relative to the population of the province of Burgos was 10.5%, while for the rest of Spain was 1.46%. The success of this project has been the topic of articles in the local and regional press (Andrés 2020; Antolín 2020; García 2020). In addition, the project won the “Innovation against the Virus” prize given by the Ministry of Education of the Autonomous Community of Castilla y León (Blanco 2020).

Finally, we want to point out that many researchers around the world have applied their skills and knowledge to support important efforts related to the COVID-19 pandemic. For instance, Duque et al. (2020) formulate and solve a stochastic optimization problem that selects daily hospitalization triggers and recommends reinstatement and relaxation of lock-down periods. Also, Mathematics Against CoronavirusFootnote 2 is a community of researchers from areas such as statistics, computer science, and operations research that have developed models ranging from system dynamics, prediction, and simulation to optimization. As stated in the INFORMS website,Footnote 3 “operations research and analytics professionals are providing valuable knowledge for policymakers and business leaders responding to the impacts of the coronavirus outbreak.” We feel honored to belong to a group that has made a modest positive contribution to meet the challenges of this world-wide crisis.