Mathematical Model and Simulated Annealing Algorithm for the Two-Dimensional Loading Heterogeneous Fixed Fleet Vehicle Routing Problem

+is paper addresses the two-dimensional loading heterogeneous fixed fleet vehicle routing problem, which is a complex and unstudied variant of the classical vehicle routing problem and has a wide range of applications in transportation and logistics fields. In this problem, each customer demands a set of rectangular two-dimensional items, and the objective is to find the minimum cost delivery routes for a limited set of vehicles with different capacities, fixed and variable operating costs, and rectangular two-dimensional loading surfaces. We formulate a mixed integer linear programming model to obtain optimal solutions for small-scale problems. To obtain solutions for large-scale problems, we develop an algorithm based on simulated annealing and local search, which uses a collection of packing heuristics to address the loading constraints, and we also propose three new heuristics. We conduct experiments on benchmark instances derived from the two-dimensional loading heterogeneous fleet vehicle routing problem. +e results indicate that the proposed model correctly describes the problem and can solve smallscale problems, that the new packing heuristics are effective in improving the collection of packing heuristics, and that the proposed simulated annealing algorithm can find good solutions to large-scale problems within an acceptable computational time. Hence, it can be used by logistic companies using a heterogeneous fixed fleet in the integrated planning of vehicle loading and routing.


Introduction
In many organizations, the management of distribution activities is a major decision-making problem. Every manufacturing system needs an efficient way to supply its products to retailers and customers. Most firms require delivery vehicles to service a network of demand locations, meaning the efficient use of a vehicle fleet is the main feature of almost all distribution problems, since transportation costs and lead time have an important impact on supply chain management [1,2]. e vehicle routing problem (VRP) is considered one of the most important factors in both logistics and freight transportation systems. It consists of proposing the most costeffective way to deliver items from depots to customers using a fleet of vehicles. e most common costs associated with the problem are related to driving distance.
ere are several variants of the VRP due to additional constraints encountered in real-world applications [3]. Regarding fleet composition, the classical problem, called the capacitated vehicle routing problem (CVRP), considers vehicles with the same limited capacity, whereas the heterogeneous fleet vehicle routing problem (HFVRP) considers vehicles with different capacities and costs [4]. Different features of these heterogeneous VRPs have been studied. For instance, the fleet size and mix vehicle routing problem (FSMVRP) deals with an unlimited number of vehicles, and in this case, in addition to routing, the problem includes the management of the fleet composition, while the heterogeneous fixed fleet vehicle routing problem (HFFVRP) considers a limited number of vehicles.
In the above VRP variants, customer demands are usually expressed simply as weights or volumes. In this case, checking the feasibility of a solution simply requires one to ensure that the sum of the customer demands assigned to each vehicle does not exceed its total loading capacity. However, a variety of real-life applications in the distribution management context involve the transportation of rectangular-shaped items, that is, items that cannot be stacked due to their weight or fragility, such as household appliances or delicate pieces of furniture. In all these cases, it is necessary to consider additional constraints to reflect the two-dimensional loading feature of the problem, because the way these items are packed into vehicles might have a significant influence over distribution costs. When dealing with rigid discrete items, their geometry may lead to infeasible solutions if the vehicle does not have sufficient capacity. In other words, it may not be possible to accommodate items in vehicles because of their geometrical characteristics, so logistic managers must simultaneously deal with routing and packing [5][6][7]. Since routing and packing are well known NP-hard problems, combining them leads to an extremely challenging optimization problem. e two-dimensional loading capacitated vehicle routing problem (2L-CVRP) is one of the first approaches to integrate vehicle routing and loading problems. In this problem, customers demand a set of rectangular two-dimensional items, and identical vehicles have a two-dimensional loading surface. In real-world problems, in contrast to the 2L-CVRP, enterprises own a diverse fleet of vehicles, which offers the flexibility to design a more profitable distribution plan. e two-dimensional loading heterogeneous fleet vehicle routing problem (2L-HFVRP) treats the 2L-CVRP with an unlimited heterogeneous fleet. Nevertheless, since most companies that have to deliver goods own a limited fleet of vehicles, it is crucial to study routing problems that involve heterogeneous fixed fleets and loading constraints. To the best of our knowledge, we are not aware of studies that have been conducted to address such a VRP, although it is a practical problem in real-world transportation and logistics enterprises.
In this paper, we thus combine the HFFVRP with twodimensional loading constraints, which is called the twodimensional loading heterogeneous fixed fleet vehicle routing problem (2L-HFFVRP). In the 2L-HFFVRP, there are limited numbers of vehicles of different types. Each vehicle has a carrying capacity, a fixed cost related to the use of the vehicle, a variable cost proportional to the distance traveled, and a rectangular loading surface with a given length and width. Customer demand is defined by a set of rectangular items of a given width, length, and weight. All items belonging to a particular customer must be assigned to the same route. e objective is to describe the most costeffective item delivery routes that start from and return to a depot. In practice, the need for different types of vehicles is determined by customer characteristics. Usually, larger vehicles are more appropriate for serving customers who require large orders, while smaller vehicles are more adequate for delivering small quantities or serving customers that have access restrictions.
Due to their structure, it happens frequently that the items may not be picked up from any side by the loading/ unloading equipment. Additionally, vehicles are generally rear-loaded, and load rearrangement at the customer site can be difficult, time-consuming, or even impossible due to the weight and size of the items or the limitations of forklift trucks; therefore, each item to be unloaded must not be blocked by other items yet to be unloaded, even partially, in the rectangular area from its loading position to the unloading side of the truck. From the viewpoint of the loading problem, there are different constraints that could be considered. Depending on these constraints, it is possible to distinguish between the following loading settings: (a) oriented loading (OL), where the rotation of items is not allowed, that is, it is assumed that all items have a fixed orientation; (b) nonoriented loading (RL), where it is allowed to rotate items by 90°during the packing process; (c) sequential loading (SL), where items are always loaded in a reverse order to the order in which customers are visited but rearrangements of the items inside the vehicle are not allowed once the route has started; and (d) unrestricted loading (UL), where items are allowed to be rearranged during the distribution process. Figure 1 illustrates the difference between unrestricted and sequential loading. e purpose of this paper is to present two approaches for solving 2L-HFFVRP. First, a mixed integer linear programming model is formulated to obtain optimal solutions for small-scale problems. To this end, we developed a fourindex formulation that simultaneously manages the routing and loading aspects and also considers the sequential loading constraints. Second, to obtain solutions for large-scale problems, we propose an algorithm based on simulated annealing (SA) and local search. Specifically, we developed a heuristic algorithm to obtain an initial feasible solution and used SA and local search to explore the routing problem search space, while the loading constraints were solved by a bundle of packing heuristics composed of five heuristics from the literature and three new heuristics. We assume that all items have a fixed orientation and consider both unrestricted and sequential loading. e remainder of this paper is organized as follows. Section 2 briefly reviews the literature on combined twodimensional loading and routing problems. In Section 3, we define the problem to be addressed. In Section 4, we present the mathematical formulation for this problem. In Section 5, the packing heuristics are described, and a hybrid SA algorithm is developed. Section 6 presents the set of instances generated and explains the computational experiments and results. Finally, Section 7 provides the managerial insights of the study and Section 8 draws concluding remarks and perspectives for future research.

