A Multi-Start Algorithm for Solving the Capacitated Vehicle Routing Problem with Two-Dimensional Loading Constraints

This work presents a multistart algorithm for solving the capacitated vehicle routing problem with 2D loading constraints (2L-CVRP) allowing for the rotation of goods. Research dedicated to graph theory and symmetry considered the vehicle routing problem as a classical application. This problem has complex aspects that stimulate the use of advanced algorithms and symmetry in graphs. The use of graph modeling of the 2L-CVRP problem by undirected graph allowed the high performance of the algorithm. The developed algorithm is based on metaheuristics, such as the Constructive Genetic Algorithm (CGA) to construct promising initial solutions; a Tabu Search (TS) to improve the initial solutions on the routing problem, and a Large Neighborhood Search (LNS) for the loading subproblem. Although each one of these algorithms allowed to solve parts of the 2L-CVRP, the combination of these three algorithms to solve this problem was unprecedented in the scientific literature. In our approach, a parallel mechanism for checking the loading feasibility of routes was implemented using multithreading programming to improve the performance. Additionally, memory structures such as hash-tables were implemented to save time by storing and querying previously evaluated results for the loading feasibility of routes. For benchmarks, tests were done on well-known instances available in the literature. The results proved that the framework matched or outperformed most of the previous approaches. As the main contribution, this work brings higher quality solutions for large-size instances of the pure CVRP. This paper involves themes related to the symmetry journal, mainly complex algorithms, graphs, search strategies, complexity, graph modeling, and genetic algorithms. In addition, the paper especially focuses on topic-related aspects of special interest to the community involved in symmetry studies, such as graph algorithms


Introduction
The vehicle routing problem (VRP) is a class of computational optimization problem that involves designing delivery routes and lower logistic cost collection to satisfy the demands of a set of customers or a group of geographically dispersed cities [1,2]. The VRP generalizes the traveling salesman problem (TSP) [3]. Although the TSP is an old problem whose origin is not well known, many studies on it are still carried out, as is the case with [4]. The VRP is addressed in studies dedicated to graph theory and algorithms [5] CGA-TS approach for the CVRP, and the LNS for the 2L-CVRP. Section 5 presents the results of the testing and a comparison with other approaches found in the literature, as well as a statistical validation for the comparison. Finally, Section 6 shows conclusions based on the study and suggestions for future works.

Literature Review
In the field of the VRP, we can find a combination of the CVRP with 2D loading constraints, which is denoted as the 2L-CVRP. The first work on the 2L-CVRP was presented by Iori et al. [14]. They proposed an exact method based on the branch-and-cut algorithm to minimize the routing costs, and a nested branch-and-bound algorithm to solve the 2D sequential oriented loading subproblem (2|SO|L). Later, Gendreau et al. presented an approach using a metaheuristic to solve the 2L-CVRP [15]. They proposed a solution based on the tabu search algorithm. They applied heuristics, lower bounds, and a branchand-bound algorithm for the loading feasibility checks, considering both sequential-and unrestricted-oriented loading, 2|SO|L and 2|UO|L, respectively.
Fuellerer et al. [16] developed a metaheuristic approach based on ant colony optimization to solve the routing problem combined with heuristics for the loading subproblem. For the first time, the results allowed for the rotation of goods by 90°(nonoriented) when considering the 2L-CVRP. Thus, two new variants of loading constraints were studied considering the nonoriented version with sequential and unrestricted loading, 2|SR|L and 2|UR|L, respectively. Zachariadis et al. [17] proposed a metaheuristic algorithm based on the combination of the tabu search with a guided local search and a collection of packing heuristics for the loading subproblem. Strodl et al. [18] developed a solution for the routing problem using a variable neighborhood search and an exact branch-and-bound algorithm for the loading. A simulated annealing algorithm with a collection of packing heuristics to solve the 2L-CVRP was presented by Leung et al. in 2010 [19]. Later, Leung et al. [20] in 2011 also proposed an approach based on the tabu search combined with an extended guided local search and a collection of packing heuristics. Duhamel et al. [21] presented a study with greedy randomized adaptive search and evolutionary local search algorithms. Zachariadis et al. [22] presented a work with an algorithm referred to as Promise Routing-Memory Packing (PRMP) to solve 2|UO|L and 2|SO|L versions of the 2L-CVRP. Dominguez et al. [23] published a paper with a multistart biased randomized algorithm (MS-BR) to solve the 2L-CVRP with two unrestricted loading configurations (2|UO|L and 2|UR|L). This was the second paper proposed in the literature to solve the 2D loading constraint allowing for the rotation of items. Wei et al. proposed a variable neighborhood search algorithm combined with a skyline heuristic to solve 2|UO|L and 2|SO|L versions of the 2L-CVRP [24]. Recently, Wei et al. [25] proposed a simulated annealing algorithm with an open space based heuristic to check the loading feasibility, which outperformed all the previous approaches on 2|SO|L, 2|UO|L, 2|SR|L, and 2|UR|L. Related to the previous 2L-CVRP work cited, the open space based heuristic for the 2D strip packing problem was presented also by Wei et al. [26].
Many other constraints related to the real-world problems of the logistics industry were studied with the 2L-CVRP. Zachariadis et al. [27] proposed an algorithm based on an extended version of their previous approach [22] to solve the 2L-CVRP with simultaneous pickups and deliveries. They also considered the results for the 2|UR|L version and obtained good results. In this paper, their algorithm is referred to as the xPRMP. Pinto et al. presented a study using variable neighborhood search algorithms to solve the 2L-CVRP with mixed linehauls and backhauls [28].
The CVRP with 3D loading constraints (3L-CVRP) is a related problem that was studied extensively in recent years. The 3L-CVRP was studied by many authors who proposed metaheuristics like a tabu search [29], a guided tabu search [30], the ant colony algorithm [31], and a GRASPxELS [32]. Bortfeldt [33] presented a hybrid algorithm for the 3L-CVRP. Koch, Bortfeldt, and Wäscher [34] presented a hybrid algorithm to solve the 3L-CVRP with backhauls and time windows. Recently, Bortfeldt and Yi [35] proposed a hybrid algorithm to solve the split delivery VRP (SDVRP) with the 3L-CVRP. We can cite some surveys on this subject area from Iori and Martello [36] and from Polaris et al. [37].
Finally, we can cite some other TSP-and VRP-related studies, such as the orienteering problem presented for the first time in [38,39], who proposed a memetic approach to solve it. The green vehicle routing problem (GVRP) is an emerging research field, and a recent systematic literature review on this field can be seen in [40].

