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

Authors: Lilian Caroline Xavier Candido; Luzia Vidal de Souza
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This 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. The results indicate that the proposed model correctly describes the problem and can solve small-scale 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.


---

## Body

## 1. 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].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 cost-effective 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 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–7]. Since routing and packing are well known NP-hard problems, combining them leads to an extremely challenging optimization problem.The 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. The 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 two-dimensional loading constraints, which is called the two-dimensional 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. The objective is to describe the most cost-effective 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 by90° 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.Figure 1
Unrestricted and sequential loading.The 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 four-index 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.The remainder of this paper is organized as follows. Section2 briefly reviews the literature on combined two-dimensional 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 2L-HFFVRP to date. This 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.The 2L-CVRP was first presented by Iori et al. [8]. They 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. 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 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. This variant was also addressed by Ferreira et al. [17] through an exact branch-and-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. 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 two-dimensional loading constraints, and this study is shown in Table 1.Table 1
Contributions of previous studies.
AuthorsLimited fleetItems rotationSequential loadingLeung et al. [18]NoNoYesDominguez et al. [19]NoYesNoZhang et al. [20]NoNoYesSabar et al. [21]NoNoNoThis paperYesNoYesOur 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. Therefore, we address the 2L-HFFVRP.


## 3. Problem Statement

The 2L-HFFVRP is defined on a complete undirected graphG=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 dijdii=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 Qk, a loading surface with length Lk and width Wk, a fixed cost Fk, and a variable cost Vk. In general, a vehicle with a larger capacity usually has higher fixed and variable costs. The traveling cost of each edge i,j∈A by vehicle k is Cijk=Vk⋅dij. Each customer ii=1,…,n demands a set of mi rectangular items, denoted by ITi, and the total weight of ITi is equal to Di. Each item Iir∈ITir=1,2,…,mi has a specific length lir and width wir.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)
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 allowedThe objective of the 2L-HFFVRP 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 2L-HFFVRP. Besides this problem version, this study also considers the sequential 2L-HFFVRP.


## 4. Mixed Integer Linear Programming Formulation

This section proposes a mathematical formulation to represent the 2L-HFFVRP, considering the characteristics described in Section3. 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 mixed-integer linear programming model depends on the sets, parameters, and variables described in Table2.Table 2
Model sets, parameters, and variables.
Indexes/setsi,j,g∈N=0,…,nSet of vertices from the depot (n=0) to customers (n>0)k∈V=1,…,vSet of vehiclest∈T=1,…,n+1Set of stages in which a customer can be served by a vehiclep,q∈I=1,…,∑i=1nmiSet of items demanded by all customers. Each customeri demands a set ITi=∑j=1i−1mj+1,…,∑j=1imjParametersnNumber of customersdijDistance from each customeri to each customer jvNumber of vehiclesQkWeight capacity of each vehiclekLkLength of the loading surface of each vehiclekWkWidth of the loading surface of each vehiclekFkFixed cost of each vehiclekVkVariable cost of each vehiclekmiNumber of items demanded by each customeriDiTotal weight demanded by each customerilpLength of each itempwpWidth of each itempVariableszijktBinary variable that indicates whether vehiclek goes directly from vertex i to vertex j in stage t (zijkt=1) or not zijkt=0′xpContinuous variable that indicates thex-coordinate of the front-left corner of item p on the vehicle loading surfaceypContinuous variable that indicates they-coordinate of the front-left corner of item p in the vehicle loading surfacespkBinary variable that indicates whether itemp is loaded into vehicle kspk=1′ or not spk=0′αpqBinary variable that indicates whether itemp is loaded behind item qαpq=1′ or not αpq=0′βpqBinary variable that indicates whether itemp is loaded in front of item qβpq=1′ or not βpq=0′γpqBinary variable that indicates whether itemp is loaded on the left side of item qγpq=1′ or not γpq=0′δpqBinary variable that indicates whether itemp is loaded on the right side of item qδpq=1′ or not δpq=0′


### 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 time-dependent 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 staget. 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. Therefore, the set of possible values for t is T=1,…,n+1.The (routing) decision variableszijkt indicate whether vehicle k goes directly from vertex i to vertex j in stage t (zijkt=1) or not (zijkt=0). It is assumed that dii=0, i,i∈A. The mathematical model for the HFFVRP is expressed as follows:(1)Minimize∑k∈V∑j∈N\0Fkz0jk1+∑i∈N∑j∈N∑k∈V∑t∈TVkdijzijkt,subject to(2)∑j∈N∑k∈V∑t∈Tzijkt=1,i∈N\0,(3)∑j∈N∑k∈V∑t∈T\1tzijkt−∑g∈N∑k∈V∑t∈Ttzgikt=1,i∈N\0,(4)∑j∈N\0z0jk1≤1,k∈V,(5)∑j∈Nzijk,t+1−∑g∈Nzgikt=0,i∈N\0,k∈V,t∈T\n+1,(6)∑i∈N\0∑j∈N∑t∈T\1Dizijkt≤Qk,k∈V,(7)zijkt∈0,1,i,j∈A,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 z0jkt=0, ∀j∈N\0,k∈V,t∈T\1; that is, the vehicles can leave the depot only in the first stage.


### 4.3. 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. 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.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 two-dimensional case and to the vehicle routing context.Let us create a setI of all items to be loaded in the vehicles as follows. The first m1 elements of I are the items of customer 1, the next m2 elements are the items of customer 2, and so on. Therefore, each item r∈ITi becomes an item p=∑j=1i−1mj+r, so each customer i demands a set of items p∈Ip=∑j=1i−1mj+1,…,∑j=1imj, that is, ITi=∑j=1i−1mj+1,…,∑j=1imj. Each item p∈I has a specific length lp and width wp. The total number of items in set I is B=∑i=1nmi. 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 coordinatesL,0 and L,W, as shown in Figure 2.Figure 2
Vehicle loading surface.The (loading) decision variables indicate coordinatesxp,yp of the front-left corner of item p in the vehicle loading surface, and the variables spk indicate whether item p is loaded into vehicle k (spk=1) or not (spk=0). Binary variables αpq, βpq, γpq, and δpq (p,q∈I,p<q) indicate the placement of items relative to each other. αpq, βpq, γpq, and δpq are equal to 1 if item p is, respectively, loaded behind, in 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.Figure 3
Relative positions of two items.Assuming thatp<q, we have αpq=1, βpq=0, γ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. Therefore, it is possible to have αpq+βpq≤1 and γpq+δpq≤1, while 1≤αpq+βpq+γpq+δpq≤2.LetM be a sufficiently large number. The 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:(8)xp+lp≤xq+1−αpqM,p,q∈I,p<q,(9)xq+lq≤xp+1−βpqM,p,q∈I,p<q,(10)yp+wp≤yq+1−γpqM,p,q∈I,p<q,(11)yq+wq≤yp+1−δpqM,p,q∈I,p<q,(12)αpq+βpq+γpq+δpq≥spk+sqk−1,p,q∈I,p<q,k∈V,(13)xp+lp≤Lk+1−spkM,p∈I,k∈V,(14)yp+wp≤Wk+1−spkM,p∈I,k∈V,(15)∑p∈ITispk=mi∑j∈N∑t∈Tzijkt,i∈N\0,k∈V,(16)xp,yp∈ℝ,p∈I,(17)spk∈0,1,p∈I,k∈V,(18)αpq,βpq,γpq,δpq∈0,1,p,q∈I,p<q.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 i must be loaded on the same vehicle k 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:(19)xq+lq≤xp+2+γpq+δpq−∑g∈Nzgjkt2+∑g∈Nzgikt1M,i,j∈N\0,i<j,p∈ITi,q∈ITj,k∈V,t1,t2∈T,t1<t2,(20)xp+lp≤xq+2+γpq+δpq−∑g∈Nzgikt2+∑g∈Nzgjkt1M,i,j∈N\0,i<j,p∈ITi,q∈ITj,k∈V,t1,t2∈T,t1<t2.If customeri is served before customer j, constraints (19) prevent the items of customer j from blocking the unloading of the items of customer i. These constraints verify whether customer i is served in stage t1 preceding stage t2 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 γ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.Figure 4
Relative positions of the items of two customers.The complete formulation for the 2L-HFFVRP with sequential loading constraints is given by the objective function (1) with the constraints (2)–(20). The 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. 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 small-scale 2L-HFFVRP cases. A metaheuristic algorithm for handling large-scale problems is presented in the next section.


## 4.1. Notations

The mixed-integer linear programming model depends on the sets, parameters, and variables described in Table2.Table 2
Model sets, parameters, and variables.
Indexes/setsi,j,g∈N=0,…,nSet of vertices from the depot (n=0) to customers (n>0)k∈V=1,…,vSet of vehiclest∈T=1,…,n+1Set of stages in which a customer can be served by a vehiclep,q∈I=1,…,∑i=1nmiSet of items demanded by all customers. Each customeri demands a set ITi=∑j=1i−1mj+1,…,∑j=1imjParametersnNumber of customersdijDistance from each customeri to each customer jvNumber of vehiclesQkWeight capacity of each vehiclekLkLength of the loading surface of each vehiclekWkWidth of the loading surface of each vehiclekFkFixed cost of each vehiclekVkVariable cost of each vehiclekmiNumber of items demanded by each customeriDiTotal weight demanded by each customerilpLength of each itempwpWidth of each itempVariableszijktBinary variable that indicates whether vehiclek goes directly from vertex i to vertex j in stage t (zijkt=1) or not zijkt=0′xpContinuous variable that indicates thex-coordinate of the front-left corner of item p on the vehicle loading surfaceypContinuous variable that indicates they-coordinate of the front-left corner of item p in the vehicle loading surfacespkBinary variable that indicates whether itemp is loaded into vehicle kspk=1′ or not spk=0′αpqBinary variable that indicates whether itemp is loaded behind item qαpq=1′ or not αpq=0′βpqBinary variable that indicates whether itemp is loaded in front of item qβpq=1′ or not βpq=0′γpqBinary variable that indicates whether itemp is loaded on the left side of item qγpq=1′ or not γpq=0′δpqBinary variable that indicates whether itemp is loaded on the right side of item qδpq=1′ or not δpq=0′


## 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 time-dependent 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 staget. 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. Therefore, the set of possible values for t is T=1,…,n+1.The (routing) decision variableszijkt indicate whether vehicle k goes directly from vertex i to vertex j in stage t (zijkt=1) or not (zijkt=0). It is assumed that dii=0, i,i∈A. The mathematical model for the HFFVRP is expressed as follows:(1)Minimize∑k∈V∑j∈N\0Fkz0jk1+∑i∈N∑j∈N∑k∈V∑t∈TVkdijzijkt,subject to(2)∑j∈N∑k∈V∑t∈Tzijkt=1,i∈N\0,(3)∑j∈N∑k∈V∑t∈T\1tzijkt−∑g∈N∑k∈V∑t∈Ttzgikt=1,i∈N\0,(4)∑j∈N\0z0jk1≤1,k∈V,(5)∑j∈Nzijk,t+1−∑g∈Nzgikt=0,i∈N\0,k∈V,t∈T\n+1,(6)∑i∈N\0∑j∈N∑t∈T\1Dizijkt≤Qk,k∈V,(7)zijkt∈0,1,i,j∈A,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 z0jkt=0, ∀j∈N\0,k∈V,t∈T\1; that is, the vehicles can leave the depot only in the first stage.


## 4.3. 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. 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.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 two-dimensional case and to the vehicle routing context.Let us create a setI of all items to be loaded in the vehicles as follows. The first m1 elements of I are the items of customer 1, the next m2 elements are the items of customer 2, and so on. Therefore, each item r∈ITi becomes an item p=∑j=1i−1mj+r, so each customer i demands a set of items p∈Ip=∑j=1i−1mj+1,…,∑j=1imj, that is, ITi=∑j=1i−1mj+1,…,∑j=1imj. Each item p∈I has a specific length lp and width wp. The total number of items in set I is B=∑i=1nmi. 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 coordinatesL,0 and L,W, as shown in Figure 2.Figure 2
Vehicle loading surface.The (loading) decision variables indicate coordinatesxp,yp of the front-left corner of item p in the vehicle loading surface, and the variables spk indicate whether item p is loaded into vehicle k (spk=1) or not (spk=0). Binary variables αpq, βpq, γpq, and δpq (p,q∈I,p<q) indicate the placement of items relative to each other. αpq, βpq, γpq, and δpq are equal to 1 if item p is, respectively, loaded behind, in 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.Figure 3
Relative positions of two items.Assuming thatp<q, we have αpq=1, βpq=0, γ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. Therefore, it is possible to have αpq+βpq≤1 and γpq+δpq≤1, while 1≤αpq+βpq+γpq+δpq≤2.LetM be a sufficiently large number. The 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:(8)xp+lp≤xq+1−αpqM,p,q∈I,p<q,(9)xq+lq≤xp+1−βpqM,p,q∈I,p<q,(10)yp+wp≤yq+1−γpqM,p,q∈I,p<q,(11)yq+wq≤yp+1−δpqM,p,q∈I,p<q,(12)αpq+βpq+γpq+δpq≥spk+sqk−1,p,q∈I,p<q,k∈V,(13)xp+lp≤Lk+1−spkM,p∈I,k∈V,(14)yp+wp≤Wk+1−spkM,p∈I,k∈V,(15)∑p∈ITispk=mi∑j∈N∑t∈Tzijkt,i∈N\0,k∈V,(16)xp,yp∈ℝ,p∈I,(17)spk∈0,1,p∈I,k∈V,(18)αpq,βpq,γpq,δpq∈0,1,p,q∈I,p<q.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 i must be loaded on the same vehicle k 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:(19)xq+lq≤xp+2+γpq+δpq−∑g∈Nzgjkt2+∑g∈Nzgikt1M,i,j∈N\0,i<j,p∈ITi,q∈ITj,k∈V,t1,t2∈T,t1<t2,(20)xp+lp≤xq+2+γpq+δpq−∑g∈Nzgikt2+∑g∈Nzgjkt1M,i,j∈N\0,i<j,p∈ITi,q∈ITj,k∈V,t1,t2∈T,t1<t2.If customeri is served before customer j, constraints (19) prevent the items of customer j from blocking the unloading of the items of customer i. These constraints verify whether customer i is served in stage t1 preceding stage t2 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 γ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.Figure 4
Relative positions of the items of two customers.The complete formulation for the 2L-HFFVRP with sequential loading constraints is given by the objective function (1) with the constraints (2)–(20). The 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. 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 small-scale 2L-HFFVRP cases. A metaheuristic algorithm for handling large-scale 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 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 heuristicsHeuri (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.Figure 5
Process of loading items.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 is0−i1−i2−i3−0. The partial loading for this route and the available loading positions are shown in Figure 6. After the loading of the items of customer i2 and before the loading of the items of customer i1, the highlighted positions are excluded from the list of available loading positions, as they would generate loadings that would contravene the sequential loading constraints.Figure 6
Infeasible positions for sequential loading.HeuristicsHeuri (i=1,…,5) are based on the study by Zachariadis et al. [10] and adopt the following criteria to determine an item loading position.Heur1: 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.Heur2: 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.Heur3: 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.Heur4: 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.Heur5: 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 (Heur6), Max Touching Perimeter Select Item (Heur7), and Max Relative Touching Perimeter Select Item (Heur8). Heur6 selects an item loading position according the following criteria.Heur6: 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. This strategy makes the packing compact.HeuristicsHeuri (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 unique visit order. If sequential loading is considered, Ord_Seq1, Ord_Seq2, and Ord_Seq3 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_Un1, Ord_Un2, and Ord_Un3 are produced simply by sorting all items by decreasing length, width, and surface area, respectively. Note that although Ord_Seq1, Ord_Seq2, and Ord_Seq3 are primarily designed to deal with the sequential problem, they succeed in producing feasible unrestricted loadings in numerous cases when Ord_Un1, Ord_Un2, and Ord_Un3 fail; therefore, they are also considered for solving unrestricted problems.Heuristics Heur7 and Heur8 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 Heur7 and Heur8 are as follows.Heur7: 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 Heur3. The combination of an item and a position that matches the largest perimeter is chosen, and the item is inserted.Heur8: 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.Figure 7
Occupied rectangular area in the heuristic Min Occupied Rectangular Area.The feasibility of a route is tested by heuristicsHeur1 to Heur8, in this sequence, such that heuristics Heuri (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.


### 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.Givenv vehicles available in the depot, it is assumed that Q1≤Q2≤⋯≤Qv, F1≤F2≤⋯≤Fv, V1≤V2≤⋯≤Vv, L1≤L2≤⋯≤Lv, and W1≤W2≤⋯≤Wv. First, for each vehicle k, a subset CustomersListk 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 CustomersListk are added to a vehicle route k,…,v according to the minimum cost insertion so that the route is feasible. The 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 CustomersListk 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. ToInsertk is the set of customers of CustomersListk that have not yet been inserted in any route. At the beginning of route construction, ToInsertk=CustomersListk, and in the course of the algorithm, this set is updated as customers are added to routes.Algorithm 1:Pseudocode for constructing an initial solution.
for each customer idok=1;while customer i cannot be served by vehicle kdok=k+1;end whileinsert customeri into set CustomersListk;end forfork=1tovdoToInsertk←CustomersListk;end forRoute construction:fork=vto 1 dowhileToInsertk is not empty doif it is possible to insert some customer of set ToInsertk into any route Rj (j=k,…,v) thenExecute the feasible insertion of customeri into route Rj, which minimizes the insertion cost;Delete customeri from the set ToInsertk;elseforj=1tovdoToInsertj←CustomersListj;Empty routeRj;end forRandomly select a vehiclej=1,…,v and a customer i∈ToInsertj;Insert customeri into route Rj;Delete customeri from the set ToInsertj;go to Route construction;end ifend whileend forreturn generated solution


### 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,NS1, NS2, or NS3, selected randomly in each loop with equal probability. Move type 1-Inter (NS1) 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.Figure 8
Interroute relocation move.Move type 2-inter (NS2) performs route exchange [25]. It swaps the positions of the two customers on different routes (Figure 9).Figure 9
Interroute exchange move.Move type 3-inter (NS3) 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.Figure 10
Interroute interchange move.


### 5.4. Local Search Mechanism

The neighborhood structures presented in Section5.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 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).Figure 11
Intraroute relocation move.Figure 12
Intraroute exchange.Figure 13
Intraroute interchange move.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. 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 NS1, NS2, and NS3 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 functionpT,Δ=exp−Δ/T, where Δ=costS′−costS and T is a parameter of SA called temperature, which decreases during the process according to Tk=0.9⋅Tk−1 to enforce the convergence of the search. At the beginning of the algorithm, T is assigned an initial value T0. As suggested in [29, 30], T0=−Δmax/lnp0 is defined to make any worst solution that leads to a cost increase Δmax accepted with a fixed probability p0. 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 ofLen 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. Therefore, 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. The algorithm ends when the current temperature is lower than 0.01. Algorithm 2 provides a framework for the proposed SA methodology.Algorithm 2:Pseudocode of the SA methodology for the 2L-HFFVRP.
Input: problem data; SA parameters: p0; LenOutput: best solution S∗ found by the algorithmConstruct the initial solutionS0;Apply the local search mechanism toS0;/∗ settingT0∗/S←S0, Δmax=0, i=1, j=1;whilei<Lenandj<LdoRandomly selectNS from NS1,NS2,NS3, and obtain a solution S′∈NSS;ifS′ is feasible thenApply the local search mechanism toS′;ifcostS′−costS>ΔmaxthenΔmax=costS′−costS;end ifi=i+1;S←S′;ifcostS<costS0thenS0←S;end ifend ifj=j+1;end whileT0=−Δmax/lnp0;/∗ SA algorithm ∗/S←S0, S∗←S0, T=T0, i=1, j=1;while stop criteria are not met dowhilei<Lenandj<LdoRandomly selectNS from NS1,NS2,NS3, and obtain a solution S′∈NSS;ifS′ is feasible thenApply the local search mechanism toS′;ifcostS′<costSthenS←S′;ifcostS<costS∗thenS∗←S;end ifelseSetS←S′ with probability p=exp−costS′−costS/Tend ifi=i+1;end ifj=j+1end whileT=0.9⋅T, i=1, j=1;end whilereturnS∗


## 5.1. 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 heuristicsHeuri (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.Figure 5
Process of loading items.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 is0−i1−i2−i3−0. The partial loading for this route and the available loading positions are shown in Figure 6. After the loading of the items of customer i2 and before the loading of the items of customer i1, the highlighted positions are excluded from the list of available loading positions, as they would generate loadings that would contravene the sequential loading constraints.Figure 6
Infeasible positions for sequential loading.HeuristicsHeuri (i=1,…,5) are based on the study by Zachariadis et al. [10] and adopt the following criteria to determine an item loading position.Heur1: 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.Heur2: 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.Heur3: 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.Heur4: 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.Heur5: 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 (Heur6), Max Touching Perimeter Select Item (Heur7), and Max Relative Touching Perimeter Select Item (Heur8). Heur6 selects an item loading position according the following criteria.Heur6: 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. This strategy makes the packing compact.HeuristicsHeuri (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 unique visit order. If sequential loading is considered, Ord_Seq1, Ord_Seq2, and Ord_Seq3 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_Un1, Ord_Un2, and Ord_Un3 are produced simply by sorting all items by decreasing length, width, and surface area, respectively. Note that although Ord_Seq1, Ord_Seq2, and Ord_Seq3 are primarily designed to deal with the sequential problem, they succeed in producing feasible unrestricted loadings in numerous cases when Ord_Un1, Ord_Un2, and Ord_Un3 fail; therefore, they are also considered for solving unrestricted problems.Heuristics Heur7 and Heur8 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 Heur7 and Heur8 are as follows.Heur7: 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 Heur3. The combination of an item and a position that matches the largest perimeter is chosen, and the item is inserted.Heur8: 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.Figure 7
Occupied rectangular area in the heuristic Min Occupied Rectangular Area.The feasibility of a route is tested by heuristicsHeur1 to Heur8, in this sequence, such that heuristics Heuri (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.


## 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.Givenv vehicles available in the depot, it is assumed that Q1≤Q2≤⋯≤Qv, F1≤F2≤⋯≤Fv, V1≤V2≤⋯≤Vv, L1≤L2≤⋯≤Lv, and W1≤W2≤⋯≤Wv. First, for each vehicle k, a subset CustomersListk 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 CustomersListk are added to a vehicle route k,…,v according to the minimum cost insertion so that the route is feasible. The 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 CustomersListk 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. ToInsertk is the set of customers of CustomersListk that have not yet been inserted in any route. At the beginning of route construction, ToInsertk=CustomersListk, and in the course of the algorithm, this set is updated as customers are added to routes.Algorithm 1:Pseudocode for constructing an initial solution.
for each customer idok=1;while customer i cannot be served by vehicle kdok=k+1;end whileinsert customeri into set CustomersListk;end forfork=1tovdoToInsertk←CustomersListk;end forRoute construction:fork=vto 1 dowhileToInsertk is not empty doif it is possible to insert some customer of set ToInsertk into any route Rj (j=k,…,v) thenExecute the feasible insertion of customeri into route Rj, which minimizes the insertion cost;Delete customeri from the set ToInsertk;elseforj=1tovdoToInsertj←CustomersListj;Empty routeRj;end forRandomly select a vehiclej=1,…,v and a customer i∈ToInsertj;Insert customeri into route Rj;Delete customeri from the set ToInsertj;go to Route construction;end ifend whileend forreturn generated solution


## 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,NS1, NS2, or NS3, selected randomly in each loop with equal probability. Move type 1-Inter (NS1) 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.Figure 8
Interroute relocation move.Move type 2-inter (NS2) performs route exchange [25]. It swaps the positions of the two customers on different routes (Figure 9).Figure 9
Interroute exchange move.Move type 3-inter (NS3) 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.Figure 10
Interroute interchange move.


## 5.4. Local Search Mechanism

The neighborhood structures presented in Section5.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 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).Figure 11
Intraroute relocation move.Figure 12
Intraroute exchange.Figure 13
Intraroute interchange move.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. 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 NS1, NS2, and NS3 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 functionpT,Δ=exp−Δ/T, where Δ=costS′−costS and T is a parameter of SA called temperature, which decreases during the process according to Tk=0.9⋅Tk−1 to enforce the convergence of the search. At the beginning of the algorithm, T is assigned an initial value T0. As suggested in [29, 30], T0=−Δmax/lnp0 is defined to make any worst solution that leads to a cost increase Δmax accepted with a fixed probability p0. 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 ofLen 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. Therefore, 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. The algorithm ends when the current temperature is lower than 0.01. Algorithm 2 provides a framework for the proposed SA methodology.Algorithm 2:Pseudocode of the SA methodology for the 2L-HFFVRP.
Input: problem data; SA parameters: p0; LenOutput: best solution S∗ found by the algorithmConstruct the initial solutionS0;Apply the local search mechanism toS0;/∗ settingT0∗/S←S0, Δmax=0, i=1, j=1;whilei<Lenandj<LdoRandomly selectNS from NS1,NS2,NS3, and obtain a solution S′∈NSS;ifS′ is feasible thenApply the local search mechanism toS′;ifcostS′−costS>ΔmaxthenΔmax=costS′−costS;end ifi=i+1;S←S′;ifcostS<costS0thenS0←S;end ifend ifj=j+1;end whileT0=−Δmax/lnp0;/∗ SA algorithm ∗/S←S0, S∗←S0, T=T0, i=1, j=1;while stop criteria are not met dowhilei<Lenandj<LdoRandomly selectNS from NS1,NS2,NS3, and obtain a solution S′∈NSS;ifS′ is feasible thenApply the local search mechanism toS′;ifcostS′<costSthenS←S′;ifcostS<costS∗thenS∗←S;end ifelseSetS←S′ with probability p=exp−costS′−costS/Tend ifi=i+1;end ifj=j+1end whileT=0.9⋅T, i=1, j=1;end whilereturnS∗


## 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 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.


### 6.1. 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. These 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 mi 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. 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 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. Table3 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.Table 3
Dataset for the 2L-HFFVRP instances.
Inst.Vehicle typeABCD10232200333035413450222604570332803139018100542110153120614130161413341505541601071713135180571914882017146212612112216812230810112404913250101215261148152716181428291716290172430391722311173626323282026331192429341332933356224636733469To validate the mathematical formulation, 8 small-scale instances, namely, P3, P4, P5, P6, P7, P8, P9, and P10, are generated based on instance 1 for the 2L-HFFVRP. 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 2L-HFFVRP, 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 probabilityp0 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 large-scale 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.Table 4
The parameter settings of SA.
ParameterTested rangeSuggested valueProbabilityp0 of accepting any solution at the initial temperature0.05–0.200.10NumberLen of feasible moves to be performed at each temperature5000–20 00010 000


### 6.3. Results and Discussion

The 2L-HFFVRP mathematical formulation of Section3 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 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 Tables5–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.Table 5
Comparison of results of CPLEX and SA on small-scale HFFVRP of class 1.
InstanceNumber of customersCPLEXSA% gapCostTime (s)CostTime (s)P33130.250.27130.2529.650.00P44188.100.27188.1037.230.00P55222.350.50222.3533.550.00P66255.960.89255.9624.460.00P77308.362.75317.5126.972.97P88345.484.36345.4829.190.00P99394.9716.70394.9731.030.00P1010416.61587.25417.4013.900.19Average76.6228.250.39Table 6
Comparison of results of CPLEX and SA on small-scale unrestricted 2L-HFFVRP of classes 2–5.
InstanceClassNumber of customersNumber of itemsCPLEXSA% gapCostTime (s)CostTime (s)P3236130.250.43130.2517.420.00336130.250.32130.2532.170.00437130.250.34130.2530.950.005312130.250.55130.2538.950.00P4248201.550.67201.5520.560.00347188.100.54188.1027.060.004411201.550.52201.5534.360.005413201.550.46201.5539.670.00P52510222.350.49222.3521.370.00359237.520.54237.5212.980.004512222.350.66222.3528.920.005517222.351.71222.3532.520.00P62611257.310.93257.3112.530.003611257.311.17257.3124.440.004616257.312.17257.3123.630.005618257.313.37257.3116.130.00P72712308.362.84317.5127.032.973712308.363.69318.8112.443.394719308.368.65318.8119.923.395723308.3613.38317.5122.642.97P82813345.482.86345.4811.630.003815371.768.27371.7613.860.004822345.4815.68345.4831.640.005826345.4846.06345.4817.400.00P92915406.736.54406.739.280.003917397.4764.76397.4726.430.004924394.9775.75394.9711.270.005930394.97996.45397.466.980.63P1021017446.91223.98446.9111.130.0031029442.883600.04453.1311.742.3141025422.042444.98430.8615.112.0951033430.863600.09439.0012.251.89Average347.7821.080.61Table 7
Occupied rectangular area in the heuristic Min Occupied Rectangular area.
InstanceClassNumber of customersNumber of itemsCPLEXSA% gapCostTime (s)CostTime (s)P3236130.250.31130.259.460.00336130.250.38130.2515.960.00437130.250.44130.2521.460.005312130.250.84130.2534.860.00P4248201.550.54201.559.280.00347188.100.58188.1015.070.004411201.551.81201.5520.670.005413201.552.21201.5539.760.00P52510222.352.34222.359.560.00359237.521.97237.5222.050.004512222.355.53222.3522.770.005517222.3515.08222.3544.130.00P62611257.319.99257.3110.080.003611257.318.60257.3118.830.004616257.3128.48257.3121.250.005618257.3134.61257.3125.920.00P72712308.3620.75317.5114.332.973712308.3639.11318.8112.733.394719308.36137.70326.8217.815.995723308.36440.71317.5112.842.97P82813345.4844.89345.4811.400.003815371.76146.19371.7612.430.004822345.48653.90345.4815.170.005826345.482207.99345.4820.630.00P92915406.73473.62406.738.420.003917397.47986.25397.479.440.004924394.973600.24∗394.911.430.005930429.883600.26406.7311.30−5.39P1021017––∗446.910.46–31029453.433600.20453.139.81−0.0741025477.253600.33430.8614.39−9.7251033––430.8617.55–Average655.5317.230.004The results for HFFVRP small-scale problems of class 1 are presented in Tables5, 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 (Table5) 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.Tables5–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.The results provided by SA are very close to optimality. According to Table5, 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. 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 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. Therefore, the SA algorithm is effective in solving problems with larger numbers of customers and items.The large-scale problems described in Section5.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. 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 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.Table 8
Results of SA on large-scale instances.
Class 1Classes 2-5HFFVRPUnrestricted 2L-HFFVRPSequential 2L-HFFVRPInstanceCostTime (s)CostTime (s)CostTime (s)1683.7810.43715.7710.48719.5022.272716.7513.47758.2411.28760.1721.613800.2415.83877.2912.02881.6026.544728.5913.10791.6810.01797.6224.145782.5214.39931.0211.761009.9124.326890.9720.78968.1713.65968.2329.8774012.2227.666109.1721.026233.7657.4084105.1726.727267.2723.017269.3051.1991130.4114.431198.1311.951206.9524.24105824.7356.618454.7824.588456.4054.55116168.4856.199530.8429.819591.7061.07121806.8213.981822.9510.961824.9722.291321 551.3875.0230 529.2646.2931 009.94116.031410 804.3619.8111 529.6813.0911 544.3129.051510 529.4426.9412 541.3917.8712 816.6639.12161454.3112.411486.7410.351506.3122.80172166.9011.992232.449.922234.9821.06183464.5860.796350.5529.786652.6582.37192275.3932.504735.3424.374905.9660.60202387.6047.246578.1322.656793.6350.48214462.1137.2310 038.3126.0910 549.4375.65225678.0929.6610 997.6025.5811 487.4268.85235166.6735.0810 549.3424.8210 968.8550.65243803.3523.255422.3819.925769.5054.52255509.0943.1114 162.4728.2114 683.9191.52267658.3055.0815 482.9220.9816 300.1561.13273745.4526.716799.6521.616999.7953.24285053.2274.7425 484.5029.3327 079.76121.112911 004.28113.2327 967.0946.6528 522.72202.55308688.9244.0119 487.9329.1120 362.33116.05319799.8833.4323 957.0332.6824 363.94293.833211 043.7722.6124 775.0727.1725 988.46157.523311 065.0144.1825 983.1935.3326 493.88118.10347258.4919.3316 023.9831.7216 455.94143.01355367.1468.9011 148.2054.8811 596.82598.67363557.0326.266226.4543.116366.40316.60Average5309.6035.2010 275.4123.9510 588.1693.45In this sense, the impact of the three new packing heuristics, Heur6, Heur7, and Heur8, 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 6, Heur7, and Heur8 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 Heur7. 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 Heur7 stands out regarding managing more complex loadings in both unrestricted and sequential problems. For sequential problems, heuristic Heur6 manages to obtain new feasible loadings. In fact, heuristic Heur6 tends to form rectangular blocks of items, which favors sequential loading of blocks of items of the same customer. Although heuristic Heur8 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.Table 9
Percentage of the feasible loadings obtained by heuristicsHeur6, Heur7, and Heur8
HeuristicHeur6Heur7Heur8Unrestricted 2L-HFFVRP0.053.110.33Sequential 2L-HFFVRP1.281.960.45


## 6.1. 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. These 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 mi 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. 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 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. Table3 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.Table 3
Dataset for the 2L-HFFVRP instances.
Inst.Vehicle typeABCD10232200333035413450222604570332803139018100542110153120614130161413341505541601071713135180571914882017146212612112216812230810112404913250101215261148152716181428291716290172430391722311173626323282026331192429341332933356224636733469To validate the mathematical formulation, 8 small-scale instances, namely, P3, P4, P5, P6, P7, P8, P9, and P10, are generated based on instance 1 for the 2L-HFFVRP. 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 2L-HFFVRP, 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 probabilityp0 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 large-scale 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.Table 4
The parameter settings of SA.
ParameterTested rangeSuggested valueProbabilityp0 of accepting any solution at the initial temperature0.05–0.200.10NumberLen of feasible moves to be performed at each temperature5000–20 00010 000


## 6.3. Results and Discussion

The 2L-HFFVRP mathematical formulation of Section3 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 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 Tables5–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.Table 5
Comparison of results of CPLEX and SA on small-scale HFFVRP of class 1.
InstanceNumber of customersCPLEXSA% gapCostTime (s)CostTime (s)P33130.250.27130.2529.650.00P44188.100.27188.1037.230.00P55222.350.50222.3533.550.00P66255.960.89255.9624.460.00P77308.362.75317.5126.972.97P88345.484.36345.4829.190.00P99394.9716.70394.9731.030.00P1010416.61587.25417.4013.900.19Average76.6228.250.39Table 6
Comparison of results of CPLEX and SA on small-scale unrestricted 2L-HFFVRP of classes 2–5.
InstanceClassNumber of customersNumber of itemsCPLEXSA% gapCostTime (s)CostTime (s)P3236130.250.43130.2517.420.00336130.250.32130.2532.170.00437130.250.34130.2530.950.005312130.250.55130.2538.950.00P4248201.550.67201.5520.560.00347188.100.54188.1027.060.004411201.550.52201.5534.360.005413201.550.46201.5539.670.00P52510222.350.49222.3521.370.00359237.520.54237.5212.980.004512222.350.66222.3528.920.005517222.351.71222.3532.520.00P62611257.310.93257.3112.530.003611257.311.17257.3124.440.004616257.312.17257.3123.630.005618257.313.37257.3116.130.00P72712308.362.84317.5127.032.973712308.363.69318.8112.443.394719308.368.65318.8119.923.395723308.3613.38317.5122.642.97P82813345.482.86345.4811.630.003815371.768.27371.7613.860.004822345.4815.68345.4831.640.005826345.4846.06345.4817.400.00P92915406.736.54406.739.280.003917397.4764.76397.4726.430.004924394.9775.75394.9711.270.005930394.97996.45397.466.980.63P1021017446.91223.98446.9111.130.0031029442.883600.04453.1311.742.3141025422.042444.98430.8615.112.0951033430.863600.09439.0012.251.89Average347.7821.080.61Table 7
Occupied rectangular area in the heuristic Min Occupied Rectangular area.
InstanceClassNumber of customersNumber of itemsCPLEXSA% gapCostTime (s)CostTime (s)P3236130.250.31130.259.460.00336130.250.38130.2515.960.00437130.250.44130.2521.460.005312130.250.84130.2534.860.00P4248201.550.54201.559.280.00347188.100.58188.1015.070.004411201.551.81201.5520.670.005413201.552.21201.5539.760.00P52510222.352.34222.359.560.00359237.521.97237.5222.050.004512222.355.53222.3522.770.005517222.3515.08222.3544.130.00P62611257.319.99257.3110.080.003611257.318.60257.3118.830.004616257.3128.48257.3121.250.005618257.3134.61257.3125.920.00P72712308.3620.75317.5114.332.973712308.3639.11318.8112.733.394719308.36137.70326.8217.815.995723308.36440.71317.5112.842.97P82813345.4844.89345.4811.400.003815371.76146.19371.7612.430.004822345.48653.90345.4815.170.005826345.482207.99345.4820.630.00P92915406.73473.62406.738.420.003917397.47986.25397.479.440.004924394.973600.24∗394.911.430.005930429.883600.26406.7311.30−5.39P1021017––∗446.910.46–31029453.433600.20453.139.81−0.0741025477.253600.33430.8614.39−9.7251033––430.8617.55–Average655.5317.230.004The results for HFFVRP small-scale problems of class 1 are presented in Tables5, 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 (Table5) 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.Tables5–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.The results provided by SA are very close to optimality. According to Table5, 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. 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 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. Therefore, the SA algorithm is effective in solving problems with larger numbers of customers and items.The large-scale problems described in Section5.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. 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 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.Table 8
Results of SA on large-scale instances.
Class 1Classes 2-5HFFVRPUnrestricted 2L-HFFVRPSequential 2L-HFFVRPInstanceCostTime (s)CostTime (s)CostTime (s)1683.7810.43715.7710.48719.5022.272716.7513.47758.2411.28760.1721.613800.2415.83877.2912.02881.6026.544728.5913.10791.6810.01797.6224.145782.5214.39931.0211.761009.9124.326890.9720.78968.1713.65968.2329.8774012.2227.666109.1721.026233.7657.4084105.1726.727267.2723.017269.3051.1991130.4114.431198.1311.951206.9524.24105824.7356.618454.7824.588456.4054.55116168.4856.199530.8429.819591.7061.07121806.8213.981822.9510.961824.9722.291321 551.3875.0230 529.2646.2931 009.94116.031410 804.3619.8111 529.6813.0911 544.3129.051510 529.4426.9412 541.3917.8712 816.6639.12161454.3112.411486.7410.351506.3122.80172166.9011.992232.449.922234.9821.06183464.5860.796350.5529.786652.6582.37192275.3932.504735.3424.374905.9660.60202387.6047.246578.1322.656793.6350.48214462.1137.2310 038.3126.0910 549.4375.65225678.0929.6610 997.6025.5811 487.4268.85235166.6735.0810 549.3424.8210 968.8550.65243803.3523.255422.3819.925769.5054.52255509.0943.1114 162.4728.2114 683.9191.52267658.3055.0815 482.9220.9816 300.1561.13273745.4526.716799.6521.616999.7953.24285053.2274.7425 484.5029.3327 079.76121.112911 004.28113.2327 967.0946.6528 522.72202.55308688.9244.0119 487.9329.1120 362.33116.05319799.8833.4323 957.0332.6824 363.94293.833211 043.7722.6124 775.0727.1725 988.46157.523311 065.0144.1825 983.1935.3326 493.88118.10347258.4919.3316 023.9831.7216 455.94143.01355367.1468.9011 148.2054.8811 596.82598.67363557.0326.266226.4543.116366.40316.60Average5309.6035.2010 275.4123.9510 588.1693.45In this sense, the impact of the three new packing heuristics, Heur6, Heur7, and Heur8, 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 6, Heur7, and Heur8 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 Heur7. 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 Heur7 stands out regarding managing more complex loadings in both unrestricted and sequential problems. For sequential problems, heuristic Heur6 manages to obtain new feasible loadings. In fact, heuristic Heur6 tends to form rectangular blocks of items, which favors sequential loading of blocks of items of the same customer. Although heuristic Heur8 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.Table 9
Percentage of the feasible loadings obtained by heuristicsHeur6, Heur7, and Heur8
HeuristicHeur6Heur7Heur8Unrestricted 2L-HFFVRP0.053.110.33Sequential 2L-HFFVRP1.281.960.45


## 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 2L-HFFVRP. The results of this method can be effectively used in decision-making at the strategic and tactical levels, leading to better use of the fleet of vehicles and cost reduction.There is a trade-off between solution cost and computational time. Better solutions can be obtained by increasing the values of parametersp0 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. 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 time-consuming. 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 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. This 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. 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 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. 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 two-dimensional 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 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.



---

*Source: 6012105-2022-01-28.xml*