The multi-depot k-traveling repairman problem

In this paper, we study the multi-depot k-traveling repairman problem. This problem extends the traditional traveling repairman problem to the multi-depot case. Its objective, similar to the single depot variant, is the minimization of the sum of the arrival times to customers. We propose two distinct formulations to model the problem, obtained on layered graphs. In order to find feasible solutions for the largest instances, we propose a hybrid genetic algorithm where initial solutions are built using a splitting heuristic and a local search is embedded into the genetic algorithm. The efficiency of the mathematical formulations and of the solution approach are investigated through computational experiments. The proposed models are scalable enough to solve instances up to 240 customers.


Introduction
In recent years, among practitioners and researchers, there has been a growing concern about the investigation of different problems arising in the class of Traveling Repairman Problem (TRP). Originally proposed in the repair and maintenance service context, this problem class has recently received attention from practitioners and researchers, stimulated by the wide range of applications of the problem in fields where the minimization of the sum of arrival times to customers is more important than other objective functions. Routing problems for emergency operations and postdisaster relief activities [4], home delivery services [23] and disk-head scheduling [3] and other customer-centered problems are examples of elective applicative fields.
The Multi-Depot k-TRP (MDk-TRP) generalizes the traditional single depot TRP to the multi-depot case and consists in minimizing the sum of the arrival times to the customers visited by a homogeneous fleet of k vehicles, housed at different depots. The introduction of cumulative costs into the objective function (i.e. the arrival time instead of the travel time as common in most routing problems) poses severe computational challenges [28] and the presence of multiple depots makes the problem even more complicated.
There is a lack of contributions on this problem variant, when compared to the standard multi-depot vehicle routing problem [25], even though it better reflects many real-world logistics systems. In the field of city logistics, for instance, widespread customer positions coupled with the sustainability challenge compelled the use of multiple distribution centers (also called satellites) [30]. This two-tier system is advantageous since each customer can be visited by vehicles from the closest satellite, reducing traffic jams and decreasing the level of emissions of associated pollutants.
The characteristics of spread customer positions, and multiple distribution centers is not only related to the last mile distribution system. Also in the health-care context, especially from the perspective of the design and optimization of emergency logistic systems in the aftermath of a disaster, relevant problems have been often modeled as a multi-depot routing problem [40,43]. Disasters can be devastating to communities and governments, and decision-making, in this customer-centric setting, needs to be quick and efficient since it becomes crucial to plan life-saving relief plans. There are numerous areas of disaster relief management for which the benefits of using multiple distributions centers are valuable. Timely delivery of medical supplies or reliefs, post-disaster assessment and recovery, are all linked to effective, customeroriented multi-depot routing problems [45].
In this paper, we present two mathematical formulations for the MDk-TRP and we develop a genetic algorithm. Extensive computational experiments, run on a set of benchmark instances, show the efficiency of the proposed approach. The remainder of this paper is organized as follows. In Sect. 2, we provide a detailed review of the relevant literature. In Sect. 3, we describe the problem and present two new mathematical formulations using the concept of layered network. In Sect. 4, we explain the proposed genetic algorithm in detail. Section 5 presents the computational results carried out on a set of benchmark instances. Finally, our major conclusions and future research directions are presented in Sect. 6.

