A solution approach for multi-trip vehicle routing problems with time windows, fleet sizing, and depot location

We present a solution approach for a multi-trip vehicle routing problem with time windows in which the locations of a prescribed number of depots and the fleet sizes must also be optimized. Given the complexity of the task, we divide the problem into subproblems that are solved sequentially. First, we address strategic decisions, which are solved once and remain constant thereafter. Depots are allocated by solving a p -median problem and fleet sizes are determined by identifying the vehicle requirements of several worst-case demand instances. Then, we address the operational planning aspect: optimizing the vehicle routes on a daily basis to satisfy the fluctuating customer demand. We assign customers to depots based on distance and “routing effort,” and for the routing problem we combine a tailor-made branch-and-cut algorithm with a heuristic consisting of a route construction phase and packing of routes into vehicle trips. Our strategic decision models are robust in the sense that when applied to unseen data, all customers could be visited with the allocated fleet sizes and depot locations. Our operational routing methods are both time and cost-effective. The exact method yields acceptable optimality gaps in 20 min and the heuristic runs in less than 2 min, finding optimal or near-optimal solutions for small instances. Finally, we explore the trade-off between depot and fleet costs, and routing costs to make recommendations on the optimal number of depots. Our solution approach was entered into the 12 th AIMMS-MOPTA Optimization Modeling Competition and was awarded the first prize.


INTRODUCTION
The vehicle routing problem (VRP) is among the most studied problems in the field of combinatorial optimization.It is a generalization of the well-known traveling salesman problem (TSP) where instead of a single salesman there is a fixed number of identical vehicles leaving from and returning to a given depot.The VRP has many practical applications and it can be tailored to model a wide variety of business-specific needs.Many of these special cases can be found in the literature, and a variety of exact methods and heuristic approaches have been developed to obtain good solutions in a reasonable amount of time [17].Just like the TSP, the VRP is an intriguing problem in that it is simple to explain and yet extremely challenging to solve.
Formally, the VRP can be defined as follows: let G = (V, E) be a connected, undirected graph, where V is the set of nodes and E is the set of edges.Each node represents a customer and each edge is associated with the traveling time (or distance) between two customers.There is a fleet of |K| identical vehicles of capacity Q.Each customer i ∈ V has a fixed demand of q i , and all of its demand must be satisfied by a single vehicle.The objective is to construct a single route for each vehicle such that the total distance traveled is minimized.This is commonly known as the capacitated VRP (CVRP), which is the basis of all further variations of the problem.A more detailed description of it can be found in [22] and a good overview of the work that has been done in the past years on the VRP is given in a survey paper by Laporte [17].
The modeling problem proposed for the 12 th AIMMS-MOPTA Optimization Modeling Competition [1] was a variant of the CVRP with some additional elements.(i) It is required that every customer be visited within a given time window (8:00 a.m. and 4:00 p.m.) and, additionally, vehicles can only operate within fixed time intervals (from 6:00 a.m. to 5:00 p.m.).This variant can be found in the literature as the VRP with Time Windows (VRPTW) [15].(ii) In this problem, there is not a single depot; a fixed number of depots is given as a problem parameter; the locations of these depots must be determined.These are also well studied problems and they are referred to in the literature as the multi-depot VRP and the location routing problem, respectively.Laporte et al. describe an algorithm for solving a combination of these two problems efficiently in [18].(iii) The size of the fleet assigned to each depot must be minimized and this assignment must remain constant for the entire duration of the problem.A problem of this type was first introduced by Golden et al. in [13] and significant amount of work has been done on this topic since then.Baldaci et al. introduce in [3] a new mathematical formulation for a problem of this type with a heterogeneous fleet-a generalization of our fleet sizing problem, as we only need to consider a homogeneous fleet.Lastly, (iv) vehicles are allowed to perform multiple trips, starting and terminating at the same depot.This is often referred to as the Multi-Trip VRP (MTVRP) [6].We use the term trip in the context of the multi-trip aspect of the problem: it refers to an occasion where a given vehicle leaves the depot to visit customers and return.By contrast, we use the term route to refer to the actual path that a vehicle follows during a trip.
All of the problems listed above are known to be computationally challenging, and heuristic methods have been developed as a consequence of the limitations of the exact approaches.
The data provided for this problem as part of the 12th AIMMS-MOPTA Optimization Modeling Competition [1] includes a graph of 1457 nodes representing ZIP codes in the state of Pennsylvania, USA.The given graph is not complete, but it is connected, so a path can be found between any two nodes using existing edges.
Demand of each customer is provided for every day in 2018 and 2019.Despite the large number of customers, on average only 5% must be visited on a given day.The 2018 data is used for model development while the 2019 data is reserved for testing the robustness of our solution method.In practice, demand for a given day is not known until 6:00 p.m. the day before, so all routes should be obtained quickly enough to inform the drivers on time.
We propose a method that efficiently determines the location of a given number of depots, optimizes the fleet size and produces a sensible routing strategy in a reasonable time.
This article is structured as follows: in Section 2 we present a mathematical formulation of the problem.Our methodology is outlined in Section 3. Sections 4 and 5 present an analysis of the minimum number of depots and the depot location model, respectively.In Section 6, a robust approach to the fleet sizing problem is proposed.Section 7 describes an algorithm to assign customers to depots on a daily basis.Both exact and heuristic approaches for the routing problem are described in Section 8. Solutions are presented in Section 9. We present a methodology for automating the choice of the number of depots in Section 10 and lastly, conclusions and recommendations are given in Section 11.

