A Neural Column Generation Approach to the Vehicle Routing Problem with Two-Dimensional Loading and Last-In-First-Out Constraints

1 Introduction

The vehicle routing problem (VRP) is crucial in various logistics applications such as express systems, industrial warehousing, and on-demand delivery [Zong et al., 2021]. The basic form of VRP, or the capacitated vehicle routing problem (CVRP), focuses on optimizing a fleet of vehicles to meet customer demands with the goal of minimizing the total travel costs. An important extension of CVRP is the vehicle routing problem with the two-dimensional loading constraints (2L-CVRP) [Iori et al., 2007]. In 2L-CVRP, the challenge involves transporting rectangular items, each with specific length, width, and weight. The loading constraints in 2LCVRP revolve around two main aspects: the orientation of items, which concerns whether an item can be rotated, and sequential loading, which concerns whether items belonging to successive customers along a route are allowed to be moved [Fuellerer et al., 2009]. The latter is often referred to as the last-in-first-out (LIFO) or rear-loading constraints. This problem is especially relevant in scenarios like transporting furniture or large industrial equipment where stacking items is not feasible due to safety or operational concerns. A real-world example of 2L-CVRP's application is seen at Opein, a company distributing large equipment, highlighting the model's practical importance. The practical value of 2L-CVRP has motivated numerous studies including both exact algorithms and heuristics, and solving new practical variants [Pollaris et al., 2015]. In this paper, we focus on 2L-CVRP with fixed item orientation and LIFO rule.

∗Corresponding author.

The 2L-CVRP encompasses two NP-hard problems: the CVRP and the two-dimensional bin packing problem (2DBPP) stemming from the loading constraints. Given this complexity, the 2L-CVRP is solved approximately in most of studies [Gendreau et al., 2008; Wei et al., 2018]. However, the development of exact solvers is crucial for a deeper understanding of the problem's structure and for evaluating the performance gap between approximate solutions and optimal solutions. To achieve this end, recent studies [Cˆot´e et al., 2020; Zhang et al., 2022a; Zhang et al., 2022b] have proposed several efficient exact algorithms and closed many open instances. The state-of-the-art (SOTA) exact algorithm for the 2L-CVRP [Zhang et al., 2022b] relies on the column generation (CG), which features repeatedly solving a challenging pricing problem (PP) and thus is regarded as one of the main bottlenecks blocking us from solving more open instances.

Recently, the integration of machine learning (ML) with CG has shown promise in solving combinatorial optimization problems more efficiently, while still aiming for optimal solutions [Morabit et al., 2021; Zhang et al., 2022c; Shen et al., 2022; Chi et al., 2022; Morabit et al., 2023; Yuan et al., 2023]. ML's ability to learn from data and make probabilistic predictions offers a potential acceleration in the CG process, particularly in problems where traditional methods are either too slow or impractical.

However, applying ML-based CG methods to the 2LCVRP, especially under the LIFO rule, introduces distinct challenges. Existing ML-based CG algorithms often solve the PP and generate multiple columns by traditional methods, with ML typically applied as a post-processing tool for column selection [Morabit et al., 2021; Chi et al., 2022; Yuan et al., 2023]. Some studies have attempted to bypass solving PP altogether, using ML to directly generate columns

• (a) Pipeline of the state-of-the-art method.

(b) Pipeline of our neural column generation method.

Figure 1: Comparative illustration of the SOTA method and our NCG method for the 2L-CVRP with the LIFO rule.

[Shen et al., 2022]. However, the PP in the context of 2LCVRP, which combines the elementary shortest path problem with resource constraints (ESPPRC) and the 2D-BPP with the LIFO rule, is non-trivial to be addressed through traditional methods or ML models alone. While [Morabit et al., 2023] proposed an ML approach to prune the ESPPRC pricing graph, applying this strategy to the PP of 2L-CVRP is difficult due to the added complexity of 2D-BPP. A related work to ours is [Zhang et al., 2022c], which tackled a variant of 2L-CVRP without the LIFO constraint, utilizing feedforward networks (FF) to accelerate the feasibility checking of columns generated by the pricing algorithm. Our study extends this approach by incorporating the LIFO rule, adding further complexity. We propose a novel ML model that leverages attention mechanism [Vaswani et al., 2017] for homogeneous items within the same customer and recurrence mechanism [Cho et al., 2014] for heterogeneous items across different customers. This model is used to predict the feasibility of columns generated from the ESPPRC, ensuring compliance with the 2D-BPP constraints with the LIFO rule. Our method provides a more efficient alternative to the traditional, timeconsuming feasibility checker as exemplified in [Zhang et al., 2022b], achieving a median acceleration of 29.79% and successfully solving one challenging open instance for the first time. For a visual representation of our algorithmic pipeline, see Figure 1b.