Literature Review
Integrated vehicle routing and loading problems have been increasingly studied, owing to their relevance to logistics distribution systems [5][6][7]. To the best of our knowledge, there are no studies on the 2L-HFFVRP to date. is section will focus on the previous studies on the 2L-CVRP and the 2L-HFVRP. Both these problems combine vehicle routing and two-dimensional loading constraints but consider a homogeneous and an unlimited heterogeneous fleet, respectively. As a generalization of VRP, these problems and their variants are known to be NP-hard problems, so exact methods are suitable for solving only small-scale problems, while metaheuristic algorithms are more favored for solving complex problems. e 2L-CVRP was first presented by Iori et al. [8]. ey addressed the sequential oriented version of the problem using an exact methodology that employs a branch-and-cut algorithm to deal with the routing aspects of the problem, combined with a nested branch-and-bound procedure to guarantee feasible loadings of the items into the vehicles. Subsequently, Gendreau et al. [9] proposed a metaheuristic based on tabu search for the routing problem, and checking for a feasible loading is performed via heuristics, lower bounds, and a truncated branch-and-bound procedure, considering both sequential and unrestricted oriented loading. Zachariadis et al. [10] developed a hybrid heuristic algorithm based on tabu search and guided local search and introduced a collection of packing heuristics to check the feasibility of the loading. Leung et al. [11] developed an improved algorithm based on guided tabu search and SA  and introduced a new packing heuristic to address loading constraints. Fuellerer et al. [12] developed an ant colony algorithm to solve the routing problem, combined with heuristics for the loading subproblem. e authors addressed sequential and unrestricted problems, as well as oriented and nonoriented loading, and the rotation of items was allowed for the first time. Duhamel et al. [13] proposed a method that combines greedy adaptive random search with evolutionary local search and transformed the loading constraints into a project scheduling problem. A heuristic algorithm called promise routing-memory packing, based on the compression idea, was introduced by Zachariadis et al. [14]. Wei et al. [15] proposed an SA algorithm, with a mechanism of repeatedly decreasing and increasing the temperature that uses an open space-based heuristic to deal with loading constraints, which outperformed all previous approaches on all problem versions. Recently, Ji et al. [16] proposed an enhanced neighborhood search algorithm incorporated with a based tabu search packing algorithm to solve 2L-CVRP with split delivery.
is variant was also addressed by Ferreira et al. [17] through an exact branchand-cut approach, where a tailored procedure that includes the computation of lower bounds, a constructive-based heuristic, and a constraint programming model are proposed to handle the packing problem.
Concerning 2L-HFVRP, Leung et al. [18] first addressed the problem through an SA algorithm with a heuristic local search to solve the routing problem and a collection of six packing heuristics to solve the loading constraints, considering oriented sequential and unrestricted versions of the problem. Dominguez et al. [19] developed an algorithm that combines biased-randomized versions of routing and packing heuristics to solve oriented and nonoriented unrestricted problems. e biased randomization of heuristics refers to the use of skewed probability distributions to induce nonsymmetric random behavior in a heuristic procedure and transform a deterministic heuristic into a multistart probabilistic algorithm. A hybrid swarm algorithm based on artificial bee colony and artificial immune system was proposed by Zhang et al. [20] to deal with sequential and unrestricted problems. More recently, Sabar et al. [21] proposed an adaptive memetic algorithm that integrates new multiparent crossover operators with multilocal search algorithms in an adaptive manner for solving the routing problem, while hybridization of five packing heuristics is used to perform the loading process, considering only nonoriented unrestricted loading. A comparison between the previous studies, involving heterogeneous vehicle routing problems and two-dimensional loading constraints, and this study is shown in Table 1.
Our literature review of the vehicle routing problems that considers a heterogeneous fleet and two-dimensional loading constraints reveals that other researchers only solved these problems considering an unlimited fleet. However, in practice, it is not a realistic scenario, since companies generally own a limited fleet. erefore, we address the 2L-HFFVRP.

