Arcflow Formulations and Constraint Generation Frameworks for the Two Bar Charts Packing Problem

1. Introduction

Given a set of weighted items and an unlimited number of identical capacitated bins, the bin packing problem (BPP) requires packing all the items into the minimum number of bins. The BPP was introduced several decades ago by Kantorovich (1960) and quickly became one of the most studied combinatorial optimization problems. Indeed, its relatively simple structure made it the perfect playground for the research community to develop and test the most effective solution methods, including branch-and-bound algorithms (Martello and Toth 1990), metaheuristics (Falkenauer 1996), branch-and-price algorithms (Vance 1998), pseudo-polynomial formulations (Valério de Carvalho 1999), constraint programming (Shaw 2004), and reinforcement learning (Cai et al. 2019). Besides theoretical motivations, the BPP is also studied for its applications in the cutting industry (Stadtler 1990) and for its key role when solving more complex combinatorial optimization problems such as the vehicle routing problem (Lysgaard et al. 2004) and the two-dimensional BPP (Côté et al. 2021).

Along the years, many extensions of the BPP were introduced. We mention in particular (i) the BPP with conflicts (Muritiba et al. 2010) in which pairs of items cannot be packed in the same bin, (ii) the contiguous BPP (Martello et al. 2003) in which pairs of identical items must be packed in consecutive bins, and (iii) the vector packing problem (Brandão and Pedroso 2016) in which the item weights and bin capacities are multidimensional vectors and a capacity constraint is imposed on every dimension.

Erzin et al. (2021b) introduces one of the latest BPP extensions, which is known as the bar charts packing problem (BCPP). In the BCPP, we are given a set of ordered lists of weighted items (the bar charts) that must be packed using the minimum number of identical capacitated bins. Besides the usual capacity constraints, the BCPP also requires the items of a given bar chart to be packed in consecutive bins. The BCPP originates from a class of project scheduling problems in which one must determine the starting date (the bin) of a set of projects (the bar charts) so that the total makespan (the number of bins used) is minimized. Every project runs for a certain number of periods (one item in the bar chart corresponds to one period) and consumes a given amount of resource that may vary from one period to the other (the weights of the items in the bar chart). The total amount of resource available per period (the bin capacity) is limited, and once a project is started, it cannot be interrupted (contiguity constraints). In this paper, we are interested in the special case of the problem in which every bar chart is composed of exactly two items: the 2-BCPP. We show in Figure 1 a 2-BCPP instance (on the left part) together with a feasible packing that uses seven bins (on the right part).

Even though the 2-BCPP has a limited set of practical applications, it is still an interesting problem to study because it lies at the intersection of several well-known BPP extensions. It is, therefore, reasonable to assume that some of the most effective techniques aimed at solving these extensions can be adapted to solve the 2-BCPP. However, the available literature on the topic is mostly dedicated to approximation algorithms and their worst case performance ratio (Erzin et al. 2021a, b, c), whereas lower bounding procedures and exact approaches have not been thoroughly investigated.

1.1. Our Contribution

We provide an overview of the most relevant packing problems and cluster these problems into three classes ${\mathcal{C}}_l,\hspace{0.17em}{\mathcal{C}}_u$, and ${\mathcal{C}}_e$. We show that

• Every problem $p_l\in {\mathcal{C}}_l$ is a relaxation of the 2-BCPP, and thus, algorithms aimed at solving p_l can also be used to derive valid lower bounds for our problem.

• The 2-BCPP is a relaxation of every problem $p_u\in {\mathcal{C}}_u$, and thus, algorithms aimed at solving p_u can also be used to derive valid upper bounds for our problem.

• The 2-BCPP is a special case of every member $p_e\in {\mathcal{C}}_e$, and thus, algorithms aimed at solving p_e can also be used to solve our problem, but the current state-of-the-art approaches are not well-suited for large-sized 2-BCPP instances.

We then introduce Eulerian-flow and link-flow, two new integer linear programming (ILP) models based on the well-known arcflow formulation introduced by Valério de Carvalho (1999) for the BPP. As both models use a pseudo-polynomial number of variables and an exponential number of constraints, we propose a constraint generation framework to solve each of these models to optimality.

Subsequently, we empirically evaluate the bounding procedures, the newly introduced decomposition approaches, and two methods proposed by Erzin et al. (2021b) (an ILP model and a heuristic) on a large set of random instances, and we identify instance features for which the tested approaches are particularly effective. For example, the ILP model of Erzin et al. (2021b) works very well when the maximum item weight is below a 10th of the bin capacity, whereas Eulerian-flow and link-flow are several orders of magnitude faster otherwise. We also design two new classes of hard instances for the 2-BCPP: TRIPLETS, in which the main difficulty is to find an optimal solution (i.e., find the right upper bound), and DONUTS, in which the main difficulty is to prove that the incumbent solution is optimal (i.e., find the right lower bound).

Besides making significant progress in the 2-BCPP field, the ideas we introduce are also relevant for arcflow formulations in general. In particular, to the best of our knowledge, it is the first time that an arcflow model is proposed for a BPP with a form of precedence constraints, indicating that similar approaches can be used for a larger set of packing problems. It is also the first time that an arcflow model is solved through constraint generation. This could lead to promising developments in decomposition approaches in which the master problem is a BPP (see e.g., the recent work of Côté et al. 2021 to solve the two-dimensional BPP). Indeed, arcflow formulations can now be considered to solve the master problem instead of the frequently used textbook (or “Kantorovich”) ILP model.

1.2. Related Work

The 2-BCPP is formally introduced by Erzin et al. (2021b). The authors propose several approximation algorithms with proven worst case performance ratios and compare the results obtained by their approaches to those obtained when solving a descriptive ILP model with the commercial solver CPLEX. They empirically show that the solver is more effective for instances with up to 100 bar charts, whereas the heuristics display better empirical performance otherwise. The ILP model considered by Erzin et al. (2021b) uses one variable per combination of bar chart and bin and is a direct extension of the textbook BPP model. Later on, Erzin et al. (2021a) perform additional experiments in which they use 2-BCPP instances with known optimal solutions to compute more accurate optimality gaps for the tested approaches. In a subsequent work, Erzin et al. (2021c) propose two additional approximation algorithms for special versions of the 2-BCPP and derive the associated worst case performance ratios. They study the case in which every bar chart contains at least one item with weight above half of the bin capacity and the case in which, in addition to that condition, the first item of every bar chart is at least as large as the second item.

