INTEGER FORMULATIONS FOR THE INTEGRATED VEHICLE ROUTING PROBLEM WITH TWO-DIMENSIONAL PACKING CONSTRAINTS

ABSTRACT

Two integer linear programming models are developed for the unrestricted vehicle routing problem with two-dimensional loading constraints. The first one is a complete model, and the other uses valid inequalities to guarantee that routes are connected and respect the two-dimensional loading constraints. The models are solved with a branch-and-cut algorithm. Computational experiments on benchmark instances showed the complete model has allowed optimal solutions for 5% of the instances, while the second model optimally solved 64% of the instances. Given the superior performance of the second model, we adapted it to handle the sequential variant of the problem, which is harder, and then optimal solutions were obtained for 46% of the instances within the given time limit. The second model compared with a branch-and-cut algorithm from the literature found identical or better solutions for all the instances.

Keywords:
vehicle routing problem with two-dimensional loading; multi-drop requirements; integer programming formulation; branch-and-cut

1 INTRODUCTION

This paper addresses the Vehicle Routing Problem with Two-Dimensional Loading Constraints (2L-CVRP), an NP-Hard problem (Iori et al., 2007). It integrates two other NP-hard optimization problems: Capacitated Vehicle Routing Problem (CVRP) (Dantzig & Ramser, 1959; Clarke & Wright, 1964; Mor & Speranza, 2020) and Two-Dimensional Bin Packing Problem (2BPP) (Kantorovich, 1960; Sweeney & Paternoster, 1992; Côté et al., 2021). The two-dimensional loading in the CVRP emerges when it is impossible to stack items due to fragility or size. Moreover, including loading constraints in CVRPs will yield more accurate route planning since this problem appears in many practical situations, especially in logistics. Côté et al. (2017) pointed out that solving the 2L-CVRP in a non-integrated way may increase the cost of solutions by 7% on average.

In the 2L-CVRP, we are interested in minimum-cost routes to satisfy customers’ demands for two-dimensional rectangular items. Vehicles depart from a central depot, visit one or more customers to deliver items, and return to the depot. Besides respecting the vehicles’ maximum weight and area capacities, we need to give attention to the positioning of items on the vehicles’ rectangular base. These loading subproblems require items not to overlap each other, and they must be entirely inside the base. In some cases, it is required that items must have free passage when unloading at customers. The latter refers to multi-drop requirements or the sequential variant, where re-handling items is not allowed when unloading; otherwise, we have the unrestricted variant, where items can be re-handled during unloading operations (Iori & Martello, 2010).

In this paper, two integer models are developed for the 2L-CVRP and solved with a Branch- and-Cut algorithm (B&C). Differently from the literature (e.g., Iori et al. (2007), Azevedo et al. (2009), and Hokama et al. (2016)), the first model considers the entire problem, meaning we deal simultaneously with routing and packing decisions in the model. This first model has three- index routing variables and five-index packing variables. On the other hand, the second model follows the literature and has two-index routing variables, where valid inequalities are added as demanded after running separation routines to deal with the loading subproblems.

Our contributions to the 2L-CVRP rely on: (i) proposal of a complete model; (ii) proposal of a tailored separation routine for detecting violated packing decisions; (iii) comparison between two methodologies (to insert all constraints at once and to detect/insert valid inequalities during optimization); and, (iv) computational results for a wide range of experiments, comprising both the unrestricted and sequential variants of the 2L-CVRP.

In the next section, a literature review related to the 2L-CVRP is presented. Some variants of the 2L-CVRP consider constraints, such as time windows, heterogeneous fleets, and multiple depots. Concerning the loading subproblem, additional constraints are related to cargo stability, fragility conditions, load balance, and grouping of items. Section 3 has the developed models and separation routines. Computational experiments and results are detailed in Section 4. Finally, concluding remarks and proposals for future works are given in Section 5.

2 LITERATURE RELATED TO THE 2L-CVRP

The 2L-CVRP was introduced by Iori et al. (2007). The authors described an integer linear programming model to solve the sequential variant of the 2L-CVRP. This model was solved with a B&C algorithm in which cutting planes related to sub-tours were added during the optimization process. The loading subproblem associated with each route was solved with a Branch-and-Bound algorithm (B&B) called Check-2L, where lower and upper bounds were also applied to prevent unnecessary enumeration. Computational experiments were carried out on instances adapted from Reinelt (1991) and Toth & Vigo (2001), totaling 60 instances with up to 35 customers and 114 items.

In the literature, there are few exact algorithms to deal with the 2L-CVRP and its variants. Azevedo et al. (2009) proposed a B&C to the 2L-CVRP. The authors developed new separation routines and valid inequalities based on those of the CVRP. An exact algorithm, combined with simple heuristics, was used to solve the loading subproblems. In the sequential variant, they compared their results with those of Iori et al. (2007). They were able to find optimal solutions to more instances, besides reducing the runtime required to solve many others.

Hokama et al. (2016) presented a B&C to solve the sequential 2L-CVRP, where valid inequalities were proposed following the literature in the CVRP. The authors proposed a constraint programming formulation to handle the loading subproblems. When comparing their results with those from the literature, they improved in 9% the number of optimal solutions obtained.

Côté et al. (2017) compared solutions of the 2L-CVRP with those obtained from solving each problem (i.e., the CVRP and 2BPP) independently and then combining their solutions to obtain a feasible 2L-CVRP solution. The computational experiments on benchmark instances using integer programming models showed that an improvement of approximately 7% can be achieved if the 2L-CVRP is solved directly. They also pointed that the loading subproblem may be the bottleneck of this problem.

Junqueira et al. (2013) solved the problem with three-dimensional items, presenting a complete model that was solved with a B&C algorithm. Their model combined a time-dependent formulation for the routing counterpart and a grid-based formulation for the loading subproblem. The model has four-index variables related to the routing and seven-index variables related to the loading and allowed optimal solutions to small instances within some hours of running time. The authors also studied the impact of some practical constraints in the loading subproblem, such as cargo stability, multi-drop requirements, and load-bearing strength of items. In this work, we propose a complete model for the 2L-CVRP that has fewer variables and constraints than Junqueira et al. (2013). Recently, Ferreira et al. (2021) proposed integer programming models and B&C algorithms for three variants of the 2L-CVRP: allowing split delivery, with green requirements, and combing split delivery and green requirements.

Regarding heuristics for the 2L-CVRP and its variants, Gendreau et al. (2008) proposed a tabu search-based approach, while Fuellerer et al. (2009) developed an ant colony optimization. In Zachariadis et al. (2009), there is a guided tabu search heuristic with a structure to save evaluated routes and avoid unnecessary computation. In Strodl et al. (2010), a variable neighborhood search handled routing decisions while an exact procedure dealt with loading subproblems. On the other hand, Leung et al. (2010) developed a simulated annealing able to improve many solutions reported by Gendreau et al. (2008), Zachariadis et al. (2009), and Fuellerer et al. (2009).

In Duhamel et al. (2011), the greedy randomized-adaptive search procedure (GRASP) was combined with a local evolutionary search for the 2L-CVRP. On the other hand, Leung et al. (2011) proposed an extended guided tabu search. Shen & Murata (2012) solved the 2L-CVRP with a genetic algorithm, where a heuristic based on the bottom-left was used to handle loading subproblems. In Bin et al. (2013), the 2L-CVRP was solved by an artificial bee colony algorithm combined with three heuristics for handling loading subproblems. These authors also allowed items to rotate orthogonally. Computational tests on instances with up to 100 customers and 310 items showed that such a heuristic performed better than those of Gendreau et al. (2008), Zachariadis et al. (2009), and Leung et al. (2011).

In Zachariadis et al. (2013) there is a local search method combined with a practical heuristic for loading subproblems, while Domínguez et al. (2014) proposed a multi-start biased-randomized algorithm for the 2L-CVRP in which items could orthogonally rotate. In Wei et al. (2015), there is an efficient variable neighborhood search. The loading subproblems are solved with a skyline- based heuristic, where different indicators are used to determine which items to arrange in the available positions. To save computational time, routes already evaluated are kept in a pool. Wei et al. (2018) combined simulated annealing with the open space technique to efficiently solve the 2L-CVRP. While simulated annealing drives the optimization towards good-quality solutions, the open space technique deals with loading subproblems. Recently, Ferreira & Queiroz (2022) solved a generalization of the 2L-CVRP that includes decisions related to the location of depots. The authors presented a hybrid heuristic that integrates simulated annealing and the artificial algae algorithm. The heuristic also uses the variable neighborhood descent in the local search phase and the skyline technique to obtain feasible packings.

Table 1 summarizes the literature related to the 2L-CVRP and its variants. Each row in the table presents a work, its group, i.e., exact or heuristic, and if the 2L-CVRP is the main problem; otherwise, the problem is a variant. Other features are also described, such as the variant investigated in the loading subproblem (unrestricted and/or sequential); if the problem is solved assuming that the number of routes is exactly equal to (=) or less than (≤) the number of available vehicles (|K|), or it is not imposed (≥ K); if it is allowed to have routes serving only one customer; which type of cost ci j is assumed to the edges {i, j} (integer or real); and, finally, the main method that is used to solve the problem (B&C, constraint programming, tabu search, simulated annealing, genetic algorithm, local search, among others). The entry ‘-’ indicates that no information was given by the authors.

Table 1
Literature related to the 2L-CVRP and its variants.

Group	Authors	Problem	Loading subproblem	Number of routes	Single routes	Cost of edges	Method
Exact	Iori et al. (2007)	2L-CVRP	Sequential	= |K|	No	Integer	B&C
	Azevedo et al. (2009)	2L-CVRP	Sequential and Unrestricted	= |K|	No	Integer	B&C
	Hokama et al. (2016)	2L-CVRP	Sequential	= |K|	No	Integer	B&C
	Côté et al. (2017)	2L-CVRP	Sequential	≥ |K|	Yes	Integer	B&C
Heuristics	Gendreau et al. (2008)	2L-CVRP	Sequential	= |K|	No Yes	Integer Real	Tabu Search
			Sequential and Unrestricted	≤ |K|
	Fuellerer et al. (2009)	2L-CVRP	Sequential	= |K|	No Yes	Integer Real	Ant Colony Algorithm
			Sequential and Unrestricted	≤ |K|
	Zachariadis et al. (2009)	2L-CVRP	Sequential and Unrestricted	≤ |K|	Yes	Real	Guided Tabu Search
	Strodl et al. (2010)	2L-CVRP	Unrestricted	≤ |K|	Yes	Real	Variable Neighborhood Search
	Leung et al. (2010)	2L-CVRP	Sequential and Unrestricted	≤ |K|	Yes	Real	Simulated Annealing
	Duhamel et al. (2011)	2L-CVRP	Unrestricted	≤ |K|	Yes	Real	GRASP
	Leung et al. (2011)	2L-CVRP	Sequential and Unrestricted	≤ |K|	Yes	Real	Extended Guided Tabu Search
	Shen & Murata (2012)	2L-CVRP	Sequential	≤ |K|	-	Real	Genetic algorithm
	Bin et al. (2013)	2L-CVRP	Sequential and Unrestricted	≤ |K|	Yes	Real	Artificial Bee Colony
	Zachariadis et al. (2013)	2L-CVRP	Sequential and Unrestricted	≤ |K|	Yes	Real	Local Search
	Domínguez et al. (2014)	2L-CVRP	Unrestricted	≤ |K|	Yes	Real	Multistart Biased-Randomized Algorithm
	Wei et al. (2015)	2L-CVRP	Sequential and Unrestricted	≤ |K|	Yes	Real	Variable Neighborhood Search
	Wei et al. (2018)	2L-CVRP	Sequential and Unrestricted	≤ |K|	Yes	Real	Simulated Annealing