Problem Description
Zachariadis et al. [17] describe the 2L-CVRP as follows. Let G = (V, E) be an undirected graph, where V = {0, . . . , n} is the vertices set, in which vertex 0 represents the single depot and the other n vertices correspond to the customers, and E = {(i 0 , j 0 ), . . . , (i n , j n ) : i, j ∈ V, i = j } is the set of edges. Each edge (i, j) ∈ E has an associated cost c ij that corresponds to the cost for the transport from i to j.
There are v identical vehicles, each with a weight capacity equal to Q and a rectangular loading surface of width W and length L (related (x, y) coordinates). Let A = W × L denote the loading area. The demand of customer i (i = 1, . . . , n) consists a set of m i items, denoted as IT i , of total weight q i : item I ik ∈ IT i (k = 1, . . . , m i ) has width w ik and length l ik . Let a i = ∑ m i k=1 w ik l ik denotes the total area of the items demanded by customer i. The objective of the 2L-CVRP is to determine the minimum cost set of routes that satisfy the following constraints: (a) the quantity of generated routes does not exceed the number of available vehicles; (b) all routes must start and end at the central depot; (c) each customer i can be visited only once; (d) each customer i must be served by only one vehicle; (e) the total weight of items demanded for a route does not exceed Q, and (f) all loading must be nonstackable by the set of customers covered by a route into the W L loading surface of the vehicles.
Usually, a practical constraint can be imposed, considering the convenience of unloading: each time a customer i is visited, all items I ik must be unloaded so that no items of other customers are moved, when k represents the items of a particular customer. This version is called sequential loading (also referred to as LIFO), and it can be denoted as 2|SO|L when the rotation of items is not allowed. Figure 1 depicts an example of the 2|SO|L 2L-CVRP. The version of loading without the LIFO constraint is called unrestricted. According to Fuellerer et al. [16], four combinations can be listed as follows: 2|UO|L: unrestricted, oriented loading; 2|UR|L: unrestricted, rotation allowed loading; 2|SO|L: sequential, oriented loading; and 2|SR|L: sequential, rotation allowed loading. The following describes the employed solution approach to solve the 2|UR|L 2L-CVRP. Firstly, the routing component is determined in a multistart approach. It consists of a constructive procedure to generate an initial solution. It then continues to search to find a trajectory in the search space that is the core of the optimization method to produce a final solution. This master routing algorithm utilizes a loading feasibility method that generates feasible loading patterns.

The Multistart CGA-TS-LNS Algorithm
The proposed algorithm to solve the 2|UR|L version of the 2L-CVRP in this work makes use of the combination of heuristics derived from the genetic algorithm (GA) [41], the tabu search (TS) [42], and the Large Neighborhood Search (LNS) [43].
The CGA-TS-LNS algorithm was developed with a multistart approach, wherein each global iteration a new initial solution is generated from a small sequence of iterations of the CGA-based subalgorithm Immediately after the initial solution, the TS-based subalgorithm is used intensively to improve the solutions for the CVRP. The last phase of each global iteration, the constraints for the loading subproblem, are evaluated by applying the LNS heuristic-based subalgorithm for packing the items in the routes of the best solution found by the TS. At the end of each loop, the algorithm records the feasible and unfeasible routes in terms of loading in two hashtables, respectively. Thus, with these records it avoids generating solutions with unfeasible routes found in previous loops. It also avoids reprocessing routes already evaluated as feasible that are found again in the current best solution. The three-step process is repeated until a time limit is reached or if no feasible global solution is found. Figure 2

The Algorithm Initialization
It is necessary to define parameters when beginning the execution of the algorithm. Some of them are fixed with previously established values and others are calculated from the data of the tested instances, as detailed below. The parameters are linked to the strategies and to the techniques used in the developed method. The improvement strategies used in the routing phase and the initial parameters of the algorithm are detailed in the following.

