MATH-HEURISTIC FOR THE CAPACITATED TWO-ECHELON VEHICLE ROUTING PROBLEM

ABSTRACT In this paper, we propose a math-heuristic that combines mathematical programming techniques with a heuristic approach based on the Simulated Annealing metaheuristic to solve the Two-Echelon Capacitated Vehicle Routing Problem (2E-CVRP). This problem arises in the context of a distribution network that is divided in two levels: satellite facilities that connect customers to fulfilment centers where freight originates. As it is an NP-hard problem, the proposed approach combines a cluster-first route-second math-heuristic in which approaches are more appropriate, particularly for problem sizes that are more commonly found in practice. The results of the experiments with benchmark instances show that such cluster-first route-second math-heuristic approach utilizing package solvers (CPLEX and TSP CONCORDE) can effectively help solving the CVRP for small instances when compared to an exact method. The experiments conducted on benchmark instances demonstrated the effectiveness of the proposed “cluster-first, route-second” math-heuristic approach, which utilizes package solvers such as CPLEX and TSP CONCORDE, in solving the CVRP for small instances, outperforming exact methods. This research contributes to demonstrating the potential applications of package solvers on heuristic structures for solving the CVRP. Although the presented math-heuristic has limitations, mainly due to the extensive usage of mathematical programming and the chosen characteristics of the implemented local search operators, it can quickly generate high-quality initial solutions for medium and large instances. By showcasing the “cluster-first, route-second” approach, this paper provides a framework that can be further improved by local search or embedded in other metaheuristics, such as GRASP or tabu search, and has practical implications for various industries.


INTRODUCTION
This article addresses the two-echelon capacitated vehicle routing problem (2E-CVRP).In this problem, the distribution network is divided in two levels or layers (echelons), as can be seen in Figure 1.The main depot and the intermediate facilities (satellite terminals) correspond to the main distribution structure, from where loads originate, and are usually located away from the large urban centers, configuring the first echelon or tier, denoted by -FE or 1E -with customers are serviced in the second echelon -SE or 2E -as described by Sluijk et al. (2022).The satellites, located at the edges or borders between layers, allow the consolidation and transshipment of products between the vehicles utilized on either echelon.This makes it possible to use larger vehicles in the first tier, which can transfer the loads from the depot to be delivered to the intermediate installations (satellite terminals) at times of lower traffic or that allow the traffic of larger vehicles (for example, in the evening or night), which are then distributed to the delivery points using smaller vehicles, usually more suitable for driving and parking in the most central and congested areas.
As highlighted by Crainic et al. (2009), two-tier systems are increasingly being used in large cities that face serious problems related to vehicle traffic.Examples of two-tier distribution systems include home delivery of e-commerce, distribution of supermarket products, multimodal cargo transport, among others (Sluijk et al., 2022).Solutions described in the literature present exact methods (or algorithms) and heuristic methods.State-of-the-art algorithms can solve only relatively small instances to optimality in reasonable computational time since the 2E-CVRP is of the NP-hard class (Arnold and Sorensen, 2019;Cuda et al., 2015).Heuristic approaches do not necessarily provide the optimal solution to the problem but allow obtaining near-optimal results for a NP-hard problem usually in very short times (Arnold and Sorensen, 2019).A math-3 heuristic combines heuristic techniques with operations research methods and tools to obtain effective solutions (Boschetti and Maniezzo, 2022).This paper proposes the development of a math-heuristic based on the Simulated Annealing metaheuristic for the 2E-CVRP that derives from formulation developed by Sluijk et al. (2022).
The 2E-CVRP formulation proposed by the authors allows solving both the first and the second echelons of the routing problem.In this work, we propose a model that comprises a mathheuristic for solving the 2E-CVRP by applying a "cluster-First, route-Second" approach, commonly utilized in constructive heuristics to the CVRP (Prins et. al., 2014).With this approach, customers are first allocated to the satellites thus decomposing the problem into several Capacitated VRPs and it provides the main benefit of being computationally cheaper.However, the routes obtained later can be suboptimal (Singh et al., 2023).In the proposed math-heuristic two CVRP problems are solved after the clustering process (which defines the customers and load assigned to each satellite): the capacitated vehicle routing problem from the central depot to the satellites and the capacitated routing problem within each satellite.
Our motivation stems from the importance of demonstrating possible applications of package solvers in combination with heuristic approach to solve the CVRP.Although our proposed mathheuristic has limitations due to the usage of mathematical programming and chosen specific characteristics of the implemented local search operators, it aims to illustrate how the "cluster-first, route-second" approach can be used to quickly generate good-quality initial solutions that can be further improved by local search.Additionally, our approach can be integrated into a multi-start constructive heuristic such as GRASP (Feo and Resende, 1995) or even in a metaheuristic such as tabu search (Glover, 1990) or VNS (Hansen and Mladenović, 2001).This article is organized as follows: in the next section the literature review of the problem is presented.The proposed heuristic is described in Section 3 and the results of computational experiments are reported in Section 4. Finally, the conclusions are in Section 5.

LITERATURE REVIEW
In the classic vehicle routing problem (VRP), the goal is to determine the route(s) to visit a set of service demand nodes, seeking to minimize the distance or the total cost and ensuring that each node is visited exactly once and the total demand on any route does not exceed the capacity of the respective vehicle.In other words, demands -pickups, deliveries, or service visits -must be allocated to each vehicle of an available fleet, and the routes (or sequence of stops) of each vehicle must minimize the total cost of the service, usually composed of the weighted sum of fixed costs and costs proportional to the distance traveled by vehicles and travel time, subject to constraints such as time windows (time periods in which the deliveries can be made), maximum distance or maximum duration of routes, among others (Cunha, 2006).It is a widely explored problem in the literature.Good recent references that provide a broad view of the state-of-the-art of this problem can be found, for example, in Vidal et al. (2014), Lahyani et al. (2015), Braekers et al. (2016) and Vidal et al. (2020).
The first formal definition of a 2E-VRP appeared in Crainic et al. (2009), where the authors study a rich variant of a 2E-VRP problem with various products and depots, time dependencies and vehicle synchronization, in a broader context of city logistics.More specifically, the problem studied included the selection of routes and the scheduling of the dispatched vehicles, as well as the selection of demand routes from the main terminals, via satellite terminals, to the final customers.The coordination and time-based synchronization of the vehicles at both levels were central elements of the problem.The authors investigated main aspects of the problem and proposed a mathematical formulation and a general methodology for its resolution based on a hierarchical decomposition approach.No computational experiments or application to any practical case have been reported.Perboli et al. (2011) were the first authors to formally define the term 2E-VRP and propose a mathematical model for the problem with a single deposit from which vehicles depart at the first level.Cuts derived from VRP formulations were proposed to improve the linear relaxation of the mathematical model.Two problem decomposition heuristics that seek to optimally determine the allocation of clients to satellite terminals -allowing the resolution of each of the scripting problems independently -based on the mathematical model (math-based heuristics) were also presented.Computational experiments comprised four sets of test instances with up to 50 customers and a single warehouse.The results showed that the gap of the mathematical model is small up to 32 customers, increasing in problems with 50 customers and particularly in those with four satellites.The two heuristics based on mixed-integer programming presented good performance both from the computational point of view and from the quality of the solutions obtained.Jepsen et al. (2013) proposed an arc-flow-based formulation for 2E-VRP, while Baldacci et al. ( 2013) presented a route-based formulation.Both can resolve small problem instances.Marujo et al. (2018) proposed a method to identify the economic (on costs and service levels) and environmental impacts of using motorized cargo tricycles in combination with conventional trucks -acting as mobile depots -as a solution to address the constraints imposed on access of conventional cargo vehicles in densely populated areas.The analytical method in conjunction with Monte Carlo simulation for eco-efficiency evaluation was applied to the megacity of Rio de Janeiro.The results showed that the use of cargo tricycles that replenish from parked trucks serving as mobile warehouses does not result in a significant increase in the time or distance traveled or decrease in the level of service experienced by customers, while greenhouse gas emissions and local air quality pollutants can be significantly reduced using load tricycles and mobile deposits in the last mile distribution.Li et al. (2021) addressed the problem of route optimization in a series of synchronized express service steps (line haul, pickup and delivery), in which satellite terminals are classified in two types: satellite terminals for cargo pickup, called source satellites, which serve only senders; and satellite terminals for cargo deliveries, so-called destination satellites, serve only the receivers.At the second level, a heterogeneous fleet of small-capacity trucks meets senders' requests and collects cargo for transport to a source satellite.Upon arriving at the destination satellite termi-nal -connected to a source satellite through the intercity line haul -the loads are transshipped and loaded into trucks scheduled to arrive on time to the recipients.In the first level, the loads are loaded in large capacity trucks and transported between the satellite terminals.To ensure the on time express services, synchronization constraints on the source satellites and destination satellites must be assured.The goal is to find the set of routes that minimizes the sum of costs.The authors proposed a mixed-integer programming model in which satellite bi-synchronization constraints provide an innovative method for formulating two-level network routing problems.A strategy was also developed based on the ALNS meta-heuristic ("adaptive large neighborhood search") that uses different operators of destruction and repair of solutions.The experiments involved small (three source satellite terminals, three destination and up to 9 pickup or delivery points) and large (17 source and destination satellite terminals and 120 pickup and delivery points).The exact method was able to solve only part of the small instances, whereas ALNS was able to obtain quality solutions (deviations of less than 0.5% for those instances) in reduced processing times.In a related subsequent work, Li et al. (2022) addressed the 2E-CVRP with constraints on simultaneous grouping and simultaneous pickup and delivery, a new variant of the classic 2E-VRP problem in which customers from the same administrative region are served by vehicles from the same satellite, to ensure service consistency, and the pickup and delivery are carried out simultaneously at the second level.To optimally solve the problem, a path-based formulation was proposed, and a specialized branch-and-cut-and-price algorithm was developed, which also allows to solve two related variants: the 2E-VRP with grouping constraints and the 2E-VRP with simultaneous collection and delivery.Computational experiments with three groups of test instances were conducted to evaluate the performance of the algorithm proposed in the problems.The results show that the proposed dominance rule can significantly reduce the number of labels generated and all valid inequalities have a major impact on the robustness of the path-based model, making it competitive when compared to existing exact algorithms.Liu and Jiang (2021) describe a variable neighborhood search (VNS) algorithm to solve the 2E-VRPSDP, a variant of the two-echelon vehicle routing problem with simultaneous pickups and deliveries in which both feasible and infeasible solutions can be explored.Arnold and Sorensen, 2019 demonstrate the importance of well-implemented local search and how they are sufficient to create heuristics to obtain high-quality solutions in short processing time.
Good reviews on the 2E-VRP can be found in Cuda et al. (2015) and Sluijk et al. (2022).The latter comprises 40 new articles that have been published in a relatively short time frame since the previous review; additionally, it includes variants of problems, inspired by real world applications, which have been studied to fill the gap between theoretical and canonical models and practical applications.The authors identified several challenging extensions of the 2E-VRP, including multiple warehouses, time windows, simultaneous pickup and delivery operations, heterogeneous fleet of vehicles at each of the two levels, and the inclusion of delivery alternatives such as direct warehouse deliveries, without going through any satellite, for customers located near the origin of the cargo.The synchronization of routes between the two levels and split deliveries at the first level are important aspects highlighted by the authors, especially for distribution systems for e-commerce.Finally, it is also important to highlight that the literature on the subject also comprises the multi-level distribution network design problem, which also considers decisions regarding the location of the facilities.Among them, in Winkenbach et al (2015), the total distance covered in the itineraries is made through a simplified analytical method derived from Daganzo (2004) to reduce the overall complexity of the problem.The work of Janjevic et al. ( 2021) is a good reference on the state-of-the-art models of localization and routing in multi-layered networks.

PROPOSED MODEL
The proposed model comprises a math-heuristic for solving the 2E-CVRP by applying a "clusterfirst route-second" method.With this approach, customers are first allocated to the satellites, thus decomposing the problem into several Capacitated VRPs.Two CVRP problems are solved after the clustering process (which defines the customers and load assigned to each satellite): the capacitated vehicle routing problem from the central depot to the satellites and the capacitated routing problem within each satellite.The math-heuristic steps (Figure 2) are: • Step 1: builds an initial cluster of demand nodes around the existing satellites respecting the existing capacity constraints. • Step 2: solves the first echelon CVRP based on the clusters obtained in Step 1.
• Step 3: solves the second echelon CVRP from clusters formed in Step 1.
• Step 4: and Step 5 solve the 2E-CVRP by grouping all demand nodes in each single satellite if capacity constraints allow.These two steps are intended to tackle scenarios where the optimal solution is located within a single satellite, and Step 1 fails to detect a unique cluster..
• Steps 6 and 7 aim to determine the solution with the lowest cost among the results obtained in Step 3 for each individual satellite, as well as the results obtained in Step 5 (FE and SE) for each single satellite • In Step 8, the chosen solution from Step 7 is subjected to the Simulated Annealing metaheuristic with the aim to improve it...

Math-heuristic model
The model comprises eight steps and aims to solve the 2E-CVRP using a "cluster-first, routesecond" approach.Figure 2 depicts all the steps of the math-heuristic, while its variables and parameters are defined in Table 1.x ij Binary decision variable, xij = 1 if demand node i is assigned to satellite cluster j, 0 otherwise.
Step 1 Step 2 city freighter s Number of city freighters required per each satellite cluster s∈ S.
Step 1 z j Binary decision variable, z j = 1 if satellite j has 1 or more demand nodes assigned to its cluster (satellite j is active), 0 otherwise.
Step 1 Step 2 load depot s Load assigned to depot s∈ S.
Step 1 Step 2 x ijk Binary decision variable, Step 2 Step 6 u sk Subtour elimination variable in Step 1 -w sk Load transported from main depot {0} to satellite s ∈ S by vehicle k ∈ K1.
Step 2 Step 6 y ijs Binary variable, y i js = 1 if vehicle dispatched from s ∈ S uses arc(i, j) ∈ A2, 0 otherwise.

-
Step 6 f ijs Load transported from s∈ S using arc(i, j) ∈ A2. - Step 6 Figure 2 -Description of proposed model.

Step 1 -Build initial satellite clusters
This step attempts to cluster the demand nodes around each satellite by applying a k-medoids approach (Kaufman and Rousseeuw, 2005).It will split the demand nodes into s satellite clusters, where the number s is the number of satellites.
Step 1 consists of an implementation in CPLEX solver of the following mathematical model: Model constraints: Pesquisa Operacional, Vol.43, 2023: e270829 The objective function (1) minimizes the sum of the square of the distances from each demand node to the satellites (the square of the distance is used to increase the penalty and reduce the sensibility of the formulation when applied on the solver if a node has approximately the same distance from the satellites).Constraints (2) ensure that all demand nodes must be assigned to only one satellite cluster, while constraints (3) impose that the total satellite cluster demand must respect satellite capacity.Constraints (4), ( 5) and ( 6) bound the total number of city freighters at the second echelon.Constraint (7) states that a satellite cluster j is opened if there is at least one demand node assigned to it.Expressions ( 8), ( 9) and ( 10) correspond to the non-negativity constraints.

Step 2 -Solve the CVRP for the first echelon
The objective of this step is to solve the first echelon routes -given the demand clusters around the satellites from Step 1 -so that the satellites can be loaded from the main depot with the required demand at minimum cost.The formulation of Step 2 is based on Sluijk et al. (2022), being the main difference the inclusion on the constraint (18) to allow finding a solution to the FE before the SE is solved by constraining the total demand in each satellite according to the clusters obtained in Step 1.
Step 2 consists of an implementation in CPLEX solver of the following mathematical model: Model constraints: The objective function (11) aims to minimize the route costs from the main depot to the satellites and the handling costs.Constraints ( 12) limit the number of available vehicles at the first echelon.Constraints (13) assure the flow of conservation at the first echelon.Constraints ( 14) impose that each satellite can be visited at most once by a first echelon vehicle.Constraints (15) are the subtour elimination constraints.Constraints ( 16) and ( 17) establish the first echelon vehicles capacity constraints, while constraints (18) assure that Step 1 cluster demands must be respected.Equations ( 19), ( 20) and ( 21) are the non-negativity constraints.

Step 3 -Solve the CVRP for the second echelon
The objective of this step is to solve the second echelon routes -given the demand clusters around the satellites from Step 1 -so that the demand nodes can be serviced from each satellite at minimum cost.The formulation of Step 3 is based on Sluijk et al. (2022) being the main difference the inclusion on the constraint ( 23) to allow finding a solution to the SE after the FE is solved by constraining the total demand in each satellite according to the clusters obtained in Step 1.
Step 3 consists of an implementation in CPLEX of the following mathematical model (this step solves the second echelon as a monolithic block): Model constraints: x i jk ∈ {0, 1} , ∀ (i, j) ∈ A1, ∀k ∈ K1 (33) The objective function ( 22) aims to minimize the route costs from the satellites to the demand nodes.Constraints (23) assure that the total load from each satellite must respect the satellite cluster load from Step 1. Constraints ( 24) and ( 25) impose that there is no flow between satellites.Constraints (26) assure the flow conservation in the second echelon.Constraints ( 27) limit the number of vehicles per satellite.Constraints (28) limit the number of vehicles in the second echelon.Constraints (29) assure that each demand node is visited exactly once.Constraints (30) establish that the net flow at each demand node is equal to its demand.Constraints (31) assure that the flow from each satellite must comply with the second echelon vehicle capacity constraints.Equations ( 32), ( 33) and ( 34) are the non-negativity constraints.

Steps 4 and 5 -Group demand in single cluster
In Step 4, if capacity constraints allow -second echelon freighters capacity, total number of second echelon freighters available, number of freighters allowed per satellite and satellite capacity -the demand is assigned for every single satellite.Then Step 5 applies Step 2 and Step 3 for every single cluster identified in Step 4. These two steps are designed to address the situation where the optimum solution resides within a single satellite and Step 1 does not find a single cluster.
3.1.5Steps 6 and 7 -Build solution set and obtain initial feasible solution Step 6 and Step 7 identify the minimum cost solution between the solution (FE and SE) obtained in Step 3 and each solution (FE and SE) obtained in Step 5 for each different satellite.

Step 8 -Simulated Annealing
This step executes the Simulated Annealing (SA) metaheuristic to the chosen solution in Step 7 as detailed in the pseudocode depicted in Figure 3.
S0 is the initial feasible solution obtained from Step 7. S* is the best solution and S is the incumbent solution during the SA iterations.T is the temperature parameter and Ti the initial process temperature.L(T) is the SA iterations parameter and alpha the temperature reduction rate.The neighborhood solution (S') is obtained by either two local search operators applied on the incumbent solution S: move or swap.
The move operator relocates a demand node from its cluster to another different cluster if capacity constraints allow (total number of vehicles available and vehicle capacities at the second echelon, maximum number of vehicles per satellite and satellite capacity).
The swap operator exchanges the satellites of two demand nodes in different satellites if capacity constraints allow (total number of vehicles available at the second echelon, maximum number of vehicles per satellite and satellite capacity).Operators move and swap are selected randomly with probability of 50% as described by the pseudocode in Figure 4.The local search operators (move or swap) perform as described by the pseudocode in Figure 5.

Second echelon math-heuristic variation
An alternative math-heuristic to solve the second echelon (Step 3) was also developed using the TSP Concorde algorithm libraries (Applegate et al., 2003) as described by the pseudocode in Figure 7.
In this math-heuristic, set V corresponds to the set of satellite vehicles, whose cardinality |V| is equal to the round up to the nearest integer of the division of the total satellite demand by the freighter capacity parameter.The mathematical modeling of the cluster vans step of the mathheuristic is: Model constraints: z j ∈ {0, 1} , ∀ j ∈ S (39) The objective function (35) aims to cluster the demand nodes in the minimum required set of vehicles.Constraints (36) state that each demand node must be allocated to a single second echelon vehicle.Constraints (37) assure that a route (vehicle) cluster j is opened if there is at least one demand node assigned to it.Constraints (38) impose that the total load allocated to a vehicle must be within freighter capacity.Equations ( 39) and ( 40) are the non-negativity constraints.Once the demand nodes are assigned to routes -the vehicle clusters are identified -the initial routes are obtained by applying the TSP Concorde algorithm to each vehicle cluster.Finally, within each satellite cluster a Simulated Annealing process as described in 3.2.6 is applied between the routes (or vehicles): demand nodes are moved or swapped between the same satellite cluster routes with 50% probability.

COMPUTATIONAL EXPERIMENTS
The proposed models were implemented in C++ using CPLEX 22.1.0.The second echelon mathheuristic variation (section 3. given in https://homepage.univie.ac.at/ulrich.breunig/research/.The best-known solutions (column BKS), in all the result tables presented below, marked with an asterisk are optimal solutions.Non marked best-known solutions have not been proven optimal (Breunig et al., 2016).
The chosen instances can be considered small -up to 100 customers -and medium -from 101 to 350 customers (Alesiani et al., 2022).The problem instances were solved using an Intel i7 computer with 32GB of memory running Linux operating system.Steps 1-3 and 5 correspond to mathematical models that were implemented in CPLEX as described respectively above in subsections 3.1.1,3.1.2and 3.1.3,respectively.
An exact model, solving concurrently the first and second echelon, implemented in Gurobi, using the mathematical formulation presented in Sluijk et al. (2022) was developed to certify the implemented math-heuristic model-since it uses part of the original formulation -and to check the processing time of the math-heuristic when compared to the exact monolithic approach.column runtime (running time in seconds); column 'deviation' the percentage difference from the cost obtained from the math-heuristic and the best-known solution for the instance.Instance sets named Set#2, Set#3 and Set#4 were solved using the math-heuristic for the second echelon as described in 3.2.3.
The following parameters were used in the test instances to stay within the computational time comparable to the methods analyzed in Sluijk et al. (2022): solver gap limit equal 0.005% and solver time limit = 600 seconds, local search iterations limit L(T) = 10 and simulated annealing iterations limit = 10,simulated annealing initial (Ti) temperature = 1000 and final temperature (Tf) = 0.01, temperature reduction ratio (alpha)=0.9and simulated annealing time limit = 1800 seconds.The size of the test instances contained in datasets Set#5 and Set#6 required the use of the math-heuristic described in 3.3 since no feasible quality solution was obtained after three hours of processing with the second echelon math-heuristic described in 3.2.2.The following parameters were used for instance sets Set#5 and Set#6 to ensure the computational time comparable to the methods analyzed in Sluijk et al., 2022: solver gap limit = 0.005% and solver time limit = 600 seconds, local search iterations limit = 10 and simulated annealing iterations limit = 10, simulated annealing initial temperature = 1000 and final temperature = 0.01, temperature reduction ration (alpha) = 0.9, simulated annealing process time limit = 3,600 seconds.The used intra-route parameters are local search iterations limit = 50 and simulated annealing iterations limit=50, simulated annealing initial temperature=1,000 and final temperature= 0.01, temperature reduction ration (alpha) = 0.9, simulated annealing process time limit = 300 seconds.Table 9 presents a summary of the results obtained for all the test instances.Column 'Qty Sat.' denotes the number of satellites in the instance; column 'Qty.Clt' the number of demand nodes; columns delta the difference from the cost obtained from the math-heuristic and the best-known solution for the instance.Finally, columns 't(s)' correspond to the execution time in seconds.

CONCLUDING REMARKS
The main objective of the work was to develop a math-heuristic to obtain quality solution in computational time comparable to the methods analyzed in Sluijk et al. (2022).The math-heuristic allowed to high quality, optimum or near optimum solutions with faster processing time than the exact model for 30 and 50 customers instances when the number of maximum vehicles per satellite equals the second echelon number of available vehicles and the number of available vehicles in the second echelon is between two and five.
As the number of clients and satellites grow, and the number of second echelon vehicles restriction by satellite gets stricter, the search neighborhood grows and the math-heuristic with the chosen parameters to keep the processing time comparable to the approaches in Sluijk et al. (2022) did not reach good solutions, although all under 16.25% apart from the best-known solution.This is expected since the designed model relies heavily on the mathematical program-ming model resolution for the first echelon and second echelon (reason also that has driven the approach described in item 3.3 for Set#5 and Set#6 instances).Additionally, the two local implemented search operators (move and swap) between satellites and between routes for the alternative second echelon math-heuristic are intentionally restricted for movements that maintain the capacity restrictions, which prevents temporarily accepting unfeasible solutions to explore broader neighborhoods (Liu and Jiang, 2021) specially when the solutions are very tight (trucks are almost full and demand sizes do not allow the desired movements).The small number of local implemented operators (only two) also explains the poor results for medium and large instances (Arnold and Sorensen, 2019).
The results on benchmark instances show that such "cluster-first, route-second" math-heuristic approach utilizing package solvers (CPLEX and TSP CONCORDE) can effectively help solving the CVRP for small instances when compared to an exact method.The contribution of this paper lies in demonstrating the potential applications of package solvers to solve the CVRP on heuristic structures.
While the math-heuristic presented in this paper has certain limitations in solving medium and large instances, it serves to illustrate how the 'cluster-first, route-second' approach can be utilized to develop initial construction heuristic frameworks.This is due to the simplification of the original problem and the subsequent ability to gradually replace mathematical models with heuristic structures, leading to improved scalability.This work also motivates future research aimed at enhancing the simulated annealing approach by developing more efficient local search engines, effective destroy and repair methods together with local search operators to explore several neighborhoods (e.g., pruning or temporarily accepting unfeasible solutions) or alternatively devising population-based heuristics that are all state-of-the-art approaches while taking advantage of the current existing package solvers processing capacity.
j) | i, j ∈ {0} ∪ S, i ̸ = j } = set of the first echelon arcs connecting the main depot {0} to satellite j ∈ S. (i, j) | i, j ∈ S ∪ C, i ̸ = j } = set of the second echelon arcs, connecting satellite i ∈ S with demand node j ∈ C.

Figure 6 -
Figure 6 -Obtain new solution after local search operators.

Table 1 -
Math-heuristic parameters and variables.
Table 2 and Table3present the results for the exact model.The solver time limit was set to 36,000 sec- onds (10 hours).Column 'Qty Sat' contains the number of satellites in the instance; column 'Qty Clt' the number of demand nodes; column 'Dmd' the total second echelon demand; column 'Qty Tr' the number of first echelon vehicles; column 'Truck Cap' the capacity of each first echelon vehicle; column 'Max Van/Sat' the maximum number of second echelon vehicles per satellite; columns 'Qty Vans' the total number of second echelon vehicles; column 'Van Cap' the capacity of each second echelon vehicle; column 'BKS' the cost of the best-known solution for the test instance; column 'Sol' the cost obtained from the exact model; column runtime (solver execution time in seconds); column 'Gap %' the solver gap after the model execution (or execution time limit being reached) and column deviation the percentage difference from the cost obtained from the exact model and the best-known solution for the instance.The math-heuristic test instances results are presented in: Table4(test instances from Set#2), in Table5(test instances from Set#3), Table6(test instances from Set#4), Table7(test instances from Set#5) and Table8(test instances from Set#6).The results presented are an average of five different executions per set of test instances.Column 'Qty Sat' contains the number of satellites in the instance; column 'Qty Clt' the number of demand nodes; column 'Dmd' the total second echelon demand; column 'Qty Tr' the number of first echelon vehicles; column 'Truck Cap' the capacity of each first echelon vehicle; column 'Max Van/Sat' the maximum number of second echelon vehicles per satellite; columns 'Qty Vans' the total number of second echelon vehicles; column 'Van Cap' the capacity of each second echelon vehicle; column 'BKS' the cost of the best-known solution for the test instance; column 'SA' the cost obtained from the math-heuristic;

Table 4 -
Set#2 Results obtained by the math-heuristic.

Table 5 -
Set#3 Results obtained by the math-heuristic.

Table 6 -
Set#4 results obtained by the math-heuristic.

Table 9 -
Average results obtained by the math-heuristic.