MATHEMATICAL FORMULATION
In this section, we present a mixed integer linear programming formulation for the multi-trip VRP with time windows, fleet sizing and depot location.Some assumptions were made at different stages of the problem: • Each customer must be served by exactly one vehicle.
• The demand of customers colocated with a depot is assumed to be served automatically.
• The distance between two customers is the shortest path between their locations through known routes.The Dijkstra algorithm [9] was applied to generate a distance matrix with all the shortest paths between any two nodes in the network.• No time is needed to unload the demand at a customer location or to reload at the depot.
We define an undirected graph G = (V, E) where the set of nodes V represents the customers and the set of edges E represents the shortest paths between customers.We also define a set of unknown depot locations, a set of potential vehicles to be used and a set of trips each vehicle might take.Formally, we introduce the following mathematical formulation: Sets:

Set of trips per vehicle
Parameters: Depot opening hours Decision variables: there is a depot at node i, 0 otherwise.Minimize subject to The objective function minimizes three types of costs: the variable routing costs (fuel, vehicle maintenance, tolls, etc.), the fixed depot costs (depot lease/purchase and personnel wages), and the fixed vehicle costs (lease/purchase, insurance, driver wages, etc.).
Constraints ( 2) and (3) ensure that every node i ∈ V is visited exactly once if it is not a depot and at least once otherwise.Constraints (4) stipulate that, if trip r of vehicle k enters node i, it must leave node i.We introduce constraints (5) to ensure that every vehicle leaves a depot at most once in each trip; they are necessary to avoid sub-trips within the same trip that could deceive the time window constraints by taking place simultaneously.The vehicle capacity constraints are defined in (6).If vertex i is the depot of vehicle k, then the demand of this node is excluded from the load.Constraints (7) ensure that exactly |H| depots are placed.Constraints ( 8)-( 9) assign vehicles to at most one depot.Constraints (10) ensure that if a vehicle is not used, this vehicle cannot traverse any arcs.Constraints (11) impose that a vehicle that does not complete any trip cannot be marked as used by variables w.This is needed to prevent negative demand on a node when evaluating the capacity constraints (6), which could happen if it is not costly to mark a vehicle as used.
We introduce constraints (12) to ensure that every vehicle in use always departs from its assigned depot by forbidding it to enter other depots.We also introduce symmetry breaking constraints (13) and (14) to impose a natural ordering in the selected vehicles and trips, respectively, and eliminate equivalent solutions.They ensure that vehicle k + 1 cannot be used unless k is in use; the same is true for trips r + 1 and r.
Constraints ( 15)-( 21) are our adaptation of the time window constraints for the multi-depot multi-trip VRP.(15) are used when none of the nodes considered are a depot.(16) only applies when the first node is a depot and (17), when the last node is a depot.Additional variables h kr 1 , h kr 2 are needed to define visit times to the depot, as there are two times associated with the depot for every vehicle and trip: start and end time of a trip.Because of this, we do not need the vehicle index nor the trip index in s i , as nodes that are not depots can only be visited once by one vehicle.Constraints (18) ensure that a trip r cannot start before the finish time of the previous trip completed by the same vehicle.Finally, Constraints ( 19)-( 21) enforce the time windows of customers and vehicles.Subtour elimination constraints are not necessary due to the presence of the time window constraints [15].
Using the formulation above, we can prove the complexity of the problem.By fixing the number of trips to one, the index r of any decision variables can be dropped.We can fix the depot at node 0 by forcing z 0 = 1, |H| = 1, and setting a d = a c and b d = b c .With these changes, we notice that by removing the resulting redundant constraints, we obtain the VRPTW formulation.This shows that our problem is a generalization of the VRPTW.Since the VRPTW is itself a generalization of the TSP, which belongs to the class of NP-hard problems, we conclude that our problem is also NP-hard -there is no known algorithm that can solve this problem in polynomial time.

METHODOLOGY
In this problem, there is a hierarchy of decisions that take place at different stages in the planning period: 1. Strategic level.Decisions involving capital expenditure (capex) that need to be made in advance and remain constant throughout the entire planning horizon; these are the location of depots and the fleet sizing.2. Operational level.Decisions that are made on a daily basis and have an impact on operational expenditure (opex); these are the vehicle routes.Since customer demand is only known at 6:00 p.m. the day before, these routes should be generated fast enough so that schedules can be communicated to drivers before they start their shift at 6:00 a.m.
Because of this hierarchy, there is a business need to decompose the problem into separate subproblems.We have developed a solution approach that addresses the different stages of the problem sequentially.First, the location of depots is optimized.Next, fleet sizes are determined for each depot.After these strategic decisions have been fixed, we address the operational planning aspect: on a given day, customers that require a visit are assigned to depots and vehicle routes are optimized.
All our solution methods are implemented in AIMMS, with CPLEX as the MIP solver.

MINIMUM NUMBER OF DEPOTS
Since the number of depots is a parameter and not a decision variable, we are interested in identifying the minimum number of depots that preserves feasibility.This number can be determined by solving a set covering problem (SCP).The SCP aims to place facilities at locations such that every node can be reached from at least one facility.Each depot can be considered to have a radius of coverage equal to the longest distance that a vehicle could cover in a single trip without violating the time windows.
In our problem, the radius of coverage is d max = 5.5 h × 40 mph = 220 miles.A vehicle could visit a customer located at most 5.5 h away from the depot, although this would likely result in a large fleet size with a poor vehicle utilization.Given the data, the solution to the SCP is trivial, as there are 10 ZIP codes such that the location of all customers is less than 5.5 h away from any of these 10 locations.Therefore, one depot suffices to satisfy all the demand in 2018, provided there are enough vehicles to complete the generated trips.