Improvement Strategies and their Parameters
To intensify and diversify the local search in the tabu search algorithm, two strategies were employed to obtain better results to solve the routing problem. The first is a list of better solutions, called elite solutions, which is built during the execution of the algorithm [44,45].
The list of elite solutions serves to make random choices, with a certain degree of probability. At certain times during the execution of the TS, one of the best solutions found until that point becomes the current solution. In this way, the generation of new neighbors is promoted through an elite solution. This intensifies the search in regions that can be considered promising. Elite solutions can also be evaluated for loading feasibility in the final stage of the algorithm, but only if the best solution found in the current global iteration is not feasible. The size of the list of elite solutions is identified by the sizeEL parameter.
The second strategy for improving the routing problem deals with a dynamic variation in the number of neighbors generated from the current solution. The numNeighbors variable is calculated at the beginning of the algorithm's execution, right after obtaining customer data. Stochastically, during the iterations of the TS, this variable can have its value changed within a range of values close to the value initially calculated. This strategy is pseudocoded at line 15 of Algorithm 2 in Section 4.3.
A simple strategy to facilitate loading is applied during the routing phase. The first global iteration of the algorithm considers 100% of the vehicle's loading surface area. After each global iteration, with a 50% probability, the vehicle area considered for loading can be reduced by 0.05%; note that the minimum limit is 93%. The variable used to control this restriction is called percAreaVehicle. This strategy, applied during the steps of generating the initial solution and during the routing phase, tends to facilitate the feasibility of loading in the last stage of each global iteration of the algorithm.
All percentages used in the strategies described above were chosen after an amount of tests and observations, that is, the best results were found using these values.

Other Control Parameters
The maxGen parameter indicates the maximum number of generations, that is, the maximum iterations to be performed by the genetic algorithm to build the initial solution (input data-Algorithm 1). The size of the tabu list is fixed and is identified as sizeTL.
The maxNoImproveCnt parameter signals the TS algorithm when it is time to try one of the strategies (intensification/diversification) previously mentioned in Section 4.1.1. The maxNoImproveCnt2End parameter establishes one of the exit conditions for the TS during the routing phase, indicating the maximum number of consecutive iterations without improvement to end the TS and start the next phase. MaxLNSTime, on the other hand, indicates the maximum time, in seconds, for executing the LNS algorithm in the feasibility assessment phase of loading the routes of the best solutions. Finally, maxGTime defines the maximum time, in seconds, for the execution of the main algorithm in the search for the best solution.

Solution Structures and the Initial Solution
The initial solution is generated from a variant of the genetic algorithm, called the constructive genetic algorithm (CGA) proposed by [46]. Some researchers already used genetic algorithms to generate an initial population of solutions, such as [47], who used the GA and Push-Forward Insertion Heuristic (PFIH) [48] for this purpose. We decided to use the CGA due to its constructive nature and because we did not find its use in the VRP in the literature. Before detailing the construction of the initial solution, Section 4.2.1 shows, in detail, the structures used to represent and manipulate the solutions.

Solution Structures
A matrix structure allowed to represent and manipulate the solutions in the memory during the execution of the algorithm. This matrix structure focused on the improvement of performance, mainly in stage two of the algorithm where the tabu search is performed. Figure 3 shows a sample of the solution matrix (S) for the instance 0902 from [15]. In this instance, 8 vehicles are available, and their trips are represented in rows 1 to 8 of the matrix S. The first column, identified by 0, represents the trip id in rows 1 to 8, and the first row (row id = 0) is the row of the summarized data of the solution. In the first row of S, there are 25 customers to be visited, as shown in column 1. Column 2 represents the total cost of the solution. Column 3 shows that each vehicle has 48 unit of measurement capacity, and column 4 shows that the surface area for loading is 800. As row 0 represents the summarized data of the solution, columns 5 to c max (c max = 80) are marked as empty, i.e., with value −1. Considering this example, rows 1 to 8 have the trip data. These rows contain columns 1 to 4, which have the summarized data of each trip, representing the number of customers, the trip cost, the total weight, and the total area, respectively. For the trips, the customers' IDs are stored in columns 5 to c max of the S matrix, and the order in which they are stored represents the sequence of visits to the customers. Note, the depot is excluded, which is by default the start-and end-point in this study.

Generating the Initial Solution Using Constructive Genetic Algorithm
The CGA is efficient in finding optimal or near-optimal solutions when employed in a variety of problems. Some researchers already used CGA in their works and obtained excellent results. Among them we can mention [49,50]. One of the main differences of the CGA, when compared to that of a classical genetic algorithm, is its fitness process [46].
To construct the initial solution according to the 2L-CVRP definition, the CGA's principles must be employed. This clustering problem in graphs is stated as the search for partitions on the vertex set V in a predefined number of clusters normally indicated by the quantity of available vehicles. A general overview of the CGA framework is provided by the pseudocode in the Algorithm 1.
The size of the initial population P 0 is defined as parameter initPSize and currGen counts the number of generations (line 1 and 2 of Algorithm 1). To define the P 0 of individuals, all the arcs (i, j) comprised by the set of edges have their associated costs c i,j , and are sorted by the increasing value of cost c i,j to be electable as the new routes (one for each available vehicle) (i.e., the initial edge of clusters that in some way attract the other edges which participate in the representation). Figure 4 shows an example of the ranking matrix of Euclidean distances between the customers. //initialize the current generation number 2: define initPSize; 3: Generate an initial population P 0 with size initPSize; 4: Evaluate the population P 0 ; 5: while currGen < maxGen do 6: Select P; //P is the population of routes 7: Recombination P; 8: Evaluate P; 9: Evaluate the loading feasibility of P in terms of weight and area; 10: let currGen = currGen +1; 11: end while 12: return the best feasible solution found; Using the ranking of distances, the closest customers are combined to form the initial population P 0 (line 3 of Algorithm 1) ( Figure 5). This initial population creates new generations through an iterative process of evaluation until convergence criteria are met to reach an optimal solution. For each p ∈ P 0 , the fitness is calculated based on the distance between the customers. Equation (1) calculates the fitness (line 4 of Algorithm 1).

