A space-indexed formulation of packing boxes into a larger box
AA Space-Indexed Formulation of Packing Boxes into a LargerBox
Sam D. Allen, Edmund K. Burke, Jakub Mareˇcek ∗ The research has been performed, while all authors were at the University of Nottingham, School ofComputer Science, Jubilee Campus, Nottingham, NG8 1BB, UK
Abstract
Integer programming solvers fail to decide whether 12 unit cubes can be packed into a1x1x11 box within an hour using the natural relaxation of Chen/Padberg. We presentan alternative relaxation of the problem of packing boxes into a larger box, which makesit possible to solve much larger instances.
Keywords:
Integer programming, Linear programming, Packing, Load planning,Pigeon hole principle
1. Introduction
Multi-dimensional packing problems have a number of applications in cutting, pack-ing, scheduling, and transportation. Problems in dimension three with rotations aroundcombinations of axes in multiples of 90 degrees are of particular interest in many natu-ral applications. Let us fix the order of six such allowable rotations in dimension threearbitrarily and define:The
Container Loading Problem (CLP) : Given dimensions of a large box (“con-tainer”) x, y, z > n small boxes D ∈ R n × with associated values w ∈ R n , and specification of the allowed rotations r = { , } n × , find the greatest k ∈ R such that there is a packing of small boxes I ⊆ { , , . . . , n } into the container with value k = (cid:80) i ∈ I w i . The packed small boxes I may be rotated in any of the allowed ways, mustnot overlap, and no vertex can be outside of the container.The Van Loading Problem (VLP) : Given dimensions of a large box (“van”) x, y, z >
0, maximum mass p ≥ n small boxes D ∈ R n × with associated values w ∈ R n , mass m ∈ R n , and specification of the allowed rotations r = { , } n × , find the greatest k ∈ R such that there is a packing of small boxes I ⊆ { , , . . . , n } into the container with value k = (cid:80) i ∈ I w i and mass (cid:80) i ∈ I m i ≤ p . The (cid:63) All accompanying materials are available at http://discretisation.sourceforge.net . ∗ The corresponding author is Jakub Mareˇcek.
Email address: [email protected] (Jakub Mareˇcek )
Preprint submitted to Elsevier Received: January 5, 2021 a r X i v : . [ m a t h . O C ] J a n ymbol Meaning n The number of boxes. H A fixed axis, in the set { X, Y, Z } . α An axis of a box, in the set { , , } . L αi The length of axis α of box i . l αi The length of axis α of box i halved. D H The length of axis H of the container. w i The volume of box i in the CLP. Table 1: Notation used. packed small boxes I may be rotated in any of the allowed ways, must not overlap, andno vertex can be outside of the container.Both problems are NP-Hard even to approximate [1]. For example in branch-and-price-and-cut solvers for realistic models of transportation considering both weight andvolume of vehicle load [2, 3], however, one needs to solve the van loading problem to opti-mality as the pricing subproblem. In scheduling, the problem corresponds to schedulingnon-malleable parallel jobs, requiring a certain number of clock cycles, a certain numberof processors, and a certain amount of another discrete resource. Although there area number of excellent heuristic solvers, the progress in exact solvers for the ContainerLoading Problem has been limited, so far. Chen et al. [4] proposed an integer linearprogramming (ILP) formulation, which allowed for up to six boxes to be packed intoa smaller box. Fasano [5] and Padberg [6] improved the formulation and suggested itcould cope with “about twenty boxes”. We confirm these limits using ILOG CPLEX,Gurobi, and SCIP in Section 4. This observation motivates the rest of the paper, wherewe present a new space-indexed formulation providing strong linear programming relax-ations. Throughout, we use the notation as described in Table 1.
2. Related Work
Chen et al. [4] introduced an integer linear programming formulation using the relativeplacement indicator: λ Hij = (cid:40) i precedes box j along axis H δ Hαi = (cid:40) i is rotated so that axis α is parallel to fixed H n (cid:88) i =1 (cid:88) H w i δ H i (1)s.t. (cid:88) H δ H i = (cid:88) H δ H i = (cid:88) α δ Hαi (2) L j ( i ) λ Hj ( i ) i + (cid:88) α l αi δ Hαi ≤ x Hi ≤ (cid:88) α ( D H − l αi ) δ Hαi − L j ( i ) λ Hij ( i ) (3) D H λ Hji + (cid:88) α l αi δ Hαi − (cid:88) α ( D H − l αj ) δ Hαj ≤ x Hi − x Hj (4) ≤ (cid:88) α ( D H − l αi ) δ Hαi − (cid:88) α l αj δ Hαj − D H λ Hij (5) (cid:88) H ( λ Hij + λ Hji ) ≤ (cid:88) H δ H i , (cid:88) H ( λ Hij + λ Hji ) ≤ (cid:88) H δ H j (6) (cid:88) H δ H i + (cid:88) H δ H j ≤ (cid:88) H ( λ Hij + λ Hji ) (7) n (cid:88) i =1 (cid:88) H (cid:32)(cid:89) α L αi (cid:33) δ H i ≤ (cid:89) H D H (8) δ Hαi ∈ { , } , λ Hij ∈ { , } (9) L i ≤ L i ≤ L i , j ( i ) such that L j ( i ) = max { L j } for 1 ≤ i (cid:54) = j ≤ n. (10)Padberg [6] has studied its properties. In particular, he identified the subsets ofconstraints with the integer property. Despite the interesting theoretical properties,modern integer programming solvers fail to solve instances larger than 10–20 boxes usingthis formulation, as evidenced by Table 2. This is far from satisfactory. Arguably, the formulation could be strengthened by addition of further valid con-straints. We have identified three classes of such constraints. Compound constraintsstop combinations of certain boxes being placed next to each other if they would violatethe domain constraint. For example: δ Hki + δ Hmj + λ Hij ≤ l ki + l mj > D H ∀ k, m ∈ { , , } , H ∈ { X, Y, Z } and 1 ≤ i (cid:54) = j ≤ n (11)The example for 2 boxes can be easily extended to 3 or more boxes “in a row.” One canalso attempt to break symmetries in the problem. If boxes i and j (where i < j ) areidentical (i.e. L αi = L αj ∀ α ): λ Hij = 0 and (cid:88) H δ H i ≥ (cid:88) H δ H j (12)Furthermore, if a box has two or more sides of the same length then we can limit therotations. We define function f , which provides a canonical mapping of { X, Y, Z } to3 , , } : (cid:88) H ( f ( H ) · δ Hki ) ≥ (cid:88) H ( f ( H ) · δ Hmi ) if l ki = l mi ∀ ≤ i ≤ n and 1 ≤ k < m ≤
3. A Space-Indexed Formulation
We hence propose a new formulation, where the large box is discretised into units ofspace, whose dimensions are given by the largest common denominator of the respectivedimensions of the small boxes. The small boxes are grouped by their dimensions intotypes t ∈ { , , . . . , n } . A t is the number of boxes of type t available. The formulationuses the space-indexed binary variable: µ tx,y,z = t is placed such that its lower-back-leftvertex is at coordinates x, y, z (cid:88) x,y,z,t µ txyz w t (15)s.t. (cid:88) t µ txyz ≤ ∀ x, y, z (16) µ txyz = 0 ∀ x, y, z, t where x + L t > D X oror y + L t > D Y or z + L t > D Z (17) (cid:88) x,y,z µ txyz ≤ A t ∀ t (18) (cid:88) x (cid:48) ,y (cid:48) ,z (cid:48) ,t (cid:48) ∈ f ( D,n,L,x,y,z,t ) µ t (cid:48) x (cid:48) y (cid:48) z (cid:48) ≤ ∀ x, y, z, t (19) µ txyz ∈ { , } ∀ x, y, z, t (20)where one may use Algorithm 1 in the Appendix or similar to generate the index set f in Constraint 19.The constraints are very natural: No region in space may be occupied by more thanone box type (16), boxes must be fully contained within the container (17), there maynot be more than A t boxes of type t (18), and boxes cannot overlap (19). There is one For one particularly efficient version of Algorithm 1, see class
DiscretisedModelSolver andits method addPositioningConstraints in file
DiscretisedModelSolver.cpp available at http://discretisation.sourceforge.net/spaceindexed.zip (September 30th, 2014). This code is distributedunder GNU Lesser General Public License. (cid:88) x,y,z,t µ txyz m t ≤ p (21) In order to reduce the number of regions of space, and thus the number of variables inthe formulation, a sensible space-discretisation method should be employed. Mareˇcek [7]has shown the problem of finding the best possible discretisation is ∆ p -Complete, where∆ p is the class of problems that can be solved in polynomial time using oracles from NP.There are, however, practical heuristics.The greatest common divisor (GCD) reduction can be applied on a per-axis basis,finding the greatest common divisor between the length of the container for an axis andall the valid lengths of boxes that can be aligned along that axis and scaling by theinverse of the GCD. This is trivial to do and is useful when all lengths are multiples ofeach other, which may be common in certain classes of instances. In other instances,this may not reduce the number of variables at all. In such cases, a non-linear space-discretisation approach can be applied. To do this we use the same values as before ona per-axis basis, i.e. the lengths of any box sides that can be aligned along the axis. Wethen use dynamic programming to generate all valid locations for a box to be placed.For example, given an axis of length 10 and box lengths of 3, 4 and 6, we can place boxesat positions 0, 3, 4, 6, and 7. (Positions 8, 9 and 10 are of course also possible, but nolength is small enough to still lie within the container if placed at these points.) Thishas reduced the number of regions along that axis from 10 to 5. An improvement on thisscale may not be particularly common in practice, though the approach can help wherethe GCD is 1. It is obvious that this approach can be no worse than the GCD methodat discretising the container. This also adds some implicit symmetry breaking into themodel.
4. Computational Experience
We have tested the formulations above on three sets of instances: •
3D Pigeon Hole Problem instances, Pigeon- n , where n + 1 unit cubes are to bepacked into a container of dimensions (1 + (cid:15) ) × (1 + (cid:15) ) × n . • SA and SAX datasets, which are used to test the dependence of solvers’ performanceon parameters of the instances, notably the number of boxes, heterogeneity ofthe boxes, and physical dimensions of the container. There is 1 pseudo-randomlygenerated instance for every combination of container sizes ranging from 5–100 in For n = 10 , , , ,
19, these instances are included in the MIPLIB 2010, cf. http://miplib2010.zib.de/ . • Van-loading instances, generated so as to be similar to those used in pricing sub-problems of branch-and-cut-and-price approaches to vehicle routing with both loadweight and volume considerations [2, 3]. The van-loading and Van-loading-X in-stances use containers of 10x10x30 and 10x10x50 units, respectively, representingthe approximate ratios of a small commercial van and a larger freight truck. 2, 4 or6 out of 6 pseudo-randomly chosen orientations are allowed. The load weights andknapsack values are generated pseudo-randomly and independently of the volumeof the small boxes. 10 instances were generated for each class of problem.All the tests were performed on a 64-bit computer running Linux, which was equippedwith 2 quad-core processors (Intel Xeon E5472) and 16 GB memory. The solvers testedwere IBM ILOG CPLEX 12.2.0, Gurobi Solver 4.0, and SCIP 2.0.1 [8] with CLP asthe linear programming solver. The instances are available at http://discretisation.sourceforge.net .For the 3D Pigeon Hole Problem, results obtained within one hour using the threesolvers and the Chen/Padberg formulation are shown in Table 2, while Table 3 comparesthe results of Gurobi Solver on both formulations. None of the solvers managed toprove optimality of the incumbent solution for twelve unit-cubes within an hour usingthe Chen/Padberg formulation, although the instance of linear programming (after pre-solve) had only 616 rows and 253 columns and 4268 non-zeros. Due to their surprisingdifficulty, Chen/Padberg formulation of Pigeon- n for n ∈ { , , , , } have beencontributed to and accepted into MIPLIB 2010 [9], the benchmark. In contrast, thespace-indexed formulation allowed Gurobi Solver to solve Pigeon-10000000 within anhour, where the instance of linear programming had 10000002 rows, 10000000 columns,and 30000000 non-zeros.For the SA and SAX datasets, Figures 1–2 summarise solutions obtained withinan hour using the Chen/Padberg formulation, the space-indexed formulation, and ametaheuristic approach described by Allen et al. [10, 11]. In the case of the SAX dataset,the volume utilisation is given with respect to the tightest upper bound found by any ofthe solvers. Finally, for the van-loading instances, the results of Gurobi Solver with thespace-indexed model after one hour can be seen in Table 4. As there were 10 instancestested for each class of problem, the mean number of boxes, number of unique box typesand mean gap between solution value and bound are also given.Overall, the space-indexed relaxation provides a particularly strong upper bound.The mean integrality gap, or the ratio of the difference between root linear program-ming relaxation value and optimum to optimum has been 10.49 % and 0.37 % for theChen/Padberg and the space-indexed formulation, respectively, on the SA and SAX in-stances solved to optimality within the time limit of one hour. The mean gap, or theratio of the difference between linear programming relaxation value and value of the bestsolution found to the bound, left after an hour on Van-loading instances was 1.5%. This isencouraging, albeit perhaps not particularly surprising, as similar discretised relaxations6roved to be very strong in scheduling problems corresponding to one-dimensional pack-ing [12, 13, 14] and can be shown to be asymptotically optimal for various geometricproblems both in dimensions two [15] and higher [16] dimensions. Table 2: The performance of various solvers on 3D Pigeon Hole Problem instances encoded in theChen/Padberg formulation. “–” denotes that optimality of the incumbent solution has not been provenwithin an hour.
Time (s)Gurobi 4.0 CPLEX 12.2.0 SCIP 2.0.1 + CLPPigeon-01 < < < < < < < < < < < < < < < <
5. Conclusions and Future Work
The space-indexed linear programming relaxation is strong, but employs a pseudo-polynomial numbers of variables, which presents a challenge to off-the-shelf solvers. Fu-ture work could focus on the development of solvers specific to such relaxations. Onecould either implement a matrix-free interior point method exploiting the structure ofthe relaxation, or employ a column generation schema. The Dantzig-Wolfe decomposi-tion could, perhaps, divide the problem into sub-problems pertaining to each type of box(19) and a master problem enforcing the set-packing constraint (16). Another interestingdirection for future research is automated reformulations. Just as modern integer pro-gramming solvers often fail on small instances of the Chen/Padberg formulation, theyfail on certain formulations of scheduling problems. It has been suggested [7] one couldtry to discretise such relaxations automatically, starting with the identification of whatvariables play the role of relative position in time or space, which could improve theperformance of general-purpose integer linear programming solvers considerably. Thismay help to clarify the importance of discretisation in linear programming relaxations ingeneral.
Acknowledgements.
This work has been first presented at the 7th ESICUP Meeting inBuenos Aires, Argentina, June 6-9, 2010. 7 able 3: The performance of Gurobi 4.0 on 3D Pigeon Hole Problem instances encoded in theChen/Padberg and the space-indexed formulations. “–” denotes no integer solution has been found.
Time (s)Chen/Padberg Space-indexedPigeon-01 < < < < < < < < < < < < < < < < < References [1] M. Chleb´ık, J. Chleb´ıkov´a, Hardness of approximation for orthogonal rectangle packing and coveringproblems, J. Discrete Algorithms 7 (3) (2009) 291–305. doi:10.1016/j.jda.2009.02.002 .[2] M. Gendreau, M. Iori, G. Laporte, S. Martello, A tabu search algorithm for a routing and containerloading problem, Transportation Science 40 (3) (2006) 342–350.[3] M. L¨ubbecke, Personal communication, August 27th, 2010.[4] C. S. Chen, S. M. Lee, Q. S. Shen, An analytical model for the container loading problem, EuropeanJ. Oper. Res. 80 (1) (1995) 68–76. doi:10.1016/0377-2217(94)00002-T .[5] G. Fasano, Cargo analytical integration in space engineering: A three-dimensional packing model,in: T. A. Ciriani, S. Gliozzi, E. L. Johnson, R. Tadei (Eds.), Operational Research in Industry,Purdue University Press, 1999, pp. 232–246.[6] M. Padberg, Packing small boxes into a big box, Math. Methods Oper. Res. 52 (1) (2000) 1–21.[7] J. Mareˇcek, Exploiting structure in integer programs, Ph.D. thesis, The University of Nottingham,Nottingham (2011).[8] T. Achterberg, Constraint integer programming, Ph.D. thesis, Technische Universit¨at Berlin, Berlin(2007).[9] T. Koch, T. Achterberg, E. Andersen, O. Bastert, T. Berthold, R. E. Bixby, E. Danna, G. Gamrath,A. M. Gleixner, S. Heinz, A. Lodi, H. Mittelmann, T. Ralphs, D. Salvagnin, D. E. Steffy, K. Wolter,Miplib 2010, Math. Program. Comput. 3 (2011) 103–163. doi:10.1007/s12532-011-0025-9 .[10] S. D. Allen, E. K. Burke, G. Kendall, A hybrid placement strategy for the three-dimensionalstrip packing problem, European Journal of Operational Research 209 (3) (2011) 219 – 227. doi:10.1016/j.ejor.2010.09.023 .[11] S. D. Allen, Algorithms and data structures for three-dimensional packing, Ph.D. thesis, The Uni-versity of Nottingham, Nottingham (2011).[12] J. P. Sousa, L. A. Wolsey, A time indexed formulation of non-preemptive single machine schedulingproblems, Mathematical Programming 54 (1992) 353–367. able 4: The performance of the space-indexed model on the Van-loading and Van-loading-X instancesafter one hour. Each row represents 10 pseudo-randomly generated instances. Dataset Box Types Boxes Gap(Mean) (Mean)Van-loading-01 3 78 0.10%Van-loading-02 5 92 0.31%Van-loading-03 8 81 0.42%Van-loading-04 10 73 0.32%Van-loading-05 12 73 0.09%Van-loading-06 15 70 0.42%Van-loading-07 20 68 0.50%Van-loading-08 30 71 0.99%Van-loading-09 40 69 1.26%Van-loading-10 50 73 1.10%Van-loading-X-01 3 133 0.16%Van-loading-X-02 5 155 0.16%Van-loading-X-03 8 132 0.56%Van-loading-X-04 10 119 0.87%Van-loading-X-05 12 123 0.29%Van-loading-X-06 15 115 0.36%Van-loading-X-07 20 116 0.30%Van-loading-X-08 30 122 0.29%Van-loading-X-09 40 119 0.18%Van-loading-X-10 50 121 0.28% [13] J. M. van den Akker, LP-based solution methods for single-machine scheduling problems, Ph.D.thesis, Technische Universiteit Eindhoven, Eindhoven (1994).[14] Y. Pan, L. Shi, On the equivalence of the max-min transportation lower bound and the time-indexedlower bound for single-machine scheduling problems, Math. Program. 110 (3, Ser. A) (2007) 543–559.[15] C. H. Papadimitriou, Worst-case and probabilistic analysis of a geometric location problem, SIAMJ. Comput. 10 (3) (1981) 542–557.[16] E. Zemel, Probabilistic analysis of geometric location problems, SIAM J. Algebraic Discrete Meth-ods 6 (2) (1985) 189–200. doi:10.1137/0606017 . ength of container sides N u m be r o f bo x e s Figure 1: The best solutions obtained within an hour per solver per instance from the SA dataset.Each square represents 1 pseudo-randomly generated instance. The number is the volume utilisation (inpercent) with the tight upper bound of 100.
Length of container sides N u m be r o f bo x e s Figure 2: The best solutions obtained within an hour per solver per instance from the SAX dataset.Each square represents 1 pseudo-randomly generated instance. The number displayed is 100(1 − s/b ) forsolution with value s and upper bound b . lgorithm 1 f ( D, n, L, x, y, z, t ) Input:
Discretisation as indices D ⊂ R of µ variables (14), number of boxes n ,dimensions L αi ∈ R ∀ α = x, y, z, i = 1 , . . . , n , indices ( x, y, z ) ∈ D , t ∈ { , , . . . , n } of one scalar within µ (14) Output:
Set of indices of µ variables (14) to include in a set-packing inequality (19) Set S ← { ( x, y, z, t ) } for each other unit ( x (cid:48) , y (cid:48) , z (cid:48) ) ∈ D, ( x, y, z ) (cid:54) = ( x (cid:48) , y (cid:48) , z (cid:48) ) do for each box type t (cid:48) ∈ { , , . . . , n } do if box of type t at ( x, y, z ) overlaps box of type t (cid:48) at ( x (cid:48) , y (cid:48) , z (cid:48) ), i.e., ( x ≤ x (cid:48) + L x (cid:48) t (cid:48) ≤ x + L xt ) ∧ ( y ≤ y (cid:48) + L y (cid:48) t (cid:48) ≤ y + L yt ) ∧ ( z ≤ z (cid:48) + L z (cid:48) t (cid:48) ≤ z + L zt ) then S ← S ∪ { ( x (cid:48) , y (cid:48) , z (cid:48) , t (cid:48) ) } end if end for end for returnreturn