Three-Dimensional Container Loading: A Simulated Annealing Approach

16 Pages·2017·0.83 MB·English

UUnniivveerrssiittyy ooff WWiinnddssoorr
SScchhoollaarrsshhiipp aatt UUWWiinnddssoorr
Industrial and Manufacturing Systems Department of Industrial and Manufacturing
Engineering Publications Systems Engineering
2017
TThhrreeee--DDiimmeennssiioonnaall CCoonnttaaiinneerr LLooaaddiinngg:: AA SSiimmuullaatteedd AAnnnneeaalliinngg
AApppprrooaacchh
Hanan Mostaghimi Ghomi
Universty of Windsor
Bryan Gary St. Amour
University of Windsor
Walid Abdul-Kader
University of Windsor
Follow this and additional works at: https://scholar.uwindsor.ca/industrialengpub
Part of the Industrial Engineering Commons
RReeccoommmmeennddeedd CCiittaattiioonn
Mostaghimi Ghomi, Hanan; St. Amour, Bryan Gary; and Abdul-Kader, Walid. (2017). Three-Dimensional
Container Loading: A Simulated Annealing Approach. International Journal of Applied Engineering
Research, 12 (7), 1290-1304.
https://scholar.uwindsor.ca/industrialengpub/13
This Article is brought to you for free and open access by the Department of Industrial and Manufacturing Systems
Engineering at Scholarship at UWindsor. It has been accepted for inclusion in Industrial and Manufacturing
Systems Engineering Publications by an authorized administrator of Scholarship at UWindsor. For more
information, please contact [email protected].
International Journal of Applied Engineering Research ISSN 0973-4562 Volume 12, Number 7 (2017) pp. 1290-1304
© Research India Publications. http://www.ripublication.com
Three-Dimensional Container Loading: A Simulated Annealing
Approach

1Hanan Mostaghimi, 2Bryan St. Amour and 1#Walid Abdul-Kader
1Faculty of Engineering, University of Windsor, Windsor (Ontario), Canada
2Department of Computer Science, University of Windsor, Windsor (Ontario), Canada
#ORCID: 0000-0002-2926-5059

Abstract to a compact, which results in a potential loss of the market
when compared to other competitor's vehicles (Hifi, 2002;
High utilization of cargo volume is an essential factor in the
Bortfeldt et al., 2003; and Wang et al., 2008).
success of modern enterprises in the market. Although
mathematical models have been presented for container loading The basic form of the container loading problem is packing the
problems in the literature, there is still a lack of studies that best subset of rectangular boxes (called cargo) into a large
consider practical constraints. In this paper, a Mixed Integer object (called container) to maximize the total value of the
Linear Programming is developed for the problem of packing a loaded boxes. The boxes should not be overlapped and should
subset of rectangular boxes inside a container such that the total lie entirely in the container. According to the topology
value of the packed boxes is maximized while some realistic introduced by Wascher et al. (2007), the containers can be
constraints, such as vertical stability, are considered. The either homogeneous or heterogeneous. If the boxes placed in
packing is orthogonal, and the boxes can be freely rotated into the given container are identical, it is called homogeneous;
any of the six orientations. Moreover, a sequence triple-based however, if various types of boxes are loaded, the container is
solution methodology is proposed, simulated annealing is used considered as heterogeneous. Besides the non-overlapping
as modeling technique, and the situation where some boxes are constraints, some other practical constraints should be
preplaced in the container is investigated. These preplaced considered in the real-world container loading problem, such as
boxes represent potential obstacles. Numerical experiments are cargo vertical stability, preplaced boxes, and box rotation
conducted for containers with and without obstacles. The (Junqueira et al., 2012; Bortfeldt et al., 2012). However, not
results show that the simulated annealing approach is many papers have considered these practical constraints in their
successful and can handle large number of packing instances. proposed models.
Keywords: knapsack, packing sequence, rotation, obstacles, The problem addressed in this research and as per the topology
simulated annealing presented by Dyckhoff (1990), belongs to 3/B/O/F (3: three-
dimensional, B/O: one object/bin and items selection, F: few

items of different types), while Wascher et al. (2007) classify it
INTRODUCTION as the three-dimensional single orthogonal knapsack problem.
The given problem considers the packing of rectangular items
Logistics has recently played an important role in the success
into a container to maximize the total value of the packed items
of modern enterprises. Packing boxes inside a container is an
by minimizing the amount of lost space. The value of boxes is
essential material handling activity in manufacturing and
assumed to be equal to their volume. The rotation of the boxes
transportation industries. It is also a key function for operating
is considered as well. The multidimensional knapsack Problem
supply chain efficiently. The efficient use of transportation
(MKP) is a NP-hard optimization problem that can be shown
devices, like containers and palettes, leads to significant cost
by reduction from the one-dimensional packing problem
saving. Moreover, high utilization of transportation devices
(Egeblad and Pisinger, 2009). Although technological
reduces the traffic of goods and protects natural resources.
knowledge has been enhanced, solving real knapsack problems
Therefore, optimal loading of a container decreases the
is still a challenge. Due to NP-hardness of the packing problem,
shipping cost and increases the stability of the load. Container
only heuristic methods, and a few exact algorithms have been
loading problem has practical values, and it can be applied to
presented.
various fields. Loading cars, trucks, trains, or ships can be also
considered as a container loading or a three-dimensional In this paper, a mixed integer linear model and a simulated
packing problem. Furthermore, cargo volume is an important annealing algorithm are developed to address a more
factor used by motor vehicle companies to market their comprehensive knapsack problem where practical
products as sub-compact, compact, midsize, or full size. If considerations, such as vertical stability and preplaced
inaccurate, the cargo volume found may downgrade the vehicle (obstacles) constraints, are tackled. These practical constraints
from its real size. For example, a mid-size may be downgraded and box rotation contribute significantly to a study of a more
1290
International Journal of Applied Engineering Research ISSN 0973-4562 Volume 12, Number 7 (2017) pp. 1290-1304
© Research India Publications. http://www.ripublication.com
realistic 3D knapsack problem as proposed in this research proposed for packing a batch of objects. Goncalves et al. (2012)
work. The aim of this proposed research is to provide managers propose a multi-population biased random-key genetic
and decision-makers with an adequate modeling tool that helps algorithm for the single container loading problem. Maximal-
make shipping goods more efficient and eco-friendly as fewer space representation is used to manage the container free space.
trips can be made if higher container utilization is achieved. The authors consider stability and orientation constraints;
Also, auto-makers can market the class size of their vehicles however, they do not develop a mathematical model for the
more accurately. given problem. Peng et al. (2009) present a hybrid simulated
annealing algorithm for three-dimensional container loading

problem. Firstly, a heuristic algorithm is used for encoding
LITERATURE REVIEW feasible packing solutions, and then the simulated annealing
algorithm is applied to search in the encoding neighborhood.
The focus of most of the container loading or three-dimensional
Hongtao et al. (2012) address a three-dimensional single
cutting and packing problems is on the rectangular bins. Fekete
container loading problem by using a multi-stage search based
and Schepers (2004) propose a new method for obtaining lower
simulated annealing algorithm. Wei et al. (2012) use a
bound classes for higher-dimensional packing problem. The
reference length approach to address the three-dimensional
major objective of this paper is to define good criteria for
strip packing problem. In another paper, Wei et al. (2015)
removing a candidate set of boxes. Dual feasible function is a
address the problem of multiple container loading cost
way to build conservative scales. The computational results are
minimization problem by using a new approach that combines
mainly limited to the two-dimensional packing problem. Hifi
column generation technique with a goal driven search.
(2004) proposes a dynamic algorithm and an exact depth-first
search to solve the three-dimensional cutting problem. Other research works related to design automation focus on
Orientation and guillotine constraints are considered. Optimal three-dimensional placement of circuit elements by exploring
solutions are obtained for a significant number of instances, but the layout of the integrated circuits (Obenaus and Szymansky,
not all of them. Althaus et al. (2007) consider the trunk packing 2003), and Cheng et al. (2005) address floor planning for 3-D
problem where the box dimensions are as per the SAE J1100 VLSI Design. While this technique is less known to container
standard. They propose two discretized methods. First, the loading practices, it carries some similarities and it is more
space to be packed is discretized. Then, an approximation efficient in terms of search time than other methods such as
approach is considered using linear inequalities. The space partitioning placement.
discretization causes insufficient representation of the boxes.
Models that provide information on optimal objective function
Additionally, the runtime of the enumerative algorithms is
value and bounds help to assess the solution quality of heuristic
exponential.
algorithms. Although modeling three-dimensional knapsack
Although considerable advancement has been made in the problems considers practical constraints, it is still at its
development of exact algorithms, heuristic algorithms still play beginnings. Junqueira et al. (2012) present mixed integer linear
an important role in solving three-dimensional knapsack programming models for the container loading problem.
problems. Only heuristic methods can provide reasonable Vertical and horizontal stability of the cargo, as well as cargo
solutions within an acceptable running time for instances of load bearing strength, are considered in the proposed model.
real-world sized problems. Pisinger (2002) develops a wall- However, the models are only able to handle moderate sized
building based heuristic. Both homogenous and heterogeneous problems. Table 1 compares some relevant papers and models,
instances are considered. Moreover, several ranking rules are and shows their similarities and differences.
studied to select the best layers’ depths. Bortfeldt et al. (2003)
Per the literature, not all papers consider box rotation since it
propose a parallel tabu search approach for a single container
increases the search space significantly. Bin stability
loading problem and give little consideration to heterogeneous
constraints have likewise been just considered in a few papers.
instances. Wang et al. (2008) also present a heuristic method
To the best of our knowledge, preplaced boxes (obstacles) have
for a heterogeneous container loading problem and developed
not been studied in three-dimensional knapsack problems,
a dynamic space decomposition approach based on the tertiary
although it is so essential for such problems since it is often
tree structure. Egeblad and Pisinger (2009) propose a simulated
required to place certain boxes in certain positions. These
annealing based methodology for the two and three-
constraints can also be used in the case of having a non-
dimensional knapsack problems, and a three-dimensional
rectangular container. Therefore, it is important to study more
knapsack model is presented. The authors present an iterative
practical constraints in the knapsack problem. This proposed
heuristic approach for the knapsack problem that is based on
research work aims to contribute to the literature so that a 3D
the sequence triple representation. Also, Yamazaki et al. (2000)
knapsack problem can be tackled where box rotation is
apply a variety of packing sequences including sequence-triple
considered to help finding more practical packing
in their 3-D packing solution approach. To control the heuristic
configurations. Furthermore, preplaced boxes (bin with some
method, simulated annealing is used. However, rotation of
obstacles) and vertical stability constraints are considered and
boxes and preplaced boxes (obstacles) are not considered in the
addressed. This is useful, especially when considering trunk
three-dimensional model and experiments. Wu et al. (2010)
loading for auto size classification as indicated earlier in the
consider the three-dimensional bin packing problem with
Introduction section.
variable bin heights. A mixed integer programming model is
proposed, and they also present the case when more than one
type of bin is used. A genetic algorithm-based heuristic is

1291
International Journal of Applied Engineering Research ISSN 0973-4562 Volume 12, Number 7 (2017) pp. 1290-1304
© Research India Publications. http://www.ripublication.com
Table 1: Models Parameters Comparison that the boxes can be freely rotated in six different orientations.
However, it is possible to relax this constraint and fix a box in
Paper Box Rotation Vertical Stability Preplaced Boxes athematical Model euristic Approach PArLsoboabllrelvge etmeo s altwsoaaop sbhypsesseou pctroasmiesfecfc,ie l cielodfa eistn cfh.ptd t e-hTo obratsh rohtiibe tetetio t noovhbxtmnaeaoe l tst-uit dboaoeoiern mmreho. i fetnoT hncedfehsoa iee(cnobL ahsnibcniB sdhob Be xfooblr)exoefo s codx iaor sn .l rm lenea Ieqseunbd ru so ataspx ndlhbre odeotestoipu t sliltaiaoudotncnsp dbe ,bp v deeosot hropltpeublmeaa modcecbx e kebiebd.ney so d I ixtt an hoeri ineeassr
M H integers; thus, the placement is to be done in integer steps. Let
C be a rectangular bin with width W, height H and depth D.
Junqueira The origin of the Cartesian coordinate system is located at the
et al. √ √ LBB corner of the container, w, h, and d are respectively, the
i i i
(2012) width, height and depth of box type i. and (xi, yi, zi) represent
the coordinates of the LBB corner of the box.
Goncalves
et al. √ √ √ √ A mixed integer programming formulation is presented for the
(2012) given problem. Some real-world knapsack problem constraints
are considered in the model, and to the best of our knowledge,
Wu et al. they have not been addressed previously. These constraints are
√ √ √ √
(2010) vertical stability and preplaced boxes (obstacles). Since the
three-dimensional knapsack problem is NP-hard, it is difficult
Egeblad & to solve. Additionally, the flexibility of the orientation of boxes
Pisinger Just 2D √ √ √ increases the search space significantly so that the difficulty of
(2009) finding the optimal solution is increased as well. Some exact
algorithms and heuristic methods are proposed in the published
Wang et al.
√ √ √ literature. As exact algorithms require more time to find a
(2008)
solution, heuristic approaches are more popular and can be
effective alternatives to finding an optimal or near optimal
Bortfeldt et
√ √ √ √ solution. The proposed three-dimensional solution
al. (2003)
methodology is based on sequence triple representation, which
Pisinger is defined below under Placement Algorithm. Simulated
√ √ √
(2002) annealing is used as the meta-heuristic method. As the number
of box types (or box dimensions and variety) is finite, the use
Proposed of simulated annealing is favoured by its efficiency in
√ √ √ √ √ √
Research neighborhood search.

PROBLEM DEFINITION MATHEMATICAL FORMULATION
In this study, the three-dimensional knapsack problem is A mixed-integer programming model of the 3D-knapsack
considered where there is one bin with fixed size and a set of problem is formulated in this section. The mathematical model
boxes, and each box has an associated size. The aim is to find is based on work done by Egeblad and Pisinger (2009) and Wu
an efficient solution methodology to pack rectangular boxes in et al. (2010). Some modifications are made to their models and
a single bin so that the total value of the packed boxes is include new constraints addressing vertical stability and
maximized, or equivalently the empty spaces left are obstacles, which were not considered in any previously
minimized. The boxes are assumed to be strongly published works. Constraints (1) – (4) are adapted from
heterogeneous, which means that there is relatively many Egeblad and Pisinger (2009); however, the authors did not
different types of boxes and a small number of boxes for each consider box orientation in their model. The binary position
box type (Wascher et al., 2007). Moreover, the packing is variables, which show the orientation of the boxes, are
considered feasible if each box lies entirely in the bin and the integrated in constraints (5) – (17). This makes the model more
packed boxes do not overlap. The edges of all boxes must be comprehensive. They are described below in this section.
parallel to the edges of the bin (orthogonal packing). The boxes
The following are the main assumptions considered for the mix
are assumed to be of rectangular shape; however, the bin can
integer linear model:
be considered either of rectangular or nonrectangular shape. In
the case of having preplaced boxes (obstacles), the bin is 1. The boxes are strongly heterogeneous,
assumed to be non-rectangular. 2. The boxes must be located orthogonally,
3. The boxes can freely rotate,
Some practical considerations that play an important role in
4. The box and bin dimensions are assumed to be non-
modeling more realistic knapsack problems, such as box
negative integer,
rotation and bin stability, are presented. The algorithm assumes
5. The value of a box is equal to its volume, and
1292
International Journal of Applied Engineering Research ISSN 0973-4562 Volume 12, Number 7 (2017) pp. 1290-1304
© Research India Publications. http://www.ripublication.com
6. The X, Y, and Z axes of the bin are shown in the following  Parameters:
figure.
(w, h, d): width, height, and depth of box i
i i i
Y
(W, H, D): width, height, and depth of the bin
(r, s, k): Left-Bottom Behind (LBB) coordinates of the
X preplaced boxes (obstacles)
Z
(a, b, c, d): binary orientation parameters of the preplaced boxes

P: value of box i
Figure 1: The X, Y, and Z axes of the bin i
M: large number
Notation
Pb: set of preplaced boxes
The model notation for the parameters and variables are as
follows:
 Variables:
(x,y,z): Left-Bottom Behind (LBB) coordinates of box i
i i i
Xw: 1 if width of box i is parallel to the container’s X
i
0 otherwise
Zw: 1 if width of box i is parallel to the container’s Z
i
0 otherwise
Yh: 1 if height of box i is parallel to the container’s Y
i
0 otherwise
Zd: 1 if depth of box i is parallel to the container’s Z
i
0 Otherwise
rij, lij: 1 if box i is to the right (rij) or to left (lij) of box j
0 otherwise
o , u : 1 if box i is over (o ) or under (u ) box j
ij ij ij ij
0 otherwise
b , f : 1 if box i is behind (b ) or in front (f ) of box j
ij ij ij ij
0 otherwise
s: 1 if box i is packed
i
0 otherwise
𝑦𝑎: 1 if 𝑥 ≥𝑥 (x coordinate of the LBB corner of box j is greater than or equal to x coordinate of the LBB corner of
𝑖𝑗 𝑗 𝑖
box i)
0 otherwise
𝑥𝑎: 1 if 𝑥 <𝑥′ (x coordinate of the LBB corner of box j is less than x coordinate of the Right-Bottom-Behind (RBB)
𝑖𝑗 𝑗 𝑖
corner of box i)
0 otherwise
1293
International Journal of Applied Engineering Research ISSN 0973-4562 Volume 12, Number 7 (2017) pp. 1290-1304
© Research India Publications. http://www.ripublication.com
𝑦𝑏: 1 if 𝑧 ≥𝑧 (z coordinate of the LBB corner of box j is greater than or equal to z coordinate of the LBB corner of
𝑖𝑗 𝑗 𝑖
box i)
0 otherwise
𝑥𝑏: 1 if 𝑧 <𝑧′ (z coordinate of the Left-Bottom-Front (LBB) corner of box j is less than z coordinate of the Left-
𝑖𝑗 𝑗 𝑖
Bottom-Front (LBF) corner of box i)
0 otherwise
𝑦𝑐: 1 if 𝑥′ >𝑥 (x coordinate of Right-Bottom-Behind (RBB) corner of box j is greater than x coordinate of the LBB
𝑖𝑗 𝑗 𝑖
corner of box i)
0 otherwise
𝑥𝑐: 1 if 𝑥′ ≤𝑥′ (x coordinate of RBB corner of box j is less than or equal to x coordinate of the RBB corner of box i)
𝑖𝑗 𝑗 𝑖
0 otherwise
𝑦𝑑: 1 if 𝑧′ >𝑧 (z coordinate of the LBB corner of box j is greater than z coordinate of the LBF corner of box i)
𝑖𝑗 𝑗 𝑖
0 otherwise
𝑥𝑑: 1 if 𝑧′ ≤𝑧′ (z coordinate of the LBF corner of box j is less than or equal to z coordinate of the LBF corner of box i)
𝑖𝑗 𝑗 𝑖
0 otherwise
𝑧𝑎: 1 if 𝑥 ≤𝑥 <𝑥′ (both 𝑥𝑎 and 𝑦𝑎are equal to one)
𝑖𝑗 𝑖 𝑗 𝑖 𝑖𝑗 𝑖𝑗
0 otherwise
𝑧𝑏: 1 if 𝑧 ≤𝑧 <𝑧′ (both 𝑥𝑏 and 𝑦𝑏are equal to one)
𝑖𝑗 𝑖 𝑗 𝑖 𝑖𝑗 𝑖𝑗
0 otherwise
𝑧𝑐: 1 if 𝑥 <𝑥′ ≤𝑥′ (both 𝑥𝑐 and 𝑦𝑐are equal to one)
𝑖𝑗 𝑖 𝑗 𝑖 𝑖𝑗 𝑖𝑗
0 otherwise
𝑧𝑑: 1 if 𝑧 <𝑧′ ≤𝑧′ (both 𝑥𝑑 and 𝑦𝑑are equal to one)
𝑖𝑗 𝑖 𝑗 𝑖 𝑖𝑗 𝑖𝑗
0 otherwise
𝐶𝑠 : 1 if 𝑥 ≤𝑥 <𝑥′ and 𝑧 ≤𝑧 <𝑧′ (both 𝑧𝑎 and 𝑧𝑏are equal to one)
1𝑖𝑗 𝑖 𝑗 𝑖 𝑖 𝑗 𝑖 𝑖𝑗 𝑖𝑗
0 otherwise
𝐶𝑠 : 1 if 𝑥 ≤𝑥 <𝑥′ and 𝑧 <𝑧′ ≤𝑧′ (both 𝑧𝑎 and 𝑧𝑑are equal to one)
2𝑖𝑗 𝑖 𝑗 𝑖 𝑖 𝑗 𝑖 𝑖𝑗 𝑖𝑗
0 otherwise
𝐶𝑠 : 1 if 𝑥 <𝑥′ ≤𝑥′ and 𝑧 ≤𝑧 <𝑧′ (both 𝑧𝑐 and 𝑧𝑏are equal to one)
3𝑖𝑗 𝑖 𝑗 𝑖 𝑖 𝑗 𝑖 𝑖𝑗 𝑖𝑗
0 otherwise
𝐶𝑠 : 1 if 𝑥 <𝑥′ ≤𝑥′ and 𝑧 <𝑧′ ≤𝑧′ (both 𝑧𝑐 and 𝑧𝑑are equal to one)
4𝑖𝑗 𝑖 𝑗 𝑖 𝑖 𝑗 𝑖 𝑖𝑗 𝑖𝑗
0 otherwise

1294
International Journal of Applied Engineering Research ISSN 0973-4562 Volume 12, Number 7 (2017) pp. 1290-1304
© Research India Publications. http://www.ripublication.com
𝑥′ =𝑥 +𝑤𝑋𝑤 +ℎ (𝑍𝑤 −𝑌ℎ +𝑍𝑑 )+𝑑 (1−𝑋𝑤 −𝑍𝑤 +𝑌ℎ −𝑍𝑑 )
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖
𝑧′ =𝑧 +𝑑 𝑍𝑑 +ℎ (1−𝑍𝑤 −𝑍𝑑 )+𝑤𝑍𝑤
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖
Objective Function:
The objective function is to maximize the value of the packed boxes:𝑀𝑎𝑥∑𝑛 𝑃𝑠
𝑖=1 𝑖 𝑖
Subject to:
𝑟 +𝑙 +𝑏 +𝑓 +𝑜 +𝑢 =𝑠 +𝑠 −1 ∀𝑖,𝑗,𝑖 ≠𝑗 (1)
𝑖𝑗 𝑗𝑖 𝑖𝑗 𝑗𝑖 𝑖𝑗 𝑗𝑖 𝑖 𝑗
𝑥 +𝑤𝑋𝑤 +ℎ (𝑍𝑤 −𝑌ℎ +𝑍𝑑 )+𝑑 (1−𝑋𝑤 −𝑍𝑤 +𝑌ℎ −𝑍𝑑 )≤𝑥 +𝑀(1−𝑙 ) ∀𝑖,𝑗,𝑖 ≠𝑗 (2𝑎)
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑗 𝑖𝑗
𝑥 +𝑤𝑋𝑤 +ℎ (𝑍𝑤 −𝑌ℎ +𝑍𝑑 )+𝑑 (1−𝑋𝑤 −𝑍𝑤 +𝑌ℎ −𝑍𝑑 )≤𝑥 +𝑀(1−𝑟 ) ∀𝑖,𝑗,𝑖 ≠𝑗 (2𝑏)
𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑖 𝑖𝑗
𝑧 +𝑑 𝑍𝑑 +ℎ (1−𝑍𝑤 −𝑍𝑑 )+𝑤𝑍𝑤 ≤𝑧 +𝑀(1−𝑏 ) ∀𝑖,𝑗,𝑖 ≠𝑗 (2𝑐)
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑗 𝑖𝑗
𝑧 +𝑑 𝑍𝑑 +ℎ (1−𝑍𝑤 −𝑍𝑑 )+𝑤𝑍𝑤 ≤𝑧 +𝑀(1−𝑓 ) ∀𝑖,𝑗,𝑖 ≠𝑗 (2𝑑)
𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑖 𝑖𝑗
𝑦 +ℎ 𝑌ℎ +𝑤(1−𝑋𝑤 −𝑍𝑤)+𝑑 (𝑋𝑤 +𝑍𝑤 −𝑌ℎ )≤𝑦 +𝑀(1−𝑢 ) ∀𝑖,𝑗,𝑖 ≠𝑗 (2𝑒)
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑗 𝑖𝑗
𝑦 +ℎ 𝑌ℎ +𝑤(1−𝑋𝑤 −𝑍𝑤 )+𝑑 (𝑋𝑤 +𝑍𝑤 −𝑌ℎ )≤𝑦 +𝑀(1−𝑜 ) ∀𝑖,𝑗,𝑖 ≠𝑗 (2𝑓)
𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑖 𝑖𝑗
𝑥 +𝑤𝑋𝑤 +ℎ (𝑍𝑤 −𝑌ℎ +𝑍𝑑 )+𝑑 (1−𝑋𝑤 −𝑍𝑤 +𝑌ℎ −𝑍𝑑 )≤𝑊 (3𝑎)
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖
𝑦 +ℎ 𝑌ℎ +𝑤(1−𝑋𝑤 −𝑍𝑤)+𝑑 (𝑋𝑤 +𝑍𝑤 −𝑌ℎ )≤𝐻 (3𝑏)
𝑗 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖
𝑧 +𝑑 𝑍𝑑 +ℎ (1−𝑍𝑤 −𝑍𝑑 )+𝑤𝑍𝑤 ≤𝐷 (3𝑐)
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖
𝑋𝑤 +𝑍𝑤 ≤1 (4𝑎)
𝑖 𝑖
𝑍𝑤 +𝑍𝑑 ≤1 (4𝑏)
𝑖 𝑖
0≤𝑍𝑤 −𝑌ℎ +𝑍𝑑 ≤1 (4𝑐)
𝑖 𝑖 𝑖
0≤1−𝑋𝑤 −𝑍𝑤 +𝑌ℎ −𝑍𝑑 ≤1 (4𝑑)
𝑖 𝑖 𝑖 𝑖
0≤𝑋𝑤 +𝑍𝑤 −𝑌ℎ ≤1 (4𝑒)
𝑖 𝑖 𝑖
(𝑥,𝑦,𝑧)=(𝑟,𝑠,𝑘) (5)
𝑖 𝑖 𝑖
(𝑋𝑤,𝑍𝑤,𝑍𝑑 ,𝑌ℎ )=(𝑎,𝑏,𝑐,𝑑) ∀𝑖 ∈𝑃𝑏 (6)
𝑖 𝑖 𝑖 𝑖
𝑥 −𝑥 ≤𝑀.𝑦𝑎𝑥 −𝑥 ≥𝑀(𝑦𝑎 −1) ∀𝑖,𝑗, 𝑖 ≠𝑗 (7𝑎)
𝑗 𝑖 𝑖𝑗 𝑗 𝑖 𝑖𝑗
𝑥′−𝑥 ≤𝑀.𝑥𝑎𝑥′−𝑥 ≥𝑀(𝑥𝑎 −1)+0.5 ∀𝑖,𝑗, 𝑖 ≠𝑗 (7𝑏)
𝑖 𝑗 𝑖𝑗 𝑖 𝑗 𝑖𝑗
𝑦𝑎 +𝑥𝑎 −1 𝑦𝑎 +𝑥𝑎
𝑖𝑗 𝑖𝑗 ≤𝑧𝑎 ≤ 𝑖𝑗 𝑖𝑗 ∀𝑖,𝑗,𝑖 ≠𝑗 (7𝑐)
2 𝑖𝑗 2
𝑧 −𝑧 ≤𝑀.𝑦𝑏𝑧 −𝑧 ≥𝑀(𝑦𝑏 −1) ∀𝑖,𝑗,𝑖 ≠𝑗 (8𝑎)
𝑗 𝑖 𝑖𝑗 𝑗 𝑖 𝑖𝑗
𝑧′−𝑧 ≤𝑀.𝑥𝑏𝑧′−𝑧 ≥𝑀(𝑥𝑏 −1)+0.5 ∀𝑖,𝑗,𝑖 ≠𝑗 (8𝑏)
𝑖 𝑗 𝑖𝑗 𝑖 𝑗 𝑖𝑗
𝑦𝑏 +𝑥𝑏 −1 𝑦𝑏 +𝑥𝑏
𝑖𝑗 𝑖𝑗 ≤𝑧𝑏 ≤ 𝑖𝑗 𝑖𝑗 ∀𝑖,𝑗,𝑖 ≠𝑗 (8𝑐)
2 𝑖𝑗 2
𝑥′−𝑥 ≤𝑀.𝑦𝑐𝑥′−𝑥 ≥𝑀(𝑦𝑐 −1)+0.5 ∀𝑖,𝑗,𝑖 ≠𝑗 (9𝑎)
𝑗 𝑖 𝑖𝑗 𝑗 𝑖 𝑖𝑗
𝑥′−𝑥′ ≤𝑀.𝑥𝑎𝑥′−𝑥′ ≥𝑀(𝑥𝑐 −1) ∀𝑖,𝑗,𝑖 ≠𝑗 (9𝑏)
𝑖 𝑗 𝑖𝑗 𝑖 𝑗 𝑖𝑗
𝑦𝑐 +𝑥𝑐 −1 𝑦𝑐 +𝑥𝑐
𝑖𝑗 𝑖𝑗 ≤𝑧𝑐 ≤ 𝑖𝑗 𝑖𝑗 ∀𝑖,𝑗,𝑖 ≠𝑗 (9𝑐)
2 𝑖𝑗 2
𝑧′−𝑧 ≤𝑀.𝑦𝑑𝑧′−𝑧 ≥𝑀(𝑦𝑑 −1)+0.5 ∀𝑖,𝑗,𝑖 ≠𝑗 (10𝑎)
𝑗 𝑖 𝑖𝑗 𝑗 𝑖 𝑖𝑗
𝑧′−𝑧′ ≤𝑀.𝑥𝑑𝑧′−𝑧′ ≥𝑀(𝑥𝑑 −1) ∀𝑖,𝑗,𝑖 ≠𝑗 (10𝑏)
𝑖 𝑗 𝑖𝑗 𝑖 𝑗 𝑖𝑗
𝑦𝑑 +𝑥𝑑 −1 𝑦𝑑 +𝑥𝑑
𝑖𝑗 𝑖𝑗 ≤𝑧𝑑 ≤ 𝑖𝑗 𝑖𝑗 ∀𝑖,𝑗,𝑖 ≠𝑗 (10𝑐)
2 𝑖𝑗 2
1295
International Journal of Applied Engineering Research ISSN 0973-4562 Volume 12, Number 7 (2017) pp. 1290-1304
© Research India Publications. http://www.ripublication.com
𝑧𝑎 +𝑧𝑏 −1 𝑧𝑎 +𝑧𝑏
𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗
≤𝐶𝑠 ≤ ∀𝑖,𝑗,𝑖 ≠𝑗 (11)
2 1𝑖𝑗 2
𝑧𝑎 +𝑧𝑑 −1 𝑧𝑎 +𝑧𝑑
𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗
≤𝐶𝑠 ≤ ∀𝑖,𝑗,𝑖 ≠𝑗 (12)
2 2𝑖𝑗 2
𝑧𝑐 +𝑧𝑏 −1 𝑧𝑐 +𝑧𝑏
𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗
≤𝐶𝑠 ≤ ∀𝑖,𝑗,𝑖 ≠𝑗 (13)
2 3𝑖𝑗 2
𝑧𝑐 +𝑧𝑐 −1 𝑧𝑐 +𝑧𝑑
𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗
≤𝐶𝑠 ≤ ∀𝑖,𝑗,𝑖 ≠𝑗 (14)
2 4𝑖𝑗 2
𝐶𝑠 +𝐶𝑠 +𝐶𝑠 +𝐶𝑠 =𝑢 +𝑜 ∀𝑖,𝑗,𝑖 ≠𝑗 (15)
1𝑖𝑗 2𝑖𝑗 3𝑖𝑗 4𝑖𝑗 𝑖𝑗 𝑖𝑗
𝑥′ =𝑥 +𝑤𝑋𝑤 +ℎ (𝑍𝑤 −𝑌ℎ +𝑍𝑑 )+𝑑 (1−𝑋𝑤 −𝑍𝑤 +𝑌ℎ −𝑍𝑑 ) (16)
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖
𝑧′ =𝑧 +𝑑 𝑍𝑑 +ℎ (1−𝑍𝑤 −𝑍𝑑 )+𝑤𝑍𝑤 (17)
𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖 𝑖
𝑟 ,𝑙 ,𝑜 ,𝑢 ,𝑏 ,𝑓 ∈{0,1} (18)
𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗
𝑋𝑤,𝑍𝑤,𝑍𝑑 ,𝑌ℎ ∈{0,1} (19)
𝑖 𝑖 𝑖 𝑖
𝑥𝑎,𝑥𝑏,𝑥𝑐,𝑥𝑑,𝑦𝑎,𝑦𝑏,𝑦𝑐,𝑦𝑑,𝑧𝑎,𝑧𝑏,𝑧𝑐,𝑧𝑑 ∈{0,1} (20)
𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗 𝑖𝑗
𝑠,𝐶𝑠 ,𝐶𝑠 ,𝐶𝑠 ,𝐶𝑠 ∈{0,1} (21)
𝑖 1𝑖𝑗 2𝑖𝑗 3𝑖𝑗 4𝑖𝑗
(𝑥 ,𝑦,𝑧)≥0 (22)
𝑖 𝑖 𝑖
Constraint (1) ensures that if box i and box j are packed, they first was run for an instance with 5 boxes; it reached the optimal
must be placed left (l), right (r), under (u), over (o), behind (b) solution in 53 seconds. Then the instance with 6 boxes is
or in-front (f) of each other. Constraints (2) guarantee that any considered, and the solution time is equal to 6 minutes and 14
two i and j boxes do not overlap, while considering the box seconds. However, the solution time for the instance with 7
rotation. The binary position variables (Xw, Zw, Yh, Zd) are boxes increased significantly to 4 hours and 4 minutes. The
i i i i
used to allow box rotations. Constraint set (3) ensures that all optimal results for instance with 8 boxes was obtained after 21
boxes are placed within the bin’s dimensions. Constraint set (4) hours and 39 minutes. The model was not able to reach optimal
is used to ensure that the binary variables that show the position solution for instance with 9 boxes even after 3 days; thus, the
of the boxes are controlled to represent practical positions. algorithm was terminated before reaching the solution.
Constraint (5) and (6) are used to fix the coordinates and According to the results, optimal solutions were only possible
orientations of the preplaced boxes, where Pb is a set of in a reasonable time for small size instances (up to 8 boxes).
preplaced boxes. Constraints (7) – (10) ensure vertical stability.
While this research problem is well-known as NP hard, the
These constraints compare the four corners of each newly
mathematical model as presented above has helped define and
packed box with the points that cover the top of other packed
outline this problem in more details. However, as addressed in
boxes. If one of the corners has the same x and z coordinates as
the following section, the heuristic algorithm is required to get
one of the mapped points, it means that the new box is located
solutions for larger instances in a reasonable time.
under or above that box. Constraint set (7) is used to define the
binary variable 𝑧𝑎, and it includes three parts. Constraint (7a)
𝑖𝑗
ensures that if 𝑥 ≥𝑥, then 𝑦𝑎 is equal to one; otherwise it is
𝑗 𝑖 𝑖𝑗 SOLUTION METHODOLOGY – HEURISTIC
equal to zero. Constraint (7b) makes sure that if 𝑥𝑗 <𝑥𝑖, then APPROACH
𝑥𝑎 is one; otherwise it is equal to zero. Constraint (7c)
𝑖𝑗 Three sequences are considered for the boxes to pack. These
guarantees that when 𝑦𝑎 and 𝑥𝑎 are both equal to one, then 𝑧𝑎
𝑖𝑗 𝑖𝑗 𝑖𝑗 sequences show the relative box locations. They are known as
is equal to one. Similarly, constraint sets (8), (9), and (10) are
sequence triple. Sequence triple is one of the most successful
used to define the binary variables 𝑧𝑏, 𝑧𝑐, and 𝑧𝑑. Constraints
𝑖𝑗 𝑖𝑗 𝑖𝑗 representations in the literature and defines the packing order.
(11) – (14) show whether the x and z coordinates of the corner This section will first discuss placement approach in the
of the new box are equal to x and z coordinates of the mapped subsection. Then, simulated annealing is discussed, which
points on the top of the packed boxes. Constraint (15) ensures controls the local neighbourhood search. Finally, orthogonal
that if these coordinates are the same, then the new box should rotation, preplaced boxes (obstacles), four-corner packing, and
be located on top of or under the packed box. Constraints (16) box insertion order are explained later in this section.
and (17) define 𝑥′ and 𝑧′. Finally, constraints (18) – (21)
𝑖 𝑖
represent the binary variables while constraint (22) represents
the integer variables. Placement Algorithm
The given mathematical model was coded in GAMS/Cplex, Three sequences A, B, and C represent the fully robot packable
and the computational tests were run on an Intel® Core™ i5 packing, where A, B, and C are permutations of the numbers 1
CPU @ 2.67GHz processor with 4.0 GB RAM. The model at ... n, and n is the total number of boxes to be placed in the bin.
1296
International Journal of Applied Engineering Research ISSN 0973-4562 Volume 12, Number 7 (2017) pp. 1290-1304
© Research India Publications. http://www.ripublication.com
These sequences denote the relative placement of each of the 2. Consider (x + w, y, z)
i i i i
two i and j boxes with respect to each other. Each sequence is
∀𝑗 ∈𝑃 : compute x’ and z’
defined as follows: 𝑦 j j
Where 𝑥 +1≤𝑥′ ≤𝑥 +𝑤 and 𝑧 ≤𝑧′ ≤𝑧 +𝑑 −1
· A-chain: If box i appears before box j in the A-chain, then box 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗
i is located to the left of, on top of, or in front of box j.
If (x+ w = x’ and z= z’) then
i i j i j
· B-chain: if box i appears before box j in the B-chain, then box return y+ h
j j
i is located behind, to the left of, or below, box j.
else go to 3
· C-chain: If box i appears before box j in the C-chain, then box
i is located to the right, under, or in front of box j. 3. Consider (xi, yi, zi + di)
Based on the three given sequences, it is possible to determine ∀𝑗 ∈𝑃𝑦: compute x’j and z’j
whether box i is located on the left side of, below, or behind
Where 𝑥 ≤𝑥′ ≤𝑥 +𝑤 −1 and 𝑧 +1≤𝑧′ ≤𝑧 +𝑑
box j. It is observed that box i always appears before box j in 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗
B-chain for all three given placements. Thus, the order of If (x = x’ and z+ d = z’)
i j i i j
placement of the boxes in the bin can be based on the B-chain.
return y+ h
The first box is placed at the origin, and the succeeding boxes j j
are placed according to their relative position to boxes that are else go to 4
already packed. The coordinates of each new box are calculated
4. Consider (x + w, y, z + d)
based on the following formulae: i i i i i
∀𝑗 ∈𝑃 : compute x’ and z’
𝑥 =𝑚𝑎𝑥(0,𝑚𝑎𝑥 (𝑥 +𝑤)) 𝑦 j j
𝑖 𝑗∈𝑃𝑥 𝑗 𝑗
Where 𝑥 +1≤𝑥′ ≤𝑥 +𝑤 and 𝑧 +1≤𝑧′ ≤𝑧 +𝑤
𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗
𝑦 =𝑚𝑎𝑥(0,𝑚𝑎𝑥 (𝑦 +ℎ ))
𝑖 𝑗∈𝑃𝑦 𝑗 𝑗 If (xi+ wi= x’j and zi+dj = z’j) then
𝑧 =𝑚𝑎𝑥(0,𝑚𝑎𝑥 (𝑧 +𝑑 )) return yj+ hj
𝑖 𝑗∈𝑃𝑧 𝑗 𝑗
else return 0
where P , P , and P are the subsets of packed boxes located on
x y z
the left, below, and behind the new box. return 𝑦 =𝑚𝑎𝑥(0,(𝑦 +ℎ ))
𝑖 𝑗 𝑗

The algorithm pushes each packed box downward where
To consider vertical stability and reduce the gap between the
possible such that its bottom can be supported by the floor of
boxes, some modifications have been applied to calculate the
the bin or by the top of other packed boxes.
y-coordinate of each packed box. These modifications are
explained below.
Simulated annealing
 Vertical Stability Although it is relatively not difficult to develop a simulated
annealing heuristic approach, choosing a good neighborhood
As it is assumed that the (x, y, z) coordinates of the boxes and
and cooling procedure, which itself depends on several
their dimensions are integer, it is possible to map a set of points
different parameters, is usually necessary for the algorithm to
that a certain box covers. Let (x, y, z) be the Left-Bottom
i i i work efficiently. The cooling procedure is different for various
Behind (LBB) coordinates of each box to be packed. The
types of problems, and even between instances of the same
algorithm considers four corners of the given box. If the x and
problem. Therefore, it is difficult to find out a good cooling
z coordinates of one of these corners are equal to the
procedure. In the proposed simulated annealing algorithm, the
coordinates of one of the points at the top of any packed box,
temperature is reduced when a new solution is accepted,
then the height of that box is returned. The y-coordinate of the
according to the following function:
new box would then be equal to the maximum of those values.
The proposed novel approach, which is illustrated below, 𝑡 →𝑡⁄(1+𝛽𝑡)
ensures that the vertical stability is satisfied.
where 𝛽 is the cooling parameter. Besides the cooling down
1. Consider (x, y, z) procedure, the process is allowed to heat up again whenever it
i i i
appears to be getting trapped. The heating up function is:
∀𝑗 ∈𝑃 : compute x’ and z’
𝑦 j j
𝑡 →𝑡⁄(1−𝛼𝑡)
Where 𝑥 ≤𝑥′ ≤𝑥 +𝑤 −1 and 𝑧 ≤𝑧′ ≤𝑧 +𝑑 −1
𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗 𝑗
where 𝛼 is the heating parameter. The temperature is reduced
If (xi = x’j and zi= z’j) then when the solution is accepted and increased when the solution
is rejected. Parameter 𝛼 must be smaller than 𝛽 as the number
return y+ h
j j
of acceptances is small when compared to the number of
else go to 2 rejections (Dowsland, 1993).
1297
International Journal of Applied Engineering Research ISSN 0973-4562 Volume 12, Number 7 (2017) pp. 1290-1304
© Research India Publications. http://www.ripublication.com

The neighbourhood of each solution is defined as one of these Four-corner packing
five permutations: either exchange two boxes from one of the
Four packing schemes, one for each corner, are created. First,
sequences; exchange two boxes in sequences A and B;
the coordinates of the boxes are calculated in relation to the
exchange two boxes in sequences A and C; exchange two boxes
current origin. Then, their real (x, y, z) coordinates are
in sequences C and B; or exchange two boxes in all sequences.
calculated in relation to the real origin of the bin which is its
The solutions are compared based on the bin utilization. The
LBB corner. The processing technique is as follows:
formula used for calculating the utilization percentage is as
follows: W := bin width
𝑡𝑜𝑡𝑎𝑙 𝑣𝑜𝑙𝑢𝑚𝑒 𝑜𝑓𝑝𝑎𝑐𝑘𝑒𝑑 𝑏𝑜𝑥𝑒𝑠
H := bin height
𝑢𝑡𝑖𝑙𝑖𝑧𝑎𝑡𝑖𝑜𝑛 𝑝𝑒𝑟𝑐𝑒𝑛𝑡𝑎𝑔𝑒= ×100
𝑣𝑜𝑙𝑢𝑚𝑒 𝑜𝑓 𝑏𝑖𝑛
D := bin depth

w := box width
It should be noted that the volume (or internal volume) of the
bin is a reference to an ideal solution that may not be attained h := box height
due to the discrete nature of the boxes being loaded. However,
this volume is used in the relation above to check how far the d := box depth
solution obtained by the simulated annealing algorithm is from
if (loading from front) then
the ideal solution.
// No change needed: this is the default loading method.

Orthogonal Rotation Return <x,y,z>
An extension we add to the typical A, B, C-chain representation else if (loading from rear) then
of the packing is to allow for the boxes to be rotated
return <W – x – w, y, D – z – d>
orthogonally with respect to the bin. Width, height, and depth
of all boxes are respectively parallel to x, y, and z axis, and w,
i else if (loading from left side) then
h, and d represents the width, height, and depth of box i,
i i
respectively. It is possible to obtain better packing if the boxes return <W – z – w, y, x>
are rotated in different directions. Boxes are allowed to be
else if (loading from right side) then
rotated in one of the following orientations:
WHD: Standard orientation. return <z, y, D – x – w>
WDH: Swap the height and the depth. end
HWD: Swap the width and the height.
HDW: Swap the width and the height, and then swap the height Order of box insertion
with the depth.
The order the boxes are inserted into the container is based on
DHW: Swap the depth with the width. sequence B and can either be created randomly, or based on the
volume of the boxes, which means that the boxes with larger
DWH: Swap the depth with the width, and then swap the depth
volumes are packed first (in a first fit decreasing order).
with the height.

The given rotation is applied to the simulated annealing by
adding an additional transformation to the neighbourhood NUMERICAL EXPERIMENTS
generating routine. The orientation of the boxes is generated
The proposed methodology is implemented in C++, and the
randomly at first. Thus, an additional vector R, which shows
code is tested using two different sets of boxes. The dimensions
the orientation of the boxes, is stored. The sequence triple is
of the first set of boxes are taken from the SAE J1100 Standard,
stored as well.
which includes 7 different types of boxes (after excluding the
golf bag). The dimensions of these boxes are illustrated in
Table 2a. Twelve instances are created by using the first set of
Obstacles
boxes. These instances contain 36 and 70 boxes. The maximum
If Pb is the set of rectangular obstacles (or preplaced boxes) allowed number of boxes for both types of instances is also
with known coordinates (x, y, z) and known dimensions (w, h, shown in Table 2a.
d), then the obstacles are fixed into the bin at the beginning of
The second set of boxes is generated randomly based on
the algorithm. The packing is created from the sequence triple
uniform distribution and includes 10 types of boxes. Two
(A, B, C) and those boxes that overlap with any obstacles in the
instances are created by using this set, which includes 50 boxes.
set are removed. The bin free volume is calculated as follows:
The width, height, and depth of these boxes are selected from
Bin free volume = volume of bin – total volume of obstacles the intervals [100, 250], [50, 250], [100, 300] respectively. The
1298