Fitness =
∑ distance in the route number of customers (1) From that point, until the stop condition is reached and a valid solution is found, the steps of selection, recombination (crossover operation), and evaluation are performed. This ensures that other customers are selected successively to be inserted into the routes, respecting the genetic operators necessary to maintain the adaptation characteristics acquired by previous generations.
The selections (line 6 of Algorithm 1) are performed according to the number of vehicles, which establishes the number of routes. Figure 6 shows the selection process of the CGA framework. Recombination is performed for each route, inserting the closest customers outside of the route. Figure 7 shows the recombination process of the CGA algorithm (line 7 of Algorithm 1). Then, the population is evaluated once again using the fitness Equation (1)    Since the feasibility of a solution is determined mainly by the load constraints of the problem, the proposed constructive algorithm checks the items of all the customers in each route that can be loaded into the vehicle when considering loading constraints, such as the maximum load surface and the maximum vehicle weight (line 9 of Algorithm 1). This strategy managed to successfully construct feasible initial solutions for all of the CVRP instances ( Figure 8). The stop condition is triggered at a predefined number of generations. The population increases after the initial generations, continues growing until reaching an upper limit, and decreases for higher values of the evolution parameter. The structure corresponding to the best problem solution must be kept in the process. Figure 9 depicts an example of the CGA framework steps, where the filled rectangle represents the central depot (vertex 0), the empty circles correspond to the customers, the arrows represent the routes (i, j), the dashed circles represent the vertices selected for insertion, and the dotted arrows show possible insertion edges. The vehicle loading constraint check is integrated into the insertion process.

Tabu Search Algorithm for Routing Problem
The TS is a well-known metaheuristic that is widely used to solve combinatorial analysis problems. The TS method stands out for being reasonably easy implement and for producing very satisfactory results. Several researchers [15,17,33,[51][52][53][54] successfully used this technique, either alone or combined with another method to solve some variant of the VRP.
According to Glover [42], the TS consists of a local search in the neighborhood of the current solution that was generated from movements carried out from that solution. To escape from local optima, some movements are considered prohibited. These movements are recorded in the memory. The structure that stores these movements is a finite-sized list called the tabu list. These movements are prohibited from being carried out until they cease to exist in the list. This takes place as iterations occur and the list is updated. The process is done iteratively; with each iteration, the best neighbor must be selected to become the current solution and generate a new neighborhood. This is true even if that neighbor is worse than the best solution found until that point. The process is repeated until a stop condition is reached.
In this work, the CGA and TS heuristics are applied to find better solutions to the CVRP. Algorithm 2 shows TabuSearch method that starts right after the CGA builds an initial solution for the CVRP. The initial solution becomes the first current solution of the tabu search and the best local solution as well. 10: while noImproveCnt2End < maxNoImproveCnt2End and currTime < maxTime do 11: Put currSolution in tabuList[], considering FIFO method when tabuList[] is full; 12: if noImproveCnt >= maxNoImproveCnt then 13: noImproveCnt = 0; 14: With a probability of 40% set currSolution = bestSolution, or with a probability of 30% set currSolution = one of the eliteList[] randomly or with a probability of 30% stay the same; 15: With a probability of 25% set r = random number between -n and +n, where n=(numNeighbors/2) and set numCandidates = numNeighbors + r or with a probability of 75% stay the same; 16: end if 17: Select a neighborhood structure NS randomly; 18: Generate candidateList[numCandidates] from currSolution using NS; Increase noImproveCnt and noImproveCnt2End by 1; 29: end if 30: if eliteList[] has empty entry or cost(currSolution) < worst cost(eliteList[]) then 31: Put currSolution into the eliteList[] at empty entry or replace the worst entry; 32: end if 33: end while A tabu list (tabuList) and a list of elite solutions (eliteList) are initialized to their predefined sizes at the beginning of the algorithm, as described in Section 4.1.1. The tabuList keeps the last solutions found to temporarily prohibit them from becoming the current solution. The eliteList stores the latest sizeEL best solutions found.
The eliteList has two purposes. The first is to provide, at certain times, the return of one of the best solutions previously registered as the next current solution. This purpose serves to intensify the search in promising regions of the search space [55]. The second purpose is for the best solutions found in the routing problem (CVRP) to be evaluated for loading feasibility (2L-CVRP) in the next stage of the framework. This is only carried out if the local best solution is not feasible.
A local search cycle is initiated until the stop condition is true. The main condition for stopping the tabu search is when a maximum consecutive number of iterations without improvement is achieved (noImproveCnt2End maxNoImproveCnt2End). Another stop condition is if the maximum global time minus the time for load evaluation is reached (currTime maxGTime˘− maxLNSTime).
At the beginning of each iteration, the current solution is placed on the tabu list. As the tabu list has a predefined size (sizeTL), when it is filled, the new entries follow a first in, first out method (FIFO). From the current solution, a neighborhood is generated by following a pattern of neighborhood structures as detailed in Section 4.3.1. The size of the neighborhood is initially defined by the numNeighbors variable. After the generation of the neighborhood, a local search is carried out to find the best candidate solution to become the next current solution. As already mentioned, in this choice, the neighborhood solutions cannot be on the tabuList and must have feasible routes in terms of weight and area. The best neighborhood solution is then selected and becomes the new current solution for the next iteration, even if it is worse than the most recent current solution.
To speed up the loading evaluation process, two hash table structures (HT) are used to store routes that were already evaluated in stage 3 of the algorithm by the LNS. Toffolo et al. [56], Zachariadis et al. [27], Leung et al. [19] and Wei et al. [25] used similar strategies. One structure registers the feasible routes ( f rHT) already evaluated and the other registers the infeasible routes (irHT). During the selection of the best neighbor in the tabu search, the HT structures are consulted to previously check the feasibility or unfeasibility of the solutions found.

Neighborhood Structures
Like Wei et al. [25], we used four types of neighborhood structures in our study as shown in Figure 10. When selecting one of the four types of structures, changes in the current solution can be applied to a single route or a pair of routes. The selection of the route or the pair of routes to apply the movements from a structure is random. Section 4.3.2 provides more details on how structures are used to generate neighbors for the current solution.
The first type of structure is customer relocation. In this type of movement, a customer is relocated to another position within the same route or is relocated to another route ( Figure 10a). These movements are known as intra-shift and inter-shift, respectively, in Wei et al. [26]. The second type is a customer exchange (Figure 10b), where there is an exchange of positions between two customers on the same route or the exchange of two customers between two different routes within the solution intra-swap and inter-swap according to Wei et al. [26]). The third structure is a route interchange (Figure 10c). If applied to only one route, two positions are randomly selected, and the customers between these two positions are repositioned in reverse order (intra-2opt in Wei et al. [26]). If applied to a pair of routes, where a customer on each route is randomly chosen, the initial portion of the first route to the customer's position is connected to the final portion of the second route from the position of the other customer, and vice versa. The last structure is called block exchange (Figure 10d), where there is an exchange of 2 route segments that can randomly be 1 to 3 in size.