DEPOT LOCATION
To optimize the location of a given number of depots without solving the routings for every day in the 2018 data, which would require very costly computations, we treat this as a facility location problem.We are interested in finding the best locations of a given number of depots such that if all customers are assigned to their closest depot, the weighted average distance between customers and their assigned depots will be minimum.The rationale is that minimizing this metric should result in shorter vehicle trips.
To obtain the weights for the average distance, we aggregate demand in 2018 by the number of days that each customer has positive demand.The more often a customer requires a visit, the more it will drive the location of the nearest depot toward it.We observe a strong correlation between this metric and the total annual demand of a customer in our data (the Pearson coefficient between the two metrics is  = 0.998).Therefore, choosing one metric or the other is not expected to have a significant impact on the outcome of the model.
In addition, the fixed costs of depots are also included in the objective function.Since we cannot estimate routing and vehicle costs from the average customer-depot distances, it is not straightforward to find the right balance between the depot costs and the average customer-depot distance.To overcome this difficulty, we include both terms in the objective function with a parameter  ∈ [0, 1] that controls their relative importance and we choose the best value for  by inspection.
Our formulation is a variation of the p-median problem [8].Variables z i indicate whether a depot is placed at location i ∈ V and variables y ij indicate whether node i is assigned to a depot located in j. d ij is the distance between nodes i and j, d max is the maximum distance allowed between a demand node and a depot due to the time window constraints, and Ω i is the number of days that customer i ∈ V had positive demand in 2018, namely Constraints (25) ensure each node is assigned to exactly one depot.( 26) is the p-median constraint, where the p parameter is represented by the cardinality of the set of depots, |H|.Constraints (27) stipulate that nodes can only be assigned to depots, not to customer nodes.The main difference with the p-median problem lies in constraints (28), which are needed to ensure that only nodes that can be reached from a depot can be assigned to that depot.
This p-median problem, with circa 1500 nodes or ZIP codes, is too large to solve to optimality with exact methods.Our approach consists of solving the p-median problem in two stages, considering different resolutions: • Stage 1: Group nodes that are close together into clusters, forming at most 500 clusters.Then solve the p-median problem using these clusters instead of the original nodes to identify |H| clusters that contain suitable depot locations.Since each of these |H| clusters will generally include multiple nodes (multiple ZIP codes), Stage 2 assigns a depot to a single node in each of the clusters.• Stage 2: Consider the p-median problem with the original nodes.Fix z j = 0 for all nodes j ∈ V that do not belong to the clusters identified in Stage 1.This ensures that only nodes that are in a cluster are candidate depot locations.Finally, solve the resulting p-median problem to obtain the final set of depot locations H.
Both steps run in less than a minute, while introducing only a small distortion to the spatial distribution of demand.To perform the clustering described in Stage 1, we apply k-medoids [16], an unsupervised machine learning technique similar to the well-known k-means algorithm [20].For a given number of clusters k, k-medoids minimizes the dissimilarities between points in a cluster and the center of that cluster.It can be easily implemented if the input data is a distance matrix.The code to perform the clustering of nodes was written in Python using the library PyClustering [21].
Figure 1 illustrates the trade-off between depot costs and average distance, for a range of values of  ∈ [0, 1].It can be observed that as we prioritize minimizing the total fixed cost of depots ( → 0), the average distance depot-customer grows very rapidly.Similarly, if only the average distance metric is considered ( = 1), then solutions with very high depot costs are obtained.However, by reducing  to 0.9, namely, by adding a small penalty for high depot costs, these costs are reduced between 20% and 40% while the increase in the average distance is negligible.This suggests that it is easy to find a cheaper alternative in the neighborhood of a location that is optimal under spatial considerations only.
Unlike routing costs, which are 100% opex, the fixed cost of the depots has a capex element (lease or purchase payments).This means that part of the depot costs will be amortized within a few years, whereas routing costs are running costs.Therefore, it seems sensible to invest in finding good locations for the depots at the expense of purchasing or leasing slightly more expensive facilities.For these reasons, we set  to 0.9.Using 2018 data, this scheme returned a depot configuration in 40 s on average, for different numbers of depots.
It should be noted that balancing of demand across different depots has not been included in this model.Having some variability in the volume of demand assigned to each depot is not considered to be problematic, since it is possible to assign fleets of different size to different depots.

FLEET SIZING
The next strategic decision concerns the sizing of the fleet and the distribution of this fleet across different depots.This aspect of the problem is key to obtaining robust solutions.While a poor choice of depot locations or number of depots can be compensated with a larger fleet of vehicles, the fleet sizing is critical to ensuring all demand can be routed: too small a fleet will result in some customers not being visited on certain days unless additional resources are put in place.
To obtain a fleet size that can cope with 2018 demand but also be robust against unseen demand data (e.g., 2019 demand), it is not sufficient to set the fleet size to the number of vehicles required on the busiest day in 2018, because the fleet size required depends not only on the total demand but also on the spatial distribution of the customers.To overcome this, we propose to simulate a set of demand instances based on a busy period in 2018.We then determine the number of vehicles required on those simulated days by solving a VRP on each of these instances (using the methodology outlined in Section 8) and decide the fleet size based on these results.
Simulating "busy days" first requires an understanding of the demand distribution in 2018.Let q(i, d) be the demand of customer i in day d.We assume that q(i, d) is given by the product of two random variables X(i, d) and Y(i, d).X(i, d) is a binary variable indicating whether customer i is visited on day d and Y(i, d) indicates how many demand units customer i must be served on day d.The rationale behind this assumption is that we observed that, on average, only 5% of the customers had a positive demand (q(i, j) > 0) on any given day in 2018, and in those cases where q(i, j) > 0, this demand was variable and seemed to follow a distribution centered in 10.
Through statistical analysis we identified the following properties of X and Y that can be exploited to simulate instances that resemble the actual demand in 2018:

The amount of goods required by an active customer node, Y(i, d), is independent of the customer and the day. Y(i, d)
can be replaced with Y.

The amount of goods required by an active customer node, Y, follows a
Poisson distribution with  = 9.95.The parameter  was fitted using the maximum likelihood method.The goodness of fit was tested using a chi-squared test, and it can also be visualized in Figure 2A.

X(i, d) is a Bernouilli distributed random variable that takes value 1 if customer i requires a visit on day d, which happens
with probability p id .This probability is different for every customer and day (see next point).

The seasonality behind the demand in 2018 is driven by X(i, d).
The R 2 score of a linear regression model linking the number of active customers with the total demand on a given day is 0.986, which indicates that 98.6% of the variability in the total demand is explained by the number of active customers.Busier days are so because there are more customers that have to be visited, and not because the demand per active customer is higher.The similarity between these two time series can be seen in Figure 2B.

5.
With the exception of a few days in January, the busiest period of the year is the summer season (see Figure 2B).Also, a logistic regression model showed that customers are 20.3% less likely to have demand on a weekend.6.On any given day, the demands of different customers in terms of X(i, d) are independent from each other.The Pearson correlation coefficients  that resulted from comparing the demands of all pairs of customers are very low (|| < 0.15 for 95% of all pairs of customers and max || = 0.26).As a result, we can expect that the spatial distribution of demand will be similar on different days.In other words, the relative volume of customers in each depot area will not vary significantly on different days.This is a convenient property from a fleet sizing point of view, given that vehicles must operate from the same depot every day.
Based on these observations, the following approach is used to simulate a "busy day": 1. Sample from X: each customer i is activated with a probability p i, busyday .This probability is initially estimated as the frequency with which customer i required a visit over all weekdays from June to August in 2018.In other words, we sample data from an artificial "busy day" that behaves like an average summer weekday in terms of demand volume, and the probability p i, busyday will only depend on the customer node i.Furthermore, to achieve more robust solutions, customers that were not active during this period (and therefore would have zero chance of appearing in our simulated data) are assigned a probability p equal to the lowest nonzero probability in the data.2. Sample from Y: for each customer that has been activated in the previous step, the amount of demand is sampled from Poisson(9.95).3. Calibrate demand level: to obtain robust solutions that can ensure not only that all demand is satisfied in 2018, but also cope with a possible increase in demand in the future, we introduce a calibration factor  > 1 that increases the generated demand by multiplying p i, busyday by this factor at the time of simulating data.We choose to increase p i, busyday and not the parameter  because it is variable X, and not Y, that is affected by seasonality.
In summary, the demand of customer i on a simulated "busy day" d is obtained from Equations (30)-(32): The value of  plays a crucial role in ensuring the robustness of our solutions; it will determine by how much we underestimate or overestimate the fleet size required.We tune this parameter by inspection to ensure that the average demand of the simulated instances is approximately equal to some demand level that we define in advance.In particular, we have explored  several demand levels that may be of interest to the business: the maximum in 2018, the 95 th percentile (P95) and 75 th percentile (P75).Finally, we solve the fleet sizing problem by completing the following steps: Step 1. Choose a design criterion with its associated  from Table 1.
Step 3.For each generated instance, solve a VRP and store the number of vehicles required at each depot.
Step 4. Set the final fleet size assigned to each depot to the 95 th percentile of the list of n vehicle requirements obtained in Step 3.This is different from the P95 used for simulating demand.
We choose to set the fleet size assigned to each depot to the 95 th percentile of the list of n vehicle requirements to be conservative: it is sufficient to route the demand of 95% of the generated instances and it is less sensitive to extreme values in the results than the maximum requirement.We set the number of simulations n to 50, which is high enough to obtain a stable number of vehicles with the 95 th percentile, but not too high that it takes a prohibitive amount of time to run the fleet sizing model.Lastly, we solve these 50 routing instances using our heuristic algorithm described in Section 8 with the standard savings function and no route improvements, all of which speed up the process and lead to more pessimistic and conservative results.
Results of running the fleet sizing algorithm for the four proposed demand levels in Table 1 are displayed in Figure 3.To simplify the analysis, we ran the algorithm for a single representative number of depots (four depots).The horizontal lines represent the fleet size obtained using different demand levels, whereas the time series represent the actual vehicle requirements on any day in 2018 and 2019.We can see that the same seasonality pattern is present in 2019 and there does not seem to be a positive trend in the year-over-year growth, although the maximum fleet requirement in 2019 (16 vehicles) was one unit higher than the highest requirement in 2018 (15 vehicles).
When the simulated demand is calibrated using the maximum demand in 2018, the resulting fleet size is sufficient to route all demand in 2018 and 2019, as can be seen in Figure 3.However, if the decision maker is risk-averse and does not have any contingency measures for when the existing fleet cannot cope with the upcoming demand, they should use the "maximum demand 2018 + 10%" criterion, which, in this case, gives two additional vehicles for safeguarding.
On the other hand, the P75 and P95 design criteria yielded fleet sizes that did not meet the vehicle demand in 2018 nor 2019.However, it is worth noting that the maximum demand in 2018 was significantly higher than average levels throughout the year: the maximum number of customers was 25% higher than the 95 th percentile and 45% higher than the 75 th percentile.This shows that designing for the peak will result in a low utilization of vehicles for most of the year.Therefore, if the business had the ability to rent additional vehicles for a day on short notice, there would be scope for potential savings by designing the fleet using the 75 th percentile criterion.For the case presented in Figure 3, this resulted in a fleet of 12 vehicles, which suffices to serve the demand of most days.
To estimate the impact on costs of these different approaches, we present in Table 2 an analysis of the additional rental vehicles that would be required in 2018 and 2019 using different calibration criteria.To be conservative, we have assumed that the daily rate to rent a vehicle on short notice is triple that of standard vehicles, namely, $90/day.
It is clear from these results that designing the fleet for peak demand results in a low utilization of vehicles for most of the year.If it is possible to rent vehicles on an ad hoc basis, even under our conservative assumption on rental costs, the cost of doing so using the P75 and P95 demand 2018 criteria would be negligible compared to the fixed vehicle costs.
Our approach to robust fleet sizing has the additional advantage that it allows the decision maker to manage risk and costs effectively.These criteria are intuitive and do not involve tuning a mysterious numeric parameter that is hard to interpret;  gets updated in the background.It took between 25 ("P75 demand") and 90 s ("Max demand + 10%") to run this fleet sizing model in a four-depot scenario.

