Abstract
This paper addresses the twodimensional 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 twodimensional 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 twodimensional loading surfaces. We formulate a mixed integer linear programming model to obtain optimal solutions for smallscale problems. To obtain solutions for largescale 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 twodimensional loading heterogeneous fleet vehicle routing problem. The 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 largescale 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.
1. Introduction
In many organizations, the management of distribution activities is a major decisionmaking 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].
The 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. The most common costs associated with the problem are related to driving distance. There are several variants of the VRP due to additional constraints encountered in realworld 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 reallife applications in the distribution management context involve the transportation of rectangularshaped 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 twodimensional 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–7]. Since routing and packing are well known NPhard problems, combining them leads to an extremely challenging optimization problem.
The twodimensional loading capacitated vehicle routing problem (2LCVRP) is one of the first approaches to integrate vehicle routing and loading problems. In this problem, customers demand a set of rectangular twodimensional items, and identical vehicles have a twodimensional loading surface. In realworld problems, in contrast to the 2LCVRP, enterprises own a diverse fleet of vehicles, which offers the flexibility to design a more profitable distribution plan. The twodimensional loading heterogeneous fleet vehicle routing problem (2LHFVRP) treats the 2LCVRP 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 realworld 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 (2LHFFVRP). In the 2LHFFVRP, 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. The 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 rearloaded, and load rearrangement at the customer site can be difficult, timeconsuming, 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 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.
The purpose of this paper is to present two approaches for solving 2LHFFVRP. First, a mixed integer linear programming model is formulated to obtain optimal solutions for smallscale 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 largescale 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.
The 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.
2. Literature Review
Integrated vehicle routing and loading problems have been increasingly studied, owing to their relevance to logistics distribution systems [5–7]. To the best of our knowledge, there are no studies on the 2LHFFVRP to date. This section will focus on the previous studies on the 2LCVRP and the 2LHFVRP. Both these problems combine vehicle routing and twodimensional 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 NPhard problems, so exact methods are suitable for solving only smallscale problems, while metaheuristic algorithms are more favored for solving complex problems.
The 2LCVRP was first presented by Iori et al. [8]. They addressed the sequential oriented version of the problem using an exact methodology that employs a branchandcut algorithm to deal with the routing aspects of the problem, combined with a nested branchandbound 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 branchandbound 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. The 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 routingmemory 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 spacebased 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 2LCVRP with split delivery. This variant was also addressed by Ferreira et al. [17] through an exact branchandcut approach, where a tailored procedure that includes the computation of lower bounds, a constructivebased heuristic, and a constraint programming model are proposed to handle the packing problem.
Concerning 2LHFVRP, 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 biasedrandomized versions of routing and packing heuristics to solve oriented and nonoriented unrestricted problems. The 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 twodimensional loading constraints, and this study is shown in Table 1.
Our literature review of the vehicle routing problems that considers a heterogeneous fleet and twodimensional 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. Therefore, we address the 2LHFFVRP.
3. Problem Statement
The 2LHFFVRP is defined on a complete undirected graph , where is a set of vertices, the vertex 0 denotes the depot, and the vertices correspond to the positions of the customers . is a set of edges, and each edge is associated with a distance that denotes the distance from customer to customer . A set of vehicles is available at the depot. Each vehicle has a weight capacity , a loading surface with length and width , a fixed cost , and a variable cost . In general, a vehicle with a larger capacity usually has higher fixed and variable costs. The traveling cost of each edge by vehicle is . Each customer demands a set of rectangular items, denoted by , and the total weight of is equal to . Each item has a specific length and width .
For 2LHFFVRP, 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)The 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
The objective of the 2LHFFVRP is to find a set of routes of minimum cost that fulfill the constraints and satisfy all customers’ demands. The constraints mentioned above refer to the unrestricted 2LHFFVRP. Besides this problem version, this study also considers the sequential 2LHFFVRP.
4. Mixed Integer Linear Programming Formulation
This section proposes a mathematical formulation to represent the 2LHFFVRP, considering the characteristics described in Section 3. The notation used in the model, including indexes, sets, parameters, and variables, is presented in Section 4.1. The mathematical model is presented in the following stages: initially, a model for the HFFVRP is introduced in Section 4.2; then, we present constraints for vehicle loading in Section 4.3; and subsequently, sequential loading constraints are developed in Section 4.4.
4.1. Notations
The mixedinteger linear programming model depends on the sets, parameters, and variables described in Table 2.
4.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 fourindex 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 . Moreover, it is assumed that each route starts at vertex 0 in stage and finishes at vertex 0 in stage if only one vehicle is used. Therefore, the set of possible values for is .
The (routing) decision variables indicate whether vehicle goes directly from vertex to vertex in stage () or not (). It is assumed that , . The mathematical model for the HFFVRP is expressed as follows:
subject to
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 is visited in stage , then this customer has to be the starting point for some other customer in stage . Constraints (4) ensure that each vehicle leaves the depot at most once in stage 1. Constraints (5) ensure that if vehicle travels from customer to customer in stage , then in stage , the same vehicle must travel from customer to another customer . 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 , ; that is, the vehicles can leave the depot only in the first stage.
4.3. TwoDimensional Loading Constraints
In this section, we show how twodimensional loading considerations can be embedded into formulation (1)–(7) of the last subsection. These 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.
Twodimensional 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 of all items to be loaded in the vehicles as follows. The first elements of are the items of customer 1, the next elements are the items of customer 2, and so on. Therefore, each item becomes an item , so each customer demands a set of items , that is, . Each item has a specific length and width . The total number of items in set is . 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 frontleft corner is adopted, and the vehicle door through which loading and unloading operations are performed is placed on the side between coordinates and , as shown in Figure 2.
The (loading) decision variables indicate coordinates of the frontleft corner of item in the vehicle loading surface, and the variables indicate whether item is loaded into vehicle () or not (). Binary variables , , , and () indicate the placement of items relative to each other. , , , and are equal to 1 if item is, respectively, loaded behind, in front of, on the left side of, or on the right side of item . An example of the use of these variables is shown in Figure 3, which shows two items and loaded in a vehicle.
Assuming that , we have , , , and . Note that an item is said to be loaded behind or in front of an item only if there is no intersection of its projections on the axis. Analogously, item can be on the left or on the right of item if there is no intersection of its projections on the axis. Therefore, it is possible to have and , while .
Let be a sufficiently large number. The formulation for the 2LHFFVRP 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. This 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 must be loaded on the same vehicle that visits this customer. Finally, constraints (16)–(18) are the (loading) decision variable domain constraints.
4.4. 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. This 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 is served before customer , constraints (19) prevent the items of customer from blocking the unloading of the items of customer . These constraints verify whether customer is served in stage preceding stage in which customer is served. A particular case is illustrated in Figure 4 that shows an item of customer and an item of customer . If customer is not visited before customer , item of customer and item of customer have unrestricted coordinates. However, if customer is served before customer , there are two possibilities: if item of customer is completely on the left or on the right side of item of customer , we have , and there is no intersection of their projections on the axis, so the item does not block the other one. On the other hand, if item is not completely on the left or right side of item , an item is behind or in front of the other, and the constraints impose the condition that item must be in front of item . A similar explanation is valid for constraints (20) in the case that customer is served after customer .
The complete formulation for the 2LHFFVRP with sequential loading constraints is given by the objective function (1) with the constraints (2)–(20). The fourindex formulation is essential to properly model the 2LHFFVRP when sequential loading constraints are present. Both indices and are fundamental to the modeling of the problem: index is necessary to determine in which vehicle a given item is placed, and index , in turn, makes it possible to indicate the position of a customer in the route and set the sequential loading constraints. The 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.
This formulation can be used to optimally solve smallscale 2LHFFVRP cases. A metaheuristic algorithm for handling largescale problems is presented in the next section.
5. 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. These heuristics are described, and the proposed algorithm is explained.
5.1. Heuristics for TwoDimensional 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 () are used to solve the twodimensional 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 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 . The partial loading for this route and the available loading positions are shown in Figure 6. After the loading of the items of customer and before the loading of the items of customer , 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 () are based on the study by Zachariadis et al. [10] and adopt the following criteria to determine an item loading position. : BottomLeft Fill (axis): the selected position is that with the minimum axis coordinate, where ties are broken by the minimum axis coordinate. With this heuristic, packing tends to form strips parallel to the axis. : BottomLeft Fill (axis): the selected position is that with the minimum axis coordinate, where ties are broken by the minimum axis coordinate. With this heuristic, packing tends to form strips parallel to the axis. : 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. This heuristic tends to spread items to the edges of the loading surface and later fills the inner part of it. : 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. This heuristic tends to fill the inner part of the loading surface before filling the edges. : 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. Then, we introduce three new packing heuristics, namely, the Min Occupied Rectangular Area (), Max Touching Perimeter Select Item (), and Max Relative Touching Perimeter Select Item (). selects an item loading position according the following criteria. : 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 in the given position. This strategy makes the packing compact. Heuristics () 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 unique visit order. If sequential loading is considered, , , and 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, , , and are produced simply by sorting all items by decreasing length, width, and surface area, respectively. Note that although , , and are primarily designed to deal with the sequential problem, they succeed in producing feasible unrestricted loadings in numerous cases when , , and fail; therefore, they are also considered for solving unrestricted problems. Heuristics Heur and Heur 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. The criteria employed by Heur and Heur are as follows. : 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 . The combination of an item and a position that matches the largest perimeter is chosen, and the item is inserted. : 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. The combination of item and position that matches the largest ratio is selected, and the item is inserted.
The feasibility of a route is tested by heuristics to , in this sequence, such that heuristics () 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.
5.2. Initial Solution
The initial solution is the starting point of the improvement metaheuristic proposed in the following subsections. The 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 vehicles available in the depot, it is assumed that , , , , and . First, for each vehicle , a subset of customers that must be served by vehicle or a larger vehicle is defined. Subsequently, for each vehicle , the customers of set are added to a vehicle route according to the minimum cost insertion so that the route is feasible. The cost of adding customer to vehicle route depends on whether route is empty. If route is empty, the cost is based on the fixed cost of vehicle , plus the cost of traveling from the depot to customer and back. If the route is not empty, the cost of adding customer 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 in a particular route , 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. is the set of customers of that have not yet been inserted in any route. At the beginning of route construction, , and in the course of the algorithm, this set is updated as customers are added to routes.

5.3. Neighborhood Structures
SA explores the search space by performing moves to step from a current solution to a subsequent solution. The 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.
Three types of movements are employed to disturb a solution, each of which defines a neighborhood structure, , , or , selected randomly in each loop with equal probability. Move type 1Inter () performs a customer relocation move [25]. It reassigns a customer from their current route to a position on another route (Figure 8). This move type can make some routes empty, reducing the number of vehicles used.
Move type 2inter () performs route exchange [25]. It swaps the positions of the two customers on different routes (Figure 9).
Move type 3inter () is a variant of route interchange 2opt [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.
5.4. Local Search Mechanism
The 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. Three move types, similar to those for interroute improvement, define neighborhood structures as intraroute optimization. Move type 1Intra 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 2Intra swaps two customers on the route [25] (Figure 12). Finally, move type 3Intra is defined by a 2opt edge exchange method [26] that swaps the positions of two edges of the route (Figure 13).
Move types 1intra, 2intra, and 3intra 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. This 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.
5.5. SA Algorithm
SA is an algorithmic approach for solving combinatorial optimization problems [27, 28] that has been widely applied to vehicle routing problems. The proposed SA uses neighborhood structures , , and to generate a candidate solution from a current solution ; this solution is accepted as the new current solution if it is better than .
To avoid local optima, a worse solution may be accepted subject to the acceptance probability function , where and is a parameter of SA called temperature, which decreases during the process according to to enforce the convergence of the search. At the beginning of the algorithm, is assigned an initial value . As suggested in [29, 30], is defined to make any worst solution that leads to a cost increase accepted with a fixed probability . To estimate , move types 1Inter, 2Inter, and 3Inter are randomly applied to the initial solution, and the value of is then estimated by the maximum absolute difference observed in the costs of two neighboring solutions.
For each temperature, a sequence of moves is carried out to explore the search space. If is too small, the solution space is not fully explored, whereas if is too large, the computational time required is rather long. Only feasible solutions are considered in our algorithm. Therefore, owing to the difficulty in finding feasible solutions to some instances of the problem, a maximum number of attempts 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. The algorithm ends when the current temperature is lower than 0.01. Algorithm 2 provides a framework for the proposed SA methodology.

6. Computational Experiments
This 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 2LHFFVRP, we defined a set of instances by extending to the 2LHFFVRP some 2LHFVRP instances from the literature [18]. First, a detailed discussion of the benchmark instances is provided. Subsequently, the simulation results of benchmark instances are presented.
6.1. Benchmark Instances
To verify the effectiveness of the mathematical model and of the SA algorithm, smallscale and largescale instances were used to estimate the results. These instances are described in this section.
Iori et al. [8] and Gendreau et al. [9] proposed a dataset including 36 2LCVRP 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 , a set of 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 2LCVRP instances to create 2LHFVRP instances. The data for each benchmark were generated by eliminating the limit on the number of vehicles and adding information about 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 2LHFFVRP, 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.
To validate the mathematical formulation, 8 smallscale instances, namely, P3, P4, P5, P6, P7, P8, P9, and P10, are generated based on instance 1 for the 2LHFFVRP. For classes 1–5, problems P3, P4, P5, P6, P7, P8, P9, and P10 take into account the first 3, 4, 5, 6, 7, 8, and 9 customers of instance 1 for the 2LHFFVRP, respectively. Instances of class 1 correspond to the pure HFFVRP since every route is feasible in terms of the loading constraints.
6.2. SA Parameter Setting
The proposed SA algorithm has two parameters that need to be set by the user: the probability of accepting any worst solution at the initial temperature and the number 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. The suggested parameter values, along with the tested ranges, are reported in Table 4.
6.3. Results and Discussion
The 2LHFFVRP mathematical formulation of Section 3 and the SA algorithm of Section 4 were implemented in the Visual Basic .NET programming language. The 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. The pure HFFVRP was solved for problems of class 1. For problems of classes 2–5, the 2LHFFVRP was solved regarding both unrestricted and sequential loading.
Smallscale 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. The 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.
The results for HFFVRP smallscale 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 2LHFFVRP, respectively, for smallscale 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 2LHFFVRP (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 2LHFFVRP (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.
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 45, or for any class of sequential problem P10.
The 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 2LHFFVRP instances of classes 2–5 and that the average gap is 0.61%. Table 7 shows that for 24 out of the 32 sequential 2LHFFVRP 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 2LHFFVRP are equal to the proven optimal solutions obtained by CPLEX for the unrestricted 2LHFFVRP. Therefore, 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. The average gap between the SA algorithm and CPLEX is 0.004%. The 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 2LHFFVRP, and 17.23 seconds for the sequential 2LHFFVRP, against the 76.62 seconds, 347.78 seconds, and 655.53 seconds of CPLEX, respectively. Therefore, the SA algorithm is effective in solving problems with larger numbers of customers and items.
The largescale 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 2LHFFVRP instances of classes 2–5. The 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. The computational times are longer for sequential problems. On average, the time spent by SA to solve the sequential 2LHFFVRPs is 4 times longer than the time spent solving the unrestricted 2LHFFVRPs. In fact, computational experiments reveal that the major difficulty in solving the 2LHFFVRP 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, , and , proposed in this study is analyzed. During the tests, the packing heuristic that produces each feasible loading is recorded. The ratios between feasible loadings found by heuristics Heur , , and 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 . This means that no other previous heuristic can find these feasible loadings. Therefore, all three new heuristics are able to find new feasible loadings that cannot be obtained by the other previous heuristics. Heuristic stands out regarding managing more complex loadings in both unrestricted and sequential problems. For sequential problems, heuristic manages to obtain new feasible loadings. In fact, heuristic tends to form rectangular blocks of items, which favors sequential loading of blocks of items of the same customer. Although heuristic 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.
7. 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. The results presented in the previous section have shown the efficiency and robustness of the SA algorithm in solving unrestricted and sequential 2LHFFVRP. The 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.
There is a tradeoff between solution cost and computational time. Better solutions can be obtained by increasing the values of parameters and , 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. Therefore, 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.
8. Conclusions
In this paper, the 2LHFFVRP, a new variant of the VRP, is presented. In this problem, each customer demands a set of rectangular twodimensional 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 twodimensional loading surface. Both the unrestricted and the sequential versions of the problem are handled. This problem is interesting in terms of both theoretical complexity and realworld applications. To the best of our knowledge, the 2LHFFVRP 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 smallscale instances regarding five classes of problems. The 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. Therefore, 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 2LHFFVRP.
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 largescale 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 2LHFFVRP. The 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 twodimensional loading configuration, we examined only oriented problems. Thus, future research can consider items rotation as it could substantially improve vehicle occupation. Third, although the proposed SA algorithm for solving 2LHFFVRP 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
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.