Generating Neighbors as Candidate Solutions for the Routing Problem
Each iteration of the tabu search calls for a procedure that generates new candidate solutions from the current solution. One of the four structures mentioned in Section 4.3.1 is selected randomly with an equal probability to be applied to the current solution and generate a candidate solution. The feasibility of each candidate solution is tested in terms of the total weight and total area of the load, observing the strategy of the area of the vehicle considered for loading. This is controlled by the percAreaVehicle variable. Also, to avoid wasting time, the routes of the solutions already analyzed regarding the feasibility of loading by the LNS are consulted in the hash tables irHT and frHT. If an infeasible route is generated and it is in the irHT, then it is discarded. On the other hand, if a feasible route is generated and it is in the frHT, it is then kept; even if the area of its load exceeds the maximum considered area for loading, this is controlled by the percAreaVehicle. The latter conditional is applied when the parameter percAreaVehicle < 100%. Thus, the best candidate is selected to be the next current solution and the procedure is carried out successively to generate new neighbors in search of the best solution.

The Loading Procedure
For the loading subproblem, we developed the main procedure given by Algorithm 3 (EvalSolutionLoading). This procedure is responsible for carrying out assessments regarding the feasibility of loading the best solution found by the tabu search in the routing stage, as well as elite solutions, if necessary.
To verify the feasibility of loading each route of the evaluated solutions, the Eval-SolutionLoading algorithm calls a procedure based on the LNS, denoted as Algorithm 4 (LNSpack). The code of LNSpack was partially based on the work published by Erdogan [57]. The two procedures are described in the next subsections.  5: for each route in the solution do 6: if route is in f rHT then 7: Sign the route as feasible; 8: else 9: Start new thread with procedure to evaluate route loading in parallel: 10: Set items[] = all items demanded by customers in route; 11: isFeasible = LSN pack(items[]); 12: if isFeasible then 13: Sign route as feasible and insert it into f rHT; 14: else 15: Sign route as infeasible and insert it into irHT; 16: end if 17: End thread; 18: end if 19: end for 20: Wait until all threads end; 21: if all routes are feasible then 22: if evaluated feasible solution is better than globalBestSolution or the globalBestSolution is null then 23 The EvalSolutionLoading algorithm receives two parameters: the best solution (be-stRoutingSolution) and the elite solutions (eliteList[]) found in the routing phase. Then, these solutions are put on a list to be evaluated (evalList[]). The list is ordered from lowest to highest cost so that the best solutions are evaluated first. As soon as a solution with feasible loading is found, the evaluation of the other solutions does not proceed. In the evaluation of each solution registered in evalList[], the verification of the feasibility of loading for each route is performed by the LNSpack algorithm (Algorithm 4), which is discussed in Section 4.4.2. To save time and processing, the first step in evaluating the route is to check if it already exists in the frHT, that is, to check if it was already evaluated and identified as feasible. In that case, there is no need to evaluate it again. The routes that are already identified as infeasible are no longer generated during the routing process.
For the evaluation of the routes, a multithreading mechanism was implemented to enhance the performance of the algorithm. Thus, several routes are evaluated simultaneously within a time limit that was defined in the algorithm's initial parameters. Therefore, the power of parallel processing is exploited using multiple CPU cores. Within each thread, the LNSpack heuristic is called to evaluate the route loading. The LNSpack determines whether the route is feasible or not. Then, the hash tables frHT and irHT are updated, respectively. After all the route evaluation threads are finished, the solution is flagged as feasible if all routes are feasible. If this is the case, and if the solution is better than the best global solution found or if a better global solution still does not exist, the best global solution is updated, ending the procedure and returning to the global algorithm flow. Use the first-fit-decreasing heuristic to repack the removed items; 4: while elapsedTime < maxLNSTime do 5: Perturb solution: randomly remove items from vehicle; 6: Use the first-fit-decreasing heuristic to repack the removed items; 7: if all items were packed then 8: return success; 9: if new solution is better than the best-known solution then 10: Update best-known solution; 11: else 12: Revert the best-known solution; 13: end if 14: end if 15: Set elapsedTime; 16: end while 17 In this work, a version of an LNS algorithm, which is denoted as Algorithm 4 (LNSpack), was developed to solve the 2BPP. The LNSpack requires a list of rectangular items demanded by all customers in the route as parameters. When the LNSpack is called by the procedure for evaluating the feasibility of a route, firstly, the items are sorted, and then, the First-Fit Decreasing heuristic [58] is used to try to load the vehicle. Two forms of ordering are applied: the first is by the area, and the second is by the circumference of the item. Since in this work, only the 2|UR|L was considered, the ordering of the items by area or by circumference is quite effective. Therefore, orders by the height or width of the items were not used. Then, an iteration loop is initiated where the solution is disturbed by the random removal of some loaded items. Next, the First-Fit Decreasing heuristic is performed again to try to reload them. If all items were loaded, the procedure was completed with success. Otherwise, it checks if the loaded area of the current solution is larger than the area of the best solution found and updates the best solution. If the current solution is not better, the best solution is once again the current solution for the next iteration. The condition for stopping the iterations is finding the feasible packaging route items or reaching the maximum time (maxLNSTime). If it is the case, the procedure returns as a packaging failure.