Even though the 2-BCPP literature is relatively scarce, there are a few papers focusing on closely related packing problems that are particularly relevant for our study. In the following, we provide an overview of related problems together with a brief literature review and group them into three classes: (i) ${\mathcal{C}}_l$, the problems that are relaxations of the 2-BCPP; (ii) ${\mathcal{C}}_u$, the problems for which the 2-BCPP is a relaxation; and (iii) ${\mathcal{C}}_e$, the problems for which the 2-BCPP is a special case.

1.2.1. Problem Class ${\mathcal{C}}_l$.

Every packing problem that is a valid relaxation of the 2-BCPP belongs to ${\mathcal{C}}_l$. In other words, transforming a 2-BCPP instance I into an instance $I^{\prime }$ of problem $p (p\in {\mathcal{C}}_l)$ and solving $I^{\prime }$ to optimality provides a valid lower bound for I. If we remove the constraints forcing the items of a given bar chart to be packed in consecutive bins, we obtain the BPP. The version of the BPP in which duplicate items are merged into a unique item with integer demand is called the cutting stock problem (CSP). A survey on heuristics and approximation algorithms for the BPP and the CSP is proposed by Coffman et al. (2013), whereas a survey focusing on exact approaches and mathematical models is presented by Delorme et al. (2016). Among the most effective approaches proposed in the literature, we mention branch-and-price algorithms (see Pessoa et al. 2020, Wei et al. 2020b) and pseudo-polynomial formulations (see Valério de Carvalho 1999, Brandão and Pedroso 2016, Delorme and Iori 2020). If we additionally forbid the two items of any bar chart to be packed in the same bin, we obtain the BPP with conflicts. Branch-and-price algorithms (see Sadykov and Vanderbeck 2013) and pseudo-polynomial formulations (see Brandão and Pedroso 2016) are also shown to be effective approaches for solving the BPP with conflicts. Finally, if we force the first item of a given bar chart to always be packed “before” the second item (i.e., in a bin with a lower index), we obtain the BPP with precedence constraints. So far, this problem is mostly tackled by branch-and-bound algorithms (see Dell’Amico et al. 2012, Pereira 2016) and branch-and-price algorithms (see Letelier et al. 2022).

1.2.2. Problem Class ${\mathcal{C}}_u$.

Every packing problem for which the 2-BCPP is a valid relaxation belongs to ${\mathcal{C}}_u$. In other words, transforming a 2-BCPP instance I into an instance $I^{\prime }$ of problem $p (p\in {\mathcal{C}}_u)$ and solving $I^{\prime }$ to optimality provides a valid upper bound for I. If we force the first item of every bar chart to be packed in a bin with an odd index (i.e., if the first item of every bar chart can only be packed with other first items), we obtain the two-dimensional vector packing problem (2D-VPP). This problem is extensively studied in the literature and can be efficiently solved through tailored branch-and-price algorithms (see Wei et al. 2020a) and pseudo-polynomial formulations (see Brandão and Pedroso 2016). If, instead, we merge the two items of every bar chart into a two-dimensional item (i.e., if we force the bottom left corner of the two items in a given bar chart to have the same y-coordinate), we obtain a version of the strip packing problem in which the items can either be rectangles or “L-shaped” hexagons. Algorithms aimed at solving the classic strip packing problem (e.g., the Benders’ decomposition approaches from Côté et al. 2014, Delorme et al. 2017) can be extended to also take into account L-shaped hexagons.

1.2.3. Problem Class ${\mathcal{C}}_e$.

Every packing problem for which the 2-BCPP is a special case belongs to ${\mathcal{C}}_e$. In other words, transforming a 2-BCPP instance I into an instance $I^{\prime }$ of problem $p (p\in {\mathcal{C}}_e)$ and solving $I^{\prime }$ to optimality provides an optimal solution for I. The 2-BCPP can be seen as a special case of the contiguous BPP, which is introduced by Martello et al. (2003) as a relaxation of the strip packing problem (SPP). In that relaxation, every two-dimensional item is cut into unit-height slices, and the objective is to find a packing that uses the minimum number of bins such that the slices originating from the same item are packed into contiguous bins. Even though the items involved in a contiguity constraint have the same weight in the SPP relaxation, one can also consider a version of the contiguous BPP in which this is not the case. In the literature, the contiguous BPP is tackled with ad hoc branch-and-bound algorithms (see Martello et al. 2003) and ILP models solved through commercial solvers (see Côté et al. 2014, Delorme et al. 2017). These ILP models can easily be extended to the 2-BCPP, but Erzin et al. (2021b) show that the resulting approaches are only able to solve small-sized instances to optimality.

The 2-BCPP is also a special case of the BPP with time lags recently introduced by Letelier et al. (2022). In that extension, two items involved in a time lag constraint have a minimum and a maximum time lag, which is defined as the difference between the index of the bins in which the two items are packed. In the version studied by Letelier et al. (2022), the minimum and maximum time lags can vary from one pair of items to the other and are not necessarily bounded. For example, one pair of items can be constrained to have a time lag in the range $[-\infty ;-3]$, whereas another pair can be constrained to have a time lag in the range $[-1;15]$. The 2-BCPP is, thus, a special case of the BPP with time lags in which every pair of items in a given bar chart is forced to have a time lag in the range $[1;1]$. Letelier et al. (2022) propose two ILP models for the BPP with time lags: one that can be seen as an extension of the model proposed by Erzin et al. (2021b) and one that uses an exponential number of variables where each variable decides whether a given pattern is used in a bin with a given index. The latter model is solved with a branch-and-price algorithm and can handle instances with up to a few hundred items. Even though these approaches can also be used to solve the 2-BCPP, we do not expect them to work very well in practice because they need one set of variables per bin index, which quickly becomes impractical for large-sized instances. In this paper, we introduce two new ILP models that make use of the special structure of the 2-BCPP to avoid the necessity to take the bin index into account.

