The LRPH could be represented by the following undirected graph problem: Let be a complete undirected graph, where is the set of vertices representing customers and potential depots, and is the set of edges. Vertices correspond to set of customers each with known positive demand di, and vertices correspond to the potential depots, each with capacity wi and opening costs oi. With each edge is associated a non-negative traveling cost ^_cij . At each potential depot , a set of K heterogeneous vehicles is located, each with a vehicle capacity qk. Furthermore, any vehicle that performs a route generates a fixed cos vkt. The goal of LRPH is to determine the depots to be opened, the customers to be assigned to each open depot, and the routes to be performed to fulfill the demand of the customers by considering a heterogeneous fleet. Indeed, the LRPH integrates a strategic decision (depots to be used) and an operational decision (the routes to be developed by considering a heterogeneous fleet).

The objective function of LRPH considers the minimization of the sum of the fixed cost of the depots, of the costs of the used vehicles, and of the costs related with the distance traveled by the vehicles [4]. Any LRPH solution is feasible, if the following constraints are satisfied: (i) Each route starts and finishes at the same depot, (ii) each customer must visited exactly once by one vehicle, (iii) the sum of the demand of the customers served by a vehicle must not exceed its corresponding vehicle capacity ^_qk , (iv) the sum of the demand of the customers assigned to an open depot must not exceed its corresponding depot capacity ^_wi , and (v) the flow of products between depots is not allowed.

The granular search conception (Toth and Vigo [14]), is predicated on the employment of a sparse graph (incomplete graph) containing the edges incident to the depots, the edges belonging to the best feasible solutions found so far, and the edges whose cost is smaller than a granularity threshold , where is the average cost of the edges belonging to the best solution found so far, and β is a dynamic sparsification factor which is updated during the search ([4-5] and [14]). In particular, the search starts by initializing β to a small value β ₀. After iterations the value of β is increased to the value ^_βn , and / additional iterations are performed by considering as current solution the best feasible solution found so far [7]. Conclusively, the sparsification factor β is reset to its previous value β ₀ and the search continues. β ₀, β _n, N _s, and N _r are given parameters. The main idea of the granularity is to obtain high quality solutions in based-trajectory heuristics within short computing times [6,7].

The former heuristics use well-known intra-route (moves performed in a performed route) and well-known inter-route (moves performed between two routes assigned to the same depot or to different depots) moves corresponding to five neighborhood structures : Insertion, Swap, Two-opt, Double-Insertion, and Double-Swap [4-7]. A move is performed only if all the new inserted edges in the solution belonging to the sparse graph (granular search space).

Some of the proposed approaches allow infeasible solutions respect to the depot and vehicle capacity constraints. For any feasible solution S, we calculate its objective function value as the sum of the fixed cost of the depots, of the used vehicles and of the cost of the edges traveled by the performed routes. On the other hand, if the solution S is infeasible, we calculate its objective function value , where is a penalty term obtained by multiplying the global over vehicle capacity of the solution S times a dynamically updated penalty factor , and is a penalty term obtained by multiplying the global over depot capacity 𝑆 times a dynamically updated penalty factor. In particular, and are adjusted parameters during the search, and is the objective function value of the initial solution S ₀ . If infeasible solutions with respect to the depot capacity have been not found over N _fact iterations, then the value of is set to, where. On the other hand, if feasible solutions have been not found during N _fact iterations, then the value of is set to, . A similar procedure is applied to update the value of. Note if the current solution S is feasible, (for further details see [6-7]).

Utilizing a hybrid heuristic based on a cluster approach, a good feasible initial solution 𝑆 0 is generated within short computing times. The following steps are executed until all the customers are assigned to one route and one depot:

Firstly, a giant TSP tour containing all the customers (without the consideration of the depots) is constructed by using the well-known Lin-Kernighan heuristic (LKH) ([16] and [17]). Secondly, starting from any initial customer 𝑗 ⋆ , divide the built giant TSP tour into several clusters composed of consecutive customers so that, for each cluster the vehicle capacity constraint is satisfied. In particular, we have assigned the large vehicles firstly.

Thirdly, for each depot 𝑖 and each cluster g, a TSP tour is obtained using procedure LKH to evaluate the traveling cost of the route performed starting from i and visiting all the customers belonging to g (keeping the sequence obtained by the giant TSP). Finally, the depots are assigned to the clusters by solving an ILP model for the Single Source Capacitated Facility Location problem. This step determines the depots to be opened and the clusters to be assigned to the open depots (for further details see [18] and [19]).

The former algorithm considers a standard implementation of the Simulated Annealing metaheuristic (SA) by using a granular search space (reduced graph). We have assumed the same notation used on similar problems such as LRP [4]. Let S* be the best feasible solution found so far, S the current solution (feasible or infeasible), S’a random solution obtained from the neighborhood solutions of the current solution S, α the cooling factor, and 𝑇 the current temperature. Initially, we set , and . In addition, we set the initial temperature T0 and the number of iterations . The proposed algorithm performs the following steps until the number of iterations is reached:

Decrease the current temperature every Ncool iterations (where Ncool is a given parameter) by using the following function , where 0 < α < 1; and also increase the number of iterations .