Results
The algorithm proposed in this study, called the CGA-TS-LNS, was coded in C, and the tests were performed on a virtual machine configured as an Intel 2 × Deca Core Xeon E5-2640 CPU with a clock speed of 2.40 giga-hertz and 16 gigabyte of RAM running on a Windows Server 2019 Standard edition operating system. For the benchmarks, the tests were carried out with the well-known instances of Gendreau et al. [15]. The database is composed of 180 instances that are divided into five classes according to the characteristics of the items. Class 1 is characterized as the pure CVRP, where each customer is associated with only one item with a unit of width and a unit of length, and there are no loading restrictions. In Classes 2-5, the quantity of items demanded for each customer i is generated by an uniform distribution within a given interval, according to Column 2 of Table 1. There are 3 categories in which the items can be classified according to their forms: vertical, homogeneous, and horizontal. The dimensions of the items (w × l) are uniformly distributed according to the ranges established for the categories of the items (Columns 3-8).

Initial Parameters
After carrying out experiments with the algorithm in 15 instances of the pure CVRP, with low-, medium-, and high-complexities, we observed the results from the use of values from 5-50 for the size of the tabuList[]. We chose to set sizeTL = 8, because with this value we obtained the best results in the experiments. As it is not a very high value, the performance of the TabuSearch algorithm improved. This is because in each iteration the new current solution must not be in the tabuList[], and for that, the tabuList[] must be checked lineby-line. So, with fewer positions in the list, less computational time is required. The numNeighbors variable, which dynamically defines the number of neighbors during the tabu search, is defined initially from a polynomial function of degree 3. This function was found from a data analysis as described in the following.
Firstly, tests of the algorithm were performed with instances from the pure CVRP (Class 1). Several values between 5 and 200 were randomly assigned to the textitnumNeighbors for each execution of the framework on the instances chosen for the test. The values of the variable where the best results were obtained for each instance were noted. Then, using Pearson's correlation coefficient, a correlation analysis was performed on the data of the instances in relation to the values of numNeighbors noted. The data showed that the number of customers presented a correlation coefficient of more than 80% in relation to the values of numNeighbors noted. Therefore, the number of customers together with the values of numNeighbors noted in the tests was chosen to obtain the polynomial regression function. The function y = 1.67 × 10 -5 × x 3 − 8.90 × 10 -3 × x 2 + 1.57 × 10 0 × x − 1.83 × 10 1 , where x ∈ Z | x ≥ 15 and y ∈ Z, y = numNeighbors and x is the number of customers, was obtained through polynomial regression of degree 3, as shown in the graph of Figure 11. Also, the function found for values of numNeighbors must be used for problems with more than 15 customers, otherwise negative numbers would be obtained. For problems with fewer than 15 customers, numNeighbors is limited to a minimum of 5.
The sizeEL that defines the size of the elite solution list was set to 20. The percAreaVehicle variable was initially set to 100%, as described in Section 4.3.1, and it was applied to resolve instances of Classes 2-5. This strategy proved to be more effective for routes in which the items to be loaded had very irregular dimensions. Although it can be considered a simple strategy that can go against the idea of maximizing the loading space of the vehicle used, in general it showed good results. Due to the tests carried out in our experiments, this strategy is not very efficient in the largest and most complex instances. This is because a lot of time is needed to achieve promising results. The parameters for stopping conditions of the stages of the algorithm, as well as the global stopping condition parameter, were defined as follows: the parameter maxNoImproveCnt2End, which is one of the stopping conditions of the TabuSearch in the routing phase, was set to 400,000 iterations. Namely, the tabu search algorithm implemented in this framework can perform thousands of iterations per second. The parameter maxNoImproveCnt was set to 40,000. In our tests, generally few iterations of the CGA were necessary to generate promising initial solutions, so the maxGen parameter was set to 20. The maxLNSTime parameter was set to 60 s. The maxGTime was set to 7200 s. The vehicle fleet was considered homogeneous and the dimensions of the loading area were H = 40 and W = 20. Figure 11. Polynomial regression to get function to set numNeighbors in TS algorithm.