1.3. Layout

The remainder of this paper is structured as follows. In Section 2, we formally define the 2-BCPP, and we discuss the ILP model proposed by Erzin et al. (2021b). In Section 3, we briefly introduce existing arcflow formulations for the CSP and the 2D-VPP and outline how valid lower and upper bounds for the 2-BCPP can be obtained by solving these models. In Section 4, we introduce two new ILP models for our problem. Both of these models are based on the arcflow formulation, have an exponential number of constraints, and can be solved within a constraint generation framework. In Section 5, we present two innovative ways to construct difficult instances for the problem, and we report the findings of a large number of computational experiments. Finally, conclusions are drawn in Section 6.

2. Problem Definition

In the 2-BCPP, we are given an unlimited number of bins with capacity c and a set of n bar charts, each of them composed of two bars (or items—the terms are interchangeable) that must be packed in consecutive bins: the first bar whose weight is denoted by $w_{i1}$ and the second bar whose weight is denoted by $w_{i2}$. We call $w_i=(w_{i1},w_{i2})$ the weight vector of bar chart i. For feasibility purposes, we assume that $w_{i1},w_{i2}\le c$ for every bar chart $i=1,\dots ,n$. We also assume that $w_{i1},\hspace{0.17em}{w}_{i2}$ and c are strictly positive integers and that every bar chart $i=1,\dots ,n$ has an integer demand q_i. Note that, in the original 2-BCPP definition from Erzin et al. (2021b), every bar chart had demand one and bar weights were positive real numbers. However, an instance of such a problem can be transformed to fit our definition by merging duplicate bar charts, computing the associated demands, and multiplying the bar weights and the capacity by a suitable coefficient.

Given a 2-BCPP solution, we say that bar chart i is packed into bin j if the first item of bar chart i is packed into bin j (and, thus, if the second item of bar chart i is packed into bin j + 1). In addition, we say that a solution (or a packing) is feasible if the demand, the contiguity, and the capacity constraints are satisfied. We call z the length of a feasible packing, which is defined as the number of bins in the solution containing at least one item. For a feasible solution of length z, it is assumed that at least one item is packed in every bin $j=1,2,\dots ,z$. If this is not the case, one can simply delete the empty intermediate bins and shift the solution. Finally, we say that a feasible packing with length $\overline{z}$ is optimal if there exists no solution with length z such that $z<\overline{z}$.

3. Arcflow-Based Lower and Upper Bounding Procedures

For many packing problems, researchers show that arcflow formulations provide much better results on average than “Kantorovich-like” models (see Delorme et al. 2016, Martinovic and Scheithauer 2016, Dell’Amico et al. 2019). The name “arcflow” originates from the fact that the variables in the model represent flows on individual arcs in a directed graph (see de Lima et al. 2022 for a recent survey on the topic). Arcflow formulations are known to have a very tight LP-relaxation, but this comes at the cost of not being polynomial in the input size. Indeed, both the number of variables and constraints involved in such models are usually polynomial in the capacity, whereas the latter is coded as $O({\log }_{2}c)$ bits in the input, which is why arcflow formulations are said to be pseudo-polynomial. There is a growing interest in the research community for arcflow formulations as (i) they are usually easy to implement as they do not involve advanced optimization techniques and (ii) they result in models that can often be solved to optimality with a state-of-the-art ILP solver, especially for instances with a reasonable capacity. In summary, arcflow models are seen as a good compromise between simplicity and efficiency.

Typically, an approach using an arcflow model is composed of three phases: (i) the graph-construction phase, in which the instance I of the original problem is transformed into an instance $I^{\prime }$ of a graph problem; (ii) the solving phase, in which the instance $I^{\prime }$ is solved (e.g., with an ILP solver); and (iii) the mapping phase, in which the solution of instance $I^{\prime }$ is mapped into a solution of instance I. In the following, we briefly introduce arcflow formulations for the CSP and 2D-VPP and explain how they can be used to derive valid lower and upper bounds for the 2-BCPP.

3.1. Arcflow-Based Lower Bounding Procedure

The idea behind the arcflow formulation for the CSP is introduced by Wolsey (1977) and popularized by Valério de Carvalho (1999). The CSP is transformed into a graph problem in which the objective is to minimize the flow going from a source node to a target node under demand and flow conservation constraints. In this section, n denotes the number of items in a CSP instance, w_i is the weight of item i $(i=1,\dots ,n)$, and q_i is the demand of the item.

In the graph-construction phase, we consider a graph $G=(V,A)$, where vertex set $V=\{0,1,\dots ,c\}$ and arc set $A=\{(d,e):0\le d<e\le c \text{and} e-d=w_i$ for $i=1,\dots ,n\}\cup \{(d,d+1):0\le d<c\}$. Note that symmetry-reduction techniques can be used to reduce the size of V and A (see Valério de Carvalho 1999, Côté and Iori 2018).

In the solving phase, by introducing integer decision variables f_de that indicate the number of times arc $(d,e)\in A$ is selected, the arcflow model for the CSP can be defined as follows:

The objective function (9) minimizes the flow going from zero to c. Constraints (10) ensure that z units of flow both leave node 0 and enter node c, and enforce flow conservation in the other nodes. Constraints (11) make sure that the demand is satisfied, whereas Constraints (12) limit the variable domains.

In the mapping phase, a flow decomposition algorithm is used on the optimal solution of Model (9)–(12) to obtain a set of $0-c$ paths with unitary flow corresponding to an optimal CSP solution.

3.2. Arcflow-Based Upper Bounding Procedure

The arcflow formulation for the 2D-VPP can be attributed to Brandão and Pedroso (2016), who extend the CSP version of the model to take into account additional dimensions. The 2D-VPP is again transformed into a graph problem in which the objective is to minimize the flow leaving a source node under demand and (a form of) flow conservation constraints. In this section, n designs the number of items in a 2D-VPP instance, (c₁, c₂) is the capacity of the bin, $w_i=(w_{i1},w_{i2})$ is the weight vector of two-dimensional item i $(i=1,\dots ,n)$, and q_i is the demand of the item.

