A dichotomous search-based heuristic for the three-dimensional sphere packing problem

: In this paper, the three-dimensional sphere packing problem is solved by using a dichotomous search-based heuristic. An instance of the problem is defined by a set of n unequal spheres and an object of fixed width and height and, unlimited length. Each sphere is characterized by its radius and the aim of the problem is to optimize the length of the object containing all spheres without overlapping. The proposed method is based upon beam search, in which three complementary phases are combined: (i) a greedy selection phase which determines a series of eligible search subspace, (ii) a truncated tree search, using a width-beam search, that explores some promising paths, and (iii) a dichotomous search that diversifies the search. The performance of the proposed method is evaluated on benchmark instances taken from the literature where its obtained results are compared to those reached by some recent methods of the literature. The proposed method is competitive and it yields promising results.


Introduction
Cutting and Packing problems (CP) is a family of natural combinatorial optimization problems, admitted in numerous real-world applications from industrial engineering, logistics, manufacturing, production process, automated planning, etc. One of the more recent papers addressing an optimization with a packing problem is due to Sutou and Dai (2002), where the unequal sphere problem has been used for tackling an application of the automated radiosurgical treatment planning. Wang (1999) has also considered sphere packing problems as an optimization tool for the radiosurgical treatment planning.

ABOUT THE AUTHORS
Mhand Hifi is a full professor of Computer Science and Operations Research at the University of Picardie Jules Verne, France. He is the head of the laboratory EPROAD (Eco-PRocédés, Optimisation et Aide à la Décision) of the University of Picardie Jules Verne. He is area editor for several international journals: COR, IJMOR, AOR, AJOR, etc. His research interest is NP Hard combinatorial optimisation (sequential and parallel approaches) applied to logistic and OR problems.
Labib Yousef is a PhD student at the University of Picardie Jules Verne. He received his BS in Management at the University of Techrine, Syria. He got his MS in Industrial Engineering (Management de l'Innovation Technologique) from Grenoble INP, France.

PUBLIC INTEREST STATEMENT
Cutting and Packing problems (CP) is a family of natural combinatorial optimization problems, admitted in numerous real-world applications from industrial engineering, logistics, manufacturing, production process, automated planning, etc. This paper tackles a version of the sphere packing problem (belonging to the CP family) with a dichotomous search-based method.
Such a problem can be encountered in the radiosurgical treatment where the unequal sphere problem has been used for determining an automated radiosurgical treatment planning.
Other problems of the CP family have been described and redefined in Wascher, Haussner, and Schumann (2007) where an instance of these problems can be defined by a set of predetermined items to be packed in one or many larger containers (objects) so as to minimize the unused area/space or in some cases to maximize a utility function. Furthermore, the items are bounded by their dimensions (rectangular, circular, or irregular) and the objects can be bounded (rectangular, circular, … ) or unbounded (strips/parallelepipeds, … ). This paper deals with the problem of packing spheres (called items) in a container (parallelepiped or object) with unlimited length. The problem is known as the Three-Dimensional Sphere Packing Problem (3DSPP). An instance of 3DSPP is defined by a set N of n unequal items and an object  of fixed width W and height H and, unlimited length (for the rest of the paper, the length representing the variation of P's length is denoted L). Moreover, each item i ∈ N = {1, … ,n} is characterized by its radius r i . The goal of the problem is to optimize the length L of the object  such that all items of N are packed (positioned) in the target object, without overlapping.
The 3DSPP may be formulated as follows: where the objective function (1) minimizes the length of the object  containing all the items of N, Equation (2) ensures the nonoverlap constraint of any pair of distinct items (i,j) of N; that is, the distance between the centers of both items which must be greater than or equal to the sum of their radii and, Equations (3)-(5) ensure that all items of N belong to the target object  of dimensions (L,W,H). Note that since the aim of the problem is to find the smallest length of the object containing all items of N, then it is easy to start any method trying to solve the above problem by a trivial value representing the sum of the spheres' area affected to L (Equation (6)) and a quick (feasible) solution value with a greedy algorithm affected to L.
The rest of the paper is organized as follows. Section 2 gives a literature review for the 3DSPP and its variants. The problem representation is discussed in Section 3.1. Section 3.2 exposes the principle of the greedy selection phase which serves to determine a subset of eligible positions employed for packing the predefined items in the target object. Section 3.3 presents a truncated tree search using a width beam search phase that serves to explore some promising paths. Section 3.3.2 exposes the dichotomous search where the first two phases are embedded in an interval search in order to diversify the search process. Section 4 proposes an experimental part in which the performance of the proposed algorithm is given: its provided results are compared to those reached by some best algorithms of the literature. Finally, Section 5 concludes by summarizing the contribution of this paper.
L ≤ L ≤ L http://dx.doi.org/10.1080/23311916.2014.994257 Among existing papers addressing sphere packing problems, we cite Lochmann, Oger, and Stoyan (2006) who proposed a statistical analysis for packing random spheres with variable radius distribution. Li and Ji (2013) discussed a dynamics-based collective method for random sphere packing and they tried to apply it to the problem of packing sphere into a cylinder container. In this paper, the stability of the method and its convergence were tackled. Farr (2013) studied the problem of random close packing fractions of lognormal distributions of hard spheres. The author tailored a one-directional approach in order to predict a close packing of spheres of lognormal distributions of sphere sizes.
Packing spheres into a container has been addressed by Sutou and Dai (2002) who proposed a global optimization approach. Stoyan et al. (2003) designed a mathematical model in order to pack spheres into an open container, where both height and widths are fixed whereas the length is unfixed. They use a neighboring search based upon extremum points in order to construct and improve a series of solutions. M'Hallah, Alkandari, and Mladenović (2013) proposed a heuristic based on combining VNS with nonlinear programming solver. The method iterates some moves of the current configuration and complete the partial configuration with a solver dedicated for nonlinear programs. M'Hallah and Alkandari (2012) considered the principle used in M' Hallah et al. (2013) to solve the problem of packing identical spheres into the smallest containing cube. Soontrapa and Chen (2013) tackled the problem of packing identical spheres into a smallest containing sphere by using a random search according to Monte Carlo's method. Finally, Birgin and Sobral (2008) proposed twicedifferentiable nonlinear programming models for the problem of packing both circles and spheres into different containers where the containers may be circular, rectangular, etc. In order to find a global solution for their proposed models, ALGENCAN solver was used for generating a multiple starts solutions (an extensive efficient models and methods for packing both circular and sphere problems were reviewed in Hifi & M'Hallah, 2009b).
In this paper, we propose a three phase solution procedure that combines the width-beam search and a greedy selection procedure for approximately solving the 3DSPP. The last two phases are embedded in a dichotomous procedure that is used for searching the best length L * for the object .
Note that the proposed algorithm can be viewed as an adaptation of the method proposed in Hifi et al. (2008) for solving the two-staged two-dimensional cutting problem and later in Akeb and Hifi (2008) for solving another variant of CP family. The novelty in this approach is based upon (i) a new representation used for defining a restricted space search; that is, a structured set based on target eligible positions which serves to limit the search process without degrading the quality of solutions (Sections 3.1 and 3.2) and (ii) a nice dichotomous search (Section 3.4) that tries to avoid premature convergence toward local stagnant solutions whenever the width-beam search is applied (Section 3.3).

Tackling 3DSPP with a dichotomous search
Herein, the problem representation and strategies used by the proposed algorithm are discussed. First, Section 3.1 describes how an item can be represented in a space and which eligible positions may be chosen for assigning items of N to the target object . Second, Section 3.2 describes the principle of the simple constructive procedure which tries to guaranty a feasible packing of all items of N for the 3DSPP according to the set of eligible positions. Such a simple constructive procedure is based on the minimum distance between items and the edges of the objet  that serves to create a set of eligible positions for the current item to pack. Third and last, Section 3.3 discusses the principle of the proposed algorithm and its main steps.

Representation of 3DSPP
In this section, we present the principle of the local strategy used for filtering the huge number of possible positions when an item will be assigned to the object . The local strategy is based on the simple Greedy Principle (called GP) where the minimum distance position is favored for packing a predefined item. GP is then used as an evaluation operator for finding a subset of eligible positions of the next item to pack. Furthermore, the GP can also be used either as a complete solution procedure that provides a final feasible solution for the 3DSPP or as a greedy strategy in order to complete a partial solution.
For the rest of the paper, the following notations are introduced: • The bottom-left-depth corner of  is positioned at (0,0,0) and  is characterized by a set formed with six labels (namely faces): = {left, top, right, bottom, depth, front}. Then,  is represented in the Euclidian space, as illustrated in Figure 1.
• The center of the ith item belonging to N is positioned at (x i ,y i ,z i ).
• The distance, namely i,j , between two different items i and j belonging to N is computed as follows: Observe that positioning any item i of N into  induces to define a distance between an item belonging to N and faces of . Indeed, for instance, assigning an item i ∈ N to an eligible position of  while respecting the nonoverlap of the left-face of  requires to satisfy the following distance: Table 1 reports the distance to be satisfied whenever a preselected item i of N is assigned to an eligible position of .

Defining eligible positions
It is well known that tailored heuristics are mainly based on the strategies which are able to guide well the search process. These strategies may be using some selection criteria in order to provide either partial or final solutions for the problem to solve. Herein, we consider a simple GP which is based on searching the position realizing the minimum distance position between items and feces. In fact, GP is used as a selection criterion for defining a set of eligible positions to assign to the predefined item i (not already positioned) among all eligible positions representing the set P I i .
In what follows, we assume that the (center of the) first item i = 1 of N is positioned at the position (r 1 ,r 1 ,r 1 ) and ∀i ∈ N, i ≥ 2, the following notations are considered:  • I i denotes the set of items of N already positioned in the current object .
• I i contains the items of N which are not yet assigned to .
• P I i denotes the set of distinct eligible positions for the next item i to pack given the set of packed items I i .
• An eligible position p i+1 ∈ P I i (for the item i) is determined by using three elements e 1 , e 2 and e 3 where an element is either an item of N already positioned (representing the set I i ) or one of the six faces belonging to .
represents the set composed of the three elements e 1 , e 2 and e 3 . Figure 1 illustrates the mechanism used by the GP strategy on a small example. We assume that the first four items are already positioned in the object ; then, there are three eligible positions that emerge for the next item 5 to pack. Following the above notations, I 4 = {1,2,3,4} and P I 4 = {p j 5 , j = 1, … ,3}. First, the position p 1 5 touches both items 1 and 3 and, the "bottom" face of . Second, the position p 2 5 is obtained by using the item 2 and both faces "left" and "bottom" of . Third and last, the position p 3 5 is computed by using the item 1 and both faces "left" and "depth" of . Finally, it follows that  p 1 We can observe that the item 4 is positioned around the three items 1-3 that means  4 = {1,2,3}.
For the next item 5 to pack, its coordinates are computed by using both items 1 and 3 and one of the faces (the "bottom" in this case) of  that gives  5 = {1,3,bottom}. We recall that the objective of the problem is to minimize the length of the target object . It means that the right face of  can be omitted and only the five faces can be considered when optimizing the length of the target . Hence, all eligible positions may be obtained by using the fifth faces, the positioned items and the current item to pack.
Moreover, GP works as follows: the value corresponding to the (i + 1)th item to pack, when positioned at the eligible position p k i+1 ∈ P I i , is computed as follows: where Finally, when GP is used, then it starts by positioning the first item i = 1 at the bottom-left-depth position, i.e. at the position (r 1 ,r 1 ,r 1 ), while the remaining n − 1 items are successively positioned according to the minimum distance rule (cf. Equation (8)). As illustrated in Figure 1, the item 5 will be placed at position p 1 5 because its corresponding distance realizes the minimum value (of course, this choice may also be modified along the execution of the algorithm).
Of course, one can use GP as a greedy procedure for searching a feasible solution for the 3DSPP. Indeed, by setting the starting length of the target object  to ∞ (i.e. L  = ∞) and by applying successively GP on the rest of the nonpacked items, GP terminates by building a feasible solution with length L U . The provided length L U can also be considered as an upper bound for starting the proposed algorithm.

A width-beam search heuristic for the 3DSPP
In this section, we first describe the principle of the standard beam search (Section 3.3.1). Second and last, we show how beam search can be adapted for solving the 3DSPP (Section 3.3.2), especially when using a width-beam almost of the standard beam.

A standard beam search
Beam Search (BS) has been first proposed in Ow and Morton (1988) for tackling the scheduling problem and it has since been successfully applied to many other combinatorial optimization problems (some adaptations can be found in Della Croce, Ghirardi, & Tadei, 2004;Hifi & M'Hallah, 2009a;(8) (i+1,j) Hifi et al., 2008;Hifi & Saadi, 2010;Yavuza, 2013). Such an approach can be viewed as a truncated tree search procedure where its objective is to avoid exhaustive search by performing a partial enumeration of the solution space.
At each level of the developed tree, only a subset of nodes (called the set of elite nodes) are selected for further branching and the other nodes are discarded, where no backtracking is performed. For each level, the cardinality of the elite nodes to be investigated is fixed to that is called the beam width. Generally, these selected nodes represent those having a high potential to lead the best solutions for the treated problem. Furthermore, each node is assessed via an evaluation function whose role is to provide a promising separation mechanism of the nodes of each level of the developed tree.
Algorithm 1 describes the main steps of the standard beam search which is characterized by a beam width used for filtering the set of offspring nodes B . A node corresponds to a partial feasible solution and the set B of current nodes is initialized to the root node B 0 whereas B containing the offspring nodes is initialized to the empty set. Each node of B generates a set of offspring nodes, and moves them to B . If a node of B is a leaf (i.e. no further branching is possible out of ), then its objective function value z is computed and compared to z ⋆ . If z > z ⋆ , then the incumbent solution is set to the leaf node; z ⋆ is then updated: z ⋆ = z ; and is removed from B . All nodes belonging to B are assessed using an evaluation function, and then ranked in decreasing order of their values. The first nodes of B are then considered as the elite nodes which are appended to the set B, whereas the remaining nodes of B are discarded and B is reset to the empty set. This process is iterated until no further branching is possible, i.e. until B = �.

Adaptation of the beam search for the 3DSPP
Applying BS to 3DSPP requires defining the nodes of the tree and the branching mechanism out of the nodes of B. Herein, a node i is represented by the pair of subsets: (1) The first subset I i containing all items assigned to the target object  and, (2) The complementary subset I i containing the unassigned items.
The first subset represents a partial solution realized when assigning all the items to the object  and the second subset serves to provide a complementary solution. Such a complementary solution is computed by applying the MLDP used as greedy procedure (as discussed in Section 3.2).
Moreover, branching out of a node i is equivalent to create at most |P I i | branches emanating out of the current node (related to the eligible positions as described in Section 3.2). Each resulting node corresponds to packing the subset of items I i and assigning to the current item i a favorite eligible position. Moreover, each of these created nodes will be represented by a pair of two complementary subsets of items of N.
As mentioned above, because of the huge number of feasible positions that a predefined item may generate, we preferred the use of the width-beam search almost of the standard beam search. Indeed, as we shall see in Section 3.4, the proposed algorithm uses a dichotomous search according to a series of values assigned to (the beam width). Therefore, in order to try the exploration of a maximum number of promising paths, the beam search applies the width search where all nodes belonging to the same level are simultaneously evaluated following an estimator operator and only the best ones are selected for the rest of the search. The aforementioned process is described by Algorithm 2. The adapted algorithm works according to a given node, namely , taken at the level of the developed tree. Thereafter, the initialization step (Lines 1-3) is applied for starting the set B containing the best provided nodes regarding the starting node and, the initialization of the variable feasible to false. This variable controls the (un) feasibility of the series of the solutions built. The main loop (lines 7) starts by choosing the best eligible positions for each node belonging to B. These positions are computed by using GP's selection (cf. Section 3.2). Second, all created nodes are stored in a provisional set B where the potential of each of these nodes are evaluated according to the final solution provided by iteratively applying GP as a heuristic (cf. as discussed in the last paragraph of Section 3.2). Thereafter, for each fin al solution (either feasible or unfeasible for the target object ), the potential of a node ∈ B is represented by the density of the positioned items in . Whenever one of these constructed solutions provides a feasible solution (line 10), i.e. all items are positioned in the target object , then the algorithm stops with a feasible solution (i.e. setting the variable feasible to true). Otherwise, the set B of the best nodes is updated (line 11) by the nodes which realizing the highest densities and the current level of the developed tree is incremented. Finally, the process is iterated until either B is reduced to an empty set or when the fixed runtime limit is performed.

Using a dichotomous search
The proposed algorithm applies a series of BS for a predefined interval search. First, the interval search starts by [L,L], where L denotes a lower bound for the 3DSPP and L its upper bound (in the case of a feasible solution exists, its objective value is assigned to L). Second, for each fixed interval, BS tries to construct the best feasible solution by packing all the items into the current target object; bottom-left-depth corner (in the position (r 1 ,r 1 ,r 1 )) and creates its corresponding sets I i , I i and P 1 I (as discussed in Section 3.2). At line 10, BS is called with the target value of the object (L ⋆ ,W,H) and the created sets reached at the next step. Line 11 serves to update the interval search where its upper bound is updated whenever a feasible solution i s obtained, the lower bound is updated otherwise. Thereafter, the process is iterated till the gap between both lower and upper bounds becomes closest to a certain tolerance, namely . Finally, the aforementioned process is iterated a certain number of times following the values of (line 13) and according to the runtime limit fixed.

Computational results
In this section, we investigate the effectiveness of the width beam search-based heuristic (DSBH) on two sets of benchmark instances: Set1 and Set2.
The second set "Set2" contains six instances (KBTG1, KBTG2, KBTG3, KBTG7, KBTG8, and KBTG9) taken from Kubach et al. (2011). For each instance, both dimensions W and H of the object are fixed to 10 whereas the number of the predefined items is fixed to 30 (resp. 50) for the first (resp. last) three instances. Moreover, these six instances have been already tested in Kubach et al. (2011) where they represent the six instances with unequal spheres.
The rest of the experimental part is organized as follows. First, the behavior of the proposed algorithm DSBH is evaluated on the first set of instances Set1 (Section 4.1). The results reached by DSBH are thereafter compared to those reached by four heuristics: Stoyan et al.'s (2003) method-based heuristic (noted SYS), Birgin and Sobral's (2008) method-based heuristic and both sequential and parallel heuristics proposed by Kubach et al.'s (2009Kubach et al.'s ( , 2011) (noted KBTG s and KBTG p , respectively). http://dx.doi.org/10.1080/23311916.2014.994257 Second, for the instances of Set2 (Section 4.2), the results provided by DSBH are compared to those reached by Kubach et al. (2011). Note also that the proposed algorithm was coded in C++ and tested on an Intel Core 2 Duo (2.53 GHz and with 4 GB of RAM) and the runtime limit was fixed to one hour.

Performance of DSBH vs. three heuristics: Set1
Generally, when using approximate algorithms to solve optimization problems, it is well known that different parameter settings for the approach lead to results of variable quality. As discussed in Section 3.3.2, DSBH uses two parameters: the beam width and the maximum runtime limit to fix.
Our computational study was conducted by varying in the discrete interval {5,6,7, …} and the maximum runtime limit was fixed to 3600 s (which can be considered as a standard runtime limit considered by algorithms of the literature). Of course, the upper value of depends on the limited runtime and the size of the instance. In order to show the effect of these parameters, we discuss the quality of the solutions obtained by DSBH when the beam width increases.
First, we do it on an instance extracted from Stoyan et al. (2003) by varying the value of from 1 (the minimum value assigned to in Algorithm 3) to the upper value (incrementing with one unit) till attaining the fixed runtime limit. Figure 2(a) shows the behavior of DSBH when varying the beam width for the instance SYS5. From the curve labeled Sol., one can observe that the value of the length L ⋆ of the target object  oscillates but globally its average solution value decreases and converges toward one of the best solutions (generation of a series of local optima). This is the reason why we proposed to reiterate the GP rule in order to get the final solution realizing the best local optimum. On the other hand, according to the second curve labeled Best, one can observe that the behavior of DSBH provides a series of decreased solution values; in this case, only the best solution was retained for the whole values of . For SYS5, the best local optimum reached by DSBH realizes the value closest to 10.94 and its solution configuration is illustrated by Figure 2(b). In fact, the provided solution is equal to 10.9359 (as shown later in the second part of the experimental part) and it dominates both values 11.2170 and 10.9359 reached by the sequential KBTG algorithm (Kubach et al., 2011) and its parallel version (Kubach et al., 2009) representing the best solution values of the literature. Note that = 12 seems to give the best solution value that balances between the quality of the solution and the fixed runtime since for higher values of the solution remains unchanged.
Columns from 1 to 4 show the instance label, problem size n and both height H and width W of .  the parallel KBTG version (noted L ⋆ KBTG p ) without fixing the runtime limit. Columns 9 and 10 display the solution (noted L ⋆ ) realized by the proposed algorithm DSBH for both fixed runtimes (300 and 3600 s, respectively). Finally, column 11 reports the best value of for which its best solution for the second runtime limit is performed.
All results of Table 2 are summarized in Table 3, where it tallies the percentage improvement (when it happens) yielded by DSBH when compared to the results reached by the four considered algorithms. Table 3 is composed of two parts: the first part contains the informations about the instance and the second part shows the percentage improvement realized by DSBH over the four heuristics tested (noted %SYS, %BSA, %KBTG s and %KBTG p according to the heuristics SYS, BSA, KBTG s and KBTG p , respectively) especially for the second runtime limit for DSBH.
The analysis of the results of both Tables 2 and 3 follows: (1) First, DSBH outperforms the three algorithms SYS, BSA, and KBTG s whenever the runtime limit was fixed to 300 s. Indeed, it is able to reach the best solution for all instances of Set1.
(2) Second, DSBH outperforms SYS, BSA and both KBTG s and KBTG p algorithms since for the instances of Set1 it is able to provide better results.
(3) Third, when comparing DSBH's results to SYS's ones, one can observe that the percentage of the improvement varies from 6.49% (instance SYS5) to 8.82% (instance SYS3). This percentage improvement remains interesting when comparing DSBH's results to those reached by BSA: in this case, such improvement varies between 5.96% (instance SYS1) and 7.65% (instance SYS6). Moreover, it stills interesting when comparing DSBH's results to those provided by KBTG algorithm. Indeed, the improvement varies from 0.48% (instance SYS1) to 6.26% (SYS4).  (4) Fourth and last, DSBH's results remain better than those reached by the parallel algorithm KBTG p (we recall that the parallel algorithm ran without runtime limit). Indeed, all the solution values for the instances of Set1 are improved. In this case, the percentage of improvement varies from 0.14% (instance SYS3) to 1.65% (instance SYS4), which remains very interesting for a sequential algorithm. Figure 3 shows the behavior of DSBH on the instances of Set1 where each curve represents the variation of the improvements realized according to the algorithm SYS, BSA and both KBTG s and KBTG p , respectively. Figure 4(a) illustrates GP's configuration (when applied as a heuristic for providing a starting solution for DSBG) and Figure 4(b) displays DSBH's configuration. From the starting configuration (cf. Figure 4(a)), one can observe the failure of GP when used as a heuristic: in this case, all items are concentrated on the bottom-left-depth position and all positions located along the length of  to the top are omitted. However, the configuration performed by DSBH shows the correction made when GP's process is repeated.

Performance of DSBH versus KBTG heuristic: Set2
This section compares the results reached by DSBH to those reached by KBTG s (note that, for this type of instances, KBTG s realizes the best solution values of the literature). This comparison is performed on the instances of the second group Set2 extracted from Kubach et al. (2011). For this type of instances, instead of determining the minimum length L ⋆ of the target container , Kubach et al. (2011) computed the density of all packed items in the final object . Therefore, in addition to the density, the best length L ⋆ of the final object  corresponding to that density is reported in the same table.
The results realized by the two tested methods (DSBH and KBTG s , respectively) are reported in Table 4. Columns 1-4 display all the characteristics related to the instance. Column 5 reports the solution value (expressed in term of density) realized by KBTG s 's algorithm (extracted from Kubach et al. (2009Kubach et al. ( , 2011). The next three columns report the best length L ⋆ reached by DSBH, its corresponding density and the best value of for which the best solution is reached. Finally, the last column displays the percentage improvement realized by DSBH according to the solution values reached by KBTG s . The analysis of the results of Table 4 follows.  (1) DSBH remains competitive since it improves most solutions reached by KBTG s . Indeed, it is able to improve five out of six best solutions while it matches the other solution (instance KBTG2) when compared to the results reached by KBTG s .
(3) In term of density, DSBH realizes a density of 49.932% which is much better than those realized by KBTG s (48.120%), as illustrated in Figure 5.
Finally, Figure 6(a) illustrates the starting solution reached by iterating GP till packing all the items, whereas Figure 6(b) shows the final solution when DSBH is applied. We can observe that the starting configuration leaves also an important unoccupied area along the length of the object  to the depth. However, the configuration leads a better packing when GP rule is repeated on different eligible positions of the search space.

Conclusion
In this paper, the three-dimensional sphere packing problem was solved with a dichotomous searchbased heuristic. The proposed method is an adaptation of the beam search which combines the following components: (i) a greedy selection phase for performing a local restricted space search that contains eligible positions, (ii) a standard width-beam search phase for exploring some promising paths, and (iii) a dichotomous search for diversifying the search space. Through extensive testing on different instances of the literature, the proposed method is able to improve a significant number of the solutions provided by recent methods available in the literature.