Variants	Malapert et al. (2008)	2L-PDP	Sequential and Unrestricted	-	-	-	Constraint Programming
	Munoz (2011)	2L-CVRP-H	Sequential and Unrestricted	= |K|	Yes	Real	Mixed-integer Programming
	Khebbache-Hadji et al. (2013)	2L-CVRPTW	Unrestricted	≤ |K|	No	Real	Memetic Algorithm
	Martínez & Amaya (2013)	VRPM-TW-CL	Unrestricted	= |K|	No	Real	Tabu Search
	Domínguez et al. (2016)	2L-HFVRP	Unrestricted	≤ |K|	Yes	Real	Multistart Biased-Randomized Algorithm
	Zachariadis et al. (2016)	2L-SPD	Sequential and Unrestricted	≤ |K|	-	Real	Local Search
	Côté et al. (2020)	S2L-CVRP	Sequential	= |K|	No	Integer	B&C
	Ferreira et al. (2021)	2L-SDVRP	Sequential	≤ |K|	Yes	Integer	B&C
		G2L-CVRP
		G2L-SDVRP

• 2L-PDP: Two-dimensional Pickup and Delivery Problem
• 2L-CVRP-H: Capacitated Vehicle Routing Problem with Two-dimensional Loading Constraints and Handling Costs
• S2L-CVRP: Vehicle Routing with Stochastic Two-Dimensional Items
• 2L-CVRPTW: Two-Dimensional Loading Capacitated Vehicle Routing Problem with Time Windows
• VRPM-TW-CL: Vehicle Routing Problem with Multi-Trips, Time Windows and Two Dimensional Circular Loading Constraints
• 2L-HFVRP: Two-Dimensional Loading Heterogeneous Fleet Vehicle Routing Problem
• 2L-SPD: Two-dimensional Simultaneous Pickup and Delivery Problem
• 2L-SDVRP: Capacitated Vehicle Routing Problem with Two-Dimensional Loading Constraints and Split Delivery
• G2L-CVRP: Green Vehicle Routing Problem with Two-Dimensional Loading Constraints
• G2L-SDVRP: Green Vehicle Routing Problem with Two-Dimensional Loading Constraints and Split Delivery

Table 1 shows that less than 35% of the works proposed integer linear programming models and exact methods (mainly B&C) to solve the 2L-CVRP and variants. Other works proposed heuristics to tackle large-sized instances, as well as to consider additional constraints in the loading subproblem, as allowing items to rotate (Bin et al., 2013), circular items (Martínez & Amaya, 2013), re-handling of items (Munoz, 2011), stochastic items (Côté et al., 2020), among others. Practical constraints are also included in the routing problem, such as considering pickup and delivery (Malapert et al., 2008; Zachariadis et al., 2016), a heterogeneous fleet of vehicles (Domínguez et al., 2016), time windows to visiting customers (Martínez & Amaya, 2013; Khebbache-Hadji et al., 2013), each vehicle may perform more than one route (Martínez & Amaya, 2013), splitting the delivery (Ferreira et al., 2021), and green requirements (Ferreira et al., 2021).

Concerning the sequential and unrestricted variants, approximately half of the works in Table 1 have tackled both, while less than 25% of the works focused only on the unrestricted variant. Besides that, many works (i.e., more than 65%) allowed the number of routes to be equal to or less than the number of vehicles and vehicles that can perform single routes. On the other hand, less than 35% of the works presented results of the case in which edges have integer costs, most of them when a B&C is used. It is worth mentioning that we could not find any work making a proper comparison of models for the 2L-CVRP, thus pointing their advantages and weaknesses (for example, in terms of inserting all constraints at once or cutting planes as demanded).

3 INTEGER PROGRAMMING MODELS FOR THE 2L-CVRP

The 2L-CVRP is defined on a directed graph G = (V, E), where V is the set with n + 1 nodes, corresponding to the central depot (j = 0) and customers (j = 1, 2, . . . , n), and E is the set of arcs, in which $E=\left\{\right(i,j):i,j\in V,i\ne j\}$. Each arc $(i,j)\in E$ has a non-negative cost c _ij to traverse it. Let K be the set of identical vehicles, each vehicle of capacity P and available at the depot. The dimensions of the vehicle’s rectangular base are W (width, along the x-axis) and H (height, along the y-axis), with total area A = WH, defined in the Cartesian plane with the origin at (0, 0). Each customer j demands a set M _j of rectangular items, where $M={\cup }_{j\in V\backslash \left\{0\right\}}{M}_j$. Each item $m\in M_j$ has width w _jm , height h _jm , and weight p _jm . The total area of all items belonging to the customer j is $a_j=\sum _{m\in M_j}{w}_{jm}{h}_{jm}$ and the total weight is $p_j=\sum _{m\in M_j}{p}_{jm}$. Without loss of generality, all dimensions are assumed to be positive integer numbers. The problem has the following constraints:

• The number of routes can be less than or equal to the number of vehicles;

• Each route starts and ends at the central depot and is associated with exactly one vehicle;

• Routes serving a single customer are allowed;

• Each customer must be served by exactly one vehicle, and all its items must be delivered at once (i.e., splitting deliveries is not allowed);

• Vehicle capacities must be respected, both in terms of weight and area;

• Each route produces a loading subproblem. In this subproblem, all items belonging to the customers of such a route must be orthogonally packed on the vehicle base without overlapping and fully inside the base. Moreover, items cannot be rotated (i.e., fixed orientation) and are arranged by their bottom-left corners. When items do not always have free passage during the unloading operation at a customer, the problem is called unrestricted, meaning items can be re-handled. Otherwise, when each item has a free passage in the unloading operation (re-handling items is not allowed), the problem is called sequential.

The objective of the 2L-CVRP is to obtain routes with the minimum total cost to serve all customers’ demands, respecting the above-mentioned constraints. We develop in the next subsections two integer linear programming models to solve the unrestricted 2L-CVRP. The first model has three-index variables related to routing decisions and five-index variables related to loading subproblems. It is a complete model in which all constraints are inserted at once. The second model has only two-index variables, such that valid inequalities concerning the elimination of sub-tours and infeasible loadings in each route are detected and added during the optimization process.

3.1 Complete model

The complete model can be seen as an improvement in the integer programming model of Junqueira et al. (2013) originally proposed to the 3L-CVRP. For the loading subproblem, without loss of generality, the vehicles’ rectangular base is discretized on the grid of canonical dissections proposed by Herz (1972). Items are arranged on this grid. Sets 𝒫 ^W in (1) and 𝒫 ^H in have the canonical dissection coordinates along the x- (width direction) and y-axis (height direction), respectively. Notice that variable ρ _jm is used to indicate whether item m of customer j appears in the binary linear combination $\sum _{j\in V\backslash \left\{0\right\}}\sum _{m\in M_j}{\rho }_{jm}{w}_{jm}$ to obtain the coordinate s (similarly for the coordinate t).

${\mathcal{P}}^W=\{s\in Z^{+}|s=\sum _{j\in V\backslash \left\{0\right\}}\sum _{m\in M_j}{\rho }_{jm}{w}_{jm},0\le s\le W,\text{for}{\rho }_{jm}\in \{0,1\},j\in V\backslash \left\{0\right\},m\in M_j\}$ (1)

${\mathcal{P}}^H=\{t\in Z^{+}|t=\sum _{j\in V\backslash \left\{0\right\}}\sum _{m\in M_j}{\rho }_{jm}{h}_{jm},0\le t\le H,\text{for}{\rho }_{jm}\in \{0,1\},j\in V\backslash \left\{0\right\},m\in M_j\}$ (2)

From sets (1) and (2), the respective sets ${\mathcal{P}}_{jm}^W$ and ${\mathcal{P}}_{jm}^H$ of valid coordinates for each item m of customer j are obtained. Therefore, ${\mathcal{P}}_{jm}^W$ and ${\mathcal{P}}_{jm}^H$ are the respective sets 𝒫 ^w and 𝒫 ^H , where coordinates after W − w _jm and H − h _jm are disregarded because they are no longer valid to arrange item m in the vehicle rectangular base.

The complete model uses the binary variable x _ijk , which is equal to 1 if arc $(i,j)\in E$ has vehicle $k\in K$ traversing it; otherwise, x _ijk is equal to zero. Another variable is y _j , which indicates the vehicle’s load weight at customer j. The model also has the binary variable $z_{jk}^{mst}$, which is 1 if item m of customer j is on vehicle k, having its lower-left corner positioned in the point (s,t), where $s\in {\mathcal{P}}_{jm}^W$ and $t\in {\mathcal{P}}_{jm}^H$; otherwise $z_{jk}^{mst}$ is zero. The complete model has objective function (3) and constraints (4)-(8), which are related to the routing decisions, and constraints (9)-(10), concerning the loading subproblems.

$\text{Minimize}\sum _{(i,j)\in E}\sum _{k\in K}{c}_{ij}{x}_{ijk}$ (3)

Subject to:

$\sum _{k\in K}\sum _{\{j\in V|(i,j)\in E\}}{x}_{ijk}=1,\forall i\in V\backslash \left\{0\right\}$ (4)

$\sum _{\{i\in V|(i,h)\in E\}}{x}_{ihk}=\sum _{\{j\in V|(h,j)\in E\}}{x}_{hjk},\forall h\in V\backslash \left\{0\right\},\forall k\in K$ (5)

$\sum _{j\in V\backslash \left\{0\right\}}{x}_{0jk}\le 1,\forall k\in K$ (6)

$\sum _{j\in V\backslash \left\{0\right\}}{x}_{0jk}=\sum _{i\in V\backslash \left\{0\right\}}{x}_{i0k},\forall k\in K$ (7)

$y_i-y_j+P{x}_{ijk}+(P-p_i-p_j)x_{jik}\le P-p_j,\forall i,j\in V\backslash \left\{0\right\},i\ne j,\forall k\in K$ (8)

$\sum _{s\in {\mathcal{P}}_{jm}^W}\sum _{t\in {\mathcal{P}}_{jm}^H}{z}_{jk}^{mst}=\sum _{\{i\in V|(i,j)\in E\}}{x}_{ijk},\forall j\in V\backslash \left\{0\right\},\forall m\in M_j,\forall k\in K$ (9)

$$\begin{gathered} \sum _{j\in V\backslash \left\{0\right\}}\sum _{m\in M_j}\sum _{\{s\in {\mathcal{P}}_{jm}^W|u-w_{jm}+1\le s\le u\}}\sum _{\{t\in {\mathcal{P}}_{jm}^H|v-h_{jm}+1\le t\le v\}}{z}_{jk}^{mst}\le 1, \\ \forall u\in {\mathcal{P}}^W,\forall v\in {\mathcal{P}}^H,\forall k\in K \end{gathered}$$ (10)

$p_j\le y_j\le P,\forall j\in V\backslash \left\{0\right\}$ (11)

$x_{ijk}\in \{0,1\},\forall (i,j)\in E,\forall k\in K$ (12)

$z_{jk}^{mst}\in \{0,1\},\forall j\in V\backslash \left\{0\right\},\forall m\in M_j,\forall s\in {\mathcal{P}}_{jm}^W,\forall t\in {\mathcal{P}}_{jm}^H,\forall k\in K$ (13)

The objective function (3) is related to a solution (i.e., set of routes) of minimum total cost. Constraints (4) ensure that each customer i is served by exactly one vehicle k. Constraints (5)-(7) are related to the flow conservation. In (5), each customer h has the same incoming and outgoing flows. In (6), the number of routes may be less than the total number of vehicles. In constraints (7), a vehicle must start and end at the central depot if it is used.

Constraints (8) are used to eliminate sub-tours and were proposed by Kara (2010). These constraints ensure that the vehicle’s load weight y _i when at customer i is less than y _j if vehicle k is going from i to j (i.e., x _ijk = 1). Constraints (9) link the routing decisions (i.e., variables x _ijk ) with the loading subproblems (i.e., variables $z_{jk}^{mst}$), ensuring that item m of customer j is arranged on vehicle k if k visits j from any other customer i. Therefore, constraints (10) ensure the non-overlapping between items, imposing each point of the grid, for each vehicle k, on being covered by at most one item m of a customer j that is served by k. Finally, the variables’ domain is expressed in (11)-(13), with continuous and binary variables.

We also considered constraints to remove the symmetry in the routes by forcing the first route to have the largest weight (area), the second route to have the second-largest weight (area), and so on. However, preliminary experiments with the complete model having these constraints showed an important reduction in the number of instances with a solution and the gap increased. Therefore, we decided not to use them.

3.2 Two-index model

There are different models for the CVRP in which two-index variables are used, indicating whether a vehicle traverses an arc without explicitly specifying which vehicle is performing this. Our model is based on Iori et al. (2007) and Côté et al. (2017). Let ${\overline{x}}_{ij}$ be an integer variable that is equal to 1 if a vehicle traverses arc $(i,j)\in E$; otherwise, it is equal to zero. Let 𝒮 be the set of all infeasible routes due to the loading subproblem, and let r(S) be the smallest number of vehicles needed to serve all customers in S. The exact value of r(S) can be obtained by solving an instance of the 2BPP. The two-index model for the 2L-CVRP has the objective function (14) and constraints (15)-(21).

$\text{Minimize}\sum _{(i,j)\in E}{c}_{ij}{\overline{x}}_{ij}$ (14)

Subject to:

$\sum _{j\in V\backslash \left\{0\right\}}{\overline{x}}_{0j}\le \left|K\right|$ (15)

$\sum _{j\in V\backslash \left\{0\right\}}{\overline{x}}_{0j}=\sum _{i\in V\backslash \left\{0\right\}}{\overline{x}}_{i0}$ (16)

$\sum _{\{i\in V|(i,h)\in E\}}{\overline{x}}_{ih}=1,\forall h\in V\backslash \left\{0\right\}$ (17)

$\sum _{\{i\in V|(i,h)\in E\}}{\overline{x}}_{ih}=\sum _{\{j\in V|(h,j)\in E\}}{\overline{x}}_{hj},\forall h\in V\backslash \left\{0\right\}$ (18)

$\sum _{i\in S}\sum _{\{j\in S|(i,j)\in E\}}{\overline{x}}_{ij}\le \left|S\right|-r\left(S\right),\forall S\subseteq V\backslash \left\{0\right\},S\ne \varnothing$ (19)

$\sum _{i\in S}\sum _{\{j\in S|(i,j)\in E\}}{\overline{x}}_{ij}\le \left|S\right|-1,\forall S\in \mathcal{S}$ (20)

${\overline{x}}_{ij}\in \{0,1\},\forall (i,j)\in E$ (21)

The objective function (14) minimizes the sum of the routes cost. Constraint (15) ensures that the number of routes may be smaller than the number of available vehicles, while constraints (16) guarantee the same incoming and ongoing flows to the depot. Constraints (17) impose each customer on being visited exactly once, while constraints (18) are related to the continuity of routes.

Constraints (19) ensure that routes are valid (i.e., there is no sub-tour and respect the vehicle’s weight and area capacities). Constraints (20), which are related to the loading subproblems, impose routes on having feasible packings. Both types of constraints are detected and added as cutting planes during the optimization process. Finally, the domain of variables ${\overline{x}}_{ij}$ is given in (21).

The exact value of r(S) may be obtained from solving an instance of the 2BPP, taking into consideration all items of the customers in the route S. Hence, following Lysgaard et al. (2004), r(S) is changed to r'(S) as in (22), which returns the minimum number of vehicles required to service the customers’ demands in terms of weight and area. According to Naddef & Rinaldi (2002), constraints (19) with r'(S) are called the rounded capacity inequalities. From now on, constraints (19) are redefined with r'(S) instead of r(S).

$r\text{'}\left(S\right)=\max \left\{⌈\frac{{\sum }_{i\in S}{p}_i}{P}⌉,⌈\frac{{\sum }_{i\in S}{a}_i}{A}⌉\right\}$ (22)

Because of constraints (19) and (20), separation routines are called during the optimization to detect their violation. These routines can be applied to fractional or integer solutions to eliminate infeasible solutions. Therefore, the B&C starts solving the two-index model without these constraints since there is an exponential number of them.

Rounded capacity inequalities

For separating constraints (19), where r(S) is replaced by r'(S), the CVRPSEP package of Lysgaard et al. (2004) is used, resulting in Procedure 1 and Procedure 2 below. These procedures are applied whenever a fractional or integer solution is found in the B&C tree. In the case of integer solutions, after applying the procedures, we check whether the solution is feasible after all.

Procedure 1: Let G ^* = (V, E ^* ) be a support graph that is built from the optimal (linear programming) solution ${\overline{x}}^{*}$ of the current node in the B&C tree. This graph has the original set of nodes V, while its set of edges E ^* has only the arcs $(i,j)\in E$ for which ${\overline{x}}_{ij}^{*}>0$. The cost of each arch $(i,j)\in E^{*}$ is precisely the value that the variable ${\overline{x}}_{ij}^{*}$ has. The weight p _j of each customer $j\in V$ is used as the demand of j to find sets S that violate the rounded capacity inequalities. Therefore, the CVRPSEP package’s routine to separate rounded capacity inequalities is called, and the found cuts are added as lazy constraints.

Procedure 2: It is similar to Procedure 1, but now it uses the area aj of each customer $j\in V$ as the demand of j to find sets S.

We also tested the other separation algorithms in the CVRPSEP package for the framed capacity, strengthened comb, multistar, and hypotour inequalities (Lysgaard et al., 2004). However, we have found it more effective to prioritize and use the computing time they could require to separate the loading feasibility inequalities.

Loading feasibility inequalities

The separation of constraints (20), related to loading subproblems, occurs when an integer solution ${\overline{x}}^{+}$ is found in the B&C tree. For each route ℛ, let its set of customers be S _ℛ. A Two- dimensional Orthogonal Packing Problem (2OPP) is solved for each S _ℛ. In the 2OPP, given a set of items and a rectangular bin, the objective is to find a feasible packing of all items. These items are orthogonally packed, without overlapping, and fully inside the bin. If no such packing exists, this set of items is said to be infeasible. As the 2OPP is a decision problem, any feasible packing is valid as a solution.

The 2OPP, for a set of customers S _ℛ, is solved with the integer programming model (23)-(25). In this model, let y _mpq be a binary variable that is 1 if item m of a customer in S _ℛ is packed in the point $(p,q)\in {\mathcal{P}}^W\times {\mathcal{P}}^H$; otherwise, it is zero. Constraints (23) ensure there is no overlapping between items, similarly to constraints (10). Constraints (24) guarantee that all items, from the customers in S _ℛ, are packed. The domain of variables y _mpq is given in (25).

$\sum _{j\in S_{\mathcal{R}}}\sum _{m\in M_j}\sum _{\{s\in {\mathcal{P}}_{jm}^W:u-w_{jm}+1\le s\le u\}}\sum _{\{t\in {\mathcal{P}}_{jm}^H:v-h_{jm}+1\le t\le v\}}{y}_{mst}\le 1,\forall u\in {\mathcal{P}}^W,\forall v\in {\mathcal{P}}^H$ (23)

$\sum _{s\in {\mathcal{P}}_{jm}^W}\sum _{t\in {\mathcal{P}}_{jm}^H}{y}_{mst}=1,\forall j\in S_{\mathcal{R}},\forall m\in M_j$ (24)

$y_{mst}\in \{0,1\},\forall j\in S_{\mathcal{R}},\forall m\in M_j,\forall s\in {\mathcal{P}}_{jm}^W,\forall t\in {\mathcal{P}}_{jm}^H$ (25)

As solving the integer programming model (23)-(25) may be slow, a relaxation of the 2OPP is firstly solved. For that, such an integer model is separated into two other models, relaxing the 2OPP on a problem related to the width direction (2OPP^W ) and another related to the height direction (2OPP^H ). The idea is to quickly check whether the set of items in S _ℛ is infeasible by substituting the non-overlapping constraints (23) for the relaxed ones. Therefore, the integer programming model (26)-(28), for the 2OPP^W , considers the binary variable ${\overline{y}}_{ms}$, which is 1 if item m of customer $j\in S_{\mathcal{R}}$ is packed in the coordinate $s\in {\mathcal{P}}_{jm}^W$ on the width direction; otherwise ${\overline{y}}_{ms}$ is zero.

$\sum _{j\in S_{\mathcal{R}}}\sum _{m\in M_j}\sum _{\{s\in {\mathcal{P}}_{jm}^W|u-w_{jm}+1\le s\le u\}}{h}_{jm}{\overline{y}}_{ms}\le H,u\in {\mathcal{P}}^W$ (26)

$\sum _{s\in {\mathcal{P}}_{jm}^W}{\overline{y}}_{ms}=1,\forall j\in S_{\mathcal{R}},\forall m\in M_j$ (27)

${\overline{y}}_{ms}\in \{0,1\},\forall j\in S_{\mathcal{R}},\forall m\in M_j,\forall s\in {\mathcal{P}}_{jm}^W$ (28)

Concerning the model for the 2OPP^W , constraints (26) ensure that the sum of the items’ height, for those items that cover coordinate u when they are packed in coordinate s on the width direction, must be less than or equal to H. Constraints (27) impose all items of customers $j\in S_{\mathcal{R}}$ on being packed in some coordinate s on the width direction. The domain of variables ${\overline{y}}_{ms}$ is given in (28). Similarly, the integer programming model (29)-(31) is presented for the 2OPP^H , where the binary variable ${\tilde{y}}_{mt}$ is 1 if item m of a customer $j\in S_{\mathcal{R}}$ is packed in the coordinate $t\in {\mathcal{P}}_{jm}^H$ on the height direction; otherwise, ${\tilde{y}}_{mt}$ is zero.

$\sum _{j\in S_{\mathcal{R}}}\sum _{m\in M_j}\sum _{\{t\in {\mathcal{P}}_{jm}^H|v-h_{jm}+1\le t\le v\}}{w}_{jm}{\tilde{y}}_{mt}\le W,v\in {\mathcal{P}}^H$ (29)

$\sum _{t\in {\mathcal{P}}_{jm}^H}{\tilde{y}}_{mt}=1,\forall j\in S_{\mathcal{R}},\forall m\in M_j$ (30)

${\tilde{y}}_{mt}\in \{0,1\},\forall j\in S_{\mathcal{R}},\forall m\in M_j,\forall t\in {\mathcal{P}}_{jm}^H$ (31)

In the integer programming model (29)-(31), for the 2OPP^H , constraints (29) guarantee that the sum of the items’ width must be less than or equal to W, for those items packed in the coordinate t and covering the coordinate v on the height direction. Constraints (30) impose all items of customers $j\in S_{\mathcal{R}}$ on being packed in some coordinate t on the height direction. The domain of variables ${\tilde{y}}_{mt}$ is expressed in (31). Therefore, based on the models discussed in this subsection, the procedure to separate the loading feasibility constraints, which concern the loading subproblems, is described in Procedure 3. A hash table is used to save and search for routes S _ℛ that have already been checked due to their loading subproblem, aiming to reduce computing time.

Procedure 3: Whenever an integer solution ${\overline{x}}^{+}$ is found in the B&C tree and Procedures 1 and 2 do not succeed in generating a valid inequality, we obtain the set of routes ℛ from ${\overline{x}}^{+}$. For each route ℛ, we initially check if it is in the hash table. If yes, we retrieve the saved solution; otherwise, (i) we solve the 2OPP^W , model (26)-(28). If the resulting solution is feasible, (ii) we solve the 2OPP^H , model (29)-(31). If the resulting solution is feasible, (iii) we solve the 2OPP, model (23)-(25). If the resulting solution is feasible, steps (i)-(iii) are repeated for the next route in ${\overline{x}}^{+}$ until checking all routes. If the saved/resulting solution of any of these steps is infeasible, the inequality in (20) is added for S _ℛ , as a lazy constraint, and the procedure stops. Every checked route is added to the hash table.

3.3 Sequential 2L-CVRP

We conducted preliminary computational experiments (in Section 4.1) to analyze the complete and two-index models. The results showed the two-index model performed better than the complete model. Consequently, we decided to handle the sequential 2L-CVRP only with the two- index model because this variant generally imposes more computing time on solving the loading subproblems (Côté et al., 2014). According to Junqueira et al. (2013), a complete model for the sequential 2L-CVRP needs variables and constraints that indicate the sequence on which customers are visited on each route.

We assume that items are unloaded parallel to the y-axis (along the height direction, from bottom to top) in the sequential variant. For each route ℛ whose set of customers is S _ℛ , let o _j be the order to visit the j-th customer in ℛ (for ascending order). Hence, if o _i < o _j , customer i is visited before j, meaning that all items M _i of i must be arranged so that they have free passage during the unloading operation at i. It implies that no item in M _j of customer j must block any item in M _i of i when carrying out such an operation.

Concerning the two-index model, we add constraints (32) in the model (23)-(25) of the 2OPP, following the proposals in Queiroz & Miyazawa (2013) and Nascimento et al. (2021). Once j is visited after i on the route S _ℛ , that is o _i < o _j , these constraints ensure that if item m of customer j is packed in the point $(s,t)\in {\mathcal{P}}_{jm}^W\times {\mathcal{P}}_{jm}^H$, then no item l of customer i must be packed in points (u, v) such that m blocks l from having free passage during the unloading operation at i.

$$\begin{gathered} \sum _{\{i\in S_{\mathcal{R}}|o_i<o_j\}}\sum _{l\in M_i}\sum _{\{u\in {\mathcal{P}}_{il}^W|s-w_{il}+1\le u\le s+w_{jm}-1\}}\sum _{\{v\in {\mathcal{P}}_{il}^H|v\le t-h_{il}\}}{y}_{luv}\le (1-y_{mst})\sum _{i\in S_{\mathcal{R}}}\left|M_i\right|, \\ \forall j\in S_{\mathcal{R}},\forall m\in M_j,\forall s\in {\mathcal{P}}_{jm}^W,\forall t\in {\mathcal{P}}_{jm}^H \end{gathered}$$ (32)

4 COMPUTATIONAL EXPERIMENTS

The numerical experiments consider 80 benchmark¹ instances; a subset of them was first solved by Iori et al. (2007), and the remaining ones by Gendreau et al. (2008). In such instances, organized into 16 sets, the number of customers (n) ranges from 15 to 36 and the total number of items (M) from 15 to 114. Each set contains five instances divided into classes. The number and dimensions of the items are randomly generated, considering homogeneous, horizontal, and vertical items. The traveling/edge cost c _ij is given by the real value of the Euclidean distance between i and j.

All models were coded in the C++ programming language, and the B&C framework of the Gurobi Optimizer 9.5 (with default parameters) was used. The experiments were carried out on a computer with a processor Intel Xeon E3-1245 3.50 GHz, 32 GB of RAM, and Ubuntu 16.04 LTS. To solve each instance, a time limit of 3,600 seconds was imposed on the B&C, including the time limit of 60 seconds imposed on each call when solving any of the models related to the 2OPP^W and 2OPP^H , and 120 seconds when solving the model related to the 2OPP. We first present the results of the unrestricted 2L-CVRP, discussing the performance of the developed integer programming models. Next, we present the results of the sequential 2L-CVRP obtained with the best model we have developed.

4.1 Unrestricted 2L-CVRP

Tables 2 and 3 have the results of the unrestricted 2L-CVRP. We show the value of the best solution found and the gap (in percentage) that Gurobi returned for each model. We also report the number of cuts generated by Procedures 1 and 2 (i.e., with the CVRPSEP package), the runtime (in seconds) spent with all calls to the models for the 2OPP^W and 2OPP^H , the number of cuts inserted after solving one of the models for the 2OPP^W and 2OPP^H , the runtime (in seconds) spent with all calls to the model for the 2OPP, and the number of cuts inserted after solving the model for the 2OPP. The entries marked with an “*” indicate that the imposed time limit was reached when solving at least one loading subproblem. In this case, the final solution may not be optimal.

Table 2
Results of the unrestricted 2L-CVRP for instances with 15 to 22 customers.

Instance					Complete model			Two-index model
											Cuts	Time	Cuts	Time	Cuts
Name	Class	n	|M|	|K|	Time	Sol.	GAP	Time	Sol.	GAP	CVRPSEP	2OPPW;H	2OPPW;H	2OPP	2OPP
E016-03m	1	15	15	3	268.44	278.73	0.00	0.07	278.73	0.00	136	0.00	0	0.00	0
	2	15	24	3	3,600.00	285.50	19.09	4.11	278.73	0.00	185	0.07	3	3.75	0
	3	15	31	3	3,600.00	-	-	59.96	284.52	0.00	301	1.28	12	57.82	0
	4	15	37	4	3,600.00	282.95	19.16	34.44	282.95	0.00	328	0.48	16	32.87	0
	5	15	45	4	3,600.00	280.60	21.70	155.30	278.73*	0.00	215	17.19	1	137.36	0
E016-05m	1	15	15	5	547.52	334.96	0.00	0.09	334.96	0.00	203	0.00	0	0.00	0
	2	15	25	5	3,600.00	334.96	14.86	1.04	334.96	0.00	203	0.00	0	0.84	0
	3	15	31	5	3,600.00	352.16	22.43	3.00	352.16	0.00	471	0.04	5	2.20	0
	4	15	40	5	3,600.00	334.96	16.77	3.88	334.96	0.00	203	0.03	0	3.56	0
	5	15	48	5	3,600.00	334.96	14.87	0.70	334.96	0.00	203	0.02	0	0.34	0
E021-04m	1	20	20	4	3,600.00	364.45	20.71	0.55	358.40	0.00	579	0.00	0	0.00	0
	2	20	29	5	3,600.00	387.70	32.49	79.87	387.70	0.00	2,608	0.43	264	18.93	0
	3	20	46	5	3,600.00	-	-	247.18	394.72	0.00	2,410	12.42	276	169.38	0
	4	20	44	5	3,600.00	379.14	30.57	131.20	362.41*	0.00	369	1.77	4	128.55	0
	5	20	49	5	3,600.00	377.23	30.14	7.19	358.40	0.00	340	0.08	1	5.89	0
E021-06m	1	20	20	6	3,600.00	430.88	19.81	1.30	430.88	0.00	938	0.00	0	0.01	0
	2	20	32	6	3,600.00	430.88	29.41	3.39	430.88	0.00	938	0.02	0	1.88	0
	3	20	43	6	3,600.00	445.25	29.56	8.25	430.88	0.00	938	0.06	0	6.52	0
	4	20	50	6	3,600.00	447.37	32.17	8.69	447.37	0.00	1,341	0.07	3	3.71	0
	5	20	62	6	3,600.00	440.94	31.26	4.17	430.88	0.00	938	0.02	0	2.77	0
E022-04g	1	21	21	4	3,600.00	383.52	20.23	0.11	375.28	0.00	114	0.00	0	0.00	0
	2	21	31	4	3,600.00	376.51	27.07	9.38	375.28	0.00	143	0.06	3	8.96	0
	3	21	37	4	3,600.00	402.46	30.18	5.37	381.69	0.00	121	0.14	3	4.77	0
	4	21	41	4	3,600.00	383.87	28.24	2.92	383.87	0.00	174	0.04	2	2.48	0
	5	21	57	5	3,600.00	375.28	25.15	0.70	375.28	0.00	71	0.01	0	0.47	0
E022-06m	1	21	21	6	3,600.00	498.16	29.76	0.17	495.85	0.00	263	0.00	0	0.00	0
	2	21	33	6	3,600.00	499.08	37.60	2.76	495.85	0.00	303	0.02	1	2.36	0
	3	21	40	6	3,600.00	504.68	38.99	114.45	498.16	0.00	610	0.38	5	112.72	1
	4	21	57	6	3,600.00	602.94	50.14	27.04	498.32	0.00	783	0.34	7	24.24	0
	5	21	56	6	3,600.00	495.85	37.07	1.46	495.85	0.00	263	0.01	0	0.70	0
E023-03g	1	22	22	3	1,014.34	568.56	0.00	0.19	568.56	0.00	45	0.00	0	0.01	0
	2	22	32	5	3,600.00	750.64	34.69	3,600.00	750.64*	22.39	15,424	28.22	2,722	168.16	0
	3	22	41	5	3,600.00	729.97	31.76	3,600.00	729.97*	17.31	2,286	113.10	1,284	1,836.35	0
	4	22	51	5	3,600.00	768.78	36.36	3,600.00	727.34*	15.19	1,662	154.76	1,004	3,029.68	1
	5	22	55	6	3,600.00	725.57	32.69	3,600.00	661.41*	9.95	425	1,069.23	191	2,061.02	0
E023-05s	1	22	22	5	1,034.22	568.56	0.00	0.11	568.56	0.00	45	0.00	0	0.01	0
	2	22	29	5	3,600.00	720.77	31.08	456.69	674.55	0.00	3,856	3.08	319	83.41	0
	3	22	42	5	3,600.00	771.89	37.02	3,600.00	756.04*	14.85	7,951	123.32	4,380	392.13	1
	4	22	48	5	3,600.00	730.63	31.35	3,600.00	708.47*	10.25	3,532	195.39	2,089	2,388.38	0
	5	22	52	6	3,600.00	642.22	23.57	3,600.00	612.79*	4.84	158	1,013.08	45	2,006.30	0

Table 3
Results of the unrestricted 2L-CVRP for instances with 25 to 35 customers.

Instance					Complete model			Two-index model
											Cuts	Time	Cuts	Time	Cuts
Name	Class	n	|M|	|K|	Time	Sol.	GAP	Time	Sol.	GAP	CVRPSEP	2OPPW;H	2OPPW;H	2OPP	2OPP
E026-08m	1	25	25	8	3,600.00	620.25	32.10	0.24	607.65	0.00	349	0.00	0	0.00	0
	2	25	40	8	3,600.00	624.94	35.28	2.59	607.65	0.00	349	0.00	0	2.23	0
	3	25	61	8	3,600.00	-	-	2.33	607.65	0.00	349	0.00	0	1.91	0
	4	25	63	8	3,600.00	-	-	16.31	625.10	0.00	2,222	0.07	1	5.04	0
	5	25	91	8	3,600.00	-	-	1.30	607.65	0.00	349	0.00	0	0.74	0
E030-03g	1	29	29	3	3,600.00	535.80	35.98	3.79	535.80	0.00	619	0.00	0	0.01	0
	2	29	43	6	3,600.00	-	-	3,600.00	787.50	28.46	20,917	32.21	3,028	315.91	0
	3	29	49	6	3,600.00	659.18	52.81	3,600.00	647.67*	18.83	8,724	146.34	2,269	1,229.01	1
	4	29	72	7	3,600.00	-	-	3,600.00	-	-	-	-	-	-	-
	5	29	86	7	3,600.00	762.25	60.12	3,600.00	-	-	-	-	-	-	-
E030-04s	1	29	29	4	3,600.00	505.01	37.85	0.14	505.01	0.00	75	0.00	0	0.01	0
	2	29	43	6	3,600.00	880.73	64.98	3,600.00	-	-	-	-	-	-	-
	3	29	62	7	3,600.00	-	-	3,600.00	1,020.76*	47.64	13,987	159.75	3,275	1,846.74	0
	4	29	74	7	3,600.00	-	-	3,600.00	-	-	-	-	-	-	-
	5	29	91	7	3,600.00	-	-	3,600.00	804.92*	36.09	1,255	772.81	378	2,501.07	0
E031-09h	1	30	30	9	3,600.00	632.03	40.70	530.33	610.00	0.00	10,508	0.01	0	0.01	0
	2	30	50	9	3,600.00	-	-	3,600.00	657.57	8.99	28,511	0.08	12	18.67	0
	3	30	56	9	3,600.00	645.02	43.67	801.55	610.00	0.00	12,379	0.02	1	3.98	0
	4	30	82	9	3,600.00	-	-	3,600.00	614.23	0.84	20,449	1.43	20	22.20	0
	5	30	101	9	3,600.00	-	-	567.22	610.00	0.00	10,508	2.43	0	23.54	0
E033-03n	1	32	32	3	3,600.00	2,058.86	16.68	0.77	2,006.34	0.00	252	0.01	0	0.02	0
	2	32	44	7	3,600.00	2,872.95	51.03	3,600.00	-	-	-	-	-	-	-
	3	32	56	7	3,600.00	2,967.01	52.66	3,600.00	-	-	-	-	-	-	-
	4	32	78	7	3,600.00	-	-	3,600.00	-	-	-	-	-	-	-
	5	32	102	8	3,600.00	2,727.78	48.87	3,600.00	2,727.78*	27.18	5,058	796.76	1325	2,509.95	0
E033-04g	1	32	32	4	3,600.00	872.89	56.32	1.73	837.67	0.00	612	0.01	0	0.02	0
	2	32	47	7	3,600.00	1,258.00	70.97	3,600.00	1,258.00*	33.43	22,839	36.41	1,733	639.75	1
	3	32	57	7	3,600.00	1,302.46	71.91	3,600.00	-	-	-	-	-	-	-
	4	32	65	7	3,600.00	1,194.04	69.37	3,600.00	1,194.04*	38.52	2,018	1,334.69	199	2,145.12	0
	5	32	87	8	3,600.00	1,245.73	70.69	3,600.00	1,017.84*	17.70	2,264	847.91	135	2,411.48	0
E033-05s	1	32	32	5	3,600.00	887.15	57.14	0.77	837.67	0.00	454	0.00	0	0.01	0
	2	32	48	6	3,600.00	1,248.55	70.63	3,600.00	1,248.55	33.32	23,278	49.54	1,807	422.96	0
	3	32	59	6	3,600.00	-	-	3,600.00	-	-	-	-	-	-	-
	4	32	84	8	3,600.00	-	-	3,600.00	-	-	-	-	-	-	-
	5	32	114	8	3,600.00	-	-	3,600.00	-	-	-	-	-	-	-
E036-11h	1	35	35	11	3,600.00	704.01	41.50	180.66	698.61	0.00	7,133	0.01	0	0.02	0
	2	35	56	11	3,600.00	-	-	185.81	698.61	0.00	7,133	0.01	0	6.01	0
	3	35	74	11	3,600.00	-	-	239.92	698.61	0.00	8,045	0.04	1	4.30	0
	4	35	93	11	3,600.00	-	-	405.67	703.35	0.00	10,769	0.33	5	13.33	0
	5	35	114	11	3,600.00	-	-	183.22	698.61	0.00	7,133	0.03	0	0.46	0

In the results of Tables 2 and 3, the complete model did not find a feasible solution for 21 instances within the imposed time limit, especially those with more customers and items per customer. Considering the two-index model, the number of instances for which a solution was not found was 11 out of 80 instances. In summary, concerning the solved instances:

• feasible solutions: the complete model found an optimal solution for 4 instances, while for the 55 remaining instances, the average, minimum, and maximum gaps are 37.26%, 14.86%, and 71.91%, respectively. The two-index model found an optimal solution for 51 instances, while for the 18 remaining instances, the average, minimum, and maximum gaps are 21.43%, 0.84%, and 47.64%, respectively. If considering the 54 instances for which both models found a solution, the average gap is 32.38% and 4.88% for the complete and two-index models, respectively.

• runtime: the complete model did not reach the time limit for 4 instances. If considering the 59 instances for which the complete model returned a solution, the average runtime is 3,404.48 seconds. The two-index model did not reach the time limit for 51 instances, with an average runtime of 1,004.35 seconds when considering the 69 instances with a solution. On the other hand, if considering all 80 instances, the average runtime is 3,455.81 and 1,361.25 seconds for the complete model and the two-index model, respectively.

From the comparison above, it is clear that the two-index model presented a better performance than the complete model in terms of solved instances, gap, and runtime. Thus, concerning the separation routines of the two index model, the results in Tables 2 and 3 show for the 69 instances with a solution:

• with the CVRPSEP package: the average, minimum, and maximum number of generated cuts are 4,067, 45, and 28,511, respectively;

• with the models for the 2OPP^W and 2OPP^H: the average, minimum, and maximum number of generated cuts are 389, 0, and 4,380, respectively. The number of instances for which at least one cut was generated is 39. Regarding the average, minimum, and maximum runtime, such values are 100.23, less than 0.01, and 1,334.69 seconds, respectively;

• with the model for the 2OPP: the average, minimum, and maximum number of generated cuts are 0, 0, and 1, respectively. The number of instances for which at least one cut was generated is 5. Regarding the average, minimum, and maximum spent runtime, such values are 388.74, less than 0.01, and 3,029.68 seconds, respectively;

• impact of the 2OPP^W and 2OPP^H on the 2OPP: for the 39 instances that the models of the 2OPP^W or 2OPP^H generated cuts, the model of the 2OPP needed to further generate cuts for only 5 instances.

4.2 Sequential 2L-CVRP

Table 4 reports the results of only the two-index model for the sequential 2L-CVRP since it performs better than the complete model. This table has the same information as presented in Tables 2 and 3. For this variant, with the two-index model, it was not possible to find a feasible solution for 23 out of 80 instances. This indicates that inserting multi-drop requirements demands more computing time since the loading subproblems become more complex. In summary, concerning the 57 instances solved, the two-index model has the following characteristics:

• feasible solutions: it found an optimal solution for 37 instances, while for the 20 remaining instances, the average, minimum, and maximum gaps are 13.81%, 1.24%, and 30.48%, respectively;

• runtime: it did not reach the time limit for 37 instances, with the average runtime of 1,532.73 seconds when considering the 57 instances with a solution. On the other hand, if considering all 80 instances, the average runtime is 2,127.07 seconds;

• with the CVRPSEP package: the average, minimum, and maximum number of generated cuts are 6,927, 42, and 38,459, respectively;

• with the models for the 2OPP^W and 2OPP^H: the average, minimum, and maximum number of generated cuts are 64, 0, and 1,682, respectively. The number of instances for which at least one cut was generated is 28. Regarding the average, minimum, and maximum runtime, such values are 3.09, less than 0.01, and 62.30 seconds, respectively;

• with the model for the 2OPP: the average, minimum, and maximum number of generated cuts are 7, 0, and 133, respectively. The number of instances for which at least one cut was generated is 18. Regarding the average, minimum, and maximum runtime, such values are 966.48, less than 0.01, and 3,575.91 seconds, respectively;

• impact of the 2OPP^W and 2OPP^H on the 2OPP: for the 28 instances that the models of the 2OPP^W or 2OPP^H generated cuts, the model of the 2OPP also generated cuts for 18 of them.

Table 4
Results of the sequential 2L-CVRP for all 80 instances.

Instance		Two-index model								Instance		Two-index model
					Cuts	Time	Cuts	Time	Cuts						Cuts	Time	Cuts	Time	Cuts
Name	Class	Time	Sol.	GAP	CVRPSEP	2OPPW;H	2OPPW;H	2OPP	2OPP	Name	Class	Time	Sol.	GAP	CVRPSEP	2OPPW;H	2OPPW;H	2OPP	2OPP
E016-03m	1	0.80	278.73	0.00	343	0.00	0	0.02	0	E026-08m	1	7.01	607.65	0.00	2,537	0.00	0	0.01	0
	2	3,600.00	312.69	10.68	6,017	13.29	361	3,475.23	133		2	27.51	611.49	0.00	2,399	0.00	0	20.27	0
	3	3,600.00	285.71*	7.29	1,163	7.44	54	3,244.48	1		3	372.17	613.90*	0.00	2,537	0.01	0	363.01	0
	4	1,536.53	294.25*	0.00	967	3.67	33	1,516.15	0		4	1406.15	625.10*	0.00	11,156	0.18	12	1,084.30	0
	5	1,396.37	278.73*	0.00	599	0.78	4	1,377.52	0		5	195.11	607.65	0.00	2,537	0.01	0	183.12	0
E016-05m	1	0.53	334.96	0.00	379	0.00	0	0.00	0	E030-03g	1	19.63	535.80	0.00	2,165	0.01	0	0.22	0
	2	5.77	347.73	0.00	379	0.00	0	5.06	0		2	3,600.00	-	-	-	-	-	-	-
	3	439.85	353.48	0.00	3,842	0.37	28	369.30	6		3	3,600.00	676.83*	25.28	532	6.04	170	3,517.67	1
	4	122.87	342.00*	0.00	379	0.04	0	120.29	0		4	3,600.00	-	-	-	-	-	-	-
	5	124.54	334.96*	0.00	379	0.05	0	120.42	0		5	3,600.00	-	-	-	-	-	-	-
E021-04m	1	5.73	358.40	0.00	1604	0.00	0	0.03	0	E030-04s	1	0.39	505.01	0.00	105	0.00	0	0.04	0
	2	3,600.00	-	-	-	-	-	-	-		2	3,600.00	-	-	-	-	-	-	-
	3	3,600.00	406.55*	12.23	7,129	1.61	71	1,790.8	0		3	3,,600.00	-	-	-	-	-	-	-
	4	2,260.70	368.56*	0.00	2,440	7.11	22	2,217.75	0		4	3600.00	-	-	-	-	-	-	-
	5	504.27	358.40*	0.00	1,372	0.01	0	493.14	0		5	3,600.00	-	-	-	-	-	-	-
E021-06m	1	5.20	430.88	0.00	2,297	0.00	0	0.01	0	E031-09h	1	3,600.00	648.74	8.84	36,525	0.00	0	0.03	0
	2	220.06	440.94	0.00	2,937	0.06	2	199.02	10		2	3,600.00	-	-	-	-	-	-	-
	3	246.85	440.68*	0.00	2,297	0.04	0	240.11	0		3	3,600.00	637.74	7.00	36,635	1.08	8	86.35	0
	4	146.28	447.37	0.00	3,445	0.15	6	115.97	0		4	3,600.00	675.65*	13.04	2,1675	1.42	30	2,614.28	1
	5	96.49	430.88	0.00	2,297	0.01	0	88.09	0		5	3,600.00	648.74*	8.91	33,029	0.20	0	1,002.26	0
E022-04g	1	0.17	375.28	0.00	135	0.00	0	0.01	0	E033-03n	1	1.20	2,006.34	0.00	203	0.00	0	0.12	0
	2	3,600.00	388.72*	4.22	2,486	3.31	44	3,575.91	69		2	3,600.00	-	-	-	-	-	-	-
	3	1291.53	387.82*	0.00	901	0.71	24	1,276.21	15		3	3,600.00	-	-	-	-	-	-	-
	4	1981.61	383.87*	0.00	1,127	1.93	22	1,963.77	1		4	3,600.00	-	-	-	-	-	-	-
	5	72.32	375.28	0.00	149	0.01	0	69.31	0		5	3,600.00	-	-	-	-	-	-	-
E022-06m	1	1.84	495.85	0.00	914	0.00	0	0.01	0	E033-04g	1	4.70	837.67	0.00	1315	0.00	0	0.11	0
	2	77.45	499.08	0.00	1,700	0.05	3	68.56	2		2	3,600.00	-	-	-	-	-	-	-
	3	964.51	504.68	0.00	6,975	0.82	30	871.19	12		3	3,600.00	1,199.60*	30.26	2,611	11.47	97	1,781.36	1
	4	1,715.21	498.32*	0.00	2,968	0.68	21	1,680.56	0		4	3,600.00	-	-	-	-	-	-	-
	5	111.49	495.85	0.00	914	0.01	0	105.27	0		5	3,600.00	-	-	-	-	-	-	-
E023-03g	1	0.28	568.56	0.00	42	0.00	0	0.04	0	E033-05s	1	2.17	837.67	0.00	650	0.00	0	0.08	0
	2	3,600.00	860.70*	28.04	6,481	14.06	1,682	3,289.44	103		2	3,600.00	-	-	-	-	-	-	-
	3	3,600.00	706.59*	21.61	993	24.53	467	3,452.35	4		3	3,600.00	-	-	-	-	-	-	-
	4	3,600.00	-	-	-	-	-	-	-		4	3,600.00	-	-	-	-	-	-	-
	5	3,600.00	-	-	-	-	-	-	-		5	3,600.00	-	-	-	-	-	-	-
E023-05s	1	0.29	568.56	0.00	42	0.00	0	0.04	0	E036-11h	1	3,600.00	698.61	1.24	35,540	0.01	0	0.05	0
	2	3,600.00	725.91*	30.48	1,663	3.17	204	3,370.95	6		2	3,600.00	740.49	8.31	38,459	0.04	6	98.20	13
	3	3,600.00	818.21*	29.81	454	7.64	152	3,476.83	10		3	3,600.00	698.61	1.32	34702	0.27	1	130.29	0
	4	3,600.00	-	-	-	-	-	-	-		4	3,600.00	779.49*	12.96	25,743	1.42	19	2,165.22	3
	5	3,600.00	657.49*	13.35	213	62.30	50	3,238.68	0		5	3,600.00	698.61	1.25	35363	0.04	0	29.62	0

Another comparison is carried out between the unrestricted and sequential 2L-CVRP for the results of the two-index model. From Tables 2 to 4, considering the same 56 out of 80 instances for which such a model returned a solution in both variants (notice that these instances include the ones discussed in the previous paragraph for the sequential case):

• cost of solutions: the percentage deviation (increase) in the value of solutions when comparing the unrestricted with the sequential variants is 2.22%, with the minimum and maximum of −3.20% and 14.66%, respectively;

• feasible solutions: in the unrestricted variant, it found an optimal solution for 50 instances, while for the 2 remaining instances, the average, minimum, and maximum gaps are 13.18%, 0.84%, and 22.39%, respectively. The average gap, given the 56 instances, in the unrestricted and sequential variants is 1.41% and 4.39%, respectively;

• runtime: in the unrestricted variant, it did not reach the time limit for 50 instances. The average runtime, given the 56 instances, is 464.65 and 1, 495.81 seconds in the unrestricted and sequential variants, respectively;

• with the CVRPSEP package: in the unrestricted variant, the average, minimum, and maximum number of generated cuts are 2,719, 45, and 20,449, respectively. In terms of the average number of generated cuts, the percentage increase from the unrestricted variant to the sequential one corresponds to 157.60%;

• with the models for the 2OPP^W and 2OPP^H: in the unrestricted variant, the average, minimum, and maximum number of generated cuts are 203, 0, and 4,380, respectively. The number of instances for which at least one cut was generated is 26. Regarding the average, minimum, and maximum runtime, such values are 26.18, less than 0.01, and 1,013.08 seconds, respectively. In terms of the average number of generated cuts, the percentage decrease from the unrestricted variant to the sequential one corresponds to 69.04%;

• with the model for the 2OPP: in the unrestricted variant, only 1 cut was generated for 3 instances. Regarding the average, minimum, and maximum runtime, such values are 116.31, less than 0.01, and 2,006.30 seconds, respectively.

4.3 Comparing with the literature

Results are compared with those from the Simulated Annealing (SA) of Wei et al. (2018) that reported the best overall results among all authors in Table 1 for the 2L-CVRP. They also solved the same problem (and variants) we are handling: the 2L-CVRP allowing routes with a single customer and fewer routes than the number of vehicles. The SA of Wei et al. (2018) had a time limit of 1,800 seconds (for each instance) when solving the above-mentioned 80 instances, and it ran on a different computer. Then, a comparison concerning the runtime is not taken into consideration. For the two-index model, a time limit of 1,800 seconds was imposed to the B&C, including the time limit of 30 seconds imposed on each call when solving any of the 2OPP^W and 2OPP^H , and 60 seconds when solving the 2OPP. Table 5 has the results with respect to the best solution found and the percentage deviation between the found solution Sol and the solution reported by Wei et al. (2018), which is given by $100\times \frac{Sol-Wei}{Wei}$. In summary, we observed that:

• solutions for the unrestricted variant: it found a feasible solution for 69 out of 80 instances, where 51 are optimal. Notice that in Tables 2 and 3, for the time limit of 3,600 seconds, the number of feasible solutions is 69, for which 51 are optimal;

• solutions for the sequential variant: it found a feasible solution for 55 instances, where 37 are optimal. Notice that in Table 4, for the time limit of 3,600 seconds, the number of feasible solutions is 57, for which 37 are optimal;

• solutions of the SA in the unrestricted variant: it found equal solutions for 50 instances, while for the 19 remaining instances, solutions were worse, with an average, minimum, and maximum deviation of 15.86%, 0.02%, and 44.43%, respectively;

• solutions of the SA in the sequential variant: it found equal solutions for 40 instances, while for the 15 remaining instances, solutions were worse, with an average, minimum, and maximum deviation of 6.40%, 0.23%, and 18.01%, respectively.

Table 5
Comparison between a heuristic method and the two-index model.

unrestricted 2L-CVRP										sequential 2L-CVRP
Instance		SA	Two-index		Instance		SA	Two-index		Instance		SA	Two-index		Instance		SA	Two-index
Name	Class	Sol.	Sol.	Dev.	Name	Class	Sol.	Sol.	Dev.	Name	Class	Sol.	Sol.	Dev.	Name	Class	Sol.	Sol.	Dev.
E016-03m	1	278.73	278.73	0.00	E026-08m	1	607.65	607.65	0.00	E016-03m	1	278.73	278.73	0.00	E026-08m	1	607.65	607.65	0.00
	2	278.73	278.73	0.00		2	607.65	607.65	0.00		2	290.84	-	-		2	611.49	611.49	0.00
	3	284.52	284.52	0.00		3	607.65	607.65	0.00		3	284.52	285.71	0.42		3	613.90	613.90	0.00
	4	282.95	282.95	0.00		4	625.10	625.10	0.00		4	294.25	294.25	0.00		4	625.10	625.10	0.00
	5	278.73	278.73	0.00		5	607.65	607.65	0.00		5	278.73	278.73	0.00		5	607.65	607.65	0.00
E016-05m	1	334.96	334.96	0.00	E030-03g	1	535.80	535.80	0.00	E016-05m	1	334.96	334.96	0.00	E030-03g	1	535.80	535.80	0.00
	2	334.96	334.96	0.00		2	689.68	793.41	15.04		2	347.73	347.73	0.00		2	700.20	-	-
	3	352.16	352.16	0.00		3	620.33	754.10	21.56		3	352.16	352.98	0.23		3	639.09	676.83	5.91
	4	334.96	334.96	0.00		4	710.87	-	-		4	342.00	342.00	0.00		4	715.82	-	-
	5	334.96	334.96	0.00		5	678.66	-	-		5	334.96	334.96	0.00		5	690.96	-	-
E021-04m	1	358.40	358.40	0.00	E030-04s	1	505.01	505.01	0.00	E021-04m	1	358.40	358.40	0.00	E030-04s	1	505.01	505.01	0.00
	2	387.70	387.70	0.00		2	684.21	-	-		2	403.93	-	-		2	721.54	-	-
	3	394.72	394.72	0.00		3	706.73	1,020.76	44.43		3	394.72	406.55	3.00		3	719.91	-	-
	4	362.41	362.41	0.00		4	788.54	-	-		4	368.56	368.56	0.00		4	822.61	-	-
	5	358.40	358.40	0.00		5	624.82	804.92	28.82		5	358.40	358.40	0.00		5	624.82	-	-
E021-06m	1	430.89	430.88	0.00	E031-09h	1	610.00	610.00	0.00	E021-06m	1	430.89	430.88	0.00	E031-09h	1	610.00	648.74	6.35
	2	430.89	430.88	0.00		2	610.57	674.96	10.55		2	440.94	440.94	0.00		2	619.63	-	-
	3	430.89	430.88	0.00		3	610.00	610.00	0.00		3	440.68	440.68	0.00		3	610.00	637.74	4.55
	4	447.37	447.37	0.00		4	614.24	623.39	1.49		4	447.37	447.37	0.00		4	618.23	688.82	11.42
	5	430.89	430.88	0.00		5	610.00	610.00	0.00		5	430.89	430.88	0.00		5	610.00	648.74	6.35
E022-04g	1	375.28	375.28	0.00	E033-03n	1	2,006.34	2,006.34	0.00	E022-04g	1	375.28	375.28	0.00	E033-03n	1	2,006.34	2,006.34	0.00
	2	375.28	375.28	0.00		2	2,585.72	-	-		2	388.72	388.72	0.00		2	2,669.39	-	-
	3	381.69	381.69	0.00		3	2,454.37	-	-		3	381.69	387.82	1.61		3	2,486.44	-	-
	4	383.88	383.87	0.00		4	2,548.06	-	-		4	383.88	383.87	0.00		4	2,609.36	-	-
	5	375.28	375.28	0.00		5	2,334.78	2,985.29	27.86		5	375.28	375.28	0.00		5	2,416.04	-	-
E022-06m	1	495.85	495.85	0.00	E033-04g	1	837.67	837.67	0.00	E022-06m	1	495.85	495.85	0.00	E033-04g	1	837.67	837.67	0.00
	2	495.85	495.85	0.00		2	1,038.20	1,275.49	22.86		2	499.08	499.08	0.00		2	1,092.51	-	-
	3	498.16	498.16	0.00		3	1,002.51	-	-		3	504.68	504.68	0.00		3	1,039.06	1,199.60	15.45
	4	498.32	498.32	0.00		4	981.00	1,363.46	38.99		4	498.32	498.32	0.00		4	982.25	-	-
	5	495.75	495.85	0.02		5	871.22	1,017.84	16.83		5	495.85	495.85	0.00		5	922.75	-	-
E023-03g	1	568.56	568.56	0.00	E033-05s	1	837.67	837.67	0.00	E023-03g	1	568.56	568.56	0.00	E033-05s	1	837.67	837.67	0.00
	2	725.46	854.55	17.79		2	1,017.95	1,329.40	30.60		2	734.65	866.98	18.01		2	1,041.75	-	-
	3	678.75	749.88	10.48		3	1,149.99	-	-		3	702.59	706.59	0.57		3	1,181.68	-	-
	4	700.72	716.05	2.19		4	1,184.85	-	-		4	703.49	-	-		4	1,247.39	-	-
	5	657.77	662.04	0.65		5	1,159.94	-	-		5	658.64	-	-		5	1,229.95	-	-
E023-05s	1	568.56	568.56	0.00	E036-11h	1	698.61	698.61	0.00	E023-05s	1	568.56	568.56	0.00	E036-11h	1	698.61	698.61	0.00
	2	674.55	674.55	0.00		2	698.61	698.61	0.00		2	725.91	725.91	0.00		2	698.61	742.83	6.33
	3	740.85	775.36	4.66		3	698.61	698.61	0.00		3	741.12	-	-		3	698.61	698.61	0.00
	4	692.47	735.01	6.14		4	703.35	703.35	0.00		4	697.92	-	-		4	708.20	779.49	10.07
	5	609.90	612.79	0.47		5	698.61	698.61	0.00		5	621.85	657.49	5.73		5	698.61	698.61	0.00

One last experiment was carried out with the two-index model, comparing it with the branch-and- cut algorithm (BNC) in Hokama et al. (2016). We solved the same instances and assumed the same constraints the authors imposed when solving the sequential 2L-CVRP. We also considered the same time limit of 3,600 seconds to solve each instance, including the time limit of 60 seconds imposed on each call when solving any of the 2OPP^W and 2OPP^H , and 120 seconds when solving the 2OPP. In Table 6, we present the value of the best solution found, the runtime (in seconds), and the percentage deviation between our solution and the solution reported by Hokama et al. (2016).

Table 6
Comparison with the branch-and-cut of Hokama et al. (2016) for the sequential 2L-CVRP.

Instances		Hokama et al. (2016)	Two-index model			Instances		Hokama et al. (2016)	Two-index model
Name	Class	Sol.	Time	Sol.	Dev.	Name	Class	Sol.	Time	Sol.	Dev.
E016-03m	1	273	0.65	273	0.00	E023-03g	1	558	0.44	558	0.00
	2	285	3,600.00	285*	0.00		2	724	3,600.00	724	0.00
	3	280	2658.45	280*	0.00		3	698	3,600.00	698*	0.00
	4	288	57.51	288	0.00		4	714	3,600.00	714*	0.00
	5	279	253.69	279	0.00		5	742	1,533.29	742*	0.00
E016-05m	1	329	1.03	329	0.00	E023-05s	1	657	0.44	657	0.00
	2	342	9.77	342	0.00		2	720	553.63	720*	0.00
	3	347	301.22	347	0.00		3	730	3,600.00	730*	0.00
	4	336	131.32	336*	0.00		4	701	3,600.00	701*	0.00
	5	329	248.69	329*	0.00		5	721	122.67	721*	0.00
E021-04m	1	351	1.58	351	0.00	E026-08m	1	609	0.85	609	0.00
	2	396	988.5	396	0.00		2	612	15.33	612	0.00
	3	387	3,600.00	387*	0.00		3	615	195.93	615*	0.00
	4	374	1030	374*	0.00		4	626	540.25	626*	0.00
	5	369	88.69	369	0.00		5	609	155.38	609	0.00
E021-06m	1	423	0.45	423	0.00	E030-03g	1	524	13.40	524	0.00
	2	434	276.3	434	0.00		2	-	3,600.00	687*	-
	3	432	349.34	432*	0.00		3	637	3,600.00	637*	0.00
	4	438	311.46	438	0.00		4	-	3,600.00	-	-
	5	423	197.84	423*	0.00		5	-	3,600.00	-	-
E022-04g	1	367	0.16	367	0.00	E033-03n	1	1,991	0.79	1,991	0.00
	2	380	1,956.46	380	0.00		2	-	3,600.00	-	-
	3	373	2,351.18	373*	0.00		3	2854	3,600.00	2854*	0.00
	4	377	760.45	377*	0.00		4	-	3,600.00	-	-
	5	389	44.41	389	0.00		5	2,672	3,600.00	2,672*	0.00
E022-06m	1	488	3.39	488	0.00	E036-11h	1	682	772.98	682	0.00
	2	491	64.52	491	0.00		2	682	3,600.00	682	0.00
	3	496	2,099.52	496*	0.00		3	682	937.80	682	0.00
	4	489	392.55	489*	0.00		4	691	3,600.00	691*	0.00
	5	488	62.41	488	0.00		5	-	822.03	682	-

Observing the results of Table 6, the branch-and-cut of Hokama et al. (2016) found a solution to 54 out of 60 instances within the imposed time limit. Using the two-index model, we obtained the same solution these authors obtained for these 54 instances, besides returning a solution to two more instances (E030-03g in class 2 and E036-11h in class 5). In summary, with relation to the solved instances, the two-index model has the following characteristics:

• feasible solutions: it found an optimal solution for 43 instances, although the time limit to solve any of the loading subproblems was reached for 26 instances;

• solutions of the branch-and-cut: it found equal solutions for 54 instances and better solutions for 2 instances. For the 4 remaining instances, no solution was found within the imposed time limit. The branch-and-cut of Hokama et al. (2016) did not find any solution for 6 instances;

• runtime: it did not reach the time limit for 43 instances, with the average runtime of 1,198.33 seconds when considering all the 56 instances with a solution. If considering all the 60 instances, the average runtime increases to 1,358.45 seconds.

5 CONCLUDING REMARKS

In this work, a combined vehicle routing problem with a bin packing problem (i.e., the 2L- CVRP) has its unrestricted and sequential variants handled with a branch-and-cut algorithm. Two integer linear programming models are discussed. The first one has all constraints added at once, while the second has cutting planes related to the connectivity, capacity, and loading of routes. The separation routines deal with the rounded capacity inequalities and the resolution of integer programming models for the loading subproblem of each route.

After the computational experiments conducted on 80 benchmark instances, the models’ performance could be clearly stated. Concerning the unrestricted 2L-CVRP, the two-index model is superior to the complete model, mainly in terms of optimal solutions, gap, and runtime. Although such a model obtained better results, it depends on separation routines to detect violated inequalities and guarantee a feasible solution. This dependence may be complicated when compared with a complete model that, in turn, can be solved with any solver. The complete model allowed to solve instances with up 23 customers, while the two-index model could solve instances with up 35 customers. The latter is also very competitive (even with a limitation in the runtime) when compared with a heuristic from the literature.

In the sequential variant, where multi-drop requirements are imposed on the 2L-CVRP, obtaining an optimal solution requires more computing time. The results of the two-index model for this variant show how the number of optimal solutions may decrease, and the gap may increase compared to the unrestricted variant. Regarding the comparison with a heuristic from the literature, the two-index model is still an exciting choice when a shorter time limit is imposed. Moreover, the two-index model is simpler and could be an alternative to other models from the literature since when comparing it with another branch-and-cut, it could return equal or better solutions for all the instances.

Future works could focus on practical constraints to the loading subproblem, particularly load balancing and cargo stability (Queiroz & Miyazawa, 2013; Oliveira et al., 2021a,b; Junqueira & Queiroz, 2022) since the literature in these constraints is scarce. Another interesting direction is to investigate infeasible routes further by looking for sub-routes that are still infeasible, especially concerning the loading subproblems. Thus, stronger inequalities may be derived and the computing time for solving large instances may be reduced.

Acknowledgements

The authors would like to thank the financial support of the Coordination for the Improvement of Higher Education Personnel (CAPES grant number 001), the National Council of Technological and Scientific Development (CNPq - grants 306918/2014-5, 165896/2015-9, 308312/2016-3, 308761/2018-9, 311185/2020-7) and the State of Goiás Research Foundation (FAPEG).

References

• AZEVEDO BLP, HOKAMA P, MIYAZAWA FK & XAVIER EC. 2009. A Branch-and-Cut Approach for the Vehicle Routing Problem with Two-dimensional Loading Constraints. In: XLI Simpósio Brasileiro de Pesquisa Operacional. pp. 2491-2502. Brazil.
• BIN W, HONG C & ZHI-YONG C. 2013. Artificial bee colony algorithm for two-dimensional loading capacitated vehicle routing problem. In: 2013 International Conference on Management Science and Engineering (ICMSE). pp. 406-412.
• CLARKE G & WRIGHT JW. 1964. Scheduling of Vehicles from a Central Depot to a Number of Delivery Points. Operations Research, 12(4): 568-581.
• CÔTÉ JF, GENDREAU M & POTVIN JY. 2014. An exact algorithm for the two-dimensional orthogonal packing problem with unloading constraints. Operations Research , 62(5): 1126-1141.
• CÔTÉ JF, GENDREAU M & POTVIN JY. 2020. The Vehicle Routing Problem with Stochastic Two-Dimensional Items. Transportation Science, 54(2): 453-469.
• CÔTÉ JF, GUASTAROBA G & SPERANZA MG. 2017. The value of integrating loading and routing. European Journal of Operational Research, 257(1): 89-105.
• CÔTÉ JF, HAOUARI M & IORI M. 2021. Combinatorial Benders Decomposition for the Two- Dimensional Bin Packing Problem. INFORMS Journal on Computing, 33(3): 963-978.
• DANTZIG GB & RAMSER JH. 1959. The Truck Dispatching Problem. Management Science, 6(1): 80-91.
• DOMÍNGUEZ O, JUAN AA, BARRIOS B, FAULIN J & AGUSTIN A. 2016. Using biased randomization for solving the two-dimensional loading vehicle routing problem with heterogeneous fleet. Annals of Operations Research , 236(2): 383-404.
• DOMÍNGUEZ O, JUAN AA & FAULIN J. 2014. A biased-randomized algorithm for the two- dimensional vehicle routing problem with and without item rotations. International Transactions in Operational Research, 21(3): 375-398.
• DUHAMEL C, LACOMME P, QUILLIOT A & TOUSSAINT H. 2011. A multi-start evolutionary local search for the two-dimensional loading capacitated vehicle routing problem. Computers & Operations Research , 38(3): 617-640.
• FERREIRA KM & QUEIROZ TA. 2022. A simulated annealing based heuristic for a location- routing problem with two-dimensional loading constraints. Applied Soft Computing, p. 108443.
• FERREIRA KM, QUEIROZ TA & TOLEDO FMB. 2021. An exact approach for the green vehicle routing problem with two-dimensional loading constraints and split delivery. Computers & Operations Research , 136: 105452.
• FUELLERER G, DOERNER KF, HARTL RF & IORI M. 2009. Ant colony optimization for the two-dimensional loading vehicle routing problem. Computers & Operations Research , 36(3): 655-673.
• GENDREAU M, IORI M, LAPORTE G & MARTELLO S. 2008. A Tabu search heuristic for the vehicle routing problem with two-dimensional loading constraints. Networks, 51(1): 4-18.
• HERZ JC. 1972. Recursive computational procedure for two-dimensional stock cutting. IBM Journal of Research and Development, 16(5): 462-469.
• HOKAMA P, MIYAZAWA FK & XAVIER EC. 2016. A branch-and-cut approach for the vehicle routing problem with loading constraints. Expert Systems with Applications, 47: 1-13.
• IORI M & MARTELLO S. 2010. Routing Problems with Loading Constraints. TOP, 18(1): 4-27.
• IORI M, SALAZAR-GONZÁLEZ JJ & VIGO D. 2007. An exact approach for the vehicle routing problem with two-dimensional loading constraints. Transportation Science , 41(2): 253-264.
• JUNQUEIRA L, OLIVEIRA JF, CARRAVILLA MA & MORABITO R. 2013. An optimization model for the vehicle routing problem with practical three-dimensional loading constraints. International Transactions in Operational Research , 20(5): 645-666.
• JUNQUEIRA L & QUEIROZ TA. 2022. The static stability of support factor-based rectangular packings: an assessment by regression analysis. International Transactions in Operational Research , 29(1): 574-599.
• KANTOROVICH LV. 1960. Mathematical Methods of Organizing and Planning Production. Management Science , 6(4): 366-422.
• KARA I. 2010. Tightening Bounding Constraints of the Miller-Tucker-Zemlin Based Formulation of the Capacitated Vehicle Routing Problems and Some Extensions. In: Proceedings of the 2nd International Conference on Manufacturing Engineering, Quality and Production Systems. pp. 137-142.
• KHEBBACHE-HADJI S, PRINS C, YALAOUI A & REGHIOUI M. 2013. Heuristics and memetic algorithm for the two-dimensional loading capacitated vehicle routing problem with time windows. Central European Journal of Operations Research , 21(2): 307-336.
• LEUNG SCH, ZHENG J, ZHANG D & ZHOU X. 2010. Simulated annealing for the vehicle routing problem with two-dimensional loading constraints. Flexible Services and Manufacturing Journal, 22(1-2): 61-82.
• LEUNG SCH, ZHOU X, ZHANG D & ZHENG J. 2011. Extended guided tabu search and a new packing algorithm for the two-dimensional loading vehicle routing problem. Computers & Operations Research , 38(1): 205-215.
• LYSGAARD J, LETCHFORD AN & EGLESE RW. 2004. A new branch-and-cut algorithm for the capacitated vehicle routing problem. Mathematical Programming, 100(2): 423-445.
• MALAPERT A, GUERET C, JUSSIEN N, LANGEVIN A & ROUSSEAU LM. 2008. Two-dimensional pickup and delivery routing problem with loading constraints. In: Proceedings of the First CPAIOR Workshop on Bin Packing and Placement Constraints (BPPC’08). pp. 1-46. 4. Paris, France.
• MARTÍNEZ L & AMAYA CA. 2013. A vehicle routing problem with multi-trips and time windows for circular items. Journal of the Operational Research Society, 64(11): 1630-1643.
• MOR A & SPERANZA MG. 2020. Vehicle routing problems over time: a survey. 4OR, 18: 129- 149.
• MUNOZ AL. 2011. Solución al problema de ruteo de vehículos con restricciones de capacidad y reordenamiento de carga en los sitios de demanda. Master’s thesis. Escuela de Ingenieria, Pontificia Universidad Catolica de Chile. Santiago, Chile.
• NADDEF D & RINALDI G. 2002. Branch-and-cut Algorithms for the Capacitated VRP. In: TOTH P & VIGO D (Eds.), The Vehicle Routing Problem. pp. 53-84. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics.
• NASCIMENTO OX, QUEIROZ TA & JUNQUEIRA L. 2021. Practical constraints in the container loading problem: Comprehensive formulations and exact algorithm. Computers & Operations Research , 128: 105186.
• OLIVEIRA LA, DE LIMA VL, QUEIROZ TA & MIYAZAWA FK. 2021a. Comparing a static equilibrium based method with the support factor for horizontal cargo stability in the container loading problem. Pesquisa Operacional, 41 (e240379): 1-23.
• OLIVEIRA LA, DE LIMA VL, QUEIROZ TA & MIYAZAWA FK. 2021b. The container loading problem with cargo stability: a study on support factors, mechanical equilibrium and grids. Engineering Optimization, 53: 1192-1211.
• QUEIROZ TA & MIYAZAWA FK. 2013. Two-dimensional strip packing problem with load balancing, load bearing and multi-drop constraints. International Journal of Production Economics, 145(2): 511 - 530.
• REINELT G. 1991. TSPLIB - A Traveling Salesman Problem Library. ORSA Journal on Computing, 3(4): 376-384.
• SHEN Y & MURATA T. 2012. Pick-up Scheduling of Two-dimensional Loading in Vehicle Routing Problem by using GA. In: Proceedings of the 2012 International MultiConference of Engineers and Computer Scientists (IMECS), vol. 2. pp. 1532-1537.
• STRODL J, DOERNER KF, TRICOIRE F & HARTL RF. 2010. On Index Structures in Hybrid Metaheuristics for Routing Problems with Hard Feasibility Checks: An Application to the 2- Dimensional Loading Vehicle Routing Problem. In: BLESA MJ, BLUM C, RAIDL G, ROLI A & SAMPELS M (Eds.), Proceedings of the 7th International Workshop on Hybrid Metaheuristics, vol. 6373. pp. 160-173. Springer Berlin Heidelberg.
• SWEENEY PE & PATERNOSTER ER. 1992. Cutting and Packing Problems: A Categorized, Application-Orientated Research Bibliography. Journal of the Operational Research Society , 43(7): 691-706.
• TOTH P & VIGO D. 2001. An Overview of Vehicle Routing Problems. In: The Vehicle Routing Problem. pp. 1-26. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics .
• WEI L, ZHANG Z, ZHANG D & LIM A. 2015. A variable neighborhood search for the capacitated vehicle routing problem with two-dimensional loading constraints. European Journal of Operational Research , 243(3): 798-814.
• WEI L, ZHANG Z, ZHANG D & LIM A. 2018. A simulated annealing algorithm for the capacitated vehicle routing problem with two-dimensional loading constraints. European Journal of Operational Research , 265: 843-859.
• ZACHARIADIS EE, TARANTILIS CD & KIRANOUDIS CT. 2009. A Guided Tabu Search for the Vehicle Routing Problem with two-dimensional loading constraints. European Journal of Operational Research , 195(3): 729-743.
• ZACHARIADIS EE, TARANTILIS CD & KIRANOUDIS CT. 2013. Integrated distribution and loading planning via a compact metaheuristic algorithm. European Journal of Operational Research , 228(1): 56-71.
• ZACHARIADIS EE, TARANTILIS CD & KIRANOUDIS CT. 2016. The Vehicle Routing Problem with Simultaneous Pick-ups and Deliveries and Two-Dimensional Loading Constraints. European Journal of Operational Research , 251(2): 369-386.

• 1
  http://www.or.dei.unibo.it/library

Publication Dates

• Publication in this collection
  02 May 2022
• Date of issue
  2022

History

• Received
  11 Feb 2021
• Accepted
  03 Feb 2022