In the graph-construction phase, we consider a graph $G=(V,A)$, where vertex set $V=\{(v_1,v_2):0\le v_1\le c_1,0\le v_2\le c_2\}$ and arc set $A=\{((d_1,d_2),(e_1,e_2)):0\le d_1<e_1\le c_1 \text{and} 0\le d_2<e_2\le c_2 \text{and} (e_1-d_1,e_2-d_2)=w_i$ for $i=1,\dots ,n\}$. Note that symmetry-reduction techniques can also be used to reduce the size of V and A (see Brandão and Pedroso 2016). As it is relevant for the models introduced in the next section, we provide in Algorithm 1 in Section A of the online supplement a pseudo-code to construct graph G with a reduced vertex set V and a reduced arc set A.

In the solving phase, by introducing integer decision variables $f_{(d_1,d_2),(e_1,e_2)}$ that indicate the number of times arc $((d_1,d_2),(e_1,e_2))\in A$ is selected, the arcflow model for the 2D-VPP can be defined as follows:

The objective function (13) minimizes the flow going out of node (0, 0). Constraints (14) ensure that the flow entering any node (v₁, v₂) that is not (0, 0) is never below the flow leaving that node. Constraints (15) make sure that the demand is satisfied, whereas Constraints (16) limit the variable domains.

In the mapping phase, a flow decomposition algorithm is used on the optimal solution of Model (13)–(16) to obtain a set of $(0,0)-(v_1,v_2)$ paths $((v_1,v_2)\in V)$ with unitary flow corresponding to an optimal 2D-VPP solution.

4. New Arcflow Formulations for the 2-BCPP

Even if the arcflow formulation was originally proposed for the CSP, it has been successfully extended to many other relevant packing problems in the last two decades (see Macedo et al. 2010, Nesello et al. 2018, Delorme et al. 2021). In this section, we first introduce Eulerian-flow, a nontrivial arcflow extension for the 2-BCPP in which a feasible solution is represented by a single cycle in a graph (a Eulerian cycle). We observe that, whereas forcing the arcs selected in a solution to form cycles is relatively easy, making sure that these cycles are all connected (i.e., that there are no subtours) is much harder and can only be achieved by using an exponential number of constraints. As the model size quickly becomes too large to be solved by an ILP solver, we propose a constraint generation framework that only adds the subtour elimination constraints that are necessary to find an optimal solution. We then introduce link-flow, a model inspired by Eulerian-flow in which the number of variables is dramatically reduced at the cost of an extra set of constraints.

4.1. Eulerian-Flow

Our first model builds upon the arcflow formulation for the 2D-VPP described in the previous section in which a feasible packing is represented as a set of paths starting from node (0, 0). In that formulation, one can easily modify the arc set definition to incorporate feedback arcs $((v_1,v_2),(0,0))$ for every $(v_1,v_2)\in V$ to represent the end of a bin and updating flow conservation constraints (14) to

A feasible solution f now consists of a set of arcs that can be split into a number of cycles in which each cycle contains node (0, 0). Let us define the directed multigraph $G_f=(V_f,A_f)$, also called the solution graph hereafter, in which V_f contains the set of vertices visited in solution f and A_f contains $f_{(d_1,d_2),(e_1,e_2)}$ copies of arc $((d_1,d_2),(e_1,e_2))$ for every $((d_1,d_2),(e_1,e_2))\in A$. According to Euler’s theorem (developed to solve the Königsberg bridge problem in 1736), a directed graph contains a Eulerian cycle if and only if (i) the in-degree of any vertex is equal to its out-degree and (ii) the graph is strongly connected. Solution graph G_f contains a Eulerian cycle as condition (i) is satisfied thanks to the updated flow conservation constraints (17), and condition (ii) is satisfied because the arcs in A_f can be split into a set of cycles such that node (0, 0) is included in every cycle.

This updated formulation is not yet suitable for the 2-BCPP as it fills the bins by pairs, which means that it does not allow the first item of a bar chart to be packed with the second item of another bar chart. To correct this issue, we introduce transition arcs $((v_1,v_2),(v_2,0))$ for every $(v_1,v_2)\in V$ to represent the end of a 2-BCPP bin. Intuitively, one can see a transition arc as half of a feedback arc. In fact, using two successive transition arcs $((v_1,v_2),(v_2,0))$ and $((v_2,0),(0,0))$ is the same as using one feedback arc $((v_1,v_2),(0,0))$. We point out that, after introducing transition arcs, the arcs in A_f can still be split into a set of cycles, but node (0, 0) is not necessarily included in each of these cycles anymore. This means that an additional set of constraints is required to make sure that the solution graph G_f is strongly connected (and, thus, that it contains a Eulerian cycle).

In the graph-construction phase, we consider a graph $G=(V,A)$, where vertex set $V=\{(v_1,v_2):0\le v_1\le c,0\le v_2\le c\}$ and arc set $A=\{((d_1,d_2),(e_1,e_2)):0\le d_1<e_1\le c \text{and} 0\le d_2<e_2\le c \text{and} (e_1-d_1,e_2-d_2)=(w_{i1},w_{i2})$ for $i=1,\dots ,n\}\cup \{((d_1,d_2),(d_2,0)):0\le d_1\le c \text{and} 0\le d_2\le c\}$. As the number of vertices and arcs grows quickly with c, we detail in the following how graph G can be constructed with a reduced vertex set V and a reduced arc set A. Note that, to avoid any confusion, we temporarily use different notation for the initial graph $G^{\alpha }$, the intermediate graph $G^{\beta }$, and the final graph $G^{\gamma }$.

1. Order the bar charts (e.g., from the largest bar chart to the smallest).

2. Construct graph $G^{\alpha }=(V^{\alpha },A^{\alpha })$ following Algorithm 1 in Online Section A.

3. Save in set V^h every distinct v₂ value ($v_2\ne 0$) such that $(v_1,v_2)\in V^{\alpha }$ for some v₁.

4. Construct graph $G^{\beta }=(V^{\beta },A^{\beta })$ following a modified version of Algorithm 1 in which set V is initialized to V^h instead of $\{(0,0)\}$ and create $V^{\gamma }=V^{\alpha }\cup V^{\beta }$ and $A^{\gamma }=A^{\alpha }\cup A^{\beta }$.