In summary, the main contributions of this paper are:

• We propose a neural column generation (NCG) algorithm that combines the state-of-the-art column generation for the 2L-CVRP with the newly developed ML model.
• Our novel ML model integrates attention and recurrence mechanisms, along with a symmetry-based data augmentation technique. It effectively tackles the 2D-BPP with the LIFO rule, achieving an overall accuracy of around 95%.
• The NCG algorithm, when tested on standard benchmark instances, demonstrates a significant reduction in runtime, with a median decrease of 29.79% compared to the state-of-the-art column generation algorithm.
• Additionally, the NCG algorithm is incorporated into the state-of-the-art branch-and-price-and-cut (BPC) algorithm, successfully solving an open instance for the first time.

2 RelatedWork

In this section, we review the existing methods for 2L-CVRP and recent advances in enhance CG with ML techniques.

Non-learningMethods for 2L-CVRP. The vehicle routing problem with two-dimensional loading constraints has been an active research area due to its practical relevance in logistics. The seminal work by [Iori et al., 2007] introduced a branch-and-cut (B&C) algorithm adapted for 2L-CVRP by introducing infeasible-path constraints. They also propose an exact packing algorithm as a subroutine for the B&C algorithm. The method set a benchmark for the problem, optimally solving instances with up to 35 customers. Subsequent heuristic approaches, such as the tabu search heuristic by [Gendreau et al., 2008] and the guided local search by [Zachariadis et al., 2009], improved solution times and led to several new best solutions. The most effective heuristic for the 2L-CVRP to date involves a simulated annealing heuristic coupled with a local search-based packing algorithm [Wei et al., 2018], which has shown superior solution quality. The 2L-CVRP with unloading sequential constraints has been tackled by [Pinto et al., 2013] through a column-generationbased heuristic. This approach initially omits loading constraints and then applies a heuristic to construct feasible packing solutions. The same authors later developed a more refined branch-and-price (B&P) algorithm [Pinto et al., 2016] that integrates a variable neighborhood search (VNS) algorithm [Pinto et al., 2015] into the pricing problem. Despite these advancements, exact methods still face challenges with larger instances. The B&C algorithm by [Cˆot´e et al., 2020] incorporated advanced exact packing algorithms [Cˆot´e et al., 2014b], achieving a significant breakthrough in benchmark instances. [Zhang et al., 2022a] further enhanced this approach by introducing a heuristic for separating infeasible set inequalities at fractional nodes, closing several open instances and improving dual bounds. Later on, [Zhang et al., 2022b]

propose an exact CG algorithm addressing the loading constraints with a novel data structure L-Trie. The CG algorithm leading to a successful BPC algorithm which closed many open instances for the first time and remains to be the state-ofthe-art exact column generaiton for the 2L-CVRP. For additional studies on 2L-CVRP, please refer to [Wang et al., 2009; Pollaris et al., 2015; Iori et al., 2021].

Neural Column Generation. Recent years researchers have become increasingly interested in ML to accelerate optimization tasks [Bengio et al., 2021], and several learningbased methods have been proposed for specific problems solved by CG. For example, [Morabit et al., 2021] employed a Graph Neural Network (GNN) for column selection in CG as a supervised binary classification task, aiming to mitigate degeneracy in the restricted master problem (RMP). This approach involved representing columns and constraints as a bipartite graph, with the GNN predicting whether to include or exclude each column during CG iterations. Building on this, [Chi et al., 2022; Yuan et al., 2023] modeled CG as a Markov decision process, using reinforcement learning for column selection during CG iterations. These methods, using GNN as Q-function approximator, proved to be more efficient than traditional greedy policies in problems such as the cutting stock problem. All these methods could be viewed as a post-processing step after CG, i.e., applying ML models to select from the generated columns rather than selecting greedily. [Morabit et al., 2023] proposed a supervised MLbased algorithm to prune the pricing network in the form of ESPPRC, alternating between reduced and complete graphs to accelerate the CG process. [Shen et al., 2022] utilized a support vector machine (SVM) for the graph coloring problem, directly generating columns by sampling from the SVM to enhance the CG process. The most related work to ours is [Zhang et al., 2022c], which addressed the vehicle routing problem without the LIFO rule. They utilized FF for heuristic validation of columns generated from the ESPPRC, reducing the dependence on the exact solver [Cˆot´e et al., 2014a]. However, our paper addresses a more complex variant of the problem by including additional loading constraints, specifically the LIFO rule, which significantly increases the complexity of the CG process.

3 Background

In this section, we first present the mathematical formulation of the 2L-CVRP with the LIFO rule as well as the CG process. Then, we introduce the SOTA exact method for solving the 2L-CVRP with the LIFO rule.

3.1 ProblemFormulation

The 2L-CVRP with the LIFO rule is defined on a complete undirected graph G = (V, E), where V = {0, 1, 2, ..., n, n + 1} stands for the set of vertices consisting of customers Vc = {1, 2, ..., n} and the depot 0. Vertex n + 1 represents a copy of the depot. The connections between any pair of vertices are depicted by the edge set E. ∀e ∈ E, ce represents the traveling cost associated with edge e. An alternative representation of an edge (vi, vj) is also used. Set K represents a fleet of homogeneous vehicles which are available at the depot. A loading area characterized by H and W as the length and width, respectively, is a property attached to each vehicle. The loading areas of all the vehicles are the same. Naturally, we have the total area of the loading surface of any vehicle equal to A = H ×W. Each vehicle also has a weight capacity denoted as Q.

As for the customers, ∀i ∈ Vc is characterized by a set Mi of rectangles. Any item m ∈ Mi is marked by width wi,m, length hi,m, and weight qi,m. Let νi and ci represent the total area and the total weights of all the items in customer i. In other words, νi = ∑m∈Mi wi,mhi,m and qi = ∑m∈Mi qi,m. The total number of the items in G is equal to |M|, where M is the union of the item sets of all the customers. The 2L-CVRP calls for planning routes for the fleet such that the demands of the customers are covered while respecting the following constraints:

1. Each customer has to be visited exactly once;
2. The total weight of the items to be delivered by a vehicle cannot exceed Q;
3. The carried items have to be packed in the loading area without collision;
4. Items are not allowed to be rotated over the course of packing;
5. Items delivered to subsequent customers cannot be moved when unloading the items for the current customer (the LIFO rule).

The 2L-CVRP can be formulated as a set partitioning (SP) problem:

where Ω represents the collection of feasible routes; cr represents the total traveling cost of route r; λr is a binary decision variable indicating whether route r is selected as a part of the solution; ai,r is a binary indicator, where ai,r = 1 if vertex i is visited in route r, and ai,r = 0 otherwise.

Equation (1) defines the objective function for the set partitioning problem. Constraint set (2) imposes that there should be exactly |K| routes to be selected as we assume that there is no idle vehicle in the fleet. This has been a convention when it comes to developing exact algorithms for the 2L-CVRP [Iori et al., 2007]. Constraint set (3) calls for that each customer should be visited exactly once. Constraint set (4) defines the domain of the decision variables.

In practical applications, directly solving the SP formulation is infeasible due to the need to enumerate all routes in Ω. Typically, a smaller subset of Ω is selected to create a reduced version of the problem, known as the restricted master problem. Solving the RMP yields a solution, λ∗ , that minimizes the objective value of the reduced formulation. However, λ∗

may not be optimal for the original problem. To potentially improve upon λ∗ , a sub-problem called the pricing problem is solved to identify any routes (or columns) in Ω that could enhance the solution. This iterative process, known as column generation, continues until no further improvements are found. For a detailed discussion on this approach, readers are referred to Chapter 2 in [Desaulniers et al., 2006].

The formulation of the underlying pricing problem is as follows:

where ¯ de is the reduced cost defined as ¯ di,j = ci,j - 12 πi - 12 πj, ∀(i, j) ∈ E. The binary decision variable xe indicates whether edge e in E is used. The set δ(i) denotes edges incident to vertex i, and δ(S) denotes edges connecting vertices inside S with those outside, where S is a subset of V . Σ(S) is the set of all feasible permutations of vertices in S. The route constructed by set S in order σ is represented as (S, σ), with E(S, σ) being its edge set. The dual variables πi and πf are associated with constraints 3 and 2, respectively.

The objective function (constraint set 5) aims to minimize the reduced cost. Constraint set (6) ensures proper degree constraints for vertices, and constraint set (7) addresses subtour elimination [Desrochers and Laporte, 1991]. Constraint sets (8) and (9) ensure that the vehicle's capacity and loading surface area limits are not exceeded. Finally, constraint set (10) imposes restrictions due to the loading constraints.

3.2 State-of-the-Art Pipeline

Figure 1a illustrates the pipeline of the SOTA CG-based algorithm for the 2L-CVRP with the LIFO rule, as proposed by [Zhang et al., 2022b]. The process starts with a restricted master problem which essentially enumerates a subset of set Ω. The next step involves solving the linear relaxation of this problem to obtain a dual solution, which then facilitates the establishment of the PP. The labeling algorithm, enhanced by trie [Brass, 2010], completion bounds, and ng-route relaxation [Baldacci et al., 2011], efficiently prices out improving columns without checking the loading feasibility. The last step involves filtering out infeasible columns with the exact checker [Cˆot´e et al., 2014b] and adding the feasible ones to the restricted master problem, iterating until no feasible column is found. For a more detailed understanding, please refer to [Zhang et al., 2022b].

4 Methodology

In this section, we first introduce the motivation behind our integrated NCG approach as well as the outline of our approach. Then, we elaborate on the architecture of our machine learning model, including two key mechanisms and an important data augmentation technique.

4.1 Motivation and Outline

The current SOTA pipeline employs an exact feasibility checker to evaluate each candidate column generated by the labeling algorithm, necessitating solving the 2D-BPP with LIFO rule, a known NP-hard problem. Our proposed approach, however, introduces a ML model to refine this process by reducing the dependency on solving 2D-BPP. As illustrated in Figure 1b, the ML model is positioned before the feasibility checker, categorizing candidate columns into 'feasible' and 'infeasible'. Feasible columns directly enter the restricted master problem, while infeasible ones are further checked to rectify potential false negatives and thus preventing optimal columns being discarded. To manage false positives, the loading feasibility of variables in the optimal basis is checked during master problem resolution, adding falsepositive cuts to exclude infeasible variables from the feasible solution space.

This approach leverages the insight that columns not part of the optimal basis can bypass exact checking. However, its success depends on the ML model's accuracy, since frequent erroneous predictions can increase computational iterations.

4.2 Machine LearningModel

Note that candidate columns generated by the labelling algorithm might violate the loading constraints. To address this, we frame it as a binary classification problem, developing a ML model to predict the feasibility of each candidate column. Specifically, our model predicts the probability with which each column respects the loading constraints.

Our model's architecture comprises a parallel embedding mechanism and a recursive processing strategy, designed to effectively capture both homogeneous and heterogeneous features of items within each column. For a given input candidate column, represented as a sequence of item sets [{xi,m}1≤m≤|Mi | ]1≤i≤n , each item xi,m is represented by normalized dimensions [ wi,m W , hi,m H ]. We employ an attention mechanism [Vaswani et al., 2017] for each customer i to integrate features of homogeneous items (those belonging to the same customer). To handle heterogeneous features (items across different customers), a GRU [Cho et al., 2014] is utilized to process the sequence, incorporating the order information which is crucial in adhering to the LIFO constraints. Moreover, a symmetry-based data augmentation technique is employed to incorporate permutation invariance into the model. For a visual representation of our model's architecture, please refer to Figure 2.

Item-Level Attention Mechanism

In the 2L-CVRP context, items belonging to the same customer are homogeneous and exempt from LIFO constraints, making the Multi-Head Attention (MHA) mechanism a suitable choice [Kool et al., 2018; Kwon et al., 2020]. Our model uses MHA to capture the shared features of items within the same customer. The initial embedding for the i-th customer's m-th item is formulated as h0i,m = W0xi,m + b0, where W0 is the initial projection matrix and b0 is the bias vector. The architecture includes skip-connections [He et al., 2015], feed-forward networks, and layer normalization (LN) [Ba et al., 2016] in each sublayer. The embedding for the i-th customer's m-th item is iteratively updated in the l-th layer, as depicted in the following equations:

The MHA mechanism at the core of our attention layer is defined as follows:

where j = 1, 2, ...,H and dk = dh/H. Here, H is the number of attention heads, dh is the dimension of the item embedding, and Qji,m,Kj i,m, V j i,m represent the query, key, and value vectors, respectively. WO is the projection matrix utilized to project the final MHA output. The final embedding of each item after L layers is denoted by hi,m = hLi,m.

Customer-Level Recurrence Mechanism

In the context of the 2L-CVRP with the LIFO rule, a sequential order relationship exists among customers, indicating that items from different customers are inherently heterogeneous. This order relationship dictates that for customer i visited before customer i + 1, items of customer i (denoted as xi,1, xi,2, . . . , xi,|Mi|) must be loaded into the vehicle after items of customer i + 1 (xi+1,1, xi+1,2, . . . , xi+1,|Mi+1|).

To model the recursive relationships among items of different customers, we use GRU as follows:

where ˜ht denotes the hidden state at time step t and ˜h0 = 0. In our approach, customers are processed sequentially, inputting one item ht per time step t into the GRU, with the total number of time steps T equaling the total items, T = ∑ni=1 |Mi|. Upon processing the last item of the last customer, the final state ˜hT of the GRU is transformed into a probability via a FF network and a sigmoid function:

Data Augmentation with Permutation Invariance

As discussed in Section 4.2, items belonging to the same customer are homogeneous and are not constrained by any specific order. This characteristic allows for the application of permutation invariance as a data augmentation strategy, reflecting the symmetry of combinatorial problems [Kwon et al., 2020; Kim et al., 2022]. By permuting the items of each customer, we can generate new, equivalent instances, expanding the training dataset and mitigating early overfitting.

Specifically, after applying the item-level multi-head attention mechanism to all items of each customer i, we perform a permutation πi, shuffling the sequence (1, 2, ..., |Mi|) to produce varied item orderings. This process can be represented mathematically as:

This permutation is applied independently to each customer's set of items before they are processed by the GRU.

5 Experiments

5.1 Experimental Settings

Data Sets. Our machine learning model is trained on a unique set of packing problem instances, distinct from the benchmark instances used in [Iori et al., 2007], ensuring an unbiased performance evaluation. The training dataset is derived by solving the 2L-CVRP problems as well as 2LVRPTW (the 2L-CVRP with time window constraints). In particular, we stored all the packing instances in the course of solving the routing problems. There are in total two batches of the training samples: one from 2L-VRPTW instances and another from 2L-CVRP instances. The VRP instances are randomly generated as per the description in [Zhang et al., 2022c], ensuring a diverse distribution distinct from benchmark instances. The benchmark instances are categorized into five families, as defined in [Iori et al., 2007]. For details, please refer to the Supplementary Material. We exclude family 1, as it is equivalent to the one-dimensional loading scenario, and family 5, as highlighted by [Zhang et al., 2022b], due to its focus on very small items and a specialized column generation variant, making it less relevant for testing our algorithm in standard 2L-CVRP contexts. In Section 5, terms 'family' and 'packing class (PC)' are used interchangeably.

Evaluation Metrics. For the 2L-CVRP instances, our experiment compares the NCG approach to the SOTA CG method by solving the linear relaxation of Formulation 1 - 4. Given that both algorithms are exact, their optimal objective values are expected to match. Our primary interest lies in the difference in computational time. Let TNCG represent the walltime for NCG and TSOTA for SOTA CG. The percentage gap, calculated as TSOTA-TNCG TNCG × 100%, serves to quantify the time savings and performance increase, with higher values indicating greater efficiency gains for NCG. Furthermore, we track the number of column generation iterations to evaluate the impact of the ML-induced false-positive cuts on the frequency of column generation.

Hardware. During the training phase of theMLmodel, experiments were conducted on an AMD EPYC 7V13 64-Core CPU @ 2.45GHz with an NVIDIA A100 GPU. Post training, the ML model was serialized and integrated into the NCG algorithm using Torch C++. To ensure a fair comparison, both the NCG algorithm and the SOTA algorithm [Zhang et al., 2022b] were implemented in C++ and evaluated on an Intel i5-10600KF processor @ 4.10 GHz.

5.2 Main Results

Acceleration of SOTA Algorithm

We evaluated the NCG against the SOTA on benchmark instances with up to 50 customers from families 2 - 4. As depicted in Figure 3, NCG demonstrates a significant advantage over the SOTA algorithm, with a median performance improvement of 29.79% across all tested instances. Notably, the performance gains are more pronounced in PC3 and PC4, with median percentage gaps of 44.97% and 99.22%, respectively. This improvement correlates with the higher accuracy of our ML model in PC3 and PC4, as indicated in Table 2, suggesting fewer iterations required to rectify false predictions. In contrast, the least improvement occurs with PC2, with with a median gain of 9.38%, as the ML model reaches the worst accuracy in PC2.

Figure 4 also presents an intriguing observation: an increase in the number of CG iterations due to the incorporation of false positive cuts. Interestingly, for each instance family, the increase in CG iterations inversely correlates with performance improvement. Despite this, the additional time incurred by these extra iterations is more than compensated for by the time savings from the ML model. This highlights the effectiveness of our NCG approach in enhancing overall computational efficiency despite the potential for additional CG iterations. For detailed results on each instance, please refer to the Supplementary Material.

Solving Open Instance

Our NCG algorithm was integrated into the BPC framework developed in [Zhang et al., 2022b]. This integrated approach was tested on benchmark instances, leading to the successful resolution of an open instance. As shown in Table 1, for the instance 2304, our NCG-based BPC method obtained an optimal objective value (Optimal Obj) of 1068, while the existing SOTA BPC algorithm was unable to reach an optimal solution, yielding only a lower bound of 1066.78 and a primal bound of 1069, despite running for approximately 1.4 times longer than our approach. Furthermore, the SOTA BPC algorithm suffered from memory overload, whereas our NCG-based BPC approach did not, indicating its more efficient memory usage. This efficiency is attributed to the reduced number of routes that needed to be checked by the exact packing algorithm, resulting in lower memory demands for the BPC algorithm.

Figure 5: Comparison of validation losses over epochs for the baseline and augmented models.

5.3 Ablation Study

To validate the contributions of the attention and recurrence mechanisms and the data augmentation technique described in Section 4.2, we performed an ablation study. We compared the following configurations:

1. The full NCG model that integrates both attention and recurrence mechanisms with data augmentation.
2. A variant without data augmentation.
3. A model where the attention mechanism is replaced by a multilayer perceptron (MLP).
4. A model where the recurrence mechanism is replaced by a transformer encoder with sinusoidal positional encoding [Vaswani et al., 2017].

As shown in Table 2, the complete NCG model outperforms the other variants in terms of true positive rate (TPR) and true negative rate (TNR). These results suggest that the item-level attention mechanism contributes the most, followed by the customer-level recurrence mechanism, and then the data augmentation technique. Despite the transformer encoder's capability to process sequential information, it underperforms compared to the GRU-based model, highlighting the importance of explicit modeling of sequential relationships in the 2L-CVRP with the LIFO rule.

Figure 5 shows that, compared to the baseline model which overfits sooner, the model with permutation invariance maintains a lower validation loss for a longer period, confirming the effectiveness of the data augmentation strategy in enhancing generalization capacity.

6 Conclusion

In this work, we introduced a ML-based exact algorithm to solve the 2L-CVRP with the LIFO rule. Our approach integrates neural column generation into the SOTA pipeline, formulating the feasibility checking of candidate columns as a binary classification problem. To preserve solution optimality, we implemented post-processing steps for handling false negatives and positives. The ML model in our framework utilizes an attention mechanism for capturing homogeneous item features and a recurrence mechanism for heterogeneous features. Additionally, we employed data augmentation exploiting the problem's symmetry.

Experimental evaluations on benchmark instances demonstrate that our NCG method notably accelerates the SOTA algorithm by a median of 29.79% across various problem instances. Significantly, our approach also successfully solves an open instance, marking a substantial contribution to the field. This achievement highlights the potential of integrating machine learning techniques into traditional optimization problems for enhanced performance and efficiency.