Related works
The TRP, also known as the minimum latency problem or the delivery man problem, has been extensively studied by a large number of researchers who proposed several exact and non-exact approaches. In [20] and [2] exact enumerative algorithms were developed, in which lower bounds are derived using a Lagrangian relaxation. An enumerative algorithm that makes use of lower bounds obtained from a linear integer programming formulation was proposed in [9]. A metaheuristic method which is based on a Greedy Randomized Adaptive Search Procedure (GRASP) and Variable Neighborhood Search (VNS) was proposed in [37]. A general VNS metaheuristic [24], was able to outperform the one suggested by [37].
The generalization of the TRP to the case with multiple identical vehicles (k-TRP) was first presented in [8]. In [21], the authors proposed an exact branch-and-priceand-cut algorithm for the multi-vehicle TRP with a distance constraint. An ad hoc label-setting algorithm with bi-directional search strategy is developed to solve the pricing problem. They tested 180 instances with up to 50 nodes and obtained an average gap less than 0.5%. In [28], a new and efficient mathematical formulation for the k-TRP, defined on a multi-layer network, was presented, then solved by an iterative greedy metaheuristic approach. In [1], five mathematical models for the multi-vehicle TRP are introduced. The first three models are derived from the classical and flowbased formulations, while the last two ones are generalizations of the time-dependent TRP formulation. The authors tested the formulations with instances up to 80 nodes and 16 vehicles, where the largest instance is solved in 86 s. The authors showed that the last two models perform better in terms of computational time. Interesting variants of k-TRP with service level as well as with profits or under uncertainty on travel times have been recently investigated [4,5].
A generalization of the TRP, which considers a demand associated to each node and adds capacity constraints is known as the Cumulative Capacitated Vehicle Routing Problem (CCVRP) [17]. In the last decade, a flourishing literature stream has been dedicated to the study of mathematical models [36], exact algorithms [22], as well as heuristic and metaheuristic approaches [18,29,35] for the class of CCVRPs. The investigation of heuristic and metaheuristics approaches for the CCVRP was conducted in [26] where a memetic heuristic procedure was developed; an adaptive large neighbourhood metaheuristic was proposed in [34]. In [39], the authors designed a hybrid metaheuristic algorithm for the min-max CCVRP enhanced with a two-stage adaptive variable neighbourhood search. An exact column generation approach was developed in [12], and later, the same authors proposed a heuristic approach in [13]. In [27], the authors introduced two new mathematical models for the CCVRP and proposed an iterated greedy algorithm outperforming other existing methods in the literature. Recently, in [11] was studied the CCVRP under uncertainty of demands, for which an approximation algorithm was developed.
The MDk-TRP is a variant of k-TRPs which considers the presence of multiple depots. Despite the richness of contributions for the multi-depot vehicle routing problem (interested readers are referred to [25,33,38] and references therein), a few papers deal with multi-depot variant of the TRP. We cite [19], where a three-index mathematical model is proposed for the multi-depot CCVRP. The authors also proposed a partial optimization metaheuristic that decomposes the problem into sub-problems, later solved using either an approximate or an exact approach. To overcome the computational challenge posed by large-scale multi-depot CCVRP instances, a simple and efficient heuristic was proposed in [44] able to outperform the heuristic presented in [19] in terms of solution quality, computational time, as well as stability.

Mathematical models
The MDk-TRP is defined on an undirected graph G = (V , E) where V = C ∪ D represents the vertex set composed of the set of customers C of cardinality n and the set of potential depots D of cardinality m. The set E denotes the edge set. Each edge (i, j) ∈ E has associated a travel time t i j for which the triangular inequality holds. According to the application at hand, this traveling cost can be interpreted as distance, time, fuel consumption, etc. A homogeneous fleet of k vehicles is dispatched from the m potential depots to visit all the n customers. Decisions also entail the selection of a subset of depots to host the vehicles and the determination of the number of vehicles to be sited at each active depot. The objective is to find a feasible routing plan such that the sum of the arrival times is minimized. Since the objective function is defined as the total latency (total waiting time), the travel time between the last visited customer by the vehicle and the depot to come back is not concerned. Therefore, we apply the terms tour and path (route) interchangeably throughout the paper. Strong mathematical models for the k-TRP [28] use a multi-layer network representation. The layered network is described as follows.
Let L be the set of levels L = {1, . . . , r , . . . , N }, where N = n − k + 1, and each level includes a copy of all the customers amended also with the depot in levels from 2 to N . Each path ends in the first level and starts in a copy of the depot in some level. In fact, the level number represents the position of the customer: the customer in the first level is the last one, the customer in the second level is the last but one, and so on. By using this network, it is easy to enforce the restriction that two distinct vehicles cannot visit the same customer, neither in the same level nor in different levels. In the next section, we will present two variants of the layered network enabling two new mathematical formulations for the MDk-TRP.