5. For every vertex $(v_1,v_2)\in V^{\gamma }\setminus \{(0,0)\}$, add transition arc $((v_1,v_2),(v_2,0))$ to $A^{\gamma }$.

We point out that it is not necessary to update set V^h after the fourth step because every arc $((d_1,d_2),(e_1,e_2))$ in $A^{\beta }$ has, by construction, its counterpart $((d_1-v,d_2),(e_1-v,e_2))$ in $A^{\alpha }$, where $v\in V^h$. In fact, we can use this fact to generate $V^{\gamma }$ and $A^{\gamma }$ more efficiently by duplicating $A^{\alpha }$ and “shifting” the first dimension of both the tail and the head of every arc by the values in V^h as described in Algorithm 2 in Online Section A. In the next phases, we refer to graph $G=(V,A)$ as the final graph $G^{\gamma }$ built in this phase.

In the solving phase, for the sake of compactness, we now sometimes refer to an arc as a instead of $((d_1,d_2),(e_1,e_2))$. For every vertex $(v_1,v_2)\in V$, we also define index sets ${\delta }^{+}(v_1,v_2)$ (respectively, ${\delta }^{-}(v_1,v_2)$) as the set of arcs leaving (entering) vertex (v₁, v₂). For every $i=1,\dots ,n$, we denote by A_i the set of arcs representing the packing of a bar chart of type i (i.e., the set of arcs $((d_1,d_2),(e_1,e_2))\in A$ such that $(e_1-d_1,e_2-d_2)=(w_{i1},w_{i2})$). We call $A_1\cup A_2\cup \cdots \cup A_n$ the set of bar chart arcs. We also denote by A₀ the set of transition arcs. By introducing integer decision variables f_a that indicate the number of times arc $a\in A$ is selected, the Eulerian-flow model can then be defined as follows:

The objective function (18) minimizes the number of transition arcs used (i.e., the length of the packing). Constraints (19) ensure flow conservation, whereas Constraint (20) makes sure that there is at least one unit of flow leaving node (0, 0). Furthermore, Constraints (21) enforce demand satisfaction, and Constraints (22) eliminate subtours so that the solution graph G_f is strongly connected. More specifically, big-M constraints (22) ensure that, for every potential subtour (i.e., every subset of nodes $\tilde{V}$ that does not contain (0, 0)), there can only be a flow circulating in the subtour (i.e., in an arc whose tail and head are both in $\tilde{V}$) if there is also a flow leaving the subtour (i.e., passing through an arc whose tail is in $\tilde{V}$ and whose head is not in $\tilde{V}$). A suitable big-M value for Constraints (22) is ${\sum }_{i=1}^n{q}_i+K$ as at most ${\sum }_{i=1}^n{q}_i$ units of flow circulate in bar chart arcs and at most K units of flow circulate in transition arcs in total (we recall that K is an upper bound on the minimum number of bins). Note that we can decrease big-M values by making them constraint-specific. In that case, for the constraint associated with a given subset of nodes $\tilde{V}$, we can remove the q_i values of the bar charts for which there is no arc $a\in A_i$ with both tail and head in $\tilde{V}$. Model (18)–(23) uses $O(n{c}^2)$ variables and $O(2^{(c^2)})$ constraints.

In the mapping phase, after an optimal solution f is found for Model (18)–(23), we create the solution graph G_f and we find a Eulerian cycle in the graph starting in node (0, 0). Once such a cycle is obtained, we reconstruct a 2-BCPP solution by looping over every arc a in the Eulerian cycle and either inserting the bar chart associated with a into bin j if a is an item arc or incrementing j if a is a transition arc (j is initialized to one).

4.2. Constraint Generation Framework

A common way to solve ILP models that involve an exponential number of constraints is to use a constraint generation framework. The key idea behind constraint generation is to solve the formulation with a subset of constraints $\mathcal{S}$ (e.g., with an ILP solver) and check if the returned solution is feasible for the initial problem (i.e., if the solution found also satisfies the constraints that were not included in $\mathcal{S}$). If this is the case, then the solution is optimal. Otherwise, it means that at least one of the constraints not considered in the reduced model is important and needs to be added to $\mathcal{S}$. The process is iterated until an optimal solution is found. Note that constraint generation is very well-known in the routing community (see the work of Pferschy and Staněk 2017 on constraint generation for the traveling salesman problem).

A flowchart of our constraint generation framework is presented in Figure 6. After the graph-construction phase in which graph $G=(V,A)$ is generated, we solve the Eulerian-flow model without subtour elimination constraints to obtain a solution f. We then create the solution graph G_f, split the arcs of A_f into the minimum number of disjoint cycles, and store the ones that do not contain node (0, 0) in $\mathcal{C}$. If $|\mathcal{C}|>0$, then f contains one or more subtours, and thus, cannot be mapped into a feasible 2-BCPP solution. As a result, a subtour elimination constraint is added to $\mathcal{S}$ forbidding each of the $\left|\mathcal{C}\right|$ subtours, and Eulerian-flow is solved once more with the updated subset $\mathcal{S}$. If $|\mathcal{C}|=0$, the algorithm stops as f can be mapped into a feasible 2-BCPP solution.

4.3. Link-Flow

Even though reduction techniques are applied in the graph-construction phase of Eulerian-flow, preliminary tests show that the model size grows quickly as n and c increase even when the constraint generation framework is used. This is mainly because of the fourth step of the graph-construction phase in which the initial arc set $A^{\alpha }$ is duplicated up to c – 1 times to allow any (ordered) combination of bar chart arcs to follow a transition arc. The key idea behind link-flow is to move the transition arc from the beginning of a bin to the end so that every bin now always starts from node (0, 0) and ends with a transition arc. Whereas the idea appears to be relatively simple and effective (as the duplication phase is no longer needed), its implementation is not necessarily trivial as we now need to make sure that the transition arcs “connect” to each other so that the bins form a feasible 2-BCPP solution. Intuitively, transition arcs can be seen as a set of pieces of a puzzle, and the subset we select needs to form a circle.