Construct a random solution S’ obtained by considering all the neighborhood structures from the current solution S described in section 2.3.

Compute

Generate a random number 𝑟 in the range [0, 1];

If ;

If , set S: = S’; otherwise, keep S.

Finally, the best feasible solution found so far S* is kept. The pseudocode algorithm of GSA is presented as follows:

The GVNS algorithm brings together the potentiality of the systematic changes of neighborhood structures proposed by the well-known Variable Neighborhood Search (VNS) [20], and the efficient Granular Search Space introduced by [14] and improved by [5-7]. According to [20], the VNS applies a search strategy based on the systematic change of the neighborhoods structures to elude local optima. Three main concepts are applied on a VNS: (1) All the local minimum obtained by different neighborhood structures are not necessarily equals; (2) The best local minimum (respect to the objective function) obtained from all possible neighborhood structures (described in section 2.3) is called global minimum; (3) The different local minima obtained from the neighborhood structures should be relatively close each other.

Once the initial solution is performed (S0), the VNS approach iterates through different neighborhood structures to amend the best feasible solution (S*) found so far, until a the number of iterations is reached. The algorithm starts by setting , where S is the current solution.

The Variable Neighborhood Search considers two steps: (1) selecting a random solution from the first neighborhood and (2) applying a Granular Search Space by the same exchange operator until there is no more improvement. Then, the algorithm selects another neighborhood and the search continues.

The pseudocode algorithm of GVNS is presented as follows:

The proposed algorithm is an extension of the proposed idea by [21] for the DCVRP. After the construction of the initial solution So, the pGTS algorithm iterates through different neighborhood structures (described in Section 2.3) by using a discrete probabilistic function to improve the best feasible solution (S*) found so far, until the number of iterations is reached. The algorithm starts by setting , where is the current solution (feasible or infeasible), and is the current feasible solution. The following steps then are repeated sequently. First, the former algorithm selects a neighborhood from the neighborhoods structures described in Section 2.3 by using the following function of probability , where u is the total number of neighborhoods. Second, we apply a granular tabu search (GTS) based on the idea proposed by [14] in the selected neighborhood until a local minimum S’ is found. Depending on the solution, the following choices are possible:

Increase the probability of selecting the current neighborhood (Nk) by a given factor Pinc as follow , only if any of three cases occurs: (1) S’is infeasible and is feasible and , and (3) S’ is feasible and . If (1) is performed, then , if (2) is found, set . Finally, if (3) occurs, set.

In order to preserve the probability function properties after decreasing or increasing a certain neighborhood, the values of the probability to select other neighborhoods, where, must be adjusted [21]. If the probability to select the neighborhood Nk is decreased in a Pdec value, the remaining value is distributed for the remaining neighborhoods according to their current (probability) [21]. Therefore, the new probability to select the remaining neighborhoods () is calculated as follows:

(1)

where is the previous probability of the corresponding remaining neighborhood (). If the probability to select the neighborhood N _k is increased in a P _inc value, the remaining value is distributed for the remaining neighborhoods according to their current (probability). Therefore, the new probability to select the remaining neighborhoods () is calculated as:

(2)

Finally, the best feasible solution found so far S* is kept. The algorithm explores the solution space by moving at each iteration, from a solution to the best solution in the neighborhood , even if it is infeasible. The selected move is declared as tabu. The tabu tenure is defined as a random integer value in the range, where t _min and t _max are given parameters [21].

The pseudocode algorithm of pGTS is the following:

The pseudocode of the pGTS shows a brief summary of the performance of the proposed algorithm according to the discrete probability function and the replacement of the solutions depending of the characteristics of the solution found by the neighborhood S’.

A suitable set of parameters, whose values are based on extensive computational tests on the benchmark instances, was selected for each algorithm. The parameters and values are presented in Table 1.

These values have been utilized for the comparison of the solutions obtained by the described algorithms.

The calibration of the value of each parameter was performed by a multi-objective optimization approach (considering the minimization of the objective function and the computing time). Values coming from previous works with similar algorithms ([4,5-7,21]) are used. The next step is to select a parameter and search for the given value the best result. This search is executed applying 1D function minimization. The process is carried out for each considered value. Since this process is iterative, it can be refined using more repetitions.

For each instance, the proposed algorithms are executed 5 times due to their random calculations. Tables 2 - 5 provide the detailed results of the three proposed algorithms on the four data sets respectively. The algorithms are compared based on their best and average objective function for each instance; and the CPU for obtaining the best result and the CPU time to process the five runs of each instance. As expected, the GSA algorithm is faster, in the most of the cases, than the other algorithms since a solution is selected randomly and then evaluated; while GVNS performs local search on the selected solution, and pGTS explores the granular search space. A remarkable fact of the three algorithms is that although using random numbers, the objective function value tends to converge (see Tables 2 - 5).

The results clearly show that algorithm GSA outperforms the other two algorithms for what concerns both the CPU time and the quality of the answers found. Indeed, for all the data sets, the average costs, and the values of the best results for GSA are better than the corresponding values of algorithms GVNS and pGTS. Therefore algorithm GSA is the best performing of the three described algorithms, and, it could be compared with the most effective heuristics will publish in the literature.