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 scholarship@uwindsor.ca. 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 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 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 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 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

Three-Dimensional Container Loading: A Simulated Annealing Approach

Similar Three-Dimensional Container Loading: A Simulated Annealing Approach