Problem Statement
e 2L-HFFVRP is defined on a complete undirected graph G � (N, A), where N � 0, 1, . . . , n { } is a set of n + 1 vertices, the vertex 0 denotes the depot, and the vertices 1, 2, . . . , n correspond to the positions of the customers 1, 2, . . . , n. A � (i, j); i, j ∈ N is a set of edges, and each edge (i, j) ∈ A is associated with a distance d ij (d ii � 0) that denotes the distance from customer i to customer j. A set V � 1, . . . , v { } of v vehicles is available at the depot. Each vehicle k ∈ V has a weight capacity Q k , a loading surface with length L k and width W k , a fixed cost F k , and a variable cost V k . In general, a vehicle with a larger capacity usually has higher fixed and variable costs. e traveling cost of each edge Each customer i (i � 1, . . . , n) demands a set of m i rectangular items, denoted by IT i , and the total weight of IT i is equal to D i . Each item I ir ∈ IT i (r � 1, 2, . . . , m i ) has a specific length l ir and width w ir .
For 2L-HFFVRP, a feasible solution must satisfy the following constraints: (a) Each vehicle must start and finish at the depot (b) Each customer can be served only once (c) e weight capacity, length, and width of the vehicle cannot be exceeded (d) All items of each customer must be loaded on the same vehicle (e) Each item has a fixed loading orientation and must be loaded with its sides parallel to the sides of the loading surface (f ) Overlap between items within the same vehicle is not allowed e objective of the 2L-HFFVRP is to find a set of routes of minimum cost that fulfill the constraints and satisfy all customers' demands. e constraints mentioned above refer to the unrestricted 2L-HFFVRP. Besides this problem version, this study also considers the sequential 2L-HFFVRP.

Mixed Integer Linear
Programming Formulation is section proposes a mathematical formulation to represent the 2L-HFFVRP, considering the characteristics described in Section 3. e notation used in the model, including indexes, sets, parameters, and variables, is

Authors
Limited fleet Items rotation Sequential loading Leung et al. [18] No No Yes Dominguez et al. [19] No Yes No Zhang et al. [20] No No Yes Sabar et al. [21] No

4.1.
Notations. e mixed-integer linear programming model depends on the sets, parameters, and variables described in Table 2.

Mathematical Formulation for the HFFVRP.
A number of formulations could be used to model the HFFVRP [22]. In this study, we present a formulation based on the timedependent formulation of the traveling salesman problem for integrated routing and packing problems [23]. Although this is a heavy four-index formulation, it gives us information about the position of customers in each route, which is necessary to ensure that the sequential loading constraints are fulfilled.
It is defined that a vehicle arrives at each vertex in stage t. Moreover, it is assumed that each route starts at vertex 0 in stage t � 0 and finishes at vertex 0 in stage t � n + 1 if only one vehicle is used. erefore, the set of possible values for t is T � 1, . . . , n + 1 { }. e (routing) decision variables z kt ij indicate whether vehicle k goes directly from vertex i to vertex j in stage t (z kt ij � 1) or not (z kt ij � 0). It is assumed that d ii � 0, (i, i) ∈ A. e mathematical model for the HFFVRP is expressed as follows: subject to j∈N k∈V t∈T In this formulation, the objective function (1) aims to minimize the total cost for the vehicles to visit all customers.
Constraints (2) ensure that each customer is visited exactly once. Constraints (3) ensure the connectivity of each tour; that is, they ensure that if customer i is visited in stage t, then this customer has to be the starting point for some other customer in stage t + 1. Constraints (4) ensure that each vehicle leaves the depot at most once in stage 1. Constraints (5) ensure that if vehicle k travels from customer l to customer i in stage t, then in stage t + 1, the same vehicle k must travel from customer i to another customer j. Constraints (6) ensure that the weight capacity of the vehicles is not exceeded. Constraints (7) are the (routing) decision variable domain constraints. To ensure the validity of this formulation, we assume that z kt that is, the vehicles can leave the depot only in the first stage.

Two-Dimensional Loading Constraints.
In this section, we show how two-dimensional loading considerations can be embedded into formulation (1)-(7) of the last subsection.
ese constraints address the geometrical loading aspects; that is, they ensure that items are packed completely inside the vehicles and that they do not overlap each other in each vehicle.
Two-dimensional loading constraints for routing problems have been addressed by Junqueira [23], who developed models for integrated routing and loading problems; however, the formulation does not favor the optimum use of the vehicle loading surface when sequential loading is considered since it contains a parameter that limits the number of length units allowed to exceed the frontal border of the items of customers already loaded in order to arrange the items of a customer that will be visited earlier in the route. Due to this limitation, we build on a formulation for the container loading problem [24] adapted to the twodimensional case and to the vehicle routing context.
Let us create a set I of all items to be loaded in the vehicles as follows. e first m 1 elements of I are the items of customer 1, the next m 2 elements are the items of customer 2, and so on. erefore, each item r ∈ IT i becomes an item p � i− 1 j�1 m j + r, so each customer i demands a set of items Each item p ∈ I has a specific length l p and width w p . e total number of items in set I is B � n i�1 m i . We assume that the dimensions of items and vehicles are integers, that items can be placed only orthogonally in a vehicle, and that their orientation is fixed. A Cartesian coordinate system with its origin in the vehicle's front-left corner is adopted, and the vehicle door through which loading and unloading operations are performed is placed on the side between coordinates (L, 0) and (L, W), as shown in Figure 2. e (loading) decision variables indicate coordinates (x p , y p ) of the front-left corner of item p in the vehicle loading surface, and the variables s pk indicate whether item p is loaded into vehicle k (s pk � 1) or not (s pk � 0). Binary variables α pq , β pq , c pq , and δ pq (p, q ∈ I, p < q) indicate the placement of items relative to each other. α pq , β pq , c pq , and δ pq are equal to 1 if item p is, respectively, loaded behind, in Mathematical Problems in Engineering Binary variable that indicates whether item p is loaded on the right side of item q (δ pq � 1)′ or not (δ pq � 0)′ front of, on the left side of, or on the right side of item q. An example of the use of these variables is shown in Figure 3, which shows two items p and q loaded in a vehicle.
Assuming that p < q, we have α pq � 1, β pq � 0, c pq � 0, and δ pq � 0. Note that an item p is said to be loaded behind or in front of an item q only if there is no intersection of its projections on the x-axis. Analogously, item p can be on the left or on the right of item q if there is no intersection of its projections on the y-axis. erefore, it is possible to have α pq + β pq ≤ 1 and c pq + δ pq ≤ 1, while 1 ≤ α pq + β pq + c pq +δ pq ≤ 2.
Let M be a sufficiently large number. e formulation for the 2L-HFFVRP is composed of the objective function (1) and constraints (2), (3), (4), (5), (6), and (7) presented before, plus the following constraints: Constraints (8)- (11) ensure that the items do not overlap with each other. is check for overlap is necessary only if a pair of items is placed in the same vehicle, and this is ensured by constraints (12). Constraints (13) and (14) ensure that all the items loaded in a vehicle fit within the physical dimensions of the vehicle. Constraints (15) ensure the coupling of routing and loading structures; that is, all items required by customer i must be loaded on the same vehicle k that visits this customer. Finally, constraints (16)-(18) are the (loading) decision variable domain constraints.

Sequential Loading Constraints.
In sequential loading problems, vehicle loading must take into account the delivery route of the vehicle and the sequence in which items are unloaded. In other words, items must be loaded in the reverse order in which respective customers are visited. is rule prevents unnecessary additional handling when each delivery point of the route is reached and, consequently, additional time spent unloading and reloading items of the remaining customers.
To address the sequential loading constraints and ensure that the items of a customer do not block the items of customers that will be served before them, we develop the following constraints: If customer i is served before customer j, constraints (19) prevent the items of customer j from blocking the unloading of the items of customer i. ese constraints verify whether customer i is served in stage t 1 preceding stage t 2 in which customer j is served. A particular case is illustrated in Figure 4 that shows an item p of customer i and an item q of customer j. If customer i is not visited before customer j, item p of customer i and item q of customer j have unrestricted x coordinates. However, if customer i is served before customer j, there are two possibilities: if item p of customer i is completely on the left or on the right side of item q of customer j, we have c pq + δ pq � 1, and there is no intersection of their projections on the y-axis, so the item does not block the other one. On the other hand, if item p is not completely on the left or right side of item q, an item is behind or in front of the other, and the constraints impose the condition that item p must be in front of item q. A similar explanation is valid for constraints (20) in the case that customer i is served after customer j. e complete formulation for the 2L-HFFVRP with sequential loading constraints is given by the objective function (1) with the constraints (2)- (20). e four-index formulation is essential to properly model the 2L-HFFVRP when sequential loading constraints are present. Both indices k and t are fundamental to the modeling of the problem: index k is necessary to determine in which vehicle a given item is placed, and index t, in turn, makes it possible to indicate the position of a customer in the route and set the sequential loading constraints. e way these constraints are established allows us to obtain the optimal usage of vehicle loading surfaces, even for strongly heterogeneous items, once the position of a customer is analyzed in relation to the position of all other customers, not only those that immediately precede or follow the respective customer. If these constraints considered only the pairs of customers that are served in sequence in the route, there would be no guarantee that loading constraints would be satisfied.
is formulation can be used to optimally solve smallscale 2L-HFFVRP cases. A metaheuristic algorithm for handling large-scale problems is presented in the next section.
Mathematical Problems in Engineering 7

Hybrid SA Algorithm
Most approaches for solving problems that integrate vehicle routing and loading separate the problems into the following two stages: (1) a main algorithm for routing and (2) a tool for the loading per vehicle. In this paper, we present a hybrid algorithm based on SA and local search for the routing problem that uses heuristic methods to determine the geometric loading of items into vehicles. ese heuristics are described, and the proposed algorithm is explained.

Heuristics for Two-Dimensional
Loading. Given a route, it is necessary to determine whether it is feasible, that is, whether customer demand does not exceed vehicle capacity and whether all the items ordered by customers can be loaded onto the vehicle without overlap. In addition, in sequential loading, the item loading order must respect the order of customers served. To verify the feasibility of a route, we first examine whether the vehicle capacity is exceeded and whether all the items exceed the loading area and surface dimensions. If all conditions are satisfied, then heuristic methods are subsequently applied to determine whether feasible loading exists and, if it does, to determine its configuration. Eight heuristics Heur i (i � 1, . . . , 8) are used to solve the two-dimensional loading problem; one item at a time is inserted in the most appropriate position in a list of available loading positions according to certain criteria to be described later in this section. Initially, only the left front corner of the vehicle loading surface, corresponding to the coordinate point (0,0), is available for loading an item. When an item is inserted, the position it occupies is excluded from the list of available loading positions, while at most two new positions are generated and added to this list. Figure 5 shows the mechanism of item insertion: item D is inserted in the highlighted position in the figure on the left, and the set of available positions after its insertion are shown in the figure on the right.
In sequential loading, when all items of a customer are loaded, positions that cannot be occupied by items that have not yet been loaded are excluded from the list of available loading positions Suppose that the route is 0 − i 1 − i 2 − i 3 − 0. e partial loading for this route and the available loading positions are shown in Figure 6. After the loading of the items of customer i 2 and before the loading of the items of customer i 1 , the highlighted positions are excluded from the list of available loading positions, as they would generate loadings that would contravene the sequential loading constraints.
Heuristics Heur i (i � 1, . . . , 5) are based on the study by Zachariadis et al. [10] and adopt the following criteria to determine an item loading position.
Heur 1 : Bottom-Left Fill (x-axis): the selected position is that with the minimum x-axis coordinate, where ties are broken by the minimum y-axis coordinate. With this heuristic, packing tends to form strips parallel to the x-axis. Heur 2 : Bottom-Left Fill (y-axis): the selected position is that with the minimum y-axis coordinate, where ties are broken by the minimum x-axis coordinate. With this heuristic, packing tends to form strips parallel to the y-axis. Heur 3 : Max Touching Perimeter: the selected position is that with the maximum sum of the common edges between the item to be inserted, the loaded items in the vehicle, and the vehicle loading surface. is heuristic tends to spread items to the edges of the loading surface and later fills the inner part of it.  Mathematical Problems in Engineering Heur 4 : Max Touching Perimeter no Walls: the selected position is that with the maximum sum of the common edges between the item to be inserted and the loaded items in the vehicle. is heuristic tends to fill the inner part of the loading surface before filling the edges. Heur 5 : Min Area: the selected position is that with the minimum rectangular surface, which is the rectangular area available for loading an item into this position.
A detailed explanation of these five heuristics can be found in [10]. As in the study by Zachariadis et al. [10], we first use these five packing heuristics to derive the loading strategy. en, we introduce three new packing heuristics, namely, the Min Occupied Rectangular Area (Heur 6 ), Max Touching Perimeter Select Item (Heur 7 ), and Max Relative Touching Perimeter Select Item (Heur 8 ). Heur 6 selects an item loading position according the following criteria.
Heur 6 : Min Occupied Rectangular Area: the selected position yields the minimum occupied rectangular area after item insertion. Figure 7 shows the occupied rectangular area after inserting item C in the given position. is strategy makes the packing compact.
Heuristics Heur i (i � 1, . . . , 6) presented thus far employ individual criteria to select the position to insert a predefined item; that is, the items are loaded one at a time according to a given sequence. To increase the probability of the heuristics obtaining a feasible solution, three orders are generated for sequential and unrestricted cases. In a given route, each customer has a Mathematical Problems in Engineering unique visit order. If sequential loading is considered, Ord Seq 1 , Ord Seq 2 , and Ord Seq 3 are produced by sorting all items by the reverse customer visit order, breaking ties by decreasing length, width, and surface area, respectively. For the unrestricted problem, Ord Un 1 , Ord Un 2 , and Ord Un 3 are produced simply by sorting all items by decreasing length, width, and surface area, respectively. Note that although Ord Seq 1 , Ord Seq 2 , and Ord Seq 3 are primarily designed to deal with the sequential problem, they succeed in producing feasible unrestricted loadings in numerous cases when Ord Un 1 , Ord Un 2 , and Ord Un 3 fail; therefore, they are also considered for solving unrestricted problems. Heuristics Heur 7 and Heur 8 do not order items for loading, as they analyze the loading of all possible items to be loaded in all available loading positions that are feasible and select an item and a position that best meet the stated criteria. In sequential loading, the analysis is performed for all items of the current loading customer, while in unrestricted loading, any item can be inserted at any time. e criteria employed by Heur 7 and Heur 8 are as follows.
Heur 7 : Max Touching Perimeter Select Item: for each item that can be loaded and for all the respective available loading positions that are feasible, the touching perimeter is calculated, as in Heur 3 . e combination of an item and a position that matches the largest perimeter is chosen, and the item is inserted. Heur 8 : Max Relative Touching Perimeter Select Item: for each item that can be loaded and for all the respective available loading positions that are feasible, the relative touching perimeter is calculated as the ratio between the item perimeter and its touching perimeter. e combination of item and position that matches the largest ratio is selected, and the item is inserted. e feasibility of a route is tested by heuristics Heur 1 to Heur 8 , in this sequence, such that heuristics Heur i (i � 1, . . . , 6) are executed for the orders mentioned earlier according to the sequential or unrestricted case. If any heuristic (considering any suitable order) finds a feasible loading for the route, the route is considered feasible, and the subsequent heuristics are not attempted. Otherwise, if none of the heuristics find a feasible loading for the route, the route is considered unfeasible. For each loading route feasibility check, the proposed packing heuristics are used. As these heuristics are computationally intensive, we use a data structure to avoid duplicate examinations.

Initial Solution.
e initial solution is the starting point of the improvement metaheuristic proposed in the following subsections. e limited number of vehicles available in the depot and their capacities make it relatively difficult to determine an initial feasible solution. If at least one of the following situations occurs, a vehicle cannot service a customer: (1) the total weight of all customer items exceeds vehicle capacity; (2) the total area of all customer items exceeds the loading surface area; or (3) the length/width of any customer item exceeds the length/width of the loading surface. Owing to these limitations, the initial solution is constructed by prioritizing customers that can be served only by larger vehicles.
Given v vehicles available in the depot, it is assumed that First, for each vehicle k, a subset CustomersList k of customers that must be served by vehicle k or a larger vehicle is defined. Subsequently, for each vehicle k � v, . . . , 1, the customers of set CustomersList k are added to a vehicle route k, . . . , v according to the minimum cost insertion so that the route is feasible. e cost of adding customer i to vehicle route k depends on whether route k is empty. If route k is empty, the cost is based on the fixed cost of vehicle k, plus the cost of traveling from the depot to customer i and back. If the route is not empty, the cost of adding customer i to all route positions is analyzed, and the minimum cost insertion is performed. If it is not possible to insert all the customers of some set CustomersList k in a particular route k, . . . , v, the sequence of customer insertion is perturbed by randomly selecting a customer and a vehicle that can serve this customer to initialize the first route, and the procedure is restarted. Algorithm 1 provides a pseudocode for constructing the initial solution. ToInsert k is the set of customers of CustomersList k that have not yet been inserted in any route. At the beginning of route construction, ToInsert k � CustomersList k , and in the course of the algorithm, this set is updated as customers are added to routes.

Neighborhood Structures.
SA explores the search space by performing moves to step from a current solution to a subsequent solution. e neighborhood structures are determined such that each movement causes an interroute disturbance in the solution. Intraroute improvement is accomplished by a local search procedure introduced in the next section.
ree types of movements are employed to disturb a solution, each of which defines a neighborhood structure, NS 1 , NS 2 , or NS 3 , selected randomly in each loop with equal probability. Move type 1-Inter (NS 1 ) performs a customer relocation move [25]. It reassigns a customer from their current route to a position on another route ( Figure 8). is move type can make some routes empty, reducing the number of vehicles used.
Move type 2-inter (NS 2 ) performs route exchange [25]. It swaps the positions of the two customers on different routes (Figure 9).
Move type 3-inter (NS 3 ) is a variant of route interchange 2-opt [26]. Two customers of different routes are selected; in each route, a block is created from the selected customer to the last customer, and these blocks are swapped between routes ( Figure 10). For this move type, a fictitious customer can be selected in one of the two routes to shift a block of customers from one route to another. In this case, this move type can also empty some routes and reduce the number of vehicles used.

Local Search Mechanism.
e neighborhood structures presented in Section 5.3 are employed by the SA algorithm to perform interroute improvements. In each new feasible solution found, at least one route has its set of customers modified; therefore, a local search method is used to optimize each modified route.
Given a single route, the local search performs exhaustive moves in its neighborhood until no improvement in solution can be achieved. ree move types, similar to those for interroute improvement, define neighborhood structures as intraroute optimization. Move type 1-Intra consists of the reallocation of a customer [25]; this move transfers a customer from one position to another on the same route ( Figure 11). Move type 2-Intra swaps two customers on the route [25] (Figure 12). Finally, move type 3-Intra is defined by a 2-opt edge exchange method [26] that swaps the positions of two edges of the route (Figure 13).
Move types 1-intra, 2-intra, and 3-intra are, one at a time and in this order, applied to the route to perform all possible moves until no further improvement in route cost is found.
is procedure is applied to the routes of the initial solution and every time interroute moves find a feasible solution during SA to examine other configurations of the new route.

SA Algorithm.
SA is an algorithmic approach for solving combinatorial optimization problems [27,28] that has been widely applied to vehicle routing problems. e proposed SA uses neighborhood structures NS 1 , NS 2 , and NS 3 to generate a candidate solution S ′ from a current solution S; this solution is accepted as the new current solution if it is better than S.
To avoid local optima, a worse solution may be accepted subject to the acceptance probability function p(T, Δ) � exp(− Δ/T), where Δ � cos t(S ′ ) − cos t(S) and T is a parameter of SA called temperature, which decreases during the process according to T k � 0.9 · T k− 1 to enforce the convergence of the search. At the beginning of the algorithm, T is assigned an initial value T 0 . As suggested in [29,30], T 0 � − Δ max /ln(p 0 ) is defined to make any worst solution that leads to a cost increase Δ max accepted with a fixed probability p 0 . To estimate Δ max , move types 1-Inter, 2-Inter, and 3-Inter are randomly applied to the initial solution, and the value of Δ max is then estimated by the maximum absolute difference observed in the costs of two neighboring solutions.
For each temperature, a sequence of Len moves is carried out to explore the search space. If Len is too small, the solution space is not fully explored, whereas if Len is too large, the computational time required is rather long. Only feasible solutions are considered in our algorithm. erefore, owing to the difficulty in finding feasible solutions to some instances of the problem, a maximum number of attempts L � 10 · Len is established for each temperature. For each feasible solution found, the local search procedure is employed to improve the routes that have been changed. To accelerate the speed of the algorithm, a condition is set such that if it finds a solution that is better than the best solution, that solution is replaced, and the current temperature is updated. e algorithm ends when the current temperature is lower than 0.01. Algorithm 2 provides a framework for the proposed SA methodology.

Computational Experiments
is section reports the results of computational experiments performed with the proposed mathematical formulation and the SA algorithm. Since we did not find any instances for the 2L-HFFVRP, we defined a set of instances by extending to the 2L-HFFVRP some 2L-HFVRP instances from the literature [18]. First, a detailed discussion of the benchmark instances is provided. Subsequently, the simulation results of benchmark instances are presented.

Benchmark Instances.
To verify the effectiveness of the mathematical model and of the SA algorithm, small-scale and large-scale instances were used to estimate the results. ese instances are described in this section. Iori et al. [8] and Gendreau et al. [9] proposed a dataset including 36 2L-CVRP instances, in which the number of customers varies from 15 to 255 and each instance has five classes. In class 1, each customer demands one item with unit width and length, so the problems of this class are similar to the pure CVRP. In classes 2-5, for each customer i, a set of m i items with uniform distribution in the range (1, class number) is created. Each of these items is randomly classified, with an equal probability, into one of three possible shapes, namely, vertical, homogeneous, and horizontal, and the dimensions of the items are randomly generated in a given range. Leung et al. [18] modified 2L-CVRP instances to create 2L-HFVRP instances. e data for each benchmark were generated by eliminating the limit on the number of vehicles and adding information about for each customer i do k � 1; while customer i cannot be served by vehicle k do k � k + 1; end while insert customer i into set CustomersList k ; end for for k � 1 to v do ToInsert k ←CustomersList k ; end for Route construction: for k � v to 1 do while ToInsert k is not empty do if it is possible to insert some customer of set ToInsert k into any route R j (j � k, . . . , v) then Execute the feasible insertion of customer i into route R j , which minimizes the insertion cost; Delete customer i from the set ToInsert k ; else for j � 1 to v do ToInsert j ←CustomersList j ; Empty route R j ; end for Randomly select a vehicle j � 1, . . . , v and a customer i ∈ ToInsert j ; Insert customer i into route R j ; Delete customer i from the set ToInsert j ; go to Route construction; end if end while end for return generated solution ALGORITHM 1: Pseudocode for constructing an initial solution.      capacity, loading surface, and fixed and variable costs for each of the four types of vehicles considered, namely, A, B, C, and D. For the five classes of the same instance, the vehicle types are the same. For 2L-HFFVRP, we generate benchmark data by limiting the number of vehicles of types A, B, C, and D. For the five classes of the same instance, the number of vehicles of each type was the same. Table 3 provides details on the number of vehicles of each type considered for each instance. Instances of class 1 correspond to the pure HFFVRP since every route is feasible in terms of the loading constraints.

SA Parameter Setting.
e proposed SA algorithm has two parameters that need to be set by the user: the probability p 0 of accepting any worst solution at the initial temperature and the number Len of feasible moves to be performed at each temperature. To yield experimental results with better convergence, the parameters of SA are determined by conducting experiments on a set of largescale instances using various parameter values, involving both unrestricted and sequential problem versions. e suggested parameter values, along with the tested ranges, are reported in Table 4.

Results and Discussion.
e 2L-HFFVRP mathematical formulation of Section 3 and the SA algorithm of Section 4 were implemented in the Visual Basic .NET programming language. e mathematical model was solved by CPLEX 12.6 with the default configuration parameters. All the computational tests were executed in a machine with the following specifications: a 1.4 GHz Intel Core i5 processor with 8 GB RAM memory running a Windows 10 operating system. e pure HFFVRP was solved for problems of class 1. For problems of classes 2-5, the 2L-HFFVRP was solved regarding both unrestricted and sequential loading.
Small-scale problems were solved using the mathematical model and the SA algorithm, and the results are summarized in Tables 5-7. To evaluate the model performance, the time taken by CPLEX to solve each model was limited to 3600 seconds. With respect to the quality of the solution obtained by CPLEX, there were four possible cases: (i) the optimal solution is obtained; (ii) a nonoptimal solution is obtained, with CPLEX exceeding the time limit; (iii) no solution is obtained, with CPLEX exceeding the time limit; and (iv) there is insufficient computer memory to solve the model. e last two cases are represented in the tables by the symbol "-." Concerning SA, for each instance, five replications of the algorithm were executed, and the best results were obtained. e results for HFFVRP small-scale problems of class 1 are presented in Tables 5, which shows, for each instance, the number of customers, the solution cost and the runtime (in seconds) of CPLEX and SA for solving the problem, and the percentual difference between the solution costs of CPLEX and SA. Tables 6 and 7 show these results for the unrestricted and sequential 2L-HFFVRP, respectively, for small-scale problems of classes 2-5 and include the number of items to be loaded for each instance.
According to the results, the problems become more complex as the numbers of customers and items to be loaded increase. By comparing the results for the HFFVRP (Table 5) and for the unrestricted 2L-HFFVRP (Table 6), we notice that, for the same routing problem configuration, some solution costs deteriorate due to loading constraints. In fact, many routes become infeasible as these constraints are considered. By comparing the results for the unrestricted and sequential 2L-HFFVRP (Tables 6 and 7), we find that a few solution costs deteriorate when sequential loading constraints are embedded into the problem, mainly for problems with more customers and items. Inst.
Vehicle type  A  B  C  D  1  0  2  3  2  2  0  0  3  3  3  0  3  5  4  1  3  4  5  0  2  2  2  6  0  4  5  7  0  3  3  2  8  0  3  1  3  9  0  1  8  10  0  5  4  2  11  0  1  5  3  12  0  6  14  13  0  1  6  14  1  3  3  4  15  0  5  5  4  16  0  10  7  17  1  3  13  5  18  0  5  7  19  1  4  8  8  20  1  7  14  6  21  2  6  12  11  22  1  6  8  12  23  0  8  10  11  24  0  4  9  13  25  0  10  12  15  26  1  14  8  15  27  1  6  18  14  28  2  9  17  16  29  0  17  24  30  3  9  17  22  31  1  17  36  26  32  3  28  20  26  33  1  19  24  29  34  1  33  29  33  35  6  22  46  36  7  33  46  9 14 Mathematical Problems in Engineering   Tables 5-7 reveal that when the numbers of customers and items to be loaded are small, CPLEX can quickly find an optimal solution. However, as the numbers of customers and items increase, the solution time of CPLEX increasingly lengthens. Note that CPLEX does not find the optimal solutions in the stated time for the unrestricted problem P10 of classes 3 and 5, for the sequential problem P9 of classes 4-5, or for any class of sequential problem P10. e results provided by SA are very close to optimality. According to Table 5, the SA algorithm finds the optimal solution for 8 of the 10 HFFVRP instances of class 1, and the average gap is 0.39%. Table 6 shows that SA finds the optimal solution for 24 of the 32 unrestricted 2L-HFFVRP instances of classes 2-5 and that the average gap is 0.61%. Table 7 shows that for 24 out of the 32 sequential 2L-HFFVRP instances of classes 2-5, the optimal solution is obtained by SA. For instances P9 of class 4 and P10 of class 2 (marked with an " * " in Table 7), the CPLEX results are not proven to be optimal by the software. Nevertheless, solutions obtained by SA for these instances of the sequential 2L-HFFVRP are equal to the proven optimal solutions obtained by CPLEX for the unrestricted 2L-HFFVRP. erefore, we can say that the SA solutions for sequential instances P9 of class 4 and P10 of class 2 are optimal. In addition, for sequential problem P10 of class 5, CPLEX does not manage finding a solution in the stated time, but the solution obtained by SA is equal to the CPLEX solution for the unrestricted correspondent problem; thus, for this instance, we consider the SA solution to be very good, if not optimal. For three instances (P9 of class 5 and P10 of classes 2 and 3), the solutions obtained by SA are better than solutions found by CPLEX.
e average gap between the SA algorithm and CPLEX is 0.004%. e SA computational time is shorter than the CPLEX computational time only for instances with more customers and items-namely, instance P10 of class 1; instances P9 and P10 for unrestricted problems; and instances P7, P8, P9, and P10 for sequential problems. Nevertheless, the average computational times of SA are 28.25 seconds for the HFFVRP, 21.08 seconds for the unrestricted 2L-HFFVRP, and 17.23 seconds for the sequential 2L-HFFVRP, against the 76.62 seconds, 347.78 seconds, and 655.53 seconds of CPLEX, respectively. erefore, the SA algorithm is effective in solving problems with larger numbers of customers and items. e large-scale problems described in Section 5.2 are solved by the SA algorithm. Table 8 presents the results for the HFFVRP instances of class 1 and the average results for the unrestricted and sequential 2L-HFFVRP instances of classes 2-5. e results of Table 8 reinforce the idea that loading constraints worsen the solution costs of the HFFVRP and that sequential loading produces slightly higher solution costs than the related unrestricted problem. e computational times are longer for sequential problems. On average, the time spent by SA to solve the sequential 2L-HFFVRPs is 4 times longer than the time spent solving the unrestricted 2L-HFFVRPs. In fact, computational experiments reveal that the major difficulty in solving the 2L-HFFVRP concerns finding feasible loadings. As the number of vehicles is limited, a reasonable part of the computational time is spent finding an initial feasible solution for some instances.
In this sense, the impact of the three new packing heuristics, Heur 6 , Heur 7 , and Heur 8 , proposed in this study is analyzed. During the tests, the packing heuristic that produces each feasible loading is recorded. e ratios between feasible loadings found by heuristics Heur 6 , Heur 7 , and Heur 8 and all found feasible loadings for unrestricted and sequential problems are presented in Table 9. For instance, 3.11% of all the feasible loadings found during the tests is obtained by heuristic Heur 7 . is means that no other previous heuristic can find these feasible loadings. erefore, all three new heuristics are able to find new feasible loadings that cannot be obtained by the other previous heuristics. Heuristic Heur 7 stands out regarding managing more complex loadings in both unrestricted and sequential problems. For sequential problems, heuristic Heur 6 manages to obtain new feasible loadings. In fact, heuristic Heur 6 tends to form rectangular blocks of items, which favors sequential loading of blocks of items of the same customer. Although heuristic Heur 8 is the last to be executed, the results in Table 9 indicate that this heuristic is able to find feasible solutions not found by any of the previous heuristics.

Managerial Insights
Based on our results, several managerial implications can be used to help logistic providers and managers in optimizing the routes and loading configurations of a heterogeneous fleet of vehicles. e results presented in the previous section have shown the efficiency and robustness of the SA algorithm in solving unrestricted and sequential 2L-HFFVRP. e results of this method can be effectively used in decisionmaking at the strategic and tactical levels, leading to better use of the fleet of vehicles and cost reduction.
ere is a trade-off between solution cost and computational time. Better solutions can be obtained by increasing the values of parameters p 0 and Len, but the computational time of the algorithm will also increase. According to the period of time the company usually has to plan routes, it is possible to set the algorithm parameters to obtain better solutions within this time.
We compared sequential loading to unrestricted loading and observed that, on average, for classes 2-5 in 36 instances, the total transportation cost increased by around 3% when sequential loading is considered. erefore, if rearranging other customers' items at each customer site is not a problem for the company, especially regarding the time these operations will take, unrestricted loading should be considered, as it implies lower cost solutions compared to sequential loading. Nevertheless, it is worth noticing that, in some cases, unloading and reloading operations are very timeconsuming. For cases with up to 40 customers, we observed that the percentual difference between solution costs in unrestricted and sequential problems is, on average, smaller than that for cases with more customers. In these cases, the time saving of considering sequential loading is worth analyzing.

Conclusions
In this paper, the 2L-HFFVRP, a new variant of the VRP, is presented. In this problem, each customer demands a set of rectangular two-dimensional items, and the objective is to find the minimum cost delivery routes for a limited set of vehicles with different capacities, fixed and variable operating costs, and a rectangular two-dimensional loading surface. Both the unrestricted and the sequential versions of the problem are handled.
is problem is interesting in terms of both theoretical complexity and real-world applications. To the best of our knowledge, the 2L-HFFVRP has not been previously addressed in the literature, and nor have the sequential loading constraints.
To solve this problem, a mixed integer linear programming model was formulated. Computational tests were performed with 10 small-scale instances regarding five classes of problems.
e results show that the model is consistent and properly represents the problem treated and that this approach is able to solve problems where the numbers of customers and items are relatively small. erefore, with the developed model and for small cases, it is possible to determine the optimal solutions, or in some cases a good lower bound, in order to be able to compare the results achieved with nonexact methods or have a base of comparison of their performances. Moreover, the proposed model can be useful for motivating future research exploring exact methods for solving 2L-HFFVRP.
Owing to the complexity of the problem model, a hybrid algorithm that involves SA and packing heuristics was proposed, and three new packing heuristics were developed. By testing the proposed framework on 36 large-scale instances, each within five classes, we verified that the proposed methodology can solve this difficult problem in an acceptable computational time; hence, the proposed SA algorithm can further be assessed to solve different variants of 2L-HFFVRP. e three new packing heuristics were also proven to be effective in finding feasible loadings not found by the other heuristics, thus increasing the probability of obtaining a feasible loading.
Our study has some limitations and can be extended in several aspects in future research. It would be reasonable to examine the more realistic scenarios that sometimes arise in practical applications. First, concerning the routing problem, future work can integrate practical constraints such as split delivery and time windows. Second, in terms of two-dimensional loading configuration, we examined only oriented problems.
us, future research can consider items rotation as it could substantially improve vehicle occupation. ird, although the proposed SA algorithm for solving 2L-HFFVRP has obtained good solutions, there is still room for improvement, which may be achieved by proposing more sophisticated and superior routing and packing strategies.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this article.