Hybrid approach for solving real-world bin packing problem instances using quantum annealers

Efficient packing of items into bins is a common daily task. Known as Bin Packing Problem, it has been intensively studied in the field of artificial intelligence, thanks to the wide interest from industry and logistics. Since decades, many variants have been proposed, with the three-dimensional Bin Packing Problem as the closest one to real-world use cases. We introduce a hybrid quantum-classical framework for solving real-world three-dimensional Bin Packing Problems (Q4RealBPP), considering different realistic characteristics, such as: (1) package and bin dimensions, (2) overweight restrictions, (3) affinities among item categories and (4) preferences for item ordering. Q4RealBPP permits the solving of real-world oriented instances of 3 dBPP, contemplating restrictions well appreciated by industrial and logistics sectors.


Introduction
The optimization on the packaging of products into a finite number of containers is a crucial daily task in the field of production and distribution.Depending on the characteristics of both packages and containers, multiple packaging problems can be formulated, generally known as Bin Packing Problems (BPP) 1 .Within this category, the one-dimensional BPP (1dBPP) is considered as the simplest one 2 , whose goal is to pack all items into as few containers as possible.Many variants with a variable number of constraints have been proposed to deal with real situations in logistics and industry 3 .The three-dimensional BPP (3dBPP) 4 , in which each package has three dimensions: height, width and depth, is the best-known and the most challenging variant.Highlighted in several studies [5][6][7] , 3dBPP has a practical interest in many industrial settings.In recent years, it has been formulated to possess diverse and practical applications such as pallet loading 8 , road transportation 9 , air cargo 10 , etc. Due to its complexity, 3dBPP is also recurrently employed as a benchmark for testing newly developed methods and mechanisms 11,12 .
On another front, quantum computing is still at its early stage but has gathered a lot of attention from the scientific community as it offers researchers and practitioners a revolutionary paradigm for tackling different kinds of practical optimization problems [13][14][15][16] .In particular, quantum annealers have been recently applied to a wide variety of optimization problems inspired by the fields from industry 17 , logistics 18 and economics 19 .However, the research on BPP carried out in the quantum community is still scarce, even though BPP has been widely studied classically as an optimization problem.
The pioneering work on BPP in the field of quantum computing presents a hybrid quantum-classical method for solving the 1dBPP 20 , whose solver is composed of two modules: i) a quantum subroutine with which to search a set of feasible configurations to fill one single bin and ii) a classical computational heuristic which builds complete solutions employing the subsets given by the quantum subroutine.To deepen the performance of the quantum subroutine developed, further tests were conducted against a random sampling and a random walk-based heuristic 21 .Besides those two papers, an additional study formulates an atomic energy industry related problem as a 1dBPP, solving it using the D-Wave quantum annealer 22 .Another works show quantum-inspired evolutionary computation techniques as an alternative to tackle BPP related problems [23][24][25] .Quantum-inspired techniques are a specific class of evolutionary algorithms which make use of quantum physics to define their operations and are designed to be executed on a classical computer 26 .Thus, they can not be executed on any quantum machine.
In contrast to 1dBPP, tackling 3dBPP in the quantum domain is much more challenging due to two related grounds: i) its complexity, which increases as constraints from the real-world are taken into account and ii) the incipient state of development of the current commercial quantum computers with capacities still limited by decoherence and errors, which could be an obstacle to solve highly-constrained problems.In this paper, we present a hybrid quantum-classical computing framework for the real-world oriented 3dBPP, which is coined as Quantum for Real Bin Packing Problem (Q4RealBPP).The proposed framework resorts to the the Leap Constrained Quadratic Model (CQM) Hybrid Solver (LeapCQMHybrid 27 ) of D-Wave.At the same time, Q4RealBPP is built on an existing code 28 .This reference code is an excellent starting point, which paved the way for these two main contributions developed in this this work.
• Efficiency of the code: in the reference code, the problem is formulated such that lots of variables (thus qubits) are needed.
This issue meets the problem of feasibility in the context of quantum hardware in the noisy intermediate-scale quantum (NISQ) era 29 .Therefore, the optimization of the code is crucial for dealing with complex problems.Q4RealBPP provides a step forward on this aspect by exploring how the constraints coined as Intrinsic restrictions can be polished.These restrictions are the ones needed for the BPP definition (e.g.do not place cases outside a bin), and they were previously introduced in the reference code.Thus, the novel work performed on this specific facet has involved, among others, the redefinition of some mathematical formulas that define these intrinsic restrictions, and the elimination of a redundant optimization objective.Thanks to this procedure, problems can be formulated by using fewer variables, which directly impacts the capacity of the framework to address larger instances.
• Applicability of the tool: the reference code is oriented toward solving the most basic variant of the 3dBPP by only considering the dimension of both packages and bins.This setting falls far from real client demands, where other features such as weights, load balancing or incompatibilities,among others, are likely to take part.We have elaborated on this direction by implementing a set of constraints named Real-World BPP Restrictions.All these constraints represent an added value for Q4RealBPP, having required a significant extension of the mathematical formulation of the problem.We deepen in this aspect in following sections.
Taking these factors into account, Q4RealBPP is oriented to industrial and logistics related fields, contemplating problems such as the organization of port containers, the introduction of packages in delivery vans and trucks or the placement of foodstuffs on distribution pallets, among others.With the aid of a hybrid quantum-classical method, Q4RealBPP represents a solid step forward to solve 3dBPP with the clear purpose of facing real-world focused problems well appreciated by final users and companies.Features contemplated by Q4RealBPP involve i) dimensions of packages and bins, ii) maximum weight allowed per bin, iii) positive and negative affinities among item categories and iv) preferences for package ordering (in terms of load bearing and load balancing).To demonstrate its application, we have conducted an experimentation composed of 12 different instances serving as illustrative examples.Additionally, Q4RealBPP allows users to easily build flexible and well-defined instances to adapt a plethora of real-world situations to be solved in the quantum computer.
The rest of the article is organized as follows: in Mathematical formulation the 3dBPP formulation and its corresponding notation are presented.Moreover, a detailed study of the computational resources needed for arbitrary instances is carried out.In Experimental results, the applicability of this tool is tested using a set of realistic instances as input.Finally, the conclusions led by the presented results and our future plans are given in Conclusions and future work.

Mathematical formulation
In this section, we describe in detail the mathematical formulation of the 3dBPP variant tackled in this research.First, input parameters and variables that compose the problem are shown in Table 1.

Objectives
The 3dBPP can be solved as an optimization problem where a suitable cost function to minimize must be defined.In our case, this cost function is represented as the sum of three objectives.The strength given to each objective, i.e. the relevance accounted for each one, is up to the user preferences just by multiplying each objective with a suitable weight.Thus, the problem can be stated as min ∑ 3 i=1 ω i o i with ω i the weights of each objective o i .In our study we will not consider this bias, i.e. ω i = 1 ∀i.The first and main objective minimizes the total amount of bins used to locate the packages.This can be achieved by minimizing Additionally, for ensuring that items are packed from the floor to the top of the bin, avoiding solutions with floating packages, a second objective is defined by minimizing the average height of the items for all bins Besides these two objectives reformulated from the reference code 28 , we further add a third optional objective o 3 to take into account the load balancing feature.This concern is particularly important when air cargo planes and sailings are the chosen Parameters I, J, K i , Q sets of items, bins, orientations of i-th item (see Figure 1) and relative positions between items (see Figure 2).m, n number of items and bins.l i , w i , h i , µ i length, width, height and weight of item i ∈ I. L,W, H, M length, width, height and maximum capacity (optional) of bins.A pos , A neg sets of positive and negative affinities (incompatibilities) between items of type α and β (optional).
η maximum mass ratio between two items where one is placed at the top of the other with η > 1 (optional).( L, W ) target center of mass of the resultant packings (optional).
Variables v j binary variable that represents if bin j ∈ J is used.u i, j binary variable that represents if item i is added to bin j. r i,k binary variable that represents if the orientation k ∈ K i is applied to the item i.They are used to compute the effective length, width and height 6)-( 8)).x i , y i , z i continuous variables that return the location of the back lower left corner of item i along x, y and z axes.xi , ỹi continuous variables that account the relative distance between item i and ( L, W ) along x and y axes.Both are used if ( L, W ) is defined previously.b i,k,q binary variable that returns the relative position q ∈ Q between items i, k ∈ I. See Figure 2.
Table 1.Parameters and variables used in our formulation.
conveyance 30,31 , for example.In those situations, packages should be uniformly distributed around a given xy-coordinate inside the bin.We can tackle this by computing the so-called taxicab or Manhattan distance between items and the desired center of mass for each bin.As a result, the gaps between items are also reduced.Concerning this, the third objective to be minimized is where 0 ≤ x i < nL (bins stacked horizontally) and 0 ≤ y i < W ∀i ∈ I.This objective term minimizes for each item the distance between the center of mass projection in the xy-plane and the ( L, W ) coordinate of each bin.
The objectives above defined are subject to certain restrictions, which are essential to derive realistic solutions.The whole pool of constraints is separated into two categories: the ones intrinsic to the BPP definition (Intrinsic restrictions), and the ones relevant from a real-world perspective (Real-world BPP restrictions).

Intrinsic restrictions
Item orientations: the fact that inside a bin each item must have only one orientation can be implemented by using Orientations give rise to the effective length, width, and height of the items along x, y and z axes and because of (5), only one term r i,k is nonzero in each equation.
It should be deemed that there could be items with geometrical symmetries, as with cubic ones where rotations do not apply.Redundant and non-redundant orientations are considered in the reference code 28 .In our formulation, we previously check if these symmetries exist to define K i for each item.Thanks to this, ( 6)-( 8) are simplified filtering out redundant orientations and leading to a formulation which uses less variables (thus qubits) to represent the same problem, where x y z x y z (e) x y z Table 2. Subsets of non-redundant orientations for item i for a satisfied condition following ( 6)-( 8).
variables r i,k are needed.For i ∈ I c with I c := {i ∈ I | l i = w i = h i } (cubic items), we can set r i,1 = 1 and 0 otherwise, thus satisfying (5) in advance.In Table 2, we can see the non-redundant orientation sets for an item i depending on its dimensions.This simple mechanism reduces the complexity of the problem, being favourable for the quantum hardware to implement.Non-overlapping restrictions: since we are considering rigid packages, i.e. they can not overlap, a set of restrictions need to be defined to overcome these configurations.For this purpose, at least one of these situations must occur (see Figure 2) Item i is at the right of item k (q = 4): Item i is in front of item k (q = 5): Item i is above item k (q = 6): As discussed with the orientation variable r i,k in (5), the relative position between items i and k must be unique, so Item and container allocation restrictions: the following set of restrictions guarantees an appropriate behaviour during item and bin assignment.In order to avoid packing duplicates of the same item, each item must go to exactly one bin, where The following formula verifies if items are being packed inside bins that are already in use x y z (a) x y z x y z (e) x y z Representation of b i,k,q activated for all relative positions q ∈ Q between items i and k.See ( 9)-( 14).Both are in contact but it is not mandatory.(a

4/11
so it activates v j if needed during packaging.Bins can be activated sequentially to avoid duplicated solutions ensuring that Bin boundary constraints: in order to contemplate bin boundaries, the following set of restrictions must be met where (19) guarantees that items i placed inside the bin j are not outside of the last bin (n-th bin) along the x axis, (20) ensures that item i is located inside of its corresponding bin j along the x axis (activated if n > 1), (21) confirms that item i placed inside the bin j is not outside along the y axis, while (22) ensures that item i allocated inside the bin j is not outside along the z axis.

Real-world BPP restrictions
In this subsection we introduce those restrictions related with the operative perspective of the problem, i.e. the ones that consider real-world industrial situations.All of the following constraints are optional in our formulation.
Overweight restriction: the weight of each package and the maximum capacity of containers are common contextual data to avoid exceeding the maximum weight capacity of bins, so avoid overloaded containers.We can introduce this restriction as This restriction is activated if the maximum capacity M is given.
Affinities among package categories: there are commonly preferences for separating some packages into different bins (negative affinities or incompatibilities) or, on the contrary, gathering them into the same container (positive affinities).Let us consider I α := {i ∈ I | id of i is equal to α}, i.e.I α ⊂ I is a subset of all items labelled with id equal to α.Given a set of p negative affinities A neg := {(α 1 , β 1 ), . . ., (α p , β p )}, then the restriction will be To activate this restriction, a set of incompatibilities must be given.Moreover, we can satisfy in advance ν := 6n ∑ (α,β )∈A neg |I α ||I β | non-overlapping constraints (see ( 9)-( 14)), leading to a simpler formulation.Conversely, given a set of positive affinities A − A pos as stated with A neg , then the restriction will be posed such that This restriction is activated if a set of positive affinities is given.If A pos and A neg are given, then both restrictions can be introduced using just one formula adding ( 24) and (25).
Preferences in relative positioning: relative positioning of items demands that some of them must be placed in a specific position with respect other existing items.This preference allows introducing the ordering of a set of packages according to their positions with respect to the axes.Thus, this preference assists in ordering for many real cases such as: parcel delivery (an item i that has to be delivered before item k will be preferably placed closer to the trunk door) or load bearing (no heavy package should rest over flimsy packages), among others.
Regarding this preference, we can define two different perspectives to treat relative positioning: • Positioning to avoid (P − q ): list of items (i, k) should not be in the relative position q ∈ Q specified.So, b i,k,q = 0 is expected, favouring configurations where the solver selects q ′ ∈ Q with q ′ ̸ = q for the relative positioning of items (i, k).
• Positioning to favour (P + q ): list of items (i, k) should be in a certain relative position q.Activated this preference, b i,k,q = 1 ought to hold and consequently, b i,k,q ′ = 0 ∀q ′ ̸ = q.
Formally, these preferences are written as

5/11
These preferences could be also treated as compulsory pre-selections.In such case, the number of variables needed would be reduced, so would the search space.If we let p − = ∑ q∈Q |P − q | and p + = ∑ q∈Q |P + q | with P − q ∩ P + q ′ = ∅, based on (15), the amount of variables reduced would be given by p − + 6p + .Moreover, n(p − + 5p + ) non-overlapping constraints (see ( 9)-( 14)) are satisfied directly and can be ignored, thus simplifying the problem.In this paper, for the sake of clarity, these preferences have been applied for load bearing purposes as hard constraints (HC), as explained in the upcoming Experimental results.
Load balancing: to activate this restriction, a target center of mass must be given.Global positions with respect to the bin as a whole (as described in objective o 3 in (3)), are fixed using the following constraints This feature is represented in Figure 3 for ( L, W ) = (L/2,W /2), whose red line shows the available xi and ỹi values (see ( 4)).

Complexity of the problem
Regarding the complexity of the 3dBPP proposed in this research, the total amount of variables needed to tackle an arbitrary instance is given in Table 3, where our formulation scales as O[m 2 + nm] in terms of variables.Additionally, the total amount of constraints required is provided in Table 4, whose quantity grows quadratically as Table 3. Amount of mandatory (related to Intrinsic restrictions) and optional (related to Real-world BPP restrictions) used variables depending on the number of items and bins.See Table 1.

Experimental results
In this section, we conduct an experimentation to demonstrate the applicability of Q4RealBPP, where the problem has been modelled as a CQM and then solved by the LeapCQMHybrid provided by D-Wave 27 .Initially, it should be made clear that CQM refers to the mathematical model that uses quadratic objective functions and quadratic constraints on binary and integer variables.As this concept was introduced by D-Wave Systems, this company developed the hybrid solver LeapCQMHybrid, which also supports the definition of equality and inequality requirements.This feature brings an advantage compared to Quadratic Unconstrained Binary Optimization (QUBO), which is the native formulation for the QPUs.The CQM model and Table 5. Brief description of the real-world oriented restrictions activated and the amount of variables and constraints used for each instance (see Tables 3 and 4).In brackets, number of items considered.Here, OW: overweight, PA: positive affinities, INC: incompatibilities, LB: load bearing, CM: center of mass.
the LeapCQMHybrid introduce some interesting shortcuts for a more friendly use, allowing the user to provide a problem in a more intuitive and descriptive way, avoiding the translation into a mathematical formulation in the shape of a QUBO matrix.More specifically, the LeapCQMHybrid workflow is as follows: first, a classical front end receives the CQM formulation of the problem as input.Then, it runs a number of parallel computation threads, where each of them are composed of two parts: a heuristic module (HM) and a quantum module (QM).On the one hand, HM tries to solve the problem by using a state of the art heuristic solver.On the other hand, QM aids this resolution by formulating different quantum queries aiming to guide the HM toward promising regions of the search space, or finding improvements to already existing solutions.This QM performs its operations by acceding the latest available quantum computer of D-Wave.At the time of writing this paper, the most updated architecture is the Advantage_system, composed of 5616 qubits organized in a Pegasus topology.For more information, we refer interested readers to the related report provided by D-Wave 32 .
Having said that, for demonstrating the applicability of Q4RealBPP, we have built an ad-hoc benchmark composed of 12 different instances of the 3dBPP.In order to analyze the impact of every restriction deeply, each instance is devoted to evaluate a specific feature of the problem.Also, we have generated two specific instances that bring together all the restrictions of our modelled 3dBPP.We describe in Table 5 the main characteristics of the 12 used instances.
The whole dataset has been generated employing an own-developed Python script (coined as Q4RealBPP-DataGen).In order this paper to be as self-contained as possible, we briefly explain how Q4RealBPP-DataGen works.This script performs two steps to generate Q4RealBPP-compliant instances: firstly Q4RealBPP-DataGen randomly generates a defined number of items, following the package distribution and dimensions established in Ref. 8 .Then, the attributes regarding overweight restriction, affinities among item categories, load bearing and load balancing are completed by means of Q4RealBPP-DataGen.These constraints are randomly configured by the generator, except for the last two (load bearing and load balancing).These last features, as well as the bin dimensions, have been empirically set, in the search of a realistic scenario.It should be clarified that Q4RealBPP is a flexible framework, letting users to build their own setting not only configuring the instance of the problem but also by activating or deactivating the real-world oriented restrictions.For more information, the complete benchmark of instances and Q4RealBPP-DataGen are openly available in 33 .Furthermore, a deep explanation about how the dataset has been generated as well as the format of each instances is provided in 34 .
In our specific use case, the preferences in relative positioning (see Real-world BPP restrictions) are tested as HC for load bearing.Accounting fragility issues, one could apply the rule of choosing pairs of packages to decide on what height to place each of them based on a mass ratio η (assuming that weight is related to fragility).Thus, defining , this instantiation avoid configurations where items whose mass are more than η times the mass of other ones are placed above of them.
Figure 4 represents the results provided by Q4RealBPP for each of the instances described in Table 5. Regarding the running time, we have empirically determined the time it takes for the LeapCQMHybrid to resolve these instances to be 30 s, presenting a maximum Quantum Processing Unit (QPU) access time of 0.032 s per execution.Lastly, for interested readers, all obtained results are freely available 33 .
Besides that, we have conducted additional experimentation with the goal of further understanding the performance of 7/11 0 0  Q4RealBPP.For this purpose, we have solved all the instances described in Table 5 under different time limits.Thus, each instance has been run 10 times with 5, 10, 30 and 60 seconds as time limits.In the top of Figure 5, the obtained results are shown using the mean µ and and standard deviation σ of the energy value {e i } 10 i=1 , being e i the i-th energy returned by the solver upon the assigned time limit for each instance.Moreover, the mean QPU access time in nanoseconds is also attached.As a complementary analysis, in the botton part of Figure 5, since the returned energies vary considerably depending on the nature 8/11 for illustrating more clearly the stability of the solution.This additional analysis is depicted in the bottom part of Figure 5. Three main conclusions can be drawn from Figure 5: 1) generally, the longer the time limit given, the lower the energy returned (thus the better solution quality); 2) the deviation around the mean values of the objective function remains stable against the different time limits given as well as it does not show a remarkable dependency with the nature of the instance studied; and 3) although different time limits have been given to the solver, the QPU access time did not vary significantly for each one, indicating that the optimization process on the LeapCQMHybrid workflow mainly relies on the heuristic module of the solver.Despite the clear behaviour exhibited by our results, further work should be done for having a deeper understanding of how the presented framework performs.
To prove the applicability of Q4RealBPP, we have tested it over 12 instances of different nature, with the main intention of showcasing the capacity of the method to encompass real-world constraints.As depicted in Figure 4, Q4RealBPP has successfully tackled all the generated instances, contemplating different real-world situations.Particularly noteworthy are the last two instances, 3dBPP_11 and 3dBPP_12 (Figures 4(l) and 4(m), respectively), where all the constraints are activated.
The future work comprises 3 main interests: i) to develop a more advanced version of the Q4RealBBP to generalise the framework to other BPP variants, and considering further features such as multi-class categorisation of items; ii) to exploit this framework in conjunction to complementary Artificial Intelligence algorithms to enhance its potential real applications; and iii) to further analyse the results by comparing performance with classical solvers.

Figure 4 .
Figure 4. (a) Colour palette used in the solved instances.(b)-(m) Brief description of instances given in Table 5 and solutions provided by Q4RealBPP.Red lines show the bin boundaries.The activated restrictions work as expected.

Figure 5 .
Figure 5. Top: mean + std for the 10 runs belonging to each instance.From left to right on each instance, the time limits given to the solver are 5 s, 10 s, 30 s and 60 s.Numbers placed above and below the shaded areas show the mean QPU access time for every time limit within each instance (in units of nanoseconds).Bottom: using(28), deviation around the mean energy of the results depicted above.For both plots, shaded area shows the region where the minimum and maximum values fell during the study.