Results Comparison
This work solved the pure CVRP (Class 1) and the 2|UR|L version of the 2L-CVRP (Classes 2-5), where the rotation of the items is allowed, and the loading is unrestricted. We compared our results for Class 1 with 4 of the best previous approaches: PRMP [22], VNS [24], SA [25], and BR-LNS [59]. Importantly, after an extensive research, we found only these four published studies that involve 2|UR|L-CVRP, which is the focus of our work, and none of them implemented parallelism in their approaches. To solve the 2|UR|L, we found only four approaches with which to compare the results: ACO by Fuellerer et al. [16], MS-BR by Dominguez et al. [23], xPRMP by Zachariadis et al. [27], and SA by Wei et al. [25]. The framework was run 10 times for each instance with 1-10 random seeds, as was done in several previous works, and the best solutions were compared to that of the previous approaches.
The running time to find best cost values for all instances can be considered compatible with the previous works. The running time for each one of the 180 instances remains less than the maximum global time limit to run the algorithm, i.e., less than 7200 s. For the pure CVRP, our method proved to be more effective and the running time to find best solutions in this class was less than the time informed by the other previous approaches. For the 2L-CVRP, our algorithm takes a little more time running, but it was able to find some unpublished results. Table A1 shows the results for the Class 1 instances, where the best costs are compared with the other four previous frameworks. The CGA-TS-LNS found better and higher quality solutions for the 6 most complex instances (30)(31)(32)(33)(34)(35)(36), and matched the BKS for 28 instances.

Results for Class 1 Instances
In addition, Figure 12 shows the average cost of the solutions found by the CGA-TS-LNS for the 36 instances of Class 1 is lower than all previous approaches, and lower than the average cost of all the BKS. This demonstrates the effectiveness of our algorithm to solve the CVRP.

Results for 2|UR|L Instances
When comparing our results for the 2|UR|L in Classes 2-5, we find that our average cost in each class is better than two of the only four algorithms proposed to solve this version found in the literature. In total, the CGA-TS-LNS surpassed the BKS in 7 of the 144 instances and matched in another 66. The CGA-TS-LNS only loses to the ACO and MS-BR in 5 instances of 144, with 4 instances of Class 2, and 1 instance of Class 4. For the 2L-CVRP, the SA approach is evidently still the best, achieving the lowest average cost. Table A2 compares the best 2|UR|L results for the Class 2 instances. In comparison with that of the other four approaches and the BKS of the 36 instances, the CGA-TS-LNS found better solutions for 2 instances and matched the other 15 better solutions. Considering the average cost, xPRMP is the best framework to solve the Class 2 instances. Table A3 compares the best results from the 2|UR|L over the 36 instances of Class 3, where our CGA-TS-LNS found better higher quality solutions for 3 instances and matched 18 other better solutions. Considering the average cost for 36 instances, SA is the best approach for this Class. Table A4 compares the best results of the 2|UR|L for the 36 instances of Class 4. The CGA-TS-LNS found better, higher quality solutions for 1 instance and matched the other 14 best solutions. Again, the SA algorithm is the best to solve Class 4 instances. Table A5 compares the best results of the 2|UR|L for the 36 instances of Class 5. Our algorithm found better, higher quality solutions for 1 instance and matched 19 other better solutions. One more time, the SA algorithm is the best approach to solve the Class 5 instances.