ASSIGNMENT OF CUSTOMERS TO DEPOTS
We now enter the operational planning stage of the problem, which aims to determine the routes to be completed by each vehicle in each trip such that all the demand of a given day is served.
The mathematical model presented in Section 2 would simultaneously assign customers to depots and determine the optimal routes.However, this would require to solve a large multi-depot VRP.To overcome this, we assign customers to depots a priori.This allows us to treat the routing problem as a set of independent, smaller VRP's, one for each depot.
A reasonable approach to achieve this would be to assign customers to their closest depot for any demand realization.However, this does not take into account the spatial distribution of the demand on a given day, for example, it might be "easier" to visit a customer by a vehicle completing a trip nearby, even if this vehicle is based on a depot that is not the closest to the customer.
For this reason, in this section we explore the effect of factoring in the "routing effort" required to visit a customer.Jarrah and Bard [14] designed a metric based on the intra-customer travel time to address the problem of geographic clustering of customers.They estimate the travel time from customer i to the next customer based on the distance to the k nearest customers in the same cluster.
We have explored this idea of the distance to the k nearest customers as a proxy of the routing effort in the problem of assigning customers to depots on a given day.The approach we propose is as follows: starting from an initial assignment of customers to their closest depots, we iterate over the set of customers and we assign each customer i ∈ V to a depot  ∈ H such that the average distance from i to the k nearest neighbors currently assigned to depot  plus the distance from i to the depot is minimum.This metric is depot of customer i ≡  * = arg min where V k  is the set k customers currently assigned to depot  that are the closest to customer i.Because every time we change the depot assigned to a customer this might change the k nearest neighbors of customers that we have already scanned, we repeat this process iteratively until the assignments are stable, that is, until there is no change in the assignment of any customer in one round.
Including the distance to the depot in Equation ( 33) is necessary to ensure that, in a stable solution, the clusters of customers are still centered around the depots.Without this term, the iterative assignment process could converge in |H| clusters of customers with no connection with the depot locations.This was not considered in [14], since this study was only concerned with cluster construction.
It only remains to choose a value of k. Figure 4 shows routing costs on randomly selected instances and the impact of parameter k.Including the nearest neighbors in the assignment always results in savings in routing costs for the values of k tested in this experiment.With k = 2, there is a cost reduction of circa $45/day with respect to assigning customers to the closest depot (k = 0).Although this value is small compared to the daily routing costs, it adds up to an annual $18k saved on routing costs and it can be achieved with very little computational effort.

VEHICLE ROUTING
At this point, the original routing problem has been divided into |H| subproblems (one for each depot) of a more manageable size.These subproblems belong to the class of Multi-Trip Vehicle Routing Problems with Time Windows (MTVRPTW) [2,5].
There is now a single depot in each subproblem and a known fleet size denoted by |K|.The mathematical model in Section 2 can be employed to solve some of these instances, with some considerations: 1. Decision variables z i , which determine the location of depots, are now known and can be treated as problem data (z i = 1 if customer i is the depot of a given subproblem, 0 otherwise).2. Constraint (7), which fixes the number of depots to |H|, is no longer needed.3. Constraints (12), which ensure that the same vehicle cannot complete two trips in two different depots, can also be removed.4. We set the cardinality of the set of vehicles K to the number of vehicles assigned to the depot to ensure the number of vehicles used does not exceed the fleet size.Even though the daily fleet costs are fixed regardless of the daily vehicle utilization, we keep them in the objective function to encourage solutions with high fleet utilization.This avoids solutions where many or all the available vehicles are assigned to one or two short trips, which would not be desirable from a driver rostering perspective.Given the low cost of the vehicles, it is safe to add fleet costs in the objective function without significantly affecting the routing costs, which are the main cost component to minimize at this stage.
The computational experiments carried out to test the limits of this mathematical model, presented in Section 8.3 below, revealed the need for a more sophisticated exact method for smaller instances and a heuristic approach that can handle the larger instances in reasonable time at the expense of generating suboptimal routings in some cases.Based on these results, we adopt a hybrid approach: on any day, instances (or depots) with fewer than 15 customers are solved using the branch-and-cut algorithm described in Section 8.1 with a time limit of 20 min; for larger instances we resort to the heuristic method described in Section 8.2.This time limit of 20 min is high enough to obtain acceptable solutions for small instances but low enough that routes can be generated in an acceptable amount of time after the demand data becomes available at 6:00 p.m.

Exact method
In an effort to improve the efficiency of the exact method implemented by CPLEX, we propose a VRP-specific branch-and-cut approach.For this purpose, we define a new set of valid inequalities that are added sequentially in the form of cuts in the branching tree.We use the dual simplex method to solve the linear programming (LP) relaxation of the original MIP problem at every node, disable all of the built-in cuts used by CPLEX, apply a coefficient reduction technique and use our heuristic solution to warm start the MIP.
We ran a series of experiments to decide on the most appropriate method to apply the cuts throughout our branching tree, while balancing the time needed for separation, the number of additional constraints and the rate of improvement.As a result of our experiments, we perform a separation procedure in 10 consecutive nodes in the branching tree, remove all possible violations, and repeat this separation once in every 500 inspected nodes.This way, when the algorithm enters the 1500 th node, only 30 nodes will have been inspected for violations, and the risk of adding an overwhelming number of cuts that may slow down the solution process is lowered.All the cuts added are globally valid, and the results obtained are presented in Section 8.3.

Valid inequalities
If the depot nodes are excluded from the graph of our problem, a feasible solution should only consist of open paths.Therefore, any cycle in this graph represents a violation of a constraint of the problem and should be eliminated.This logic leads to the most common valid inequalities for the TSP, known as subtour elimination constraints.Even though in our case the time window constraints ( 15)-( 21) prevent subtours in any integer solution, the fractional solutions resulting from solving the LP relaxation of the problem can easily contain subtours.By implementing an adequate separation procedure, one is able to identify these violations and apply an appropriate constraint that is globally valid, thus making this solution infeasible in all future nodes.The constraints known as lifted subtour inequalities [12] are facet defining for the TSP for |S| ≥ 3, S ⊆ V.These inequalities have shown good results when applied to the VRPTW [4] and will be used to tighten the polyhedron of the LP relaxation of our problem.Due to the form of our degree constraints (2) and ( 3) and the fact that we are allowing multiple trips per vehicle, we need to implement the lifted subtour inequalities in a different form: The subtours of size S = {i, j} (|S|=2) will be cut off using x kr ji ≤ 1 (36)

Separation procedure
Once a feasible solution for the LP relaxation has been found, the edges from the set E in the graph G(V, E), where V is the set of original customers, are assigned a weight w(i, j) = ∑ k∈K ∑ r∈R x kr ij .To separate the violated subtour inequalities (34)-(36), we use a procedure called shrinking.A detailed description of this procedure applied to the TSP can be found in Lawler [19], and similar procedures are commonly used to separate this kind of inequalities.The separation procedure consists of creating additional nodes in the graph, while simultaneously removing both end nodes of an edge with positive weight.Each new node created in the graph is a child node, and the nodes that are removed during this process are its predecessors.That way, each child node has a line of predecessors that can be backtracked to the set of original nodes that generated it.Let n = |V| and assume a feasible solution has been obtained for the LP relaxation.The procedure is as follows: Step 0. For each node i ∈ V initiate an empty set of predecessors.
Step 1. Find an edge (i, j) such that w(i, j) > 0.5 and w(j, i) > 0.5.If such edge exists, go to Step 5. Otherwise, go to Step 2.
Step 2. Find an edge (i, j) such that w(i, j) > 0. If such edge exists, go to Step 3. Else, Stop.
Step 3. Create child node  := n + 1, and update the weights w(s, ) := w(s, i) + w(s, j) and w(, s) := w(i, s) + w(j, s), for every s ∈ V, s ≠ i, j.Update n := n + 1. Step 4. Remove 's predecessors i and j from V (V := V \ {i, j}).Update the set of 's predecessors with nodes i and j and all of their own predecessors.Go to Step 1.
This procedure is illustrated in Figure 5. From F, we can see that the set S of original nodes contains nodes i = 6 and j = 4 with their predecessors: S = {1, 2, 3, 4}.Note that node i = 6 is not in S as it is not an original node.

Routing heuristic
For instances with more than 15 customer nodes, we have developed a heuristic algorithm.This heuristic has three stages: first, the multi-trip characteristic is ignored and routes are constructed considering only vehicle capacities and time windows.Next, solution improvement techniques are employed to improve the routes constructed in the previous step.Lastly, a packing algorithm assigns routes to trips of vehicles.