First formulation
Hereafter, we present the first mathematical formulation for the MDk-TRP. To extend the network to the multi-depot context, we add copies of the depots from levels 2 to N . Figure 1 illustrates this approach.
In Fig. 1, customers are displayed by squares and depots by circles. Node 0 denotes an auxiliary depot to specify the start of each path. As it can be seen, for each level, all the copies of the depots in the set D are linked to the corresponding starting node 0 at the same level. In other words, when the node 0 is active, one of the depots must be associated to it indicating the start of the path. The variables are defined as follows. Let x r i be a binary variable that takes value 1 iff customer i ∈ C is visited at level r ∈ L. If x r i = 1, customer i is active at level r and there are r − 1 customers to be visited after in the same route. Let y r i j be another binary variable that is set to 1 iff j ∈ C is visited right after i ∈ C ∪ {0} and there are exactly r customers to be visited after node i. Finally, w r i j represents a binary variable that takes value 1 iff j ∈ C is the first customer to be visited (at level r , r ≥ 2) by a vehicle dispatched from the depot i ∈ D and there are r customers to be visited after. In other words, w r i j = 1 only if y r 0 j = 1, and 0 otherwise. The first mathematical model is expressed as follows.
r ∈L The objective function (1) minimizes the sum of arrival times to customers; here the coefficient r counts the number of times the travel time of a link is considered in the evaluation of the total latency. Constraints (2) guarantee that each customer is visited exactly once. Constraints (3) and (4) ensure the generation of exactly k routes. The constraints in (5)-(7) are the connectivity constraints and show the relation between the binary variables x r i and y r i j . Constraints in (5) require that any customer i visited at the upper level r + 1 should be connected to exactly one customer (let say j) at level r by traversing edge (i, j). Constraints (6) impose that any customer j visited at level r should be linked to exactly one recently visited customer (let say i) by traversing edge (i, j) or linked directly to the auxiliary node by traversing edge (0, j) in the same level. Constraints (7) require that customers visited at level N should be visited right after the node 0 by the link (0, j) in the same level. Constraints (8) show the relation between binary variables w r i j and y r 0 j ; if y r 0 j = 1, customer j is the first visited node at level r , hence it should follow the depot (let say i) in the same level r . Finally, constraints (9)-(12) establish the nature of the variables. This formulation has n[N (m + n + 1) − (n − 1)] binary variables.

Second formulation
Another variant of the multi-depot layered network can be obtained by replacing each copy of node 0 in levels 2 to N , with m nodes, indexed from 1 to m, representing the multiple depots. In this way, levels from 2 to N include a copy of each customer and of each depot. In the first level, only the customer copies are present. Each tour in the network is represented by a route that ends in a first level customer and starts in a copy of the depot. Figure 2 illustrates this approach.
The way in which the layered network is defined implies that the linkage variables y r 0 j are replaced with binary variables y i j r with iin D and j ∈ C. Clearly, in this case, variables w r i j are no longer needed. The second mathematical formulation is presented as follows.
r ∈L r ∈L i∈D j∈C x i r ∈ {0, 1} i ∈ C, r ∈ L (20) Similar to (1), the objective function (13) is the waiting time of all customers and it is composed of two terms: the first one accounts for links originating from depots and the second one for the other links. In a similar way, the set of constraints (2)-(7) are equivalent to the set of constraints in (14) We should mention that for the special case of |D| = 1 (single depot), both models reduce to the classical k-TRP which is known to be NP-hard [16].
vide the outline of the metaheuristic approach an then detail the algorithm elements and the search mechanism.

Outline of the genetic algorithm
Genetic algorithms are evolutionary algorithms inspired by the theory of evolution, that use concepts such as natural selection, reproduction, genetic heritage and mutation [15]. Genetic algorithms may be enhanced with local search and/or other heuristic schemes. In this case, they are called hybrid genetic algorithms. Several routing problems have been successfully addressed by genetic algorithms [31,41,42].
Let P andP denote the old and the new population composed of p s individuals (solutions/chromosomes); e s represents the size of elite solutions, μ r shows the mutation rate and p * denotes the best solution with the fitness value f ( p * ) . The pseudo-code of the proposed metaheuristic is presented in Algorithm 1. First, in Lines 2-6, the population set P is filled with randomly generated initial solutions { p 1 , . . . , p i , . . . , p p s }. In Line 7, the individuals are evaluated and sorted in nondecreasing order with respect to the fitness measure (see Sect. 4.4); also the best solution p * is updated (Line 8).
A pre-specified portion e s ∈ [0.1, 1] of the fittest individuals, let say the elite solutions, are selected to directly go into the next populationP (Lines 10-14). Then, a loop is run for p s − e s iterations where at each iteration, a new individual ofP is generated (Lines 15-30). Line 16 shows the tournament selection, a procedure which chooses a pair of parents from the population to be recombined using the crossover operator (Lines 17-18). The offsprings generated within the crossover phase are denoted withp (.) . In a local search procedure, the offsprings are educated (Line 19) and only the offspring with the lowest objective value (p) is passed into the next generation (Lines 20-29). In Line 28, the mutation rate μ r is updated appropriately (see Mutation rate update in Sect. 4.9). After, the new population setP is sorted (Line 31) and the sets P,P are updated. Lines 34-43 refer to the random restart mechanism which is, in fact, a diversification strategy that modifies the non-elite individuals in the population and updates the mutation rate, if necessary (Sect. 4.9). The algorithm ends at Line 45 and the best solution p * is returned. Lines 9-44, are executed until the stopping criterion is met. In particular, we set a certain number of iterations proportional to the problem size.