For the sake of clarity, we describe the link-flow model using two distinct graphs G¹ and G² (even though we will observe in the solving phase that building G² is not strictly necessary). The primary graph G¹ contains item arcs and link arcs (the counterpart of transition arcs in link-flow), whereas the secondary graph G² contains link and loss arcs. Constraints in G¹ make sure that the bins created are feasible and item requirements are satisfied. Constraints in G² make sure that link arcs (and, thus, bins) can be combined together to form a feasible 2-BCPP solution (possibly with the help of loss arcs). This is done by ensuring that the secondary solution graph $G_f^2$ contains a Eulerian cycle. Supplementary constraints ensure that the link arcs selected in G¹ and in G² are the same.

In the graph-construction phase, the primary graph $G^1=(V^1,A^1)$ is initially constructed as the arcflow graph for the 2D-VPP (see Section 3.2 and Algorithm 1 in Online Section A). Link arcs $((v_1,v_2),(0,0))$ are then added to A¹ for every node $(v_1,v_2)\in V^1$. The secondary graph $G^2=(V^2,A^2)$ has vertex set $V^2=\{0,1,\dots ,c\}$ and arc set $A^2=\{(c-v_1,v_2,lk):(v_1,v_2)\in V^1\}\cup \{(d,d+1,ls):0\le d<c\}$ in which the third index separates the link arcs (identified by “lk”) from the loss arcs (identified by “ls”).

In the solving phase, for every primary vertex $(v_1,v_2)\in V^1$, we still use ${\delta }^{+}(v_1,v_2)$ as the set of primary arcs leaving (v₁, v₂) and ${\delta }^{-}(v_1,v_2)$ as the set of primary arcs entering (v₁, v₂). Analogously, we define index sets ${\delta }^{+}(v)$ and ${\delta }^{-}(v)$ for every secondary vertex $v\in V^2$ as the set of secondary arcs leaving and entering v. For every $i=1,\dots ,n$, we denote by $A_i^1$ the set of primary arcs representing the packing of a bar chart of type i. We call $A_1^1\cup A_2^1\cup \cdots \cup A_n^1$ the set of bar chart arcs, and we denote by $A_0^1$ the set of link arcs in G¹. By introducing integer decision variables $f_a^1$ ($f_a^2$) that indicate the number of times primary (secondary) arc $a\in A^1$ ($a\in A^2$) is selected, the link-flow model can then be defined as follows:

The objective function (24) minimizes the length of the packing. Constraints (25) and (26) enforce flow conservation in the primary and secondary graphs, whereas Constraint (27) makes sure that at least one unit of flow leaves secondary vertex 0. Furthermore, Constraints (28) ensure that the link arcs selected in G¹ are the same as the link arcs selected in G², Constraints (29) enforce demand satisfaction, and Constraints (30) eliminate subtours so that the solution graph $G_f^2$ is strongly connected. A suitable big-M value for Constraints (30) is K + Kc as at most K link arcs are selected in G¹ (and, thus, in G²) and at most c loss arcs are required in G² per link arc. Note that we can decrease big-M values by making them constraint-specific. In that case, for the constraint associated with a given subset of nodes ${\tilde{V}}^2$, we can use the actual number of loss arcs having both tail and head in ${\tilde{V}}^2$ instead of the upper bound c. Model (24)–(32) uses $O(n{c}^2)$ variables and $O(2^c)$ constraints. Note that we can reduce the model size even further by replacing variables $f_{d,e,lk}^2$ associated with secondary link arcs $(d,e,lk)\in A^2$ by their counterpart $f_{(c-d,e),(0,0)}^1$ in the primary graph and by removing Constraints (28). We can also use the constraint generation framework introduced previously to handle subtour elimination Constraints (30).

In the mapping phase, after an optimal solution $f=(f^1,f^2)$ is found for Model (24)–(32), we use a flow-decomposition algorithm to obtain a set of $(0,0)-(0,0)$ cycles with unitary flow in the primary graph. Each of these cycles corresponds to a bin, and we now need to order these bins. To do so, in the secondary graph, we create solution graph $G_f^2$, and we find a Eulerian cycle starting from node 0. Note that $G_f^2$ is constructed as described in Eulerian-flow: $G_f^2=(V_f^2,A_f^2)$, where $V_f^2$ contains the set of secondary vertices visited in solution f and $A_f^2$ contains $f_{d,e,l}^2$ copies of arc (d, e, l) for every $(d,e,l)\in A^2$. After initializing j to one, we reconstruct a 2-BCPP solution by looping over every arc a in the Eulerian cycle and inserting in bin j either one unit of loss space if a is a loss arc or the bar charts of the $(0,0)-(0,0)$ cycle that includes the counterpart of a in the primary graph otherwise. Index j is incremented every time a is a link arc.

5. Computational Experiments

We empirically evaluate the performance of each of the proposed approaches, namely,

• CSP, the arcflow Model (9)–(12) described in Section 3.1 to obtain a lower bound.

• EUL-REL, the Eulerian-flow Model (18)–(23) without subtour elimination Constraints (22) described in Section 4.1 to obtain a lower bound.

• LINK-REL, the link-flow Model (24)–(32) without subtour elimination Constraints (30) described in Section 4.3 to obtain a lower bound.

• 2DV, the arcflow Model (13)–(16) described in Section 3.2 to obtain an upper bound.

• ERZIN, Model (1)–(8) described in Section 2 introduced by Erzin et al. (2021b).

• EUL-FLOW, the Eulerian-flow Model (18)–(23) in which subtour elimination constraints are added through the constraint generation framework described in Section 4.2.

• LINK-FLOW, the link-flow Model (24)–(32) in which subtour elimination constraints are added through the constraint generation framework described in Section 4.2.