Route construction heuristic
Our route construction approach is inspired by the adaptive routing heuristic proposed by Batarra et al. in [5].The authors first determine the period of time during which each route is strongly active.A route is strongly active at some time interval if it is guaranteed that the route will be in progress at that time, regardless of when it starts, see Figure 6.
Since the number of routes that are strongly active at the same time define a lower bound on the number of vehicles required, Batarra et al. [5] penalize routes that are strongly active at the busiest time, the critical time, and adjust the insertion costs in their route construction heuristic to include these penalties.The expected result is that routes constructed using these penalties will avoid the critical times and will therefore be easier to pack in the packing stage of the algorithm.They repeat this two-stage process iteratively (constructing routes -calculating penalties) until some stopping condition is met.
In the VRPTW considered in this work, all customers can be visited between 8:00 a.m. and 4:00 p.m. and therefore the critical time will always occur around the center of this time window.This is not the case in [5], where different customers can have different time windows and therefore critical times may occur at any time of the day.We take advantage of this property of our problem to introduce two substantial simplifications to the approach in [5]: first, the penalty for the critical time interval is a fixed value that we set in advance; second, routes are only constructed once, including the penalties.
We propose a parametric savings heuristic based on the general parallel savings heuristic developed by Clarke and Wright [7].We define the savings function as follows: The first term is the well-known savings function, namely, s routingtime (i, j) = t 0i + t j0 − t ij .It represents the savings in traveling time when routing a vehicle from depot 0 to blackcustomer i to customer j and back to the depot as opposed to visiting i and j alone in different trips.We introduce a second term, s stronglyactive (i, j), which represents the benefit of joining route finishing in customer i, with route starting with customer j in terms of how "easy" the resulting route would be to pack compared to how "easy" it would be to pack the individual routes.s stronglyactive (i, j) is defined in Equation (38): s stronglyactive (i, j) = f correcteSA (route with i and j) − f correcteSA (route of i) − f correcteSA (route of j). ( The term f correctedSA represents a function that penalizes routes that are strongly active at some time.We introduce here an additional modification with respect to the approach in [5].We seek to penalize the amount of strongly active time only up to a certain point.As the duration of a route starts to approach the vehicle time window length, that is, 11 h, it becomes more convenient to keep expanding this route, because the remaining free time of the vehicle would become more and more difficult to fill by another route.In other words, it is the routes of medium duration which can make the packing of routes in vehicles inefficient.Using this logic, we define f correctedSA as follows: where f SA is simply the amount of time that a route is strongly active.Based on this definition, f correctedSA = 0 when a route is not strongly active at any time (f SA = 0), it increases linearly reaching a maximum when f SA = 5.5 h, at which point it decreases FIGURE 6 Graphical definition of the strongly active interval of a route until f SA = 11 h, where f correctedSA is equal to 0 again.This means that no penalty is applied to routes that are not strongly active at any time nor routes that are 11 h long.Lastly, the term  in (37) controls the weight of the term s stronglyactive (i, j) in the savings function and it has been set to 0.05 through computational experiments.
The rest of the algorithm follows the same structure as the parallel savings algorithm by Clarke and Wright [7].Pairs of customers (links) are sorted in decreasing order by their savings.Starting from the first link on the list, links are used to join smaller routes into larger routes as long as they are feasible in terms of vehicle capacity and time window restrictions.The algorithm stops when the list contains no more feasible links.
Including the term s stronglyactive to the savings function does add a small amount of computational time, as this term varies with the duration of the routes, and the savings list needs to be updated multiple times.Nevertheless, the total time in the worst case scenario, a single depot on a "busy day," was under 6 s.

Solution improvement
After constructing the initial routes, two solution improvement techniques are employed to possibly improve the current routes by reducing the routing costs.First, a TSP with time windows (TSPTW) [11] is solved to optimality for each of the constructed routes to ensure customers are visited in the shortest possible path within the route.For the sake of brevity, the MILP formulation of this TSPTW is not presented here.Last, we apply an exchange heuristic where we swap two customers visited in different routes with the aim to reduce the routing distances further [23].

Packing heuristic
The heuristic we propose is a modification of the well-known first fit decreasing (FFD) heuristic [10], which sorts items in descending order of size, takes one item at a time and either assigns it to the first bin where it fits or opens a new bin if it did not fit in any of the open bins.
We propose a variation of the FFD heuristic in which the items (the routes) are also sorted in decreasing order by their size (route duration) and they are assigned to the first vehicle in which they fit, at the earliest possible start time, such that the time windows are not violated.The main difference with the FFD is that, when we assign a route to one trip to a vehicle, we pay attention to where in the current schedule of that vehicle it gets inserted.
We will insert the route at the start of the schedule if doing so reduces the difference between the depot opening time and the earliest start time of the first trip currently scheduled on that vehicle.If that is not the case, we insert it at the end if the route has a longer travel time from the last customer to the depot than the current last trip.Otherwise, we insert it in any intermediate position.The rationale for this modification of the FFD algorithm is that we want to maximize the use of the 2-h and 1-h gaps at the start and end of the day, respectively, where the depot is open but customers cannot be served.
Furthermore, we observe that most routes can be operated in either direction, namely, customers can be visited in reversed order.This is true because the time windows of all customers are identical and the graph is undirected.The only case in which a route cannot be reversed is illustrated in Figure 7B.In light of this observation, we allow routes to be inserted in reverse order if it is feasible and beneficial.A detailed description of this heuristic is presented in Algorithm 1 in the Appendix.Lastly, we define the key parameters of this packing heuristic below, which are used in Algorithm 1.These definitions are accompanied by a graphical description in Figure 7A.

Comparison of routing methods
Table 3 shows the results obtained using (i) the branch-and-cut algorithm with an initial heuristic solution and our subtour elimination constraints as the only cuts (34)-( 36), (ii) the CPLEX default settings, which automatically decide between dynamic search and branch-and-cut with default CPLEX MIP cuts, depending on the model, and (iii) our heuristic method.The number of cuts found by our separation procedure is also shown for every instance.In both exact methods, the optimality gap tolerance was set to 0.001.Other variants of our branch-and-cut algorithm were also explored, namely, a combination of subtour elimination cuts and CPLEX cuts (with and without an initial heuristic solution) and subtour elimination constraints without warm start.However, we have not included these results in Table 3, as these algorithms were outperformed by our final exact method.Our heuristic method obtained optimal routing costs in all instances with less than 15 customers and it often outperformed both exact methods with a 20-min time limit for larger instances.The number of vehicles obtained is different in some cases, but there does not seem to be a method that consistently requires more vehicles.
The comparison between our branch-and-cut algorithm and CPLEX with default settings looks promising.Even though this comparison gives CPLEX an advantage, our branch-and-cut algorithm outperformed CPLEX in many instances and when CPLEX's results were better, this difference was relatively small.Furthermore, we ran several experiments comparing our Notes: Results compare our branch-and-cut algorithm with subtour elimination cuts and warm started with a heuristic solution, the default strategy in CPLEX and our routing heuristic.The time limit for the exact methods was set to 1200 s, as per our proposed routing method, and the optimality gap tolerance was set to 0.001.Results in bold are the lowest routing costs of the two exact methods.branch-and-cut algorithm to a branch-and-bound with no cuts.We observed that ours performed significantly better, highlighting the efficiency of our adapted subtour elimination cuts.Another observation is that our algorithm seems to reduce the peak memory used.Based on several experiments run in AIMMS for instances of up to 25 customers, we observed that after 20 min the peak memory was over 60% lower with our subtour elimination cuts with an initial heuristic solution compared to CPLEX's default.
Lastly, Table 4 presents a separate computational experiment run on larger instances with 25-150 customers (the highest number of customers in a single day), solved using the routing heuristic.The average run time is 30 s and all instances were solved in under 2 min.

SOLUTIONS
Solutions have been generated for each day in the provided data sets and every number of depots from 1 to 10.The resulting costs are summarized in Figure 8. Due to the long computational times that would have been required to run our hybrid routing algorithm on 730 days, these routing problems have been solved using the heuristic method only.They represent a worst case scenario, as using an exact method for small instances can only improve these results.
There is a clear trade-off between depot costs and fleet and routing costs.As expected, fleet costs contribute the least to the total daily cost due to their low price, and they remain almost constant from five depots onward.Routing costs are the only cost component that varies on a daily basis and as such, average routing costs are different in 2018 and in 2019.The routing costs in 2019 are slightly lower, although the difference is not significant.Two conclusions can be drawn from this: 2019 demand did not change substantially with respect to 2018 and our methods produce solutions that are robust against unseen data.
Lastly, it is clear that a single depot is not cost-effective and the lowest total costs are obtained for two to five depots.This may be explained by the fact that there are two major cities in the State of Pennsylvania: Pittsburgh and Philadelphia, as illustrated in Figure A1a.These urban areas concentrate a significant fraction of the demand such that once a depot is placed in each of these city areas, the marginal benefit of additional depots is greatly reduced.

OPTIMIZING THE NUMBER OF DEPOTS
In this section, we develop a framework for automating the choice of the number of depots making use of our existing models and performing an enumeration.To obtain a solution in a reasonable amount of time, we simplified our methodology by (i) halving the number of simulated days needed to determine the fleet size; (ii) designing the fleet for the "P75 demand in 2018"; and (iii) using for the routing stage the same fast routing approach employed in the fleet sizing model.Furthermore, instead of considering all days in 2018, we selected 15 at random to obtain an average daily routing cost.

h kr 1 Time the depot is visited by its assigned vehicle k at the start of trip r h kr 2 Time
the depot is visited by its assigned vehicle k at the end of trip r

FIGURE 1
FIGURE 1 Visualization of the trade-off between optimizing for minimum cost of locating depots and minimum average distance customer-depot.Optimal solutions are obtained for combinations of 11 different, evenly spaced, values of  ∈ [0, 1] and four representative numbers of depots

FIGURE 2 (
FIGURE 2 (A) Comparison between actual distribution of Y (units of demand per active customer) in 2018 and the Poisson distribution fitted to this data.(B) Total demand and number of active customers per day in 2018 [Color figure can be viewed at wileyonlinelibrary.com]

FIGURE 3
FIGURE 3 Minimum number of vehicles required to serve the demand everyday in 2018 and 2019, compared to the fleet size calculated using different design criteria.All the VRP instances considered in this analysis have been solved using the heuristic method from Section 8 to speed up the solution [Color figure can be viewed at wileyonlinelibrary.com]

FIGURE 4
FIGURE4 Impact on routing costs of including the average distance to the k nearest neighbors when assigning customers to depots.The results for the different values of k represent the average routing costs from 50 randomly selected instances with different numbers of depots

FIGURE 7 (
FIGURE 7 (A) Illustrative schedule with two trips displaying problem parameters.(B) Example of a case where inverting the direction of the trip results in an infeasible schedule because the vehicle returns to the depot after b D

TABLE 1
Estimated values of  to generate a set of instances with an average demand (in number of customers) approximately equal to four different levels of demand from 2018 data

TABLE 2
Comparison between different design criteria in terms of fleet size, fleet fixed costs, additional vehicles required when the existing fleet is insufficient to meet demand, associated costs of these extra vehicles and total vehicle costs, for each year Earliest time that trip r could be scheduled to start, namely, max{a D , a C − r.first} k.gapTime gap between a D and the earliest start time (r.earliest) of the first trip for vehicle k

TABLE 3
Computational results of several real instances in January 2018 ranging from 6 to 22 customers

TABLE 4
Run times of 50 routing instances (10 for each range of number of customers) in 2018 ranging from 25 to 150 customers solved using our routing heuristic Average daily costs in 2018 and 2019 for 1-10 depots, broken down into depot, fleet and routing costs.Routing problems are solved for all days in the data sets provided using the heuristic method to speed up the solutions [Color figure can be viewed at wileyonlinelibrary.com]