Statistical Validation for Algorithm Comparison
According to [60] inside the field of inferential statistics, hypothesis testing can be employed to draw inferences about one or more populations of given samples (results). To perform that, two hypotheses, the null hypothesis H 0 and the alternative hypothesis H 1 , are defined. The null hypothesis is a statement of no effect or no difference, whereas the alternative hypothesis represents the presence of an effect or a difference (in our case, significant differences between algorithms). When applying a statistical procedure to reject a hypothesis, a level of significance α is used to determine at which level the hypothesis may be rejected. The Sign test for multiple comparisons described in [60] allows us to highlight those algorithms whose performances are statistically different when compared to that of the CGA-TS-LNS algorithm.
As defined in [60], when using a level of significance α = 0.10 and setting our hypotheses to be H 0 :M j ≥ M 1 and H 1 :M j < M 1 and m = 4 (m = k − 1) and n = 36, that is, our CGA-TS-LNS algorithm performs significantly better than the remaining algorithms. Also, according to Table A.21 of [60], for m = 4 (m = k − 1) and n = 36 reveals that the critical value of R j is 10.
Tables A6-A10 present the number of times the CGA-TS-LNS algorithm was superior to its competitors, as well as the number of occasions in which the CGA-TS-LNS algorithm obtained the optimal solution. In Tables A1-A5, the best-known solution from the literature (BKS) is the optimal solution to the problem and, therefore, it is impossible to obtain better solutions. Thus, we add the number of times the GGA-TS-LNS algorithm was better to the competitors to the number of times the BKS was obtained. Table A6 shows the performance comparison of GGA-TS-LNS with PRMP, VNS, BR-LNS, and SA algorithms for the instances of Class 1 (pure CVRP). Since the number of minuses in the pairwise comparison between the GGA-TS-LNS algorithm and all others is equal to 36, we can conclude that GGA-TS-LNS has a significantly better performance than all of other competitors. Table A7 shows the performance comparison of GGA-TS-LNS with ACO, MS-BR, SA, and xPRMP algorithms for the 2|UR|L instances of Class 2. Since the number of minuses in the pairwise comparison between the GGA-TS-LNS algorithm and ACO and MS-BR algorithms is equal to 34 and 25, respectively, we can conclude that GGA-TS-LNS performs significantly better than these two competitors. Table A8 shows the performance comparison of GGA-TS-LNS with ACO, MS-BR, SA, and xPRMP algorithms for the 2|UR|L instances of Class 3. Since the number of minuses in the pairwise comparison between the GGA-TS-LNS algorithm and ACO and MS-BR algorithms is equal to 35 and 27, respectively, we can conclude that GGA-TS-LNS performs significantly better than these two competitors. Table A9 shows the performance comparison of GGA-TS-LNS with ACO, MS-BR, SA, and xPRMP algorithms for the 2|UR|L instances of Class 4. Since the number of minuses in the pairwise comparison between the GGA-TS-LNS algorithm and ACO and MS-BR algorithms is equal to 36 and 30, respectively, we can conclude that GGA-TS-LNS performs significantly better than them. Also, in the pairwise comparison between the GGA-TS-LNS algorithm and XPRMP algorithm, with a level of significance α = 0.05, we can conclude that GGA-TS-LNS performs significantly better than it as well. Table A10 shows the performance comparison of GGA-TS-LNS with ACO, MS-BR, SA, and xPRMP algorithms for the 2|UR|L instances of Class 5. Since the number of minuses in the pairwise comparison between the GGA-TS-LNS algorithm and ACO and MS-BR algorithms is equal to 36 and 34, respectively, we can conclude that GGA-TS-LNS performs significantly better than these two competitors.

Conclusions
This paper presents a new hybrid algorithm to solve the 2L-CVRP. The combination of algorithms based on the CGA, the TS, and the LNS proved to be competitive when applied to the 2L-CVRP. Our results outperformed two of the only four algorithms proposed, so far, to solve version 2|UR|L of this problem. As stated by Wei et al. [25], it is not efficient to simply combine conventional algorithms for the CVRP and the 2BPP, so some loading strategies need to be employed to facilitate the solution of these integrated problems. In our case, a simple strategy was to consider the loading area dynamically during the iterations of the algorithm. This was combined with the diversification and intensification strategies, plus the use of multithreading, to evaluate the loading of the routes. This facilitated the process to find appropriate solutions for the 2L-CVRP, but it was not enough to outperform the algorithms proposed by Wei et al. [25] and by Zachariadis et al. [27].
When analyzing the results for the pure CVRP, the framework proved to be extremely efficient, surpassing all other approaches that were applied to the instances proposed by Gendreau et al. [15]. This study brings, as a reference for future works, new best-known solutions for the most complex instances of the pure CVRP and a better average of solutions for the 36 instances of this Class.
The algorithm proposed in this work has great potential, but it can still be improved. For future works, we can think about the use of other strategies to improve solutions regarding the loading subproblem. The ability to solve types 2|OU|L, 2|SO|L, and 2|SR|L can be added to the algorithm. The algorithm can also be adapted to solve other problems, such as the 3L-CVRP. Additionally, it can be tested in other databases. Also, to provide a more practical information measure, it is possible to use average values to compare the results with that of other research.
The main application would be to solve the CVRP and to compare the results and the effectiveness. In addition, we believe that the algorithm can serve as a basis for other studies to obtain improvements in logistical and transportation processes. Its application can contribute to reducing costs in the routing problems of capacitated vehicles in the real world. In this sense, the integration of this algorithm with intelligent models for logistics management [61] will improve the routes processing, and consequently, the safety of transportation. The algorithm also will allow the implementation of ubiquitous intelligent services for vehicular users [62]. Finally, future work will integrate the algorithm with strategies used to treat Context Histories [63][64][65][66] such as pattern analysis [67], context prediction [68], similarity analysis [69], cryptography [70], and IoT challenges in smart environments [71,72]. This integration will support better solutions; specifically, to routing problems and in general for problems involving context histories. Acknowledgments: This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brasil (CAPES)-Finance Code 001. The work described in this paper was supported by Master's in Systems and Industrial Processes of the University of Santa Cruz do Sul. This work was supported by national funds through the Fundação para a Ciência e a Tecnologia, I.P. (Portuguese Foundation for Science and Technology) by the project UIDB/05064/2020 (VALORIZA-Research Centre for Endogenous Resource Valorization).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: