# Algorithms for the Bin Packing Problem with Scenarios

Yulle G. F. Borges111Corresponding Author, Email Address:
<EMAIL_ADDRESS>222Institute of Computing, University of Campinas,
Campinas–SP, Brazil. Vinícius L. de Lima††footnotemark:  Flávio K.
Miyazawa††footnotemark:  Lehilton L. C. Pedrosa††footnotemark:  Thiago A. de
Queiroz333Institute of Mathematics and Technology, Federal University of
Catalão, Catalão–GO, Brazil. Rafael C. S. Schouery††footnotemark:

###### Abstract

This paper presents theoretical and practical results for the bin packing
problem with scenarios, a generalization of the classical bin packing problem
which considers the presence of uncertain scenarios, of which only one is
realized. For this problem, we propose an absolute approximation algorithm
whose ratio is bounded by the square root of the number of scenarios times the
approximation ratio for an algorithm for the vector bin packing problem. We
also show how an asymptotic polynomial-time approximation scheme is derived
when the number of scenarios is constant. As a practical study of the problem,
we present a branch-and-price algorithm to solve an exponential model and a
variable neighborhood search heuristic. To speed up the convergence of the
exact algorithm, we also consider lower bounds based on dual feasible
functions. Results of these algorithms show the competence of the branch-and-
price in obtaining optimal solutions for about 59% of the instances
considered, while the combined heuristic and branch-and-price optimally solved
62% of the instances considered.

#### Keywords

Bin Packing Problem; Scenarios; Approximation Algorithm; Variable Neighborhood
Search; Branch-and-Price Algorithm;

##  1 Introduction

Cutting and packing problems have been widely studied in the context of
Operations Research, mainly because of their properties and real-world
applicability (Wäscher et al., 2007). Among these problems, the _Bin Packing
Problem_ (BPP) asks to pack a set of one-dimensional items, each of a given
size, into the least possible number of bins of given identical capacities,
where the total size of items in a bin does not exceed the bin capacity.
Several variants of this problem have been studied in the literature, either
due to their theoretical interest, or their high applicability in different
industries (Sweeney and Paternoster, 1992; Dahmani et al., 2013).

The BPP and its variants have been extensively investigated since the thirties
(Kantorovich, 1960). Such problems are among the most studied in approximation
contexts. A recent survey by Coffman et al. (2013) presented over 180
references related to approximation results for the BPP and its variants. In
particular, the BPP is shown to be APX-hard, and indeed no algorithm has an
approximation factor smaller than $3/2$, unless $\mbox{P}=\mbox{NP}$
(Vazirani, 2001). Regarding practical techniques for the BPP, a survey by
Delorme et al. (2016a) reviewed models and solution methods, and
experimentally compared the available software tools to solve the problem.
Recent practical contributions for the BPP include the generalized arc-flow
formulation by Brandão and Pedroso (2016), the cooperative parallel genetic
algorithm by Kucukyilmaz and Kiziloz (2018), the reflect formulation by
Delorme and Iori (2019), the branch-and-cut-and-price algorithm by Wei et al.
(2020b) and the framework of de Lima et al. (2023).

Several variants of the BPP have also been extensively investigated in the
literature. For example, in the variant with fragile objects, the capacity of
a bin is directly influenced by its most fragile item (Clautiaux et al.,
2014). In the BPP with precedence constraints, a precedent item cannot be
packed in a bin later than its successors (Pereira, 2016). The variant with
overlapping items assumes that some subsets of items, when packed in the same
bin, occupy less capacity than the sum of their individual size (Grange et
al., 2018). In the BPP with item fragmentation, items can be split and
fragmentally packed (Bertazzi et al., 2019). In the generalized BPP, bins may
have different cost and capacity and items are associated to a profit and are
divided by compulsory (i.e., mandatory to load) and non-compulsory. The aim is
to minimize the overall cost based on the bins cost and items profit (Baldi et
al., 2019). In the online BPP, items arrive sequentially and each must be
packed before the arrival of the next one, with no information about the next
items (Balogh et al., 2019). The temporal variant considers one additional
dimension in the BPP, and the bin capacity should not be violated at any unit
of a discretized time horizon (Dell’Amico et al., 2020).

One of the current challenges to solving practical logistic problems is to
deal with the uncertainty arising from real-world applications (Saint-Guillain
et al., 2021; Xu et al., 2016; Juan et al., 2018; Baghalian et al., 2013). In
particular, one popular way to deal with this issue is by describing
scenarios. A scenario is defined as a possible outcome that may arise
depending on the problem’s uncertain variables. Some authors use scenarios to
consider the problem description as combinatorial rather than stochastic, as,
e.g., Feuerstein et al. (2014). Considering the scarcity of works that apply
this concept of scenarios to deal with uncertainty in bin packing problems and
the existence of practical applications, Bódis and Balogh (2018) recently
introduced the _Bin Packing Problem with Scenarios_ (BPPS).

The BPPS is a generalization of the BPP where each item belongs to a subset of
scenarios, and a packing must respect the capacity constraints of the bins for
each scenario individually. Whilst the set of scenarios is known in advance,
only one of the scenarios will be realized. This introduces the possibility of
packing in a single bin items whose combined sizes surpasses the bin capacity,
as long as the bin capacity is respected in each individual scenario. Although
Bódis and Balogh (2018) introduced three objective functions for the BPPS, we
are concerned with the objective of minimizing the number of bins of the
worst-case scenario, which is the most challenging one. In this way, the BPP
is the particular case of the BPPS where there is a single scenario.

A generalization of the BPP that is similar to the BPPS is the _Vector Bin
Packing Problem_ (VBPP), where bins and items have multiple dimensions and the
bin capacity must be respected in all of its dimensions. The objective is the
same as in the BPP, to minimize the number of bins. Caprara and Toth (2001)
presented lower bounds, heuristics, and a branch-and-bound (B&B) algorithm for
the VBPP with two dimensions. Alves et al. (2014) investigated dual feasible
functions for this problem. Buljubašić and Vasquez (2016) proposed a local
search algorithms where a Tabu search and descent search with add and drop
moves were used to explore the search space. In Hu et al. (2017), a set-
covering model is solved by a branch-and-price (B&P) algorithm. In Heßler et
al. (2018), there is a stabilized branch-and-price algorithm with dual optimal
cuts. Recently, Wei et al. (2020a) developed a branch-and-price algorithm for
the vector packing with two dimensions, where a goal cut approach is used to
obtain lower bounds and a dynamic programming with branch-and-bound solves the
pricing problem, improving previous results of the literature.

To the best of our knowledge, the recent work by Bódis and Balogh (2018) is
the only one in the literature concerned with the BPPS. Motivated by its
interesting theoretical aspects, its general applicability and its relation to
well-studied problems, this paper proposes theoretical and practical results
for the BPPS. The contributions of the current work are the following: (i) an
absolute approximation algorithm for the BPPS based on the VBPP; (ii) an
_asymptotic polynomial time approximation scheme_ (APTAS) for the version of
the problem with a constant number of scenarios; (iii) an exact method for the
problem based on a branch-and-price algorithm for an exponential _Integer
Linear Programming_ (ILP) model; and (iv) a _Variable Neighborhood Search_
(VNS) heuristic.

This paper is organized as follows. In Section 2 we provide formal definitions
to be used to contextualize the contributions of the paper. In Section 3, we
show that an approximate solution of the BPPS can be obtained by solving an
instance of the VBPP, leading to an absolute approximation algorithm for the
BPPS. In Section 4, we present an APTAS for the BPPS, when the number of
scenarios is constant. As for practical results, Section 5 describes a VNS
algorithm for the BPPS and Section 6 presents a branch-and-price algorithm to
solve an exponential model for the BPPS. The evaluation of the proposed VNS
and branch-and-price algorithms is performed in Section 7, from experiments
based on the solution of randomly generated instances. Finally, Section 8
presents the conclusions and directions for future research.

##  2 Formal Definitions

In this section, we provide some definitions to contextualize the
contributions of this work. For simplicity, throughout the paper, we denote,
for any natural $n$, the set $\\{1,\dots,n\\}$ simply by $[n]$.

#### Bin Packing Problem with Scenarios (BPPS)

In the BPPS, we are given a number $d$ of scenarios, a set $\mathcal{I}=[n]$
of $n$ items, with each $i\in\mathcal{I}$ having a size
$s_{i}\in\mathbb{Q}_{+}$ and a set of scenarios $\mathcal{K}_{i}\subseteq[d]$,
and an unlimited number of identical bins of capacity $W\in\mathbb{Q}_{+}$.
For each $k\in[d]$, the set of items in scenario $k$ is denoted by
$S_{k}=\\{i\in\mathcal{I}\mid k\in\mathcal{K}_{i}\\}$. A solution
$\mathcal{B}$ of the BPPS is a partition of $\mathcal{I}$ such that for each
part $B\in\mathcal{B}$ and for each scenario $k\in[d]$, $\sum_{i\in B\cap
S_{k}}s_{i}\leq W$. The objective of the BPPS is to find a solution that
minimizes

$\textit{V}_{\textit{BPPS}}(\mathcal{B})=\max_{k\in[d]}|\\{B\in\mathcal{B}:B\cap
S_{k}\neq\emptyset\\}|.$

The objective corresponds to minimize the number of bins of the worst-case
scenario, i.e., the scenario with the largest number of bins. A closely
related problem is the VBPP.

#### Vector Bin Packing Problem (VBPP)

In the VBPP, we are given a number $d$ of resources, a set $\mathcal{I}=[n]$
of $n$ items, each $i\in\mathcal{I}$ having a vector
${s_{i}\in\mathbb{Q}_{+}^{d}}$ of resource consumption, and an unlimited
number of identical bins of resource capacity given as a $d$-dimensional
vector ${W=(w_{1},w_{2},\dots,w_{d})}\in\mathbb{Q}_{+}^{d}$. A solution
$\mathcal{B}$ is a partition of $\mathcal{I}$ such that for each part
$B\in\mathcal{B}$ and each resource $k\in[d]$, $\sum_{i\in B}s_{ik}\leq
w_{k}$. The objective is to find a solution that minimizes

$\textit{V}_{\textit{VBPP}}(\mathcal{B})=|\mathcal{B}|.$

The VBPP is a generalization of the BPP, and thus it is also APX-hard.
Woeginger (1997) showed that there is no APTAS for the two-dimensional version
of this problem unless $\mbox{P}=\mbox{NP}$. Chekuri and Khanna (2004) gave an
asymptotic $O(\ln d)$-approximation for the case where $d$ is a fixed
constant. Later on, this result was improved to $1+\ln d$ by Bansal et al.
(2009). Also, as noted by Christensen et al. (2017), the APTAS of Fernandez de
la Vega and Lueker (1981) implies a $(d+\varepsilon)$-approximation for the
VBPP.

For the theoretical results of this paper in Sections 3 and 4, we assume,
without loss of generality, that the instances of the BPPS and the VBPP are
normalized by a proper scaling of the items, so that the capacity of the bins
correspond $1$ on each dimension.

An important concept in both theoretical and practical results for the BPP and
related variants is the one of _cutting pattern_ : a feasible combination of
items in a single bin. For the BPPS, we represent a cutting pattern as a
vector $p=(a_{1p},\dots,a_{np},b_{1p},\dots,b_{dp})$ of binary coefficients,
such that: each coefficient $a_{ip}$ is equal to $1$ if and only if item $i$
is in pattern $p$, and each coefficient $b_{kp}$ is equal to $1$ if and only
if some item of the scenario $k$ is in pattern $p$. We denote by $\mathcal{P}$
the set of all feasible patterns for the BPPS.

##  3 An Absolute Approximation Algorithm

In the following, we consider a mapping between the set of instances of the
BPPS and a subset of instances of the VBPP. Starting with a normalized
instance ${I=(d,\mathcal{I},\mathcal{K},s)}$ of the BPPS, we construct a
normalized instance ${I^{\prime}=(d,\mathcal{I},s^{\prime})}$ of the VBPP. For
that, let, for each $i\in\mathcal{I}$ and $k\in[d]$,

$s_{ik}=\begin{cases}s_{i}&\mbox{if }i\in S_{k},\\\
0&\mbox{otherwise.}\end{cases}$

It is simple to verify that a solution $\mathcal{B}$ for $I$ is a solution for
$I^{\prime}$, and vice-versa.

######  Definition 1.

A solution $\mathcal{B}$ is minimal for the BPPS if, for any two bins
${A,B\in\mathcal{B}}$, there exists a scenario $k$, such that

$\sum_{i\in A\cap S_{k}}s_{i}+\sum_{i\in B\cap S_{k}}s_{i}>1.$

Notice that the value of $\mathcal{B}$ for the BPPS is not larger than the
value of $\mathcal{B}$ for the VBPP. In the opposite direction, we have
Theorem 1.

######  Theorem 1.

If $\mathcal{B}$ is a minimal solution for the BPPS, then
$\textit{V}_{\textit{VBPP}}(\mathcal{B})\leq\sqrt{d}\textit{V}_{\textit{BPPS}}(\mathcal{B})$.

###### Proof.

We say that two bins $A$ and $B$ are incompatible in a scenario $k$ if both
bins, $A$ and $B$, are non-empty in scenario $k$. Let $G$ be a multigraph
where each bin represents a vertex, and for each pair of bins $A$ and $B$ and
scenario $k$, there is one edge between $A$ and $B$ (labeled as $k$) if they
are incompatible in $k$.

For each scenario $k\in[d]$, we define $h_{k}=|\\{B\in\mathcal{B}:B\cap
S_{k}\neq\emptyset\\}|$, i.e., $h_{k}$ is the number of bins used by
$\mathcal{B}$ in $k$. Let $h=\max_{k\in[d]}h_{k}$, and notice that, since, for
each scenario, there is an edge for each two used bins, the number of edges in
$G$ is

$|E(G)|=\sum_{k=1}^{d}\frac{h_{k}(h_{k}-1)}{2}\leq d\frac{h(h-1)}{2}.$

Let $r$ be the number of bins. Since $\mathcal{B}$ is minimal, any two bins
are incompatible in at least one scenario and $G$ must have at least one edge
between any two vertices, i.e., $|E(G)|\geq\frac{r(r-1)}{2}$. Thus,

$\frac{r(r-1)}{2}\leq d\frac{h(h-1)}{2},$

and, since $r\leq dh$, as each bin is non-empty in at least one scenario, we
have that

$r^{2}\leq dh^{2}+r-dh\leq dh^{2}.$

The result follows as $r=\textit{V}_{\textit{VBPP}}(\mathcal{B})$ and
$h=\textit{V}_{\textit{BPPS}}(\mathcal{B})$.∎

######  Corollary 1.

Suppose there exists an $\alpha$-approximation for the VBPP, then there exists
an $\alpha\sqrt{d}$-approximation for the BPPS.

###### Proof.

Consider an instance $I$ of BPPS and the instance $I^{\prime}$ of VBPP
obtained as described above. Also, let $\mathcal{B}$ be the solution found by
the $\alpha$-approximation for instance $I^{\prime}$, $\mathcal{B}_{V}^{*}$ be
an optimal solution for VBPP with instance $I^{\prime}$ and
$\mathcal{B}_{S}^{*}$ be a minimal optimal solution for BPPS with instance
$I$. Then, we have

$\textit{V}_{\textit{BPPS}}(\mathcal{B})\leq\textit{V}_{\textit{VBPP}}(\mathcal{B})\leq\alpha\textit{V}_{\textit{VBPP}}(\mathcal{B}^{*}_{V})\leq\alpha\textit{V}_{\textit{VBPP}}(\mathcal{B}^{*}_{S})\leq\alpha\sqrt{d}\textit{V}_{\textit{BPPS}}(\mathcal{B}^{*}_{S}).\qed$

The bound given by this strategy cannot improve the dependency on $d$, by the
following lemma.

######  Theorem 2.

For each $d\geq 1$, there exists an instance $I$ of the BPPS with $d$
scenarios and a solution $\mathcal{B}$ for $I$, which is minimal for the BPPS,
such that
${\textit{V}_{\textit{VBPP}}(\mathcal{B})\geq\frac{\sqrt{d}}{2}\textit{V}_{\textit{BPPS}}(\mathcal{B})}$.

###### Proof.

Let $r=\left\lceil\sqrt{d}\right\rceil$. We construct an instance
$I=(d,\mathcal{I},\mathcal{K},s)$ that contains $r$ items.

We consider a complete graph $G$ with $r$ vertices. For each vertex $i$ of
$G$, create an item $i$ with size $s_{i}=1$, and for each edge $\\{i,j\\}$ of
$G$, create a scenario $k$ with items $S_{k}=\\{i,j\\}$. Notice that each item
is in exactly $r-1$ scenarios, and that the number of scenarios is

$|E(G)|=\frac{r(r-1)}{2}=\frac{\lceil\sqrt{d}\rceil(\lceil\sqrt{d}\rceil-1)}{2}\leq\frac{(\sqrt{d}+1)\sqrt{d}}{2}\leq
d.$

Now, we create a solution $\mathcal{B}$ with $r$ bins, such that, for each
item $i$, there is one bin $B_{i}=\\{i\\}$. We claim that $\mathcal{B}$ is
minimal for the BPPS. Indeed, for any two bins $B_{i}$ and $B_{j}$, items $i$
and $j$ are contained in a common scenario $k$, with $S_{k}=\\{i,j\\}$, and
$s_{i}+s_{j}>1$.

Notice that $\textit{V}_{\textit{VBPP}}(\mathcal{B})=r$. Let
$h=\textit{V}_{\textit{BPPS}}(\mathcal{B})$, so $h$ is the maximum number of
non-empty bins in a given scenario. Since each scenario contains exactly $2$
items, then $h=2$. Therefore,

$\displaystyle
r=\lceil\sqrt{d}\rceil\geq\frac{\sqrt{d}}{2}2=\frac{\sqrt{d}}{2}h.$ ∎

Finally, using the $1+\ln d$-asymptotic approximation algorithm of Bansal et
al. (2009) for the VBPP (for constant $d$) and the fact that the APTAS of
Fernandez de la Vega and Lueker (1981) implies a
$(d+\varepsilon)$-approximation for the VBPP Christensen et al. (2017), we
have the following result.

######  Corollary 2.

There exists a $(d+\varepsilon)\sqrt{d}$-approximation for the BPPS. Also, for
any constant $d$, there exists a $(1+\ln d)\sqrt{d}$-approximation for the
BPPS.

##  4 An Asymptotic Polynomial Time Approximation Scheme

We present an APTAS for the BPPS when the number of scenarios $d$ is a
constant. First, we observe that, if all items are large and the number of
distinct item sizes is bounded by a constant, then the number of valid
patterns of items in a bin is also bounded by a constant. It implies that for
this restricted case, a polynomial algorithm can be readily obtained by simply
enumerating the frequencies of all patterns.

To obtain such a restricted instance, we use the _linear grouping_ technique,
following Fernandez de la Vega and Lueker (1981), by grouping items of similar
sizes and creating a map from smaller to larger items. A naive approach does
not work, however, since one must take into account the scenarios. Thus, we
group items according to their scenarios separately, but we do the mapping
simultaneously.

Finally, to solve a general instance, we combine small items into artificial
large items. Again, this is done on a per-scenario basis, so as only
compatible items are combined. If the small items do not fit exactly into the
area of the artificial items, then this change may cause bins to be slightly
augmented. The surpassing items are then relocated to new bins greedily,
increasing the number of bins by just a fraction of the optimal value.

In the following, we use the notion of type. The _type_ of item $i$ is the set
of scenarios which contain $i$. The set of all types is denoted by
$\mathcal{T}$. Also, the value of an optimal solution for an instance $I$ is
denoted by $\operatorname{\mbox{OPT}}(I)$.

###  4.1 Restricted instances

In this subsection, we consider instances of the BPPS in which the items are
large and the number of distinct sizes is bounded by a constant. Precisely,
let $V$ be a set of positive numbers and consider a constant $\varepsilon$,
with ${0<\varepsilon\leq 1/4}$ and $1/\varepsilon$ is an integer. An instance
$(d,\mathcal{I},\mathcal{K},s)$ of the BPPS is said to be _$V$ -restricted_ if
$\\{s_{i}:i\in\mathcal{I}\\}=V$ and $s_{i}\geq\varepsilon^{2}$ for every
$i\in\mathcal{I}$.

Recall that a bin $B\subseteq\mathcal{I}$ is feasible if, for any scenario
$k\in[d]$, ${\sum_{i\in B\cap S_{k}}s_{i}\leq 1}$. For a type
$t\in\mathcal{T}$ and $v\in V$, let $B_{tv}$ be the subset of items in $B$
with type $t$ and size $v$. Observe that the family of all sets $B_{tv}$ form
a partition of $B$. The _pattern_ of $B$ is the vector
${(|B_{tv}|)}_{t\in\mathcal{T},v\in V}$.

######  Lemma 1.

If $|V|$ is constant, the number of distinct patterns for a $V$-restricted
instance is bounded by a constant.

###### Proof.

Consider a feasible bin $B$. Since $s_{i}\geq\varepsilon^{2}$, for any $i\in
B$,

$|B|\varepsilon^{2}\leq\sum_{i\in B}s_{i}\leq\sum_{k=1}^{d}\sum_{i\in B\cap
S_{k}}s_{i}\leq\sum_{k=1}^{d}1=d,$

where the last inequality follows because $B$ is feasible. It means that the
number of items in any bin is at most $d/\varepsilon^{2}$. Since a vector
corresponding to a pattern has $|\mathcal{T}||V|$ elements, the number of
patterns is at most ${(1+d/\varepsilon^{2})}^{|\mathcal{T}||V|}$, which is
constant since $d$ is constant.∎

######  Lemma 2.

If $|V|$ is constant, then an optimal solution of a $V$-restricted instance
can be computed in polynomial time.

###### Proof.

Consider a $V$-restricted instance $I$, and let $M$ be the set of distinct
patterns for $I$. Given a solution $\mathcal{B}$ and a pattern $p\in M$, the
number of bins in $\mathcal{B}$, which has (a non-empty) pattern $p$, is
denoted by $n_{p}$. Clearly $n_{p}\leq|\mathcal{I}|$ for every $p\in M$, as
there are only $|\mathcal{I}|$ items, and each bin has at least one item. The
configuration of $\mathcal{B}$ is the vector ${(n_{p})}_{p\in M}$. Thus, the
number of possible configurations is bounded by $|\mathcal{I}|^{|M|}$, which
is polynomial since, by Lemma 1, $|M|$ is constant.

Notice that, given a vector ${(n_{p})}_{p\in M}$, one can either find a
solution $\mathcal{B}$ with this configuration, or decide there is no such
solution. Indeed, it is sufficient to count the total number of items of each
type and size over all the patterns. Therefore, one can list all
configurations in polynomial time and return the one with the minimum value.∎

###  4.2 An APTAS for large items

Suppose now that we have an instance $I=(d,\mathcal{I},\mathcal{K},s)$, such
that all items are large, i.e., for each $i\in\mathcal{I}$, we have
$s_{i}\geq\varepsilon^{2}$, but the number of distinct sizes is not
necessarily bounded by a constant.

######  Definition 2.

Given instances $I=(d,\mathcal{I},\mathcal{K},s)$ and
$\bar{I}=(d,\bar{\mathcal{I}},\bar{\mathcal{K}},\bar{s})$, we say that $I$
dominates $\bar{I}$ if there exists a subset
$\mathcal{I}^{\prime}\subseteq\mathcal{I}$ and a bijection
$f:\mathcal{I}^{\prime}\to\bar{\mathcal{I}}$, such that, for every
$i\in\mathcal{I}^{\prime}$, items $i$ and $f(i)$ have the same type and
$s_{i}\geq\bar{s}_{f(i)}$. In this case, we write $I\succcurlyeq\bar{I}$.

######  Lemma 3.

If $I\succcurlyeq\bar{I}$, then
$\operatorname{\mbox{OPT}}(I)\geq\operatorname{\mbox{OPT}}(\bar{I})$.∎

In the following, we create an instance
$\bar{I}=(d,\bar{\mathcal{I}},\bar{\mathcal{K}},\bar{s})$ with
$I\succcurlyeq\bar{I}$, such that all items of $\bar{I}$ are large and the
number of distinct sizes is bounded by a constant.

First, for each $t\in\mathcal{T}$, let $\mathcal{I}_{t}$ be the set of items
of type $t$ and let $m=\frac{2^{d}}{\varepsilon^{3}}-1$. We consider the set
$\mathcal{I}_{t}$ sorted in decreasing order of size, and build $m+1$ groups
of consecutive items, obtaining a partition
$\\{G^{0}_{t},G^{1}_{t},\dots,G^{m}_{t}\\}$, where the first $m$ groups have
size $\lceil|\mathcal{I}_{t}|/(m+1)\rceil$, and the last group either has the
same size, or is smaller.

Now, for each $t\in\mathcal{T}$ and $\ell\in\\{1,2,\dots,m\\}$, we create a
set $\bar{G}^{\ell}_{t}$ as follows: for each item $i\in G^{\ell}_{t}$, add an
item to $\bar{G}^{\ell}_{t}$ with the same set of scenarios $t$, and with size
equal to the size of the largest item in $G^{\ell}_{t}$. Let
$\bar{\mathcal{I}}$ be the union of all sets $\bar{G}^{\ell}_{t}$, and
$\bar{I}$ be the instance induced by $\bar{\mathcal{I}}$.

######  Lemma 4.

$I\succcurlyeq\bar{I}$.

###### Proof.

For $\ell\in\\{1,2,\dots,m\\}$, as
$|G_{t}^{\ell-1}|\geq|G_{t}^{\ell}|=|\bar{G}_{t}^{\ell}|$, we can find a
bijection $f$ between a subset of $G_{t}^{\ell-1}$ and $\bar{G}_{t}^{\ell}$
for $\ell\in\\{1,2,\dots,m\\}$ such that $i$ and $f(i)$ has the same type.
Finally, as $s_{i}\geq s_{j}$ for $i\in G_{t}^{\ell-1}$ and $j\in
G_{t}^{\ell}$ and $\bar{s}_{f(i)}$ is the maximum size of an item in
$G_{t}^{\ell}$, we have that $s_{i}\geq\bar{s}_{f(i)}$.∎

For a type $t\in\mathcal{T}$, the group of the largest items is $G^{0}_{t}$.
For each item $i$ in this group, we pack it into a separate bin $\\{i\\}$. By
joining all these bins, we obtain a packing $\mathcal{P}^{0}$.

######  Lemma 5.

$\textit{V}_{\textit{BPPS}}(\mathcal{P}^{0})\leq\varepsilon\operatorname{\mbox{OPT}}(I)+2^{d}$.

###### Proof.

Let $t\in\mathcal{T}$ and fix an arbitrary scenario $k$ of $t$. We consider an
optimal solution for $I$, and let $\mathcal{P}_{k}$ be the set of bins used in
scenario $k$ in this solution. Therefore,
$|\mathcal{P}_{k}|\leq\operatorname{\mbox{OPT}}(I)$. Observe that, since $k$
belongs to $t$, the set of items included in bins of $\mathcal{P}_{k}$ must
contain $\mathcal{I}_{t}$. As, for each $i\in\mathcal{I}_{t}$,
$s_{i}\geq\varepsilon^{2}$,

$\varepsilon^{2}|\mathcal{I}_{t}|\leq\sum_{i\in\mathcal{I}_{t}}s_{i}\leq\sum_{B\in\mathcal{P}_{k}}\sum_{i\in
B\cap
S_{k}}s_{i}\leq\sum_{B\in\mathcal{P}_{k}}1=|\mathcal{P}_{k}|\leq\operatorname{\mbox{OPT}}(I).$

This implies that the value of $\mathcal{P}^{0}$ can be bounded by

$\sum_{t\in\mathcal{T}}|G^{0}_{t}|=\sum_{t\in\mathcal{T}}\left\lceil\frac{\varepsilon^{3}}{2^{d}}\,|\mathcal{I}_{t}|\right\rceil\leq
2^{d}+\sum_{t\in\mathcal{T}}\frac{\varepsilon\operatorname{\mbox{OPT}}(I)}{2^{d}}\leq\varepsilon\operatorname{\mbox{OPT}}(I)+2^{d}.$

To obtain a packing of the remaining items, we solve instance $\bar{I}$. For
this, let $V$ be the set of distinct sizes in $\bar{I}$. For each
$t\in\mathcal{T}$, we created $m$ groups, where each group has items of the
same size. Therefore, ${|V|\leq 2^{d}m}$. This means that $\bar{I}$ is a
$V$-restricted instance for $|V|$ bounded by a constant. Using Lemma 2, we
obtain a solution $\bar{\mathcal{P}}^{1}$ for $\bar{I}$ in polynomial time.
Since $I\succcurlyeq\bar{I}$, we can obtain a solution
$\mathcal{P}^{0}\cup\mathcal{P}^{1}$ for $I$, where $\mathcal{P}^{1}$ is
obtained from $\bar{\mathcal{P}}^{1}$ by replacing items of
$\bar{G}^{\ell}_{t}$ by items of $G^{\ell}_{t}$. Let $I^{1}$ be the instance
obtained from the items of $\mathcal{P}^{1}$. Observe that $\mathcal{P}^{1}$
is feasible, since we replaced items of larger size which appear in the same
set of scenarios, thus for each scenario, the capacity of each bin remains
respected.

######  Lemma 6.

$\textit{V}_{\textit{BPPS}}(\mathcal{P}^{1})\leq\operatorname{\mbox{OPT}}(I)$.

###### Proof.

Note that the number of unused bins in each scenario is the same for
$\mathcal{P}^{1}$ and $\bar{\mathcal{P}}^{1}$, thus
$\textit{V}_{\textit{BPPS}}(\mathcal{P}^{1})=\textit{V}_{\textit{BPPS}}(\bar{\mathcal{P}}^{1})$.
Since, $\bar{\mathcal{P}}^{1}$ is optimal for $\bar{I}$, we have
$\textit{V}_{\textit{BPPS}}(\bar{\mathcal{P}}^{1})=\operatorname{\mbox{OPT}}(I^{1})$.
Now, using Lemma 3 and the fact that $I\succcurlyeq\bar{I}$, we have
$\operatorname{\mbox{OPT}}(I^{1})\leq\operatorname{\mbox{OPT}}(I)$.∎

Combining the previous results, one obtains a APTAS for instances with large
items only.

######  Lemma 7.

Suppose for every $i\in\mathcal{I}$, $s_{i}\geq\varepsilon^{2}$. Then one can
find in polynomial time a packing $\mathcal{P}$ of $\mathcal{I}$ with
$\textit{V}_{\textit{BPPS}}(\mathcal{P})\leq(1+\varepsilon)\operatorname{\mbox{OPT}}(I)+2^{d}$.

###### Proof.

Define $\mathcal{P}=\mathcal{P}^{0}\cup\mathcal{P}^{1}$. Notice that
$\mathcal{P}$ contains every item of $\mathcal{I}$, and thus it is feasible.
Using lemmas 5 and 6,

$\textit{V}_{\textit{BPPS}}(\mathcal{P})\leq\textit{V}_{\textit{BPPS}}(\mathcal{P}^{0})+\textit{V}_{\textit{BPPS}}(\mathcal{P}^{1})\leq\varepsilon\operatorname{\mbox{OPT}}(I)+2^{d}+\operatorname{\mbox{OPT}}(I).\qed$

###  4.3 An APTAS for the general case

In the general case of the BPPS, an instance $I=(d,\mathcal{I},\mathcal{K},s)$
may contain large and small items. Let
$\mathcal{L}=\\{i\in\mathcal{I}:s_{i}\geq\varepsilon^{2}\\}$ be the set of
large items, and $\mathcal{S}=\mathcal{I}\setminus\mathcal{L}$ be the set of
small items.

To obtain a packing of $\mathcal{L}\cup\mathcal{S}$, we first replace
$\mathcal{S}$ by the set of large items, $\hat{\mathcal{S}}$, which we define
as follows. For each type $t\in\mathcal{T}$, let $\mathcal{S}_{t}$ be the set
of small items with type $t$. Now, let $\hat{\mathcal{S}}_{t}$ be a set of
$\lceil\sum_{i\in\mathcal{S}_{t}}s_{i}/(\varepsilon-\varepsilon^{2})\rceil$
items, such that each $j\in\hat{\mathcal{S}}_{t}$ has size
$\hat{s}_{j}=\varepsilon$ and is of type $t$. The set $\hat{\mathcal{S}}$ is
the union of sets $\hat{\mathcal{S}}_{t}$ for every type $t$.

Thus, we create an instance
$\hat{I}=(d,\hat{\mathcal{I}},\hat{\mathcal{K}},\hat{s})$ whose set of items
is $\hat{\mathcal{I}}=\mathcal{L}\cup\hat{\mathcal{S}}$, each of which has
size at least $\varepsilon^{2}$. Lemma 8 compares the optimal value of $I$ and
$\hat{I}$.

######  Lemma 8.

$\operatorname{\mbox{OPT}}(\hat{I})\leq(1+2^{d+1}(d+1)\varepsilon)\operatorname{\mbox{OPT}}(I)+1$.

###### Proof.

Consider an optimal solution $\mathcal{P}^{*}$ for $I$. For each bin
$B\in\mathcal{P}^{*}$ and type $t\in\mathcal{T}$, we define
$B_{t}=B\cap\mathcal{S}_{t}$, i.e., $B_{t}$ is the set of small items with
type $t$. Then, we replace the items in $B_{t}$ by $\lceil\sum_{i\in
B_{t}}s_{i}/(\varepsilon-\varepsilon^{2})\rceil$ new items of size
$\varepsilon$, and with the same scenarios of $t$. We call the modified
packing by $\mathcal{P}^{\prime}$.

Observe that $\mathcal{P}^{\prime}$ may be infeasible, since new items may
surpass the bins’ capacity. We can obtain a feasible solution
$\mathcal{P}^{\prime\prime}$ by relocating some created items from
$\mathcal{P}^{\prime}$ to new bins. For each type $t$ and each bin
$B\in\mathcal{P}^{\prime}$, the number of items needed to be picked is the
additional area divided by $\varepsilon$ rounded up, that is,

$\displaystyle\left\lceil\frac{\left\lceil\sum_{i\in
B_{t}}\frac{s_{i}}{\varepsilon-\varepsilon^{2}}\right\rceil\varepsilon-\sum_{i\in
B_{t}}s_{i}}{\varepsilon}\right\rceil$ $\displaystyle\leq
1+\frac{\varepsilon+\varepsilon\sum_{i\in
B_{t}}\frac{s_{i}}{\varepsilon-\varepsilon^{2}}-\sum_{i\in
B_{t}}s_{i}}{\varepsilon}$ $\displaystyle=2+\frac{1}{1-\varepsilon}\sum_{i\in
B_{t}}s_{i}$ $\displaystyle\leq 2+2\sum_{i\in B_{t}}s_{i}$

where the last inequality follows from the fact that $\varepsilon\leq 1/2$.
Let $R$ be the set of all picked items. Observe that ${|\mathcal{P}^{*}|\leq
d\operatorname{\mbox{OPT}}(I)}$, since each bin is not empty for at least one
scenario. Thus, as there are $2^{d}$ types and
${\sum_{B\in\mathcal{P}^{*}}\sum_{i\in
B_{t}}s_{i}\leq\operatorname{\mbox{OPT}}(I)}$, the number of items in $R$ is
at most

$\displaystyle\sum_{t\in\mathcal{T}}\sum_{B\in\mathcal{P}^{*}}\left(2+2\sum_{i\in
B_{t}}s_{i}\right)$
$\displaystyle=2^{d+1}|\mathcal{P}^{*}|+2^{d+1}\operatorname{\mbox{OPT}}(I)$
$\displaystyle=2^{d+1}d\operatorname{\mbox{OPT}}(I)+2^{d+1}\operatorname{\mbox{OPT}}(I)$
$\displaystyle=2^{d+1}(d+1)\operatorname{\mbox{OPT}}(I).$

Since $1/\varepsilon$ is an integer, one can rearrange all items of $R$ into
at most

$\left\lceil\frac{2^{d+1}(d+1)\operatorname{\mbox{OPT}}(I)}{1/{\varepsilon}}\right\rceil$

new bins of unit capacity. Let $\mathcal{R}$ be this set of new bins, and
obtain a feasible packing $\mathcal{P}^{\prime\prime}$ of $I$, by making a
copy of $\mathcal{P}^{\prime}$, removing the items of $R$, and adding the bins
of $\mathcal{R}$.

Notice that
$\textit{V}_{\textit{BPPS}}(\mathcal{P}^{\prime})=\textit{V}_{\textit{BPPS}}(\mathcal{P}^{*})=\operatorname{\mbox{OPT}}(I)$.
Using these facts, the value of $\mathcal{P}^{\prime\prime}$ can be estimated
as

$\displaystyle\textit{V}_{\textit{BPPS}}(\mathcal{P}^{\prime\prime})$
$\displaystyle\leq\textit{V}_{\textit{BPPS}}(\mathcal{P}^{\prime})+|\mathcal{R}|$
$\displaystyle=\operatorname{\mbox{OPT}}(I)+\left\lceil\frac{2^{d+1}(d+1)\operatorname{\mbox{OPT}}(I)}{1/{\varepsilon}}\right\rceil$
$\displaystyle=\operatorname{\mbox{OPT}}(I)+2^{d+1}(d+1)\varepsilon\operatorname{\mbox{OPT}}(I)+1.$

We notice that $\mathcal{P}^{\prime\prime}$ contains all large items
$\mathcal{L}$ and, for each type $t$, a certain number of items with scenarios
$t$ and size $\varepsilon$. Thus, to obtain a packing for $\hat{\mathcal{I}}$
from $\mathcal{P}^{\prime\prime}$, it is sufficient to show that for each $t$,
the number of created items corresponding to $t$ is not smaller than
$|\hat{\mathcal{S}}_{t}|$. Indeed, for a given $t$, the number of created
items is

$\sum_{B\in\mathcal{P}^{*}}\left\lceil\sum_{i\in
B_{t}}s_{i}/\varepsilon\right\rceil\geq\left\lceil\sum_{i\in\mathcal{S}_{t}}s_{i}/\varepsilon\right\rceil=|\hat{\mathcal{S}}_{t}|.$

To conclude the proof, we create a packing $\hat{\mathcal{P}}$ for
$\hat{\mathcal{I}}$ by removing the exceeding items of
$\mathcal{P}^{\prime\prime}$. Since removing such items can only decrease the
objective value,
${\textit{V}_{\textit{BPPS}}(\hat{\mathcal{P}})\leq\textit{V}_{\textit{BPPS}}(\mathcal{P}^{\prime\prime})}$.∎

Finally, to obtain a packing for $I$, one can first obtain a packing of
$\hat{I}$, then replace items from $\hat{\mathcal{S}}_{t}$ by items in
$\mathcal{S}_{t}$. This leads to Theorem 3.

######  Theorem 3.

For every constant $\varepsilon^{\prime}>0$, one can find in polynomial time a
packing $\mathcal{P}$ of $\mathcal{I}$ such that
$\textit{V}_{\textit{BPPS}}(\mathcal{P})\leq(1+\varepsilon^{\prime})\operatorname{\mbox{OPT}}(I)+2+2^{d}$.

###### Proof.

Define $r=2^{d+1}(d+1)$ and define $\varepsilon$ as the largest number such
that $\varepsilon\leq\min(\frac{\varepsilon^{\prime}}{3r},1/4)$ and
$1/\varepsilon$ is an integer. Let $\hat{I}$ be the instance obtained from $I$
by replacing, for each type, the items smaller than $\varepsilon^{2}$ with
items of size $\varepsilon$, as previously mentioned. Since each item of
$\hat{I}$ has size at least $\varepsilon^{2}$, by Lemma 7, we obtain a packing
$\hat{\mathcal{P}}$ for $\hat{I}$ with value
$\textit{V}_{\textit{BPPS}}(\hat{\mathcal{P}})\leq{(1+\varepsilon)\operatorname{\mbox{OPT}}(\hat{I})}+2^{d}$
in polynomial time.

For each $t\in\mathcal{T}$, we place all items of $\mathcal{S}_{t}$ over the
items of $\hat{\mathcal{S}}_{t}$. Precisely, we do the following: start by
making a copy $\mathcal{P}$ of $\hat{\mathcal{P}}$; select a set $A$ of
unpacked items in $\mathcal{S}_{t}$ whose total size is at least
$\varepsilon-\varepsilon^{2}$ and at most $\varepsilon$. If there is not such
a set $A$, then let $A$ be the remaining set of unpacked items in
$\mathcal{S}_{t}$. Find a bin ${B\in\mathcal{P}}$ with an item
$j\in\hat{\mathcal{S}}_{t}\cap B$ and replace $j$ by the items of $A$. Since
there are
$\lceil\sum_{i\in\mathcal{S}_{t}}s_{i}/(\varepsilon-\varepsilon^{2})\rceil$
items in $\hat{\mathcal{S}}_{t}$, and each group $A$, except perhaps one of
them, has size at least $\varepsilon-\varepsilon^{2}$, we always find an
unused item $j\in\hat{\mathcal{S}}_{t}$; repeat these steps while there are
unpacked items in $\mathcal{S}_{t}$ for some $t$.

The resulting packing $\mathcal{P}$ contains all items in $\mathcal{I}$. As no
bin was added, $\mathcal{P}$ has the same value of $\hat{\mathcal{P}}$. Using
Lemma 8, we finally obtain

$\displaystyle\textit{V}_{\textit{BPPS}}(P)$
$\displaystyle=\textit{V}_{\textit{BPPS}}(\hat{P})\leq(1+\varepsilon)\operatorname{\mbox{OPT}}(\hat{I})+2^{d}$
$\displaystyle\leq(1+\varepsilon)[(1+2^{d+1}(d+1)\varepsilon)\operatorname{\mbox{OPT}}(I)+1]+2^{d}$
$\displaystyle=(1+r\varepsilon^{2}+r\varepsilon+\varepsilon)\operatorname{\mbox{OPT}}(I)+2^{d}+1+\varepsilon$
$\displaystyle\leq(1+3r\varepsilon)\operatorname{\mbox{OPT}}(I)+2^{d}+1+\varepsilon$
$\displaystyle\leq(1+\varepsilon^{\prime})\operatorname{\mbox{OPT}}(I)+2+2^{d}.$
∎

##  5 A Variable Neighborhood Search Algorithm

The Variable Neighborhood Search is a metaheuristic that searches in different
neighborhood structures and has systematic change of neighborhood to find
better solutions and escape from local optima. Mladenović and Hansen (1997)
proposed the VNS not only for combinatorial optimization problems, but for
optimization in general. Due to its high applicability to similar problems
(Puchinger and Raidl, 2008; Fleszar and Hindi, 2002; Beltrán et al., 2004), we
apply the VNS to solve the BPPS.

For successfully implement it, we need to clarify three key characteristics of
the problem: the local search neighborhood structures and how to navigate the
search space; sensitive and computationally viable objective functions; and
the stopping criterion. In the following, we detail each of these and describe
the resulting algorithm.

###  5.1 Neighborhood structures

A neighborhood structure specifies a well-defined way to get from any given
solution into another solution that is “close” to the first one. We refer to
these neighborhood structures as a “movement” that takes one solution as input
and, depending on some structure-dependent input parameters, it returns a
different but “close” solution. We define four neighborhood structures for the
BPPS in the following order:

  * •

$N_{1}$: Given items $i$ and $j$ packed into different bins, move $i$ to the
bin currently containing $j$, then move $j$ to the bin that contained $i$;

  * •

$N_{2}$: Given an item $i$ and a bin $B$ that currently does not contain $i$,
move $i$ to $B$;

  * •

$N_{3}$: Given bins $B_{1},B_{2}$, and $B_{3}$ and items $i$ and $j$ packed
into $B_{1}$, move these items to $B_{2}$ and $B_{3}$, respectively;

  * •

$N_{4}$: Given a bin $B$, remove it from the solution and repack the items in
the first bin that can accommodate each of them, respecting the order that
they appeared in $B$.

Neighborhoods $N_{1},N_{2}$, and $N_{3}$ are prioritized according to their
complexity. Although the neighborhood structure $N_{4}$ is the smallest, it
can take comparatively more computational effort to compute since we have to
repack all items of a given bin, so we consider it as the last one.

###  5.2 Objective function

The objective function for the BPPS may not be sensitive enough to rank two
neighbor solutions. Thus, to guide the VNS, we define a fitness function that
penalizes solutions in which the bin occupation per scenario is small.

Let $\mathcal{B}$ be a solution for the BPPS and $\mathcal{B}_{k}$ be the
subset of bins that contain items in scenario $k$, where $|\mathcal{B}_{k}|$
is the number of bins used in $k$. Furthermore, let $\text{ocp}(k,B)$ be the
total size of the items in bin $B$ from scenario $k$. Given a solution
$\mathcal{B}$, the VNS moves to a neighbor solution $\mathcal{B}^{\prime}$
such that $f(\mathcal{B})>f(\mathcal{B}^{\prime})$, for $f$ defined as:

$\displaystyle f(\mathcal{B})$
$\displaystyle=\textit{V}_{\textit{BPPS}}(\mathcal{B})-\sum_{k\in\mathcal{K}}\sum_{B\in\mathcal{B}_{k}}{\left(\frac{\text{ocp}(k,B)}{|\mathcal{B}_{k}|W}\right)}^{2}$

###  5.3 Stopping criterion

Since the VNS cannot tell whether an optimal solution is reached, we use two
stopping criteria to decide when the algorithm stops. The first criterion is
the _timeout_ , where a preset CPU time limit ($t_{\max}$) is imposed, so the
VNS stops as soon as the time limit is exceeded. The second criterion is the
_local-convergence_ , where a preset number of iterations ($c_{\max}$) is
imposed, so the VNS stops if the best known solution is not updated within
this number of iterations.

###  5.4 Algorithm

The VNS takes as input an initial solution $x$, a number of neighborhood
structures $N_{\max}=4$, a time limit $t_{\max}=1800$ seconds and a
convergence limit $c_{\max}=500$ iterations (these values were achieved after
preliminary computational tests). The initial solution $x$ is obtained as
follows: items are sorted in non-increasing order of their sizes and then
packed in the first open bin where it fits, considering the capacity
constraints in every scenario and opening a new bin whenever necessary.

The overall structure of the proposed VNS is given in the following algorithm:
$\proc{VNS}(x,N_{\max},t_{\max},c_{\max})$ $t\leftarrow 0$, $c\leftarrow 0$
$t<t_{\max}$ and $c<c_{\max}$ $\kappa\leftarrow 1$
$improvement\leftarrow\text{False}$ Assume there will be no improvement
$x^{\prime}\leftarrow\proc{Shake}(x,\kappa)$ Get a random solution in
$N_{\kappa}(x)$ $x^{\prime\prime}\leftarrow\proc{Local-
Search}(x^{\prime},\mathcal{N}_{\max})$ Perform VND to improve $x^{\prime}$
$f(x^{\prime\prime})<f(x)$ There has been an improvement $x\leftarrow
x^{\prime\prime}$ $\kappa\leftarrow 1$ $improvement\leftarrow\text{True}$
Local optima, change neighborhood $\kappa\leftarrow\kappa+1$ $\kappa\isequal
N_{\max}$ $improvement$ There has been an improvement $c\leftarrow 0$ No
improvement in this iteration $c\leftarrow c+1$ $t\leftarrow\proc{CPU-Time}()$
Update processing time.

The VNS uses a Shake procedure that simply returns a solution $x^{\prime}$
randomly selected from the neighborhood $N_{\kappa}$ of $x$. The Local-Search
procedure is based on the Variable Neighborhood Descent (VND) approach, which
is the deterministic descent-only version of the VNS. Our local search takes
as input a solution $x$ and the number of neighborhood structures
$N_{\max}=4$, resulting in the following algorithm:  $\proc{Local-
Search}(x,N_{\max})$ $\kappa\leftarrow 1$ $x^{\prime}\leftarrow\arg\min_{y\in
N_{\kappa}(x)}f(y)$ Find the first best neighbor of $x$ $f(x^{\prime})<f(x)$
There has been an improvement $x\leftarrow x^{\prime}$ $\kappa\leftarrow 1$
Local optima, change neighborhood $\kappa\leftarrow\kappa+1$ $\kappa\isequal
N_{\max}$

##  6 A Branch-and-Price Algorithm

We propose an exact method for the BPPS based on column generation and branch-
and-bound, that is, a branch-and-price algorithm. This method consists of
modelling the problem as an ILP with a very large (usually exponential) number
of variables, which are generated and added to the model only as they are
needed, in each node of the branch-and-bound tree. Branch-and-price has been
used successfully to solve practical instances of the BPP and several of its
variants (Delorme et al., 2016b; Sadykov and Vanderbeck, 2013; Dell’Amico et
al., 2020; Borges et al., 2020; Bettinelli et al., 2010).

One of the most successful formulations for the classical BPP and variants is
based on associating a variable to each cutting pattern (see, e.g., Delorme et
al. (2016b)). In practice, such kinds of formulations lead to models with an
exponential number of variables, as the number of cutting patterns is usually
exponential with respect to the number of items. Due to this, to solve
practical instances with such models, one has to rely on column generation to
solve the corresponding linear relaxation.

The _column generation method_ , proposed by Ford and Fulkerson (1958) and
generalized by Dantzig and Wolfe (1961), solves the linear relaxation of
models with a large number of variables. Gilmore and Gomory (1961, 1963) were
the first to present practical experiments based on column generation, by
solving the linear relaxation of an exponential model for the BPP. The column
generation method solves the linear relaxation by iteratively solving the
_restricted master problem_ (RMP), which considers only a subset of variables
(columns) of the original model. At each iteration, the method has to
determine the existence of non-basic columns which are candidate to improve
the current solution of the RMP, i.e., columns with negative reduced cost (on
minimization problems). If the method determines that no such column exists,
then the solution of the RMP is an optimal solution for the original linear
relaxation, providing a bound for the integer problem.

###  6.1 A pattern-based model

Recalling that $\mathcal{P}$ is the set of all cutting patterns for the BPPS,
the following exponential ILP model solves the BPPS:

$\displaystyle\mathrm{minimize}\ \ \ $ $\displaystyle\mathcal{F},$ (1)
$\displaystyle\mathrm{s.t.:}\ \ \ $
$\displaystyle\displaystyle\sum_{p\in\mathcal{P}}a_{ip}X_{p}\geq 1,$
$\displaystyle\forall i\in\mathcal{I};$ (2)
$\displaystyle\displaystyle\sum_{p\in\mathcal{P}}b_{kp}X_{p}\leq\mathcal{F},$
$\displaystyle\forall k\in\mathcal{K};$ (3) $\displaystyle X_{p}\in\\{0,1\\},$
$\displaystyle\forall p\in\mathcal{P};$ (4)
$\displaystyle\mathcal{F}\in\mathbb{Z}_{+}.$ (5)

For each pattern $p\in\mathcal{P}$, binary variable $X_{p}$ is equal to $1$ if
and only if $p$ is used in the solution. The integer variable $\mathcal{F}$
represents an upper bound on the number of bins (patterns) that is used by
each individual scenario. The objective function (1) is to minimize
$\mathcal{F}$. Constraints (2) guarantees that every item belongs to at least
one cutting pattern of the solution. Constraints (3) ensures that, for every
scenario $k\in\mathcal{K}$, the number of bins that packs at least one item
from $k$ is at most $\mathcal{F}$. The domain of variables is in (4) and (5).

Next, we describe a column generation algorithm to solve the linear relaxation
of (1)–(5). As the initial set of columns for the RMP, we consider each
cutting pattern with a single item. The dual solution of the RMP consists of a
value $\alpha_{i}\geq 0$ for each item $i$ associated with a constraint in
(2), and a value $\beta_{k}\leq 0$ for each scenario $k$ associated with a
constraint in (3). The dual solution is used to solve the pricing problem,
which determines a column with minimum reduced cost, that is

$\min_{p\in\mathcal{P}}\left\\{c_{p}-\left(\sum_{i\in\mathcal{I}}\alpha_{i}a_{ip}+\sum_{k\in
K}\beta_{k}b_{kp}\right)\right\\},$

which is equivalent to

$\max_{p\in\mathcal{P}}\left\\{\sum_{i\in\mathcal{I}}\alpha_{i}a_{ip}+\sum_{k\in
K}\beta_{k}b_{kp}\right\\}.$ (6)

###  6.2 The pricing problem

In order to generate a new column, we need to find a feasible pattern that
optimizes (6). Since it includes a term for each item and a term for each
scenario used in the pattern, we can describe the pricing problem as the
following Knapsack Problem with Scenarios.

#### Knapsack Problem with Scenarios (KPS)

In the KPS, we are given: a number $d$ of scenarios; a set $\mathcal{I}=[n]$
of $n$ items, with each item $i\in\mathcal{I}$ having a size $s_{i}$, item
value $v_{i}\in\mathbb{Q}$, and a set of scenarios
$\mathcal{K}_{i}\subseteq[d]$; a vector $u\in\mathbb{Q}^{d}$ of scenario
values, and a knapsack of capacity $W\in\mathbb{Q}_{+}$. For simplicity, for
each item $i\in\mathcal{I}$, define:

$s_{i}^{k}=\begin{cases}s_{i},\text{ if item }i\text{ is in the scenario
}k;\\\ 0,\text{ otherwise}.\end{cases}$

A solution is a subset $B$ of $\mathcal{I}$ such that for each scenario
$k\in\mathcal{K}$, $\sum_{i\in B}s_{i}^{k}\leq w_{k}$. The objective is to
find a solution that maximizes

$\sum_{i\in B}v_{i}+\sum_{k:s_{i}^{k}>0,\,i\in B}u_{k}.$

######  Theorem 4.

KPS is NP-hard in the strong sense.

###### Proof.

Let us recall the Maximum Independent Set problem, which given a graph $G$
consists of finding a subset of maximum cardinality of its vertices such that
no two of its elements are adjacent in $G$. This problem is known to be
strongly NP-hard (Garey and Johnson, 1978), so all we need to do is present a
polynomial reduction from this problem into the KPS.

Let $I=G(V,E)$ be an instance of the Maximum Independent Set problem. We
construct an instance $I^{\prime}=(d,\mathcal{I},s,v,\mathcal{K},u,W)$ of the
KPS such that $|\mathcal{I}|=|V|$, $d=|V|$. For each item $i\in\mathcal{I}$
there is a scenario $k_{i}\in[d]$ associated with it. We define
$\mathcal{K}_{i}=\\{k_{i}\\}\cup\\{k_{j}\colon(i,j)\in E\\}$ and we let
$s_{i}^{k}=1$ if $k\in\mathcal{K}_{i}$. We define $v_{i}=1$ for all
$i\in\mathcal{I}$, $u_{k}=0$ for all $k\in[d]$ and $W=1$.

Now notice that there is a bijection between solutions of $I$ and $I^{\prime}$
which preserves the value of the objective function two items corresponding to
two adjacent vertices of $V$ cannot be in the same solution of $I^{\prime}$,
the items are unit-valued and the scenarios are zero-valued.∎∎

As a consequence of Theorem 4, our pricing subproblem does not admit a pseudo-
polynomial time algorithm, unless $P=\mbox{NP}$, thus, we propose the
following ILP model to solve it.

We can use the dual solution $(\alpha,\beta)$ of the RMP as item-values and
scenario-values $(u,v)$ to describe the following model for the pricing
subproblem

$\displaystyle\mathrm{maximize}\ \ \ $
$\displaystyle\sum_{i\in\mathcal{I}}\alpha_{i}A_{i}+\sum_{k\in\mathcal{K}}\beta_{k}B_{k}$
(7) $\displaystyle\mathrm{s.t.:}\ \ \ $
$\displaystyle\displaystyle\sum_{i\in\mathcal{I}}s_{i}^{k}A_{i}\leq WB_{k},$
$\displaystyle\forall k\in\mathcal{K};$ (8) $\displaystyle
A_{i},B_{k}\in\\{0,1\\},$ $\displaystyle\forall
i\in\mathcal{I},k\in\mathcal{K}.$ (9)

In model (7)-(9), variables $A_{i}$ are related to the coefficients $a_{ip}$,
and variables $B_{k}$ are related to the coefficients $b_{kp}$. Hence,
variables $A$ and $B$ determine the cutting pattern that is generated. The
objective function (7) is related to the pricing formula in (6). Constraints
(8) assures that the solution of the model represents a valid cutting pattern
for the BPPS. Notice that $W$ is the capacity of a bin, then if a scenario
$k\in\mathcal{K}$ is used ($B_{k}=1$), the sum of the sizes of all items that
belong to the scenario $k$ must be at most $W$, and if this scenario is not
used ($B_{k}=0$), then the sum of the sizes of all items that belong to $k$
must be equal to $0$.

###  6.3 The branch-and-price algorithm

To solve the model (1)-(5), we propose a branch-and-price algorithm. The
branch-and-price is based on an enumeration of the fractional solution of the
column generation algorithm. For this end, we use the branching scheme of Ryan
and Foster (1981). Given a fractional solution $X^{*}$, we know, from Vance et
al. (1994), that there exists rows $l$ and $m$ such that

$0<\sum_{p:a_{lp}=1,a_{mp}=1}X_{p}<1.$

We use this pair of rows to create two branches for the current node. On one
side we enforce that items $l$ and $m$ must be packed in the same pattern, and
on the other side we enforce that they must be packed in different patterns.
This is achieved by the following branching constraints

$\displaystyle\sum_{p:a_{lp}=1,a_{mp}=1}X_{p}\geq 1$ (10)
$\displaystyle\sum_{p:a_{lp}=1,a_{mp}=1}X_{p}\leq 1$ (11)

Rather than explicitly adding these constraints to the RMP, they are enforced
by setting the infeasible columns’ upper bound to zero. On the branch of (10),
this is equivalent to combine items $l$ and $m$. For the branch of (11), this
is equivalent to making rows $l$ and $m$ disjoint. However, this is not
enough, since we must also avoid infeasible columns to be re-generated on each
branch, which can be achieved by modifying the underlining pricing subproblem
as follows:

  * •

For the branches in which items $l$ and $m$ are combined, it suffices to
modify the set of items to
$\mathcal{\bar{I}}=\\{1,\dots,n\\}\setminus\\{l,m\\}\cup\\{n+1\\}$, with
$\alpha_{n+1}=\alpha_{l}+\alpha_{m}$,

$\displaystyle s_{n+1}^{k}=\begin{cases}s_{l}^{k},&\text{ if
}s_{l}^{k}>0\text{ and }s_{m}^{k}=0,\\\ s_{m}^{k},&\text{ if
}s_{l}^{k}=0\text{ and }s_{m}^{k}>0,\\\ s_{l}^{k}+s_{m}^{k},&\text{ if
}s_{l}^{k}>0\text{ and }s_{m}^{k}>0,\\\ 0,&\text{ otherwise;}\end{cases}$ (12)

  * •

For branches in which items $l$ and $m$ must be packed in different patterns,
we simply add the following disjunction constraint

$A_{l}+A_{m}\leq 1.$ (13)

These changes can be easily accommodated into the model (7)-(9) without much
effort.

To obtain an upper bound at each node, we solve the current restricted master
problem with the integrality constraint. In other words, we solve an ILP model
that consists only of the columns that have been already generated. The use of
this method at each node can be computationally expensive, but it usually
provides tight upper bounds.

Some instances of the problem may be too complex to be solved by the column
generation algorithm in practical time. In this manner, we would need faster
methods for obtaining lower bounds for such instances. Two methods are
proposed to get a lower bound at the root node of the branch-and-price
algorithm. One of them is based on the continuous lower bound for the bin
packing problem, and it is given by the following equation:

$\textit{LB}_{\textit{CON}}=\max_{k\in\mathcal{K}}\left\\{\left\lceil\sum_{i\in
S_{k}}s_{i}^{k}/W\right\rceil\right\\}.$ (14)

The idea behind the continuous lower bound (14) is to break items into
unitary-sized items to fill empty spaces in a packing, which could not be
filled otherwise. The $LB_{CON}$ computes the continuous lower bound for every
scenario, and returns the largest value between them.

The other lower bound that we use is based on dual feasible functions. A _dual
feasible function_ (DFF) is a function $f$ such that, for any finite set $X$
of real numbers, if $\sum_{x\in X}x\leq 1$, then $\sum_{x\in X}f(x)\leq 1$.
For a survey on DFFs, we refer to Clautiaux et al. (2010). Given a DFF $f$, we
can construct a lower bound based on it, similarly to the continuous lower
bound, as given in the following equation:

$\textit{LB}_{\textit{DFF}}=\max_{k\in K}\\{\lceil\sum_{i\in
S_{k}}f(s_{i}^{k}/W)\rceil\\}$ (15)

In the proposed branch-and-price algorithm, given
$\lambda\in\\{1,\ldots,W/2\\}$, we consider the DFF proposed by Fekete and
Schepers (2001) and described in (16):

$f(s)=\begin{cases}W,\text{ if }s>W-\lambda;\\\ 0,\text{ if }s\leq\lambda;\\\
s,\text{ otherwise.}\end{cases}$ (16)

##  7 Numerical Experiments

The efficiency of the proposed VNS and B&P algorithm is evaluated with
computational experiments. Both solutions methods were implemented in the C++
programming language, with the B&P algorithm using the Gurobi Solver version
8.1 Gurobi Optimization (2023) to solve the linear programming models. The
experiments were carried out in a computer with an Intel(R) Xeon(R) CPU
E5–2630 v4 processor at 2.20 GHz, with 64 GB of RAM and under the Ubuntu 18.04
operating system.

A total of 120 instances were generated, divided into 12 different classes of
10 instances each. The classes are combinations of the number of items
$n\in\\{10,50,100,200\\}$ and the number of scenarios $d\in\\{0.5n,n,2n\\}$,
which depends on the number of items of the class. For every instance, the
size of the items is randomly generated over the set $\\{1,\ldots,99\\}$ under
a uniform probability distribution, whereas the size of the bin is fixed to
$100$. The scenarios were generated randomly, similarly to Bódis and Balogh
(2018), such that an item $i\in\mathcal{I}$ belongs to the scenario
$k\in\mathcal{K}$ with probability $0.5$.

A summary of the computational results is presented in Table 1. Columns $n$
and $d$ present, respectively, the number of items and the number of scenarios
of each class of instances. We consider three sets of columns, each
representing one of the tested algorithms, with the last one, VNS+B&P, being
the B&P algorithm using the solution obtained from the VNS as a warm-start.
The columns $gap$ represent the average optimality gap that each algorithm
achieved in each class, computed as $\frac{UB-LB}{UB}$, where $UB$ and $LB$
are respectively the upper and lower bounds obtained. Notice that since the
VNS does not produce a lower bound on its own, we used the lower bound from
the VNS+B&P to compute its optimality gap. Columns $time$ represent the
average time in seconds that each algorithm spent to solve each instance of
the corresponding class, being $1800$ seconds the time limit for the VNS and
VNS + B&P algorithms, and a time limit of $3600$ for B&P. Columns $|OPT|$
represent the number of instances from each class that were solved to
optimality (out of 10) by each algorithm. Columns $cols$ and $nodes$ represent
the average number of columns generated and the number of explored nodes,
respectively, considering the B&P algorithm.

|  |  | VNS |  | B&P |  | VNS + B&P
---|---|---|---|---|---|---|---
$n$ | $d$ |  | $gap$ | $time$ | $|\text{OPT}|$ |  | $gap$ | $nodes$ | $cols$ | $time$ | $|\text{OPT}|$ |  | $gap$ | $nodes$ | $cols$ | $time$ | $|\text{OPT}|$
10 | 5 |  | 0.00% | 0.00 | 10 |  | 0.00% | 1.0 | 6.9 | 0.02 | 10 |  | 0.00% | 0 | 0 | 0.00 | 10
10 | 10 |  | 0.00% | 0.00 | 10 |  | 0.00% | 1.6 | 14.0 | 0.03 | 10 |  | 0.00% | 0.3 | 5.4 | 0.00 | 10
10 | 20 |  | 1.67% | 0.00 | 9 |  | 0.00% | 1.2 | 11.5 | 0.03 | 10 |  | 0.00% | 0.3 | 3.4 | 0.00 | 10
50 | 25 |  | 7.11% | 14.80 | 1 |  | 0.00% | 6.4 | 156.1 | 11.21 | 10 |  | 1.01% | 955.8 | 98256.4 | 364.99 | 8
50 | 50 |  | 5.24% | 28.80 | 1 |  | 0.00% | 15.7 | 289.8 | 30.08 | 10 |  | 0.50% | 321.8 | 34501.3 | 187.00 | 9
50 | 100 |  | 7.26% | 28.80 | 0 |  | 0.00% | 270.2 | 2401.1 | 569.41 | 10 |  | 0.93% | 129.8 | 10871.7 | 370.18 | 8
100 | 50 |  | 9.34% | 134.30 | 0 |  | 1.44% | 80.8 | 2381.6 | 1832.76 | 5 |  | 0.56% | 6.6 | 4575.5 | 516.59 | 8
100 | 100 |  | 8.25% | 181.20 | 0 |  | 2.21% | 85.3 | 2399.6 | 3019.09 | 2 |  | 1.38% | 29.2 | 12419.1 | 1193.98 | 5
100 | 200 |  | 7.78% | 240.40 | 0 |  | 1.51% | 41.4 | 1028.1 | 2710.84 | 4 |  | 1.72% | 5.8 | 3577.8 | 1052.04 | 6
200 | 100 |  | 16.32% | 1048.80 | 0 |  | 13.74% | 1.0 | 532.7 | 3814.57 | 0 |  | 20.13% | 1 | 3323.1 | 1817.98 | 0
200 | 200 |  | 18.56% | 1407.70 | 0 |  | 17.52% | 1.0 | 460.3 | 3925.27 | 0 |  | 20.07% | 1 | 4801.9 | 1828.71 | 0
200 | 400 |  | 18.31% | 1595.30 | 0 |  | 20.16% | 1.0 | 362.7 | 3866.15 | 0 |  | 20.29% | 1 | 4832.4 | 1817.93 | 0
Average |  | 8.32% | 390.01 | 2.58 |  | 4.71% | 42.21 | 837.0 | 1648.28 | 5.91 |  | 5.55% | 121.05 | 14764 | 762.45 | 6.16
Table 1: Summary of the computational results for all the algorithms tested.

First, we discuss the performance of the VNS algorithm. From Table 1, we can
notice that most instances with $10$ items are easily solved by the heuristic.
However, for instances of size $50$ and $100$ we can notice a significant
larger gap, even though the algorithm time never reaches its time limit. This
indicates that the VNS quickly found a local-optima and was not able to escape
it, thus prematurely converging. Something similar can be observed for the
instances of size $200$, although the running time gets closer to the time
limit in these cases. For this algorithm, only two instances of size $50,100$,
and $200$ are solved to optimality, achieving a total average gap of $8.32$%
in an average time of $390.01$ seconds.

Observing the results for algorithm B&P, we can notice a significant
improvement. All instances of $10$ and $50$ items are solved to optimality
very quickly (100 seconds on average). We start to notice a drop in
performance for the instances of size $100$ and $200$. For instances of size
$100$, the algorithm was able to solve $11$ instances to optimality, out of
$30$, averaging over $1800$ seconds for all number of scenarios. We can also
observe that the average gap for these instances is very small (around $2$%
for all number of scenarios). The algorithm times out for all instances of
size $200$, with an average gap of 14%. Interestingly, the average gaps
achieved by the B&P for these instances are not far from the ones obtained by
the VNS, with the latter having half the time limit of the former. The total
average gap for this algorithm is $4.71$% in an average time of $1648.28$
seconds.

Finally, we analyze the VNS + B&P algorithm, that is, B&P using the VNS
solution as a warm-start. Once again, all instances of size $10$ are
completely solved in negligible time. Differing from the pure B&P algorithm
though, this one could not find optimal solutions for a few instances of size
$50$. These $5$ instances brought all the averages for this group of instances
up, since they were tried until the time limit of $1800$ seconds, totaling
thousands of columns and hundreds of nodes more than their non-warm-started
counterpart. This might happen because the solution of the VNS algorithm might
impact the branching step of the B&P, making it harder for the algorithm to
converge depending on the branch it takes initially. However, we can notice
that for instances of size $100$, the VNS + B&P solved $19$ instances to
optimality, $8$ more than the pure B&P and with a significant smaller average
time. This algorithm also timed out for all instances of size $200$, while
presenting a greater gap in average compared to the pure B&P because of the
smaller time limit. When compared to the pure B&P algorithm, the VNS + B&P
gives an average gap greater by $0.16$% while using less than half the time in
average.

The data from Table 1 shows us that although the VNS algorithm struggles to
avoid local-optima, its results serve as good upper bounds to be fed into a
more elaborated model such as the B&P presented. By using the VNS’s result as
a warm-start, the B&P algorithm was able to solve more instances in less time,
this was not enough to tackle instances with $200$ items.

##  8 Concluding Remarks

This work deals with a variant of the bin packing problem where items occupy
the bin capacity when they belong to the same scenario. For this problem, a
variety of solution methods are proposed, including an absolute approximation
algorithm, an asymptotic polynomial time approximation scheme, a variable
neighborhood search based heuristic, and a branch-and-price that solves a set-
cover-based formulation. Results of the heuristic and the branch-and-price
illustrate how effective is the latter to solve small and medium size
instances, having its efficiency improved when the heuristic provides a start
solutions.

While the approximation algorithms depend either on an approximation algorithm
for the vector bin packing problem, for the absolute approximation algorithm,
or the number of scenarios, for the APTAS, the heuristic and branch-and-price
can be used to handle the problem in practical contexts. The heuristic is not
so competitive with the branch-and-price algorithm, since the former only
optimally solved 26% of the instances, while the latter obtained optimal
solutions for 59% of the instances. If solutions of the heuristic are used as
a warm-start for the branch-and-price, the number of optimal solutions
increases to about 62%. In total, we were able to prove optimality for $79$
out of the $120$ proposed instances.

Future research can focus on new approximation algorithms that explore the
influence of the scenarios. Improvements in the variable neighborhood search
are also expected, including new neighborhood structures and local searches.
Regarding the branch-and-price, we will work on valid cuts and the possibility
of dynamic programming algorithms for the pricing problem.

## Compliance with Ethical Standards

The authors acknowledge the financial support of the National Counsel of
Technological and Scientific Development (CNPq grants 144257/2019–0,
312186/2020–7, 311185/2020–7, 311039/2020–0, 405369/2021–2, 313146/2022–5),
the State of Goiás Research Foundation (FAPEG), and the State of São Paulo
Research Foundation (FAPESP grant 2017/11831–1). This study was financed in
part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -
Brasil (CAPES) - Finance Code 001.

## Ethical approval

This article does not contain any studies with human participants or animals
performed by any of the authors.

## Data availability

The datasets generated during and/or analysed during the current study are
available from the corresponding author on reasonable request.

## References

  * Alves et al. (2014) Alves, C., Valério de Carvalho, J., Clautiaux, F., and Rietz, J. (2014).  Multidimensional dual-feasible functions and fast lower bounds for the vector packing problem.  European Journal of Operational Research, 233(1):43–63.
  * Baghalian et al. (2013) Baghalian, A., Rezapour, S., and Farahani, R. Z. (2013).  Robust supply chain network design with service level against disruptions and demand uncertainties: A real-life case.  European journal of operational research, 227(1):199–215.
  * Baldi et al. (2019) Baldi, M. M., Manerba, D., Perboli, G., and Tadei, R. (2019).  A generalized bin packing problem for parcel delivery in last-mile logistics.  European Journal of Operational Research, 274(3):990–999.
  * Balogh et al. (2019) Balogh, J., Békési, J., Dósa, G., Sgall, J., and van Stee, R. (2019).  The optimal absolute ratio for online bin packing.  Journal of Computer and System Sciences, 102:1–17.
  * Bansal et al. (2009) Bansal, N., Caprara, A., and Sviridenko, M. (2009).  A new approximation method for set covering problems, with applications to multidimensional bin packing.  SIAM Journal on Computing, 39(4):1256–1278.
  * Beltrán et al. (2004) Beltrán, J. D., Calderón, J. E., Cabrera, R. J., Moreno-Pérez, J. A., and Moreno-Vega, J. M. (2004).  Grasp-vns hybrid for the strip packing problem.  Hybrid metaheuristics, 2004:79–90.
  * Bertazzi et al. (2019) Bertazzi, L., Golden, B., and Wang, X. (2019).  The bin packing problem with item fragmentation:a worst-case analysis.  Discrete Applied Mathematics, 261:63–77.
  * Bettinelli et al. (2010) Bettinelli, A., Ceselli, A., and Righini, G. (2010).  A branch-and-price algorithm for the variable size bin packing problem with minimum filling constraint.  Annals of Operations Research, 179(1):221–241.
  * Bódis and Balogh (2018) Bódis, A. and Balogh, J. (2018).  Bin packing problem with scenarios.  Central European Journal of Operations Research, pages 1–19.
  * Borges et al. (2020) Borges, Y. G., Miyazawa, F. K., Schouery, R. C., and Xavier, E. C. (2020).  Exact algorithms for class-constrained packing problems.  Computers & Industrial Engineering, 144:106455.
  * Brandão and Pedroso (2016) Brandão, F. and Pedroso, J. P. (2016).  Bin packing and related problems: General arc-flow formulation with graph compression.  Computers & Operations Research, 69:56–67.
  * Buljubašić and Vasquez (2016) Buljubašić, M. and Vasquez, M. (2016).  Consistent neighborhood search for one-dimensional bin packing and two-dimensional vector packing.  Computers & Operations Research, 76:12–21.
  * Caprara and Toth (2001) Caprara, A. and Toth, P. (2001).  Lower bounds and algorithms for the 2-dimensional vector packing problem.  Discrete Applied Mathematics, 111(3):231–262.
  * Chekuri and Khanna (2004) Chekuri, C. and Khanna, S. (2004).  On multidimensional packing problems.  SIAM journal on computing, 33(4):837–851.
  * Christensen et al. (2017) Christensen, H. I., Khan, A., Pokutta, S., and Tetali, P. (2017).  Approximation and online algorithms for multidimensional bin packing: A survey.  Computer Science Review, 24:63–79.
  * Clautiaux et al. (2010) Clautiaux, F., Alves, C., and Valério de Carvalho, J. (2010).  A survey of dual-feasible and superadditive functions.  Annals of Operations Research, 179(1):317–342.
  * Clautiaux et al. (2014) Clautiaux, F., Dell’Amico, M., Iori, M., and Khanafer, A. (2014).  Lower and upper bounds for the bin packing problem with fragile objects.  Discrete Applied Mathematics, 163:73–86.
  * Coffman et al. (2013) Coffman, E. G., Csirik, J., Galambos, G., Martello, S., and Vigo, D. (2013).  Bin packing approximation algorithms: survey and classification.  In Handbook of combinatorial optimization, pages 455–531. Springer.
  * Dahmani et al. (2013) Dahmani, N., Clautiaux, F., Krichen, S., and Talbi, E.-G. (2013).  Iterative approaches for solving a multi-objective 2-dimensional vector packing problem.  Computers & Industrial Engineering, 66(1):158–170.
  * Dantzig and Wolfe (1961) Dantzig, G. and Wolfe, P. (1961).  The decomposition algorithm for linear programs.  Econometrica, 29(4):767–778.
  * de Lima et al. (2023) de Lima, V. L., Iori, M., and Miyazawa, F. K. (2023).  Exact solution of network flow models with strong relaxations.  Mathematical Programming, 197(2):813–846.
  * Dell’Amico et al. (2020) Dell’Amico, M., Furini, F., and Iori, M. (2020).  A branch-and-price algorithm for the temporal bin packing problem.  Computers & Operations Research, 114:104825.
  * Dell’Amico et al. (2020) Dell’Amico, M., Furini, F., and Iori, M. (2020).  A branch-and-price algorithm for the temporal bin packing problem.  Computers & Operations Research, 114:104825.
  * Delorme and Iori (2019) Delorme, M. and Iori, M. (2019).  Enhanced pseudo-polynomial formulations for bin packing and cutting stock problems.  INFORMS Journal on Computing, 32(1):101–119.
  * Delorme et al. (2016a) Delorme, M., Iori, M., and Martello, S. (2016a).  Bin packing and cutting stock problems: Mathematical models and exact algorithms.  European Journal of Operational Research, 255(1):1–20.
  * Delorme et al. (2016b) Delorme, M., Iori, M., and Martello, S. (2016b).  Bin packing and cutting stock problems: Mathematical models and exact algorithms.  European Journal of Operational Research, 255(1):1–20.
  * Fekete and Schepers (2001) Fekete, S. P. and Schepers, J. (2001).  New classes of fast lower bounds for bin packing problems.  Mathematical Programming, 91(1):11–31.
  * Fernandez de la Vega and Lueker (1981) Fernandez de la Vega, W. and Lueker, G. S. (1981).  Bin packing can be solved within $1+\epsilon$ in linear time.  Combinatorica, 1(4):349–355.
  * Feuerstein et al. (2014) Feuerstein, E., Marchetti-Spaccamela, A., Schalekamp, F., Sitters, R., van der Ster, S., Stougie, L., and van Zuylen, A. (2014).  Scheduling over scenarios on two machines.  In International Computing and Combinatorics Conference, pages 559–571.
  * Fleszar and Hindi (2002) Fleszar, K. and Hindi, K. S. (2002).  New heuristics for one-dimensional bin-packing.  Computers & operations research, 29(7):821–839.
  * Ford and Fulkerson (1958) Ford, L. and Fulkerson, D. (1958).  A suggested computation for maximal multi-commodity network flows.  Management Science, 5(1):97–101.
  * Garey and Johnson (1978) Garey, M. R. and Johnson, D. S. (1978).  ”strong” np-completeness results: Motivation, examples, and implications.  Journal of the ACM (JACM), 25(3):499–508.
  * Gilmore and Gomory (1961) Gilmore, P. and Gomory, R. (1961).  A linear programming approach to the cutting stock problem.  Operations Research, 9:849–859.
  * Gilmore and Gomory (1963) Gilmore, P. and Gomory, R. (1963).  A linear programming approach to the cutting stock problem - part II.  Operations Research, 11:863–888.
  * Grange et al. (2018) Grange, A., Kacem, I., and Martin, S. (2018).  Algorithms for the bin packing problem with overlapping items.  Computers & Industrial Engineering, 115:331–341.
  * Gurobi Optimization (2023) Gurobi Optimization, L. (2023).  Gurobi optimizer reference manual.
  * Heßler et al. (2018) Heßler, K., Gschwind, T., and Irnich, S. (2018).  Stabilized branch-and-price algorithms for vector packing problems.  European Journal of Operational Research, 271(2):401–419.
  * Hu et al. (2017) Hu, Q., Zhu, W., Qin, H., and Lim, A. (2017).  A branch-and-price algorithm for the two-dimensional vector packing problem with piecewise linear cost function.  European Journal of Operational Research, 260(1):70–80.
  * Juan et al. (2018) Juan, A. A., Kelton, W. D., Currie, C. S., and Faulin, J. (2018).  Simheuristics applications: dealing with uncertainty in logistics, transportation, and other supply chain areas.  In 2018 winter simulation conference (WSC), pages 3048–3059.
  * Kantorovich (1960) Kantorovich, L. V. (1960).  Mathematical methods of organizing and planning prodution.  Management Science, 6:363–422.  (in Russian 1939).
  * Kucukyilmaz and Kiziloz (2018) Kucukyilmaz, T. and Kiziloz, H. E. (2018).  Cooperative parallel grouping genetic algorithm for the one-dimensional bin packing problem.  Computers & Industrial Engineering, 125:157–170.
  * Mladenović and Hansen (1997) Mladenović, N. and Hansen, P. (1997).  Variable neighborhood search.  Computers & operations research, 24(11):1097–1100.
  * Pereira (2016) Pereira, J. (2016).  Procedures for the bin packing problem with precedence constraints.  European Journal of Operational Research, 250(3):794–806.
  * Puchinger and Raidl (2008) Puchinger, J. and Raidl, G. R. (2008).  Bringing order into the neighborhoods: relaxation guided variable neighborhood search.  Journal of Heuristics, 14(5):457–472.
  * Ryan and Foster (1981) Ryan, D. M. and Foster, B. A. (1981).  An integer programming approach to scheduling.  Computer scheduling of public transport urban passenger vehicle and crew scheduling, pages 269–280.
  * Sadykov and Vanderbeck (2013) Sadykov, R. and Vanderbeck, F. (2013).  Bin packing with conflicts: a generic branch-and-price algorithm.  INFORMS Journal on Computing, 25(2):244–255.
  * Saint-Guillain et al. (2021) Saint-Guillain, M., Paquay, C., and Limbourg, S. (2021).  Time-dependent stochastic vehicle routing problem with random requests: Application to online police patrol management in brussels.  European Journal of Operational Research, 292(3):869–885.
  * Sweeney and Paternoster (1992) Sweeney, P. E. and Paternoster, E. R. (1992).  Cutting and packing problems: a categorized, application-orientated research bibliography.  Journal of the Operational Research Society, 43(7):691–706.
  * Vance et al. (1994) Vance, P. H., Barnhart, C., Johnson, E. L., and Nemhauser, G. L. (1994).  Solving binary cutting stock problems by column generation and branch-and-bound.  Computational optimization and applications, 3(2):111–130.
  * Vazirani (2001) Vazirani, V. (2001).  Approximation Algorithms.  Springer-Verlag.
  * Wäscher et al. (2007) Wäscher, G., Haussner, H., and Schumann, H. (2007).  An improved typology of cutting and packing problems.  European Journal of Operational Research, 183(3):1109–1130.
  * Wei et al. (2020a) Wei, L., Lai, M., Lim, A., and Hu, Q. (2020a).  A branch-and-price algorithm for the two-dimensional vector packing problem.  European Journal of Operational Research, 281(1):25–35.
  * Wei et al. (2020b) Wei, L., Luo, Z., Baldacci, R., and Lim, A. (2020b).  A new branch-and-price-and-cut algorithm for one-dimensional bin-packing problems.  INFORMS Journal on Computing, 32(2):428–443.
  * Woeginger (1997) Woeginger, G. J. (1997).  There is no asymptotic ptas for two-dimensional vector packing.  Information Processing Letters, 64(6):293–297.
  * Xu et al. (2016) Xu, X.-y., Liu, J., Li, H.-y., and Jiang, M. (2016).  Capacity-oriented passenger flow control under uncertain demand: Algorithm development and real-world case study.  Transportation Research Part E: Logistics and Transportation Review, 87:130–148.


[a,b]Benjamin Jäger

# Exploring the large-$N_{c}$ limit with one quark flavour

Michele Della Morte  Sofie Martins  Francesco Sannino  J. Tobias Tsang  Felix
P. G. Ziegler

###### Abstract

We use one-flavour QCD ($N_{c}=3$) as a proxy to understand $\mathcal{N}=1$
SYM. For our simulations, we use tree-level improved Wilson fermions and a
Symanzik improved gauge action. The hadron spectrum is obtained by using LapH
smearing for different masses and simulation volumes. We also report on our
efforts to increase the number of colours in our simulations, where we find
that the simulations show increasing topological freezing for larger $N_{c}$.

20(15.0,1.70) CERN-TH-2022-209

##  1 Introduction

We follow the Corrigan and Ramond large-$N_{c}$ expansion [1] to check the
predictions from super-symmetric Yang-Mills (SYM) theories. The large $N_{c}$
dynamics of this limit was investigated in [2] with crucial differences with
respect to the ordinary large $N_{c}$ dynamics. Additionally at leading order
in this limit, it emerges a connection with the spectrum and dynamics of super
Yang Mills [3, 4, 5] with nontrivial consequences also for the thermodynamics
of the system [6]. Assuming that $N_{c}=3$ is large, we check whether such
predictions are realised by the data. A previous study [7] computed the
mesonic spectrum. We improve upon their findings by considering a finer
lattice spacing, larger volumes, and tree-level improved actions. In the near
future, we plan to extend this to $N_{c}=4,5$ and $6$.

##  2 Lattice Setup

For three numbers of colours ($N_{c}=3$) the two-index anti-symmetric
representation coincides with the conjugate representation, which allows us to
employ "standard" lattice QCD simulations and methods. For our lattice
simulations we use one flavour of tree-level improved Wilson fermions
($c_{sw}=1$) at a fixed gauge coupling of $\beta=4.5$. For the gluonic degrees
of freedom, we utilise the Symanzik-improved gauge action [8]. We use the
openQCD software package [9] and, since we only simulate a single fermion
species, we rely solely on the RHMC algorithm [10]. For the approximation
needed in the RHMC, we use a Zolotarev functional with degree of tenth order.
The integration scheme in the HMC is divided into three levels with two
Omelyan 4th and one 2nd order integrators [11]. The number of integration
steps is tuned to achieve a high acceptance, which we find to be at least
$84\%$ for all our ensembles. To suppress auto-correlation effects we separate
configurations by at least 64 MD units.

$L/a$ | 12 | 12 | 12 | 12 | 16 | 16 | 16 | 16 | 16 | 16
---|---|---|---|---|---|---|---|---|---|---
$\kappa$ | 0.135 | 0.137 | 0.139 | 0.140 | 0.135 | 0.137 | 0.139 | 0.140 | 0.1405 | 0.141
$N_{\mathrm{cnfg}}$ | 7007 | 6210 | 5840 | 5380 | 12059 | 4285 | 9458 | 8465 | 6184 | 7879
$L/a$ | 20 | 20 | 20 | 24 | 24 | 24 | 24 | 32 | 32 |
$\kappa$ | 0.135 | 0.137 | 0.139 | 0.135 | 0.139 | 0.1405 | 0.141 | 0.139 | 0.140 |
$N_{\mathrm{cnfg}}$ | 4046 | 1435 | 2747 | 2844 | 2540 | 2242 | 2000 | 1429 | 428 |
Table 1: Overview of lattice ensembles.

Table 1 shows an overview of the ensembles generated for this study, which
spans a wide range of volumes and hopping parameters $\kappa$. To compare this
setup with other lattice simulations, we apply the Wilson flow [12] to the
Yang-Mills action density and obtain an indicative lattice spacing of
$a=0.06\,$fm. For this we assumed $\sqrt{8\,t_{0}}=0.45\,$fm, which is between
the values for $N_{f}=0$ [12] and $N_{f}=2$ [13].

Simulations of one flavour of Wilson fermions may suffer a sign problem. We
discussed this issue in our previous Lattice contribution [14], where we
presented results for the sign of the Wilson Dirac operator on our ensembles.

##  3 One Flavour for $N_{c}=3$

Figure 1: Example of a simultaneous fit to three vector-vector correlation
functions with different LapH smearings (indicated by the differently coloured
markers). The spectrum is extracted using a three-state fit ansatz. We show
the effective masses for the three correlation functions (data points). We
superimpose the ground state fit result (magenta horizontal band) and the
approach to the plateau (faintly coloured bands).

To check predictions from $\mathcal{N}=1$ SYM we study the spectrum of mesons.
Since only a single fermion species is present, disconnected diagrams appear
in all correlation functions. To include disconnected diagrams we use the
Laplacian Heaviside (LapH) method [15, 16]. For comparison, we also compute
the purely connected correlation function as done in standard lattice
calculations. We use the connected correlation function of the pseudo-scalar
to define the chiral point, i.e. where the quark mass vanishes. As the
connected pseudo-scalar is identified in standard lattice simulations with the
pion, we name the corresponding ground state _fake pion_. This nomenclature is
used to emphasise that this state is unphysical in our setup.

The LapH method has the additional advantage that we can change the level of
smearing by changing the number of eigenvalues used in the approximation of
the Laplacian Heaviside kernel. Due to the multiple smearings, we have access
to more data, with the same spectrum but different overlap coefficients. As
can be seen from the effective mass plot in Figure 1, different smearing
choices result in different approaches to the ground state mass.

We perform a combined correlated fit to multiple correlation functions, which
differ in their LapH smearings. The fit ansatz includes three states, giving
access to the ground state and two excited states. However, we only use the
two lower masses in our subsequent analysis as the second excited state is
prone to fit systematics from the choice of fitrange. The bands in Figure 1
show the effective mass behaviour of the fit result for the correlation
functions that enter the fit. The horizontal magenta band indicates the ground
state mass obtained from the fit.

(a)  (b)

Figure 2: Results for the pseudo-scalar meson as a function for the bare quark
mass (left) and the fake pion mass (right).

Figure 2 shows the extracted pseudo-scalar masses as a function of the bare
quark mass $m_{0}=\tfrac{1}{2\kappa}-4$ on the left and of the connected
pseudo-scalar mass ($m_{\pi}^{\mathrm{fake}}$) on the right. Using the fake
pion mass instead of the bare quark mass, we obtain a significantly smoother
and less volume-sensitive behaviour. Furthermore, this choice facilitates the
definition of the chiral limit, as it eliminates the need to determine the
value $\kappa_{\mathrm{crit}}$ of the hopping parameter at which the quark
mass vanishes.

(a)  (b)

Figure 3: Results of the scalar meson as a function of volume (left) and the
mass of the fake pion (right).

Some of the extracted energies strongly depend on the volume, whereas others
display a clear quark mass dependence. This behaviour is most visible for the
scalar-glue correlation function, which is shown in Figure 3. By carefully
studying the quark-mass and volume dependencies of the meson masses and
overlap coefficients we disentangle mass-independent (glueballs and torelons)
from mesonic states. One conclusion from this study is that the smaller
volumes are affected by sizeable finite-volume effects, especially for the
smallest quark masses used here. For the chiral extrapolations, we explore
linear and quadratic forms in terms of the fake pion mass with and without an
extrapolation to infinite volume. Examples of these chiral extrapolation fits
for the pseudo-scalar, vector, and scalar mesons are shown in Figure 4. We
note that in the chiral limit the mass of the pseudo-scalar meson extrapolates
to a non-vanishing value. This is in agreement with expectations due to the
disconnected contributions.

(a) Pseudo-scalar (b) Vector (c) Scalar

Figure 4: Examples of the chiral extrapolation of the pseudo-scalar (top
left), vector (top right) and scalar-glue (bottom) meson masses to vanishing
quark masses.

##  4 One Flavour for $N_{c}>3$

In order to simulate larger numbers of colours on the lattice, we changed our
simulation software to Grid [17]. As an initial test, we explore quenched
calculations, for which we employ the Symanzik gauge action [8]. The
large-$N_{c}$ limit in pure gauge theory was proposed in the seminal work by
t’Hooft [18]. Accordingly, we rescale $\beta\propto N_{c}^{2}$. Our choices
are listed in table 2.

For all setups, we use a trajectory length of 2.0 to minimise auto-correlation
effects [19].

Name | $N_{c}$ | $L/a$ | $\beta$ | MD steps | Ncnfg | Acceptance
---|---|---|---|---|---|---
pgnc4L12 | 4 | 12 | 8.0 | 30 | 500 | 97.4%
pgnc5L12 | 5 | 12 | 12.5 | 38 | 500 | 97.0%
pgnc6L12 | 6 | 12 | 18.0 | 38 | 500 | 98.4%
pgnc4L16 | 4 | 16 | 8.0 | 30 | 500 | 95.8%
pgnc5L16 | 5 | 16 | 12.5 | 25 | 500 | 86.4%
pgnc6L16 | 6 | 16 | 18.0 | 28 | 500 | 82.8%
pgnc4L24 | 4 | 24 | 8.0 | 30 | 500 | 93.4%
pgnc5L24 | 5 | 24 | 12.5 | 35 | 500 | 91.4%
pgnc6L24 | 6 | 24 | 18.0 | 35 | 500 | 85.2%
Table 2: Overview of the pure gauge simulations for $N_{c}=4,5$ and $6$.

We use a single integrator level with step sizes tuned for acceptance above
80%. We observe that the 2nd-order Omelyan integrator implemented in Grid
scales unsatisfactorily with the numbers of colours and therefore use a force-
gradient integrator [20]. The Hamiltonian violations of this integrator type
are smaller and are less impacted by increases in the lattice volume, see
[21], and number of colours, which lead to a significant performance
improvement in Grid for $N_{c}>3$ in our preliminary examinations. Figure 5
shows the plaquette for three spatial extents ($L/a=12,16,24$) and fixed time-
extent of $T/a=48$ after a conservative estimate of 500 thermalisation
trajectories.

$500$$600$$700$$800$$900$$1{,}000$$0.6070$$0.6075$$0.6080$$0.6085$$0.6090$$500$$600$$700$$800$$900$$1{,}000$$0.5935$$0.5940$$0.5945$$0.5950$$500$$600$$700$$800$$900$$1{,}000$$0.5855$$0.5860$$0.5865$Trajectory
Number$L/a$121624$N_{c}=4$$N_{c}=5$$N_{c}=6$Plaquette Figure 5: Values of the
plaquette as function of MD time for $N_{c}=4,5,6$ (top, middle and bottom
panel, respectively) for three different spatial volumes.

The evolution of the topological charge is shown in Figure 6. As expected, we
observe that fluctuations in the topological charge increase as the volume
increases. However, histograms in the right column of this Figure show that we
have not yet generated sufficiently many configurations to draw clear
conclusions about the exploration of different topological regimes. We observe
indications that topological freezing is exacerbated at larger $N_{c}$, due to
the need to simulate at values of $\beta$ as large as $18$.

$500$$600$$700$$800$$900$$1{,}000$$-10$$0$$10$$500$$600$$700$$800$$900$$1{,}000$$-10$$0$$10$Topological
Charge$500$$600$$700$$800$$900$$1{,}000$$-10$$0$$10$Trajectory
Number$L/a$121624$N_{c}=4$$N_{c}=5$$N_{c}=6$04080Frequency Figure 6:
Topological charge for different number of colours and spatial volumes.

##  5 Conclusions and Outlook

We have presented an update on our current efforts to verify the predictions
from supersymmetric Yang-Mills theories by studying the mesonic spectrum of
QCD with one flavour at $N_{c}=3$. Additionally, we have shown that extending
our simulations to larger $N_{c}$ becomes more costly, in particular as the
topological charge shows significantly reduced fluctuations.

## Acknowledgements

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement № 813942 and 894103. M.D.M. and J.T.T. are partially supported by
DFF Research project 1. Grant n. 8021-00122B. The computations were performed
on UCloud interactive HPC system and ABACUS2.0 supercomputer, which are
managed by the eScience Center at the University of Southern Denmark and LUMI
supercomputer, hosted by the LUMI consortium, as part of the LUMI-G pilot
program.

## References

  * Corrigan and Ramond [1979] E. Corrigan and P. Ramond, Phys. Lett. B 87, 73 (1979).
  * Sannino and Schechter [2007] F. Sannino and J. Schechter, Phys. Rev. D 76, 014014 (2007), arXiv:0704.0602 [hep-ph] .
  * Armoni  _et al._ [2003a] A. Armoni, M. Shifman, and G. Veneziano, Nucl. Phys. B 667, 170 (2003a), arXiv:hep-th/0302163 .
  * Armoni  _et al._ [2003b] A. Armoni, M. Shifman, and G. Veneziano, Phys. Rev. Lett. 91, 191601 (2003b), arXiv:hep-th/0307097 .
  * Sannino and Shifman [2004] F. Sannino and M. Shifman, Phys. Rev. D 69, 125004 (2004), arXiv:hep-th/0309252 .
  * Sannino [2005] F. Sannino, Phys. Rev. D 72, 125006 (2005), arXiv:hep-th/0507251 .
  * Farchioni  _et al._ [2007] F. Farchioni, I. Montvay, G. Munster, E. E. Scholz, T. Sudmann, and J. Wuilloud, Eur. Phys. J. C 52, 305 (2007), arXiv:0706.1131 [hep-lat] .
  * Iwasaki [1983] Y. Iwasaki, (1983), arXiv:1111.7054 [hep-lat] .
  * [9] Lüscher, M. and Schäfer, S., “OpenQCD 1.6,” https://luscher.web.cern.ch/luscher/openQCD/.
  * Clark and Kennedy [2007] M. A. Clark and A. D. Kennedy, Phys. Rev. Lett. 98, 051601 (2007), arXiv:hep-lat/0608015 .
  * I.P. Omelyan and I.M. Mryglod and R. Folk [2003] I.P. Omelyan and I.M. Mryglod and R. Folk, Computer Physics Communications 151, 272 (2003).
  * Lüscher [2010] M. Lüscher, JHEP 08, 071 (2010), [Erratum: JHEP 03, 092 (2014)], arXiv:1006.4518 [hep-lat] .
  * Bruno and Sommer [2014] M. Bruno and R. Sommer (ALPHA), PoS LATTICE2013, 321 (2014), arXiv:1311.5585 [hep-lat] .
  * Ziegler  _et al._ [2022] F. P. G. Ziegler, M. Della Morte, B. Jäger, F. Sannino, and J. T. Tsang, PoS LATTICE2021, 225 (2022), arXiv:2111.12695 [hep-lat] .
  * Morningstar  _et al._ [2011] C. Morningstar, J. Bulava, J. Foley, K. J. Juge, D. Lenkner, M. Peardon, and C. H. Wong, Phys. Rev. D 83, 114505 (2011), arXiv:1104.3870 [hep-lat] .
  * Peardon  _et al._ [2009] M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek, R. G. Edwards, B. Joo, H.-W. Lin, D. G. Richards, and K. J. Juge (Hadron Spectrum), Phys. Rev. D 80, 054506 (2009), arXiv:0905.2160 [hep-lat] .
  * Boyle  _et al._ [2015] P. Boyle, A. Yamaguchi, G. Cossu, and A. Portelli, (2015), arXiv:1512.03487 [hep-lat] .
  * ’t Hooft [1974] G. ’t Hooft, Nucl. Phys. B 72, 461 (1974).
  * Meyer  _et al._ [2007] H. B. Meyer, H. Simma, R. Sommer, M. Della Morte, O. Witzel, and U. Wolff, Comput. Phys. Commun. 176, 91 (2007), arXiv:hep-lat/0606004 .
  * Chin and Kidwell [2000] S. A. Chin and D. W. Kidwell, Phys. Rev. E 62, 8746 (2000), arXiv:physics/0006082 .
  * Kennedy  _et al._ [2013] A. D. Kennedy, P. J. Silva, and M. A. Clark, Phys. Rev. D 87, 034511 (2013), arXiv:1210.6600 [hep-lat] .