Solution representation
Each solution in our algorithm is uniquely encoded by two ordered lists: one, called service pattern, with size n representing a giant tour, without trip delimiters. The order of the customers follows the visit order in the k routes. To specify the starting point of each route in the giant tour, we consider another list, called depot assignment, with size k, whose q th element contains the index of the first customer visited in the q th route and the index of the depot assigned to the route q. This TSP-like encoding style simplifies the implementation of crossover and mutation operators since it is easy to retrieve information from the solution [41]. if (no improvement after a pre-specified number of iterations) then 35:

Algorithm 1 The proposed metaheuristic algorithm
Mutation.rate.U pdate(μ r ) 36: for (e s < i ≤ p s ) do 37:  Figure 3 displays the chromosome structure for a solution with three depots and ten customers:

Initialization
The initialization procedure generates a pool of p s feasible solutions and it is composed of two phases, including the giant tour generation and the splitting heuristic.
In the giant tour generation phase, an ordered sequence of customers is built in a greedy fashion. To do this, first, a non-visited customer is selected randomly. The giant tour is completed on the basis of the least-cost criterion until all the customers are inserted.
After, we apply a polynomial splitting heuristic proposed in [6,32] to find the optimal splitting points over the giant tour. This specifies the first customer to be visited in each route and also the optimal depot assigned to it. Obviously, once the starting point of a route is determined, we can retrieve it easily from the giant tour. The splitting can be done by solving a shortest-path problem on an auxiliary acyclic graph defined as A = (N , D) in which N = {0, 1, 2, . . . , n} and D is the set of possible subsequences of customers in the giant tour. The shortest path from 0 to n is calculated using Bellman's algorithm for directed acyclic graph (for more information, see [6,32]).

Fitness evaluation: biased fitness criterion
The performance of any population-based algorithm is affected by the fitness measure that enables to evaluate and order the individuals. We use the fitness criterion to select individuals for mating and to determine offsprings to keep in the next generation. In particular, we apply the biased fitness function criterion, f ( p i ) which is defined, for the individual p i , as follows [41,42].
where f it( p i ) refers to the rank of p i in the population, with respect to the objective value (total latency criterion) and dc( p i ) is the rank of p i with respect to the diversity contribution ( P( p i )) which is the average distance to its N C closest neighbours as given in (24).
σ H ( p i , p j ) is the normalized Hamming distance which quantifies the differences between two solutions p i and p j in terms of the visited customers and the depot assignment patterns. In general, N C is set to nc × |P| where nc is an input parameter [41]. The inclusion of a diversity criterion into the fitness evaluation helps to avoid early convergence to local optimal solutions.

Parent selection
The parent selection procedure (Parent.Selection(P)) returns two different parents from the population set to generate new offsprings. The first parent is selected randomly (with uniform probability) from the top 50% solutions, evaluated and ranked based on the biased fitness criterion while the second parent is chosen randomly from the whole population.

Crossover operators
The genetic search progress is obtained by two essential genetic operations, including crossover (recombination) and mutation operators. As the name indicates, a crossover operator recombines the information of two parents based on some operator-specific rules and generates one or more offsprings (in our case, two offsprings), representing new solutions. Generally, the crossover operator exploits a better solution while the mutation operator explores a wider search space. In particular, we propose three crossover operators as follows: -One-point crossover This operator picks a unique point randomly chosen in the range of φ ∈ (1, n), known as "crossover point" which divides each parent's giant tour into two segments with lengths of φ and n − φ. The left segment of the first parent's giant tour and the right segment of the second parent's giant tour are directly copied into the first offspring. By changing the order of parents and with the same crossover point, the giant tour of the second offspring is built. The same process is repeated for the depot assignment lists of both parents where the cutting point, here, is picked from (1, k), to be within the depot assignment list size. -Two-point crossover This operator is similar to One-point crossover in which two random "crossover points" φ 1 , φ 2 are selected to divide each parent's giant tour into three segments, including φ 1 , φ 2 − φ 1 , and n − φ 2 customers. The middle segment of the first (second) parent's giant tour is copied exactly into the first (second) offspring as the middle part of its giant tour and the first and the third segments in the second (first) parent's giant tour are copied into the first and third segments of the offspring's customer service. In order to preserve the structure of the routes, the "crossover points" φ 1 , φ 2 are selected in a biased way and pointing to the start of routes. The same procedure is also followed for depot assignment lists. -Uniform crossover In the uniform crossover, each gene is selected randomly from one of the corre-sponding genes of the parents with equal probability. The same rule is followed for building the depot assignment list.

Mutation
The mutation operator modifies the offsprings and acts as a diversification mechanism to avoid getting trapped into local optimal solutions. We propose four mutation operators that modify only the giant tours of the offsprings.
-Intratour 3-swap This operator randomly selects one route. If it visits more than three customers, a 3-swap intratour operator is run swapping the position of three randomly selected customers. Otherwise, three random positions in the giant tour are selected and their contents are swapped. The latter case, in general, can be regarded as a 3-swap intertour move.

-Intertour 5-swap
This operator selects five randomly positions in the giant tour and swaps their contents.

-Scramble
The scramble operator selects randomly a segment of the giant tour, and rearranges its content in a random fashion.

-Inversion
Similar to the Scramble, the Inversion operator picks randomly a segment and reverses the order of the elements.
It is worth noting that after the mutation and the crossover, a repair mechanism is run to retain the feasibility of the solution, if needed. In general, the feasibility is checked with respect to set of visited customers and the number of dispatched vehicles.
If some of the customers have not been visited, they are inserted in the position leading to the least increase into the objective function. We should note that since the giant tour contains exactly n elements, in the latter case, there exist at least another customer visited more than once. Hence, whenever a non visited customer is inserted, one of the repeated customers is removed from the tour. Moreover, in the case that more than k vehicles have been dispatched, the extra tours with the highest latency are found and discarded.

Local search
The local search procedure explores the 2-swap neighbourhood and updates the current solution whenever an improving move is detected (based on the first improving strategy). The neighbourhood is searched iteratively and the solution is updated until no further improvement is possible. The solution is modified to keep feasibility (from an open depot at least one vehicle should depart) and a right balancing among different tours (this measure is a proxy for near optimality of the solution). Each subtour is assigned to the best depot, based on the total latency criterion. Then, if it yields to an improving solution, the tour lengths are rebalanced. A tour to can be splitted, if it is beneficial in terms of total latency, if the number of dispatched vehicles is less than k.

Diversification mechanism
In order to diversify the solution pool, the algorithm applies two different diversification strategies as given below. As we mentioned earlier, a diversity criterion is also included into the fitness evaluation avoiding the premature convergence of the population.
-Random restart In case the best solution (incumbent) is not updated within a pre-specified number of iterations, a strong diversification mechanism is applied.
In particular, all the non-elite solutions in the population are replaced with new randomly generated ones. -Mutation rate update The mutation probability is dynamically modified on the basis of the search history. Whenever this restart mechanism is executed, the mutation rate is increased to diversify the search.
If the best solution is modified after the local search and the mutation rate value is greater than its initial value, the mutation rate is decreased to intensify the search.

Computational experiments
In this section, we report on the computational experiments carried out with the aim of testing the proposed models as well as the solution approach. We first assessed the efficiency of the two mathematical formulations, implemented using the AMPL optimization language [10] and solved by Gurobi 8.1.0 [14] within a time limit of 7200 s. To have a comparison, we also re-implemented (excluding the constraints on the capacity of vehicles) the models proposed for the multi-depot CCVRP, namely the one proposed in [19] and [44]. We refer to these models by the authors' names and denote our first and second mathematical formulations in Sect. 3 by F1 and F2, respectively. As test bed, we considered three sets of instances taken from the literature [7,19] 1 . The algorithm was coded in C++ language and all the experiments were executed on an Intel Core i5, with 2.4 GHz CPU and 8 GB RAM.

Results obtained considering benchmark instances
Tables 1 and 2 report the computational results obtained by Gurobi, in solving the mathematical models within the time limit. For each considered model, we report the objective value of the best solution found (O F), the computational time, in seconds, (C PU), and the optimality gap (Gap), in percentage, reported by Gurobi. The Tables also report when time limit T L is reached and a dash for instances in which no feasible solution within the time limit was found. The instance name is reported in the column with the header "Instance".       Table 1, we can observe that F1 and F2 found the optimal solutions for all the instances. The models adapted from [19] and [44] were able to find the optimal solution for the smallest instances with 10 customers. For the instances with 25 customers, after two hours, we observed considerably optimality gaps of at least 29% (Lalla's model) and 38% ( Wang's model). For all the instances with more than 50 customers, the models provided feasible solutions with optimality gap of about 53% and 99%, respectively. we should remark that, in this respect, the two models were proposed for the multi-depot CCVRP. Nevertheless, our models could be used to set good lower bound in reasonable computational time. In fact, in terms of solution time, F1 and F2 are significantly faster than the other models, and the largest Ir instance including 100 customers is solved to optimality in about 76 s. Table 2 summarizes the results for the set of p and pr instances, including up to 288 customers with a fleet size of k = 35 [7]. As shown in Table 2, the results of the proposed models F1 and F2 are superior to those obtained from other existing models. In particular, model F1 finds optimal solutions in at least 70% of cases, only in three cases no feasible solution is found and for those cases not solved to optimality, the gap is always less than 5%.
In addition, while, compared to F2, model F1 involves more binary variables, its computational time is significantly lower for all the instances solved to optimality. This behaviour is not related to the tightness of its LP relaxation bound. Tables 3 and  4 report the values of the linear relaxations of F1 and F2 and the gaps with respect to the Best Known Solution (BKS), in percentage. Both F1 and F2 present similar LP relaxation bounds, except for instances p08, p10, p18, pr01 and pr05 where the difference is less than 0.01%.
In the following, we present the results of the genetic algorithm. Specifically, we considered a population size of 25 individuals ( p s = 25); the initial mutation rate was set to 1 6 and its possible values lie in the set μ r = 1 γ , γ ∈ {1, 2, . . . , 6}. The elitism rate e s was set to 40% × p s and the parameter nc was set to 0.2. The maximum number of iterations was set to 10 4 × f p s where, depending on the instance, f ∈ {1, 2, 5, 10, 15, 20}. Finally, the random restart mechanism is activated every 100 iterations without incumbent update. The genetic algorithm was tested considering two different configurations: with one-point crossover and scramble operators used for reproduction and mutation, respectively and with operators chosen randomly (with uniform probability) among the set of operators defined in Sects. 4.6 and 4.7. Table 5 summarizes the metaheuristic results. Columns 1-4 show the instance features; Columns 5 and 6 represent the objective value of the BKS (obtained from model F1 or F2), and its corresponding solution time (he symbol "-" in the Column with header O F denotes that the solver could not obtain any feasible solution within the time limit). Columns 7 and 8 report the splitting algorithm results, including the best objective value (O F) and the relative gap, in percentage, (Gap) calculated with respect to the Gurobi O F, if available; otherwise denoted by "-". Columns 9-12 and 13-16, respectively, display the results of the genetic algorithm with and without the random selection of the operators. The column with heading T reports the percentage time fraction, calculated as ( T = C PU heu C PU Gurobi × 100) where C PU heu and C PU Gurobi denote the C PU of the metaheuristic and Gurobi, respectively. The results show that the average solution time of the heuristic over all the instances is less than 3 min, compared to Gurobi that in 12% of cases could not find any feasible solution within two hours. The effectiveness of the method is also confirmed by the average T values reported in Columns 12 and 16 which are around 15% and 22%, respectively. This means that the metaheuristic is about 7 and 5 times faster than Gurobi.
In addition, the metaheuristic finds optimal solutions for instances up to 25 customers. Even the splitting algorithm performs well since average solution gap is below 12.37% and in 3 cases (for which the Gurobi could not find any feasible solution within two hours), the splitting algorithm finds good feasible solutions. To have a clear idea about the performance of the metaheuristic compared with the splitting algorithm, we focus on the most challenging cases, for which Gurobi did not find any feasible solution. For such cases, the metaheuristics provide solutions which are on average 9.60% (with a minimum of 5.80% and a maximum of 19.74%) better than the splitting algorithm solution.
The results show a slightly better performance (about 0.20%), in terms of percentage gaps, when the random selection is present (see Columns 11 and 15). Although, in this case, the average computational time is lower, the speed up ratio T deteriorates.
In summary, for small instances optimal solutions are obtained in a negligible solution time, for medium instances good quality solutions can be obtained within acceptable computational time, and for the largest and the most challenging cases, for which Gurobi does not report any feasible solution, the algorithm is able to find solutions with gaps, calculated with respect to the lower bound, less than 23%,on average.

Results considering a reduced fleet size
With the aim of assessing the impact of the fleet size on the model and on the heuristic performance, we carried out an additional set of experiments, on the p, pr, and lr instances, considering a different number of vehicles. As a matter of fact, we observed that the optimal routes obtained considering the original datasets, contain only a few vehicles, maybe because they were originally proposed for the multi- Table 5 The metaheuristic results     Table 7 The metaheuristic results: instances with reduced fleet size  depot CCVRP. Therefore, we have considered a different number of vehicles, setting k = max m, n 10 , to make our problem more challenging. Table 6 shows the results; for the sake of completeness, we also report the solutions of both Lalla's and Wang's model imposing a time limit of two hours. The best objective function value is highlighted in bold. As confirmed by the gap values reported in Columns 7 and 10, the extant models could only solve to optimality the smallest Ir instances with 10 customers. For moderate instances, with 25 customers, the optimality gaps after two hours are very high (at least 61.62% for the Wang's model and 96.60% for the Lalla's one). When the number of customers is increased to 100, none of the existing models could find any feasible solution within the time limit. Nevertheless, both F1 and F2 solved all the set of lr instances to optimality, even if the decrease in the fleet size considerably increases the solution time. We also observe that the second formulation provides slightly better solutions, but in terms of solution time, model F1 is superior to F2.

Instance
Regarding the p and pr instances, we noticed that the existing models fail to find optimal solutions. In particular, the Lalla's model obtained feasible solutions for instances up to 100 nodes and two depots. However, for the solved instances, the average gap reported by this model is 84.47%. On the other hand, the Wang's model solved a similar number of instances (but it provided feasible solutions for instances up to 100 customers and four depots), reporting an average gap of 79.32.%. Regarding our formulations, the model F2 could solve to optimality all but eight instances, for which the reported gap is on average less than 1% and always below 12.20%. The first model F1 solves the majority of the instances to optimality. For seven cases, the model provides feasible solutions with an optimality gap below 14.02% and for three instances no feasible solution is found. Table 7 reports the metaheuristic results for the instances with reduced fleet. Clearly, the fleet size decrease makes the problem less computationally tractable as confirmed by the Gurobi average solution time (2881 s) which is almost twice the average solution time reported in Table 5. Of course, this also affects the heuristic solution time which is higher, but always less than 9 min.
Both the splitting and genetic algorithms provide optimal solutions for instances up to 25 customers. Again, the random mechanism provides solutions with slightly lower gaps (6.34% versus 6.48%). Concerning the speed up ratio, the genetic algorithm without the random mechanism is faster.
To sum up, for small instances the proposed metaheuristics provide optimal solutions in neglectable solution time and for larger instances with reasonable fleet size, we obtain quite good feasible solutions with low solution times, but still in some cases the relative gap is quite high.

Conclusions
In this study, we introduced the MDk-TRP, generalizing the k-TRP to the multi-depot case, for which we presented two effective mathematical models based on a layered network structure. We investigated the viability of solving the mathematical models by using off-the-shelf software. With respect to other two formulations adapted from the multi-depot CCVRP, the efficiency of the proposed models, in terms of both solution and time, is very satisfactory. In particular, the formulations solved to optimality instances up to 240 customers within two hours, and found feasible solutions for instances of up to 249 customers. A population-based metaheuristic was designed, including a splitting heuristic for generating good initial solutions and a local search.
The algorithm provided feasible solutions for the instances that could not be solved by any of the models in competitive computational time.
Some interesting extensions of this paper can be investigated. For instance, since we observed that the LP relaxation of the proposed formulations provide good lower bounds, this can be a good starting point to develop an exact or a matheuristic approach. Another research stream that is worth following is the incorporation of stochasticity into model parameters, such as travel time uncertainty or a random presence of customers.