For comparative purposes, we also reimplement the best performing heuristic according to the experiments of Erzin et al. (2021b), GA_LO. It is a greedy approach that fills the bins one at a time by iteratively adding the bar chart with the largest first item among those fitting in the remaining space. When no bar chart can be packed into the current bin j, because either the first bar does not fit into bin j or the second bar does not fit into bin j + 1, then GA_LO closes bin j and starts filling bin j + 1. Our algorithms are all coded in C++ and can be downloaded from https://github.com/mdelorme2/Two_Bar_Charts_Packing_Problem. The experiments were run on an Intel(R) Core(TM) i7-4700MQ, 2.40 GHz with 8 GB of memory, running under Ubuntu 20.04.4 LTS, and Gurobi 9.5.0 was used to solve the ILP models. A single core was used for the tests (i.e., parameter Threads was set to one), the barrier algorithm was used to solve the root nodes of the ILP models (i.e., parameter Method was set to two), and callbacks were used for the constraint generation (and, thus, parameter LazyConstraints was set to one). In every instance, an upper bound K on the minimum number of bins required to pack the bar charts was obtained with a trivial first fit decreasing approach. For every run, a time limit of 3,600 seconds was imposed. We tested our approaches on six newly generated data sets with integer item weights and bin capacities, namely,

• U-GEN, in which bar charts contain (a combination of) small, medium, and big items.

• U-SMA, in which bar charts only contain small items.

• U-MED, in which bar charts contain at least one medium item.

• U-BIG, in which bar charts contain at least one big item.

• TRIPLETS, in which an integer solution where every bin is entirely filled exists.

• DONUTS, in which the bound derived from EUL-REL and LINK-REL is exactly one bin away from the optimal solution value.

Our data sets can be downloaded from https://github.com/mdelorme2/Two_Bar_Charts_Packing_Problem. For each data set, we first describe the instance-generation process and its motivation. We then empirically evaluate the performance of each of the aforementioned approaches, including an overview of the model size and quality of the LP relaxation.

5.1. Data Set Generation and Experiments on U-GEN, U-SMA, U-MED, and U-BIG

5.1.1. Data Set Generation.

In data sets U-GEN and U-SMA, the item weights of every bar chart are uniformly distributed in specific ranges: $[1;c]$ in U-GEN and $\left[1;\frac{c}{10}\right]$ in U-SMA. In data sets U-MED and U-BIG, the weight of one item of every bar chart (randomly selected between the first and second items) is uniformly distributed in $[1;c]$, whereas the weight of the other item is uniformly distributed in $\left[\frac{c}{4};1\right]$ in U-MED and $\left[\frac{c}{2};1\right]$ in U-BIG. For each data set, we consider three distinct bin capacities $c\in \{50,100,500\}$ and various numbers of bar charts ${\sum }_{i=1}^n{q}_i$ (named Ω thereafter).

• For c = 50, we tried $\mathrm{\Omega }\in \{10;50;100;500;1,000;5,000;10,000;50,000;100,000\}$.

• For c = 100, we tried $\mathrm{\Omega }\in \{10;50;100;1,000\}$.

• For c = 500, we tried $\mathrm{\Omega }\in \{10;50;100\}$.

For each of the 17 combinations, we generated 10 instances resulting in 170 instances in total per data set. To create an instance, we simply generated the desired Ω bar charts according to the specified distribution, and we merged the identical ones to obtain values n and q_i $(i=1,\dots ,n)$.

5.1.2. Data Set Motivation.

Data sets similar to U-GEN and U-BIG were used in the experiments of Erzin et al. (2021a) on a standardized version of the 2-BCPP in which the bin capacity was always equal to one and the item weights were real numbers with up to eight digits after the decimal point. We decided to recreate these two data sets in order to have more flexibility over the capacity and the total number of bar charts. We also mention that there were some inconsistencies in the instance files we downloaded from the web link shared in Erzin et al. (2021a),¹ which is another motivation for generating new data sets. As far as U-MED and U-SMA are concerned, the former was created to be a more difficult version of U-BIG, whereas the latter was designed to outline the main weaknesses of our new flow formulation approaches.

5.1.3. Experiments.

Our first set of experiments aims at empirically evaluating three different rules for ordering the bar charts and creating the graphs in Eulerian-flow and link-flow. Rule “FID” sorts the bar charts according to the first item weight (the second item weight is used to break the ties), rule “SID” sorts the bar charts according to the second item weight (the first item weight is used to break the ties), and “MID” sorts the bar charts according to the sum of the first and second item weights (ties are broken randomly). In all three cases, we sorted the values in nonincreasing order as it is the most common choice in arcflow formulations to reduce the model size (see Valério de Carvalho 1999). We report in Table 1 the results obtained by EUL-FLOW and LINK-FLOW with each of the three ordering rules. The first two columns identify the approach and the ordering rule used. The three following columns give some indicators of the performance of each method: the number of optimal solutions found, average CPU time over all runs (including the ones terminated by the time limit), and the average number of cuts added in the constraint generation part. Note that instances terminated because of memory issues were counted as unsolved within the time limit. The three last columns report some details about the model size: average number of variables, constraints, and nonzero elements.

Table 1. Comparison of the Ordering Rules for Eulerian-Flow and Link-Flow on U-GEN Instances

From the table, we observe that shorter average computing times are obtained when using ordering rule MID. This can be explained by the fact that the models generated with MID are slightly smaller compared with those generated with FID and SID. We point out that, even if the research community tends to agree that there is a strong correlation between smaller models and better computing performance (all other things being equal), there are some random components within ILP solvers that could explain why EUL-FLOW solves more instances to optimality with SID than it does with MID. Nevertheless, as the model size is the main drawback of arcflow formulations, we decided to use ordering rule MID in the rest of our experiments.

We show in Table 2 the results obtained by each of the tested approaches on the 17 $(c,\mathrm{\Omega })$ combinations for data set U-GEN. The first two columns identify the combination, whereas the third column reports the average optimal solution value computed over the 10 instances of the combination. The remaining columns give some indicators of the performance of each method. For the bounding procedures (i.e., CSP, EUL-REL, LINK-REL, 2DV, and GA_LO), column “diff” indicates the average difference between the bound obtained by the approach and the optimal solution value. Columns “#opt,” “time,” and “#cuts” report the number of optimal solutions found, the average CPU time over all runs (including the ones terminated by the time limit), and the average number of cuts added in the constraint generation part, respectively. For EUL-REL and LINK-REL, an instance is considered to be solved to optimality if the solution found by the solver was optimal and contained no subtours. Finally, column “abs gap” shows the average absolute gap obtained by the exact approaches (i.e., the average difference between the upper and lower bound found by the solver). When an approach obtains particularly bad results for intermediate $(c,\mathrm{\Omega })$ combinations $($e.g., ERZIN for combination (50,100)$)$, we decided not to try larger combinations and filled the associated rows with sign “-”. When Gurobi was not able to compute the continuous relaxation before the time limit $($e.g., EUR-REL for combination (500,100)$)$, near-zero lower bound values were reported by the solver. In order to avoid abnormally large and hardly interpretable average absolute gaps in the table, we replaced those by “NI” (standing for “not interpretable”).

Table 2. Results of the Tested Approaches for U-GEN Instances

The results in the table show that

• The lower bound derived from CSP is very good, especially when the total number of bar charts increases. We mention that the average time required to solve CSP was less than five seconds for every combination.

• The lower bounds derived from EUL-REL and LINK-REL are almost always equal to the optimal solution value, but the associated models take a long time to solve (sometimes even longer than their counterpart EUL-FLOW and LINK-FLOW). The solutions found by both approaches are often free of subtours, especially for large Ω values. We also note that LINK-REL is faster than EUL-REL on average and solves more instances to optimality.

• Exact approach ERZIN is not able to solve instances with as few as 50 bar charts, confirming the results reported by Erzin et al. (2021b). Further analysis (see Table 8 in Online Section B) shows that such bad performance can be attributed to the model size (as K increases with Ω) and to the poor quality of the LP-relaxation bound (as large bar charts can conveniently be broken in an LP solution to reduce the required number of bins).

• Exact approach EUL-FLOW is particularly effective as it can solve all but seven instances in the data set. LINK-FLOW is even better as it requires less computing time on average and can solve every instance. Further analysis (see Online Table 8) shows that these better results mostly come from smaller model sizes (e.g., EUL-FLOW required 1,022,159 variables on average for instances with 1,000 bar charts and capacity 100 versus 515,538 for LINK-FLOW). We also note that more cuts are required when there are fewer bar charts for both approaches. This makes sense as subtours are more likely to happen when fewer arcs are selected in the solution. It is also interesting to notice that LINK-FLOW requires less cuts on average than EUL-FLOW. This also makes sense as subtours are more likely to happen in the graph of EUL-FLOW than in the secondary graph of LINK-FLOW because the former has up to ${(c+1)}^2$ nodes, whereas the latter has c + 1 nodes at most.

• Even though the upper bound derived from GA_LO is often several bins away from the optimal solution value, it is usually better than the one derived from 2DV, and is much faster to compute (less than 0.1 second for GA_LO versus up to 100 seconds on average for 2DV).

We show in Table 3 the results obtained by each of the tested approaches on data set U-SMA and report the most noticeable changes with respect to U-GEN hereafter. We observe that ERZIN is now able to solve instances with up to 5,000 bar charts and is barely impacted by the bin capacity. Further analysis (see Online Table 9) shows that this change in performance can be attributed to the decrease in the model size (as bar charts are smaller, K now increases at a slower pace with Ω and the number of distinct bar charts n decreases) and to the increase in the quality of the LP-relaxation bound (as breaking small bar charts is less likely to reduce the required number of bins). We also notice that the flow formulations now struggle to solve instances with more than 50 bar charts when the bin capacity is above 100. Further analysis (see Online Table 9) shows that such a decrease in performance can be attributed to the increase in the model size driven by smaller bar charts. We also observe that for c = 100, it is easier for EUL-FLOW to solve instances with 500 or 1,000 bar charts than it is to solve instances with 50 or 100 bar charts. This can be explained by the fact that integer solutions found for instances containing many bar charts are less likely to contain subtours (all other things being equal). For example, 15 cuts on average were added by EUL-FLOW for instances with 50 bar charts, whereas only 0.4 cuts on average were added for instances with 1,000 bar charts. It is also interesting to see that EUL-REL now solves more instances than EUL-FLOW. This can be explained by the fact that certain Gurobi-specific reductions and transformations are incompatible with lazy constraints and are, therefore, activated in EUL-REL but not in EUL-FLOW. Another observation is that EUL-FLOW now solves more instances than LINK-FLOW. A possible explanation could be that, even though EUL-FLOW requires more variables on average than LINK-FLOW (see Online Table 9), its structure makes it easier for Gurobi heuristics to find good quality solutions. It is also worth noticing that the upper bounds derived from 2DV are now very close to the optimal solution values even though they still take a significant amount of time to compute (e.g., more than 1,200 seconds on average for combination (100, 100)). As GA_LO clearly underperforms on this data set, we investigated further and discovered that, when solving a U-SMA instance, the algorithm tends to start by (efficiently) packing many nonincreasing bar charts together in the first bins (bar chart i is said to be nonincreasing if $w_{i1}\ge w_{i2}$ and increasing otherwise), but it struggles later on when it has to pack the increasing items together. Note that alternative strategies developed by Erzin et al. (2021b) (such as “left packing” and “right packing”) could be better suited for U-SMA instances but were not tested as they are not in the scope of our study.

Table 3. Results of the Tested Approaches for U-SMA Instances

Finally, we report in Table 4 the results obtained by each of the tested approaches on data sets U-MED and U-BIG. We observe that the comments made for U-GEN are also valid for U-MED and U-BIG: (i) the lower bounds derived from CSP are very good for U-MED (less than one bin away from the optimal solution on average) and slightly less good for U-BIG (between 0 and 35 bins away from the optimal solution), (ii) the lower bounds derived from EUL-REL and LINK-REL are almost always equal to the optimal solution value and are often free of subtours but they can be long to compute, (iii) ERZIN struggles to solve instances with more than 50 bar charts because of its model size and the poor quality of its LP-relaxation bound (see Online Table 10), (iv) EUL-FLOW and LINK-FLOW can solve every tested instance with LINK-FLOW being slightly faster on average than EUL-FLOW because of its reduced model size, and (v) the upper bounds derived from GA_LO are much better than the ones derived from 2DV.

Table 4. Results of the Tested Approaches for U-MED (Top) and U-BIG (Bottom) Instances