A local search-based method for sphere packing problems

Accepted Manuscript

A Local Search-Based Method for Sphere Packing Problems Mhand Hifi, Labib Yousef PII: DOI: Reference:

S0377-2217(18)30869-5 https://doi.org/10.1016/j.ejor.2018.10.016 EOR 15408

To appear in:

European Journal of Operational Research

Received date: Revised date: Accepted date:

17 January 2017 5 October 2018 9 October 2018

Please cite this article as: Mhand Hifi, Labib Yousef, A Local Search-Based Method for Sphere Packing Problems, European Journal of Operational Research (2018), doi: https://doi.org/10.1016/j.ejor.2018.10.016

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • A new local search method is proposed. • Two versions of the sphere packing problem are tackled. • A continuous optimization-based strategy is used.

AC

CE

PT

ED

M

AN US

CR IP T

• The method succeeds in yielding the best solutions for the most in-stances used.

1

ACCEPTED MANUSCRIPT

A Local Search-Based Method for Sphere Packing Problems Mhand Hifi∗ and

Labib Yousef

Abstract

CR IP T

EPROAD, Universit´e de Picardie Jules Verne 7 rue du Moulin Neuf, 80000 Amiens, France.

ED

M

AN US

In this paper, we study the three-dimensional sphere packing which consists in finding the greatest density of a (sub)set of predefined spheres (small items) into a three-dimensional single container (large object) of given dimensions: cuboid of fixed dimensions or cuboid of variable length. The problem with the cuboid of fixed dimensions (sizes) (called knapsack problem in W¨ ascher et al. [31]) is tackled by applying a local search-based method that combines three main features: (i) a best-local position procedure stage, (ii) an intensification stage and (iii) a diversification stage. The first stage ensures a starting feasible solution using a basic greedy local strategy. The second stage tries to solve a series of decision problems in order to place a subset of complementary spheres. The third stage tries to remove some packed items and to replace them with other spheres. The proposed method is also adapted for solving the problem of packing a set of predefined spheres (small items) into a cuboid of variable length (called open dimension problem in W¨ ascher et al. [31]). The performance of the proposed method is evaluated on a set of benchmark instances taken from the literature, where its results are compared to those reached by recent published methods. The computational results showed that the proposed method remains competitive for both treated problems.

PT

Keywords. Heuristics; cuboid container; optimization; sphere packing.

Introduction

CE

1

AC

Cutting and Packing (called CP) problems belong to the well-known family of natural combinatorial optimization problems (W¨ascher et al. [31]). These problems are admitted in numerous applications from computer science, industrial engineering, logistics, manufacturing, production processes and other real-world applications. Generally, the packing problem aims is to find the greatest density of any packing of a (sub)set of small items into predetermined large objects. The small items may be spheres (three-dimensional) and the large objects as a cuboid of fixed dimensions or a cuboid of variable length. Dealing with the three-dimensional sphere packing can be encountered on modeling granular material systems (cf. Li and Ji [19]), automated radio-surgical treatment planning (cf., Sutou [24]), the storage of cylindrical drums into containers, and other applications in medicine for radio-surgical treatment planning (cf. Wang [29]). ∗

Corresponding author (M. Hifi: [email protected]).

2

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

In this paper, two versions of the three-dimensional sphere packing problem are tackled: cuboid of fixed dimensions and cuboid of variable length problems. Following the characterization used by W¨ ascher et al. [31], the first version of the problem can be considered as the Knapsack Problem (KP), where there are an output (value) maximization (assignment), one large object with fixed dimensions (assortment of the large object) with a set of (strongly) heterogeneous small items (spheres). Further, the second version of the problem can be considered as the Open Dimension Problem (ODP), where there are an input (value) minimization (assignment), one large object with unfixed length (dimension: assortment of the large object) with a set of arbitrary small items (spheres). The dimensionality of both studied problems are three, then for the rest of the paper we denote the first version of the problem as 3DKP and 3DODP for the second version. An instance of 3DKP is characterized by a container P of fixed dimensions L × H × W and a set N of n spheres (small items) of radii ri , i ∈ N = {1, . . . , n}. The goal of the problem is to find the maximum density (occupied volume: output maximization) of the container P with the available spheres of N . A feasible packing is a configuration of packed spheres without overlapping between spheres which have to be inside the container. The second problem 3DODP is characterized by a container (large object with one unfixed dimension) of dimensions (L, H, W ), a set N of n spheres of radii ri , i ∈ N, and the goal is to pack all spheres of N into the smallest container, i.e., finding the smallest length (input minimization) of the container that is able to hold all items without overlapping. By using a simple adaptation of the model described in George et al. [7] (cf. Hifi and M’Hallah [9]), 3DKP can be formulated as the following non-linear mixed integer programming: X max α ri3 ξi (1) i∈N

s.t.

ξi ri ≤ xi ≤ ξi (L − ri ), i ∈ N,

ξi ri ≤ yi ≤ ξi (H − ri ), i ∈ N,

(2)

ED

(3)

PT

ξi ri ≤ zi ≤ ξi (W − ri ), i ∈ N, q ξi ξj (ri + rj ) ≤ (xi − xj ) + (yi − yj ) + (zi − zj ), i ∈ N, j ∈ N, i < j,

ξi ∈ {0, 1}, i ∈ N,

(4) (5) (6)

AC

CE

where the objective function (1) favors the subset of spheres maximizing the density such that α = (4π)/(3LHW ). Inequalities from (2) to (4) ensure that the spheres do not go beyond the container. Inequality (5) ensures the non-overlapping between each couple of spheres belonging to the container. Each binary variable ξi , i ∈ N (the last line (6)) takes the value 1 if the i-th sphere is placed in the container; 0 otherwise. One can observe that the above model has n(n − 1)/2 non-convex constraints and 6 × n linear constraints. On the other hand, each of these constraints is related to the binary decision variables, which makes the problem more complex to solve. The second problem 3DODP can be modeled as follows: Minimize

L

(7)

δij ≥ (ri + rj ), ∀i ∈ N, j ∈ N, i < j,

ri ≤ xi ≤ L − ri , i ∈ N,

ri ≤ yi ≤ H − ri , i ∈ N,

ri ≤ zi ≤ W − ri , i ∈ N, 3

(8) (9) (10) (11)

ACCEPTED MANUSCRIPT

2

PT

ED

M

AN US

CR IP T

p where δij = (xi − xj )2 + (yi − yj )2 + (zi − zj )2 , i ∈ N, j ∈ N, i < j. The objective function is L in the optimization problem (7). Inequalities (8) are quadratic constraints that ensure the non overlapping between any pair of distinct spheres. Finally, inequalities from (9) to (11) are the linear constraints of the problem that represent the feasibility of spheres’ positions when positioned into the container of variable length L. In this paper, we first propose to solve 3DKP with a local search-based method. Next, an adaptation of the method is considered for tackling 3DODP. The proposed method combines three main stages: (i) basic greedy procedure (noted BGP), (ii) intensification stage and (iii) diversification stage. BGP provides a starting feasible solution; that is, a constructive procedure which tries to fix, at each step, a sphere until a final feasible solution is achieved. Furthermore, in order to improve the quality of a given solution, we propose two procedures in this stage: (i) an insertion procedure is used in order to improve the packing density of starting solution, and (ii) a continuous local optimization procedure is used to repair the unfeasible solution produced during the insertion procedure. Furthermore, to jump out of the local optimum produced by the previous stage and to guide the search process towards other search spaces, we diversify the search by using two procedures: (i) drop procedure and (ii) re-build procedure. Finally, we propose an adaptation of the local search-based method for solving 3DODP. The remainder of the paper is organized as follows. Section 2 gives a literature review for the sphere packing problem. A summarized version of the local search-based method is exposed in Section 3, where three stages are used. Section 3.1 gives the basic procedure used to build a starting feasible solution. Section 3.2 details the intensification strategies which is used for improving the quality of the current feasible solution. Section 3.3 presents the strategies used in the diversification stage in order to avoid getting trapped in a local optimum. Section 3.3 provides an overview of the proposed method. Section 4 evaluates the performance of the proposed algorithm by comparing its obtained results to those reached by known algorithms available in the literature. Section 4.1 evaluates the performance of the proposed method on the instances related to 3DKP. Section 4.2 discusses the effectiveness of the method when adapted for 3DODP. Finally, Section 5 concludes by summarizing the contribution of the paper.

Related works

AC

CE

The sphere packing problem arises in several branches of the industry, like, modeling random geometric objects used for representing the structure of liquids and glassy materials (cf., Stoyan and Chugay [25]). Such a problem can also be used for studying phenomena such as electrical conductivity, stress distribution, medicine applications for the radio-surgical treatment planning (cf. Wang [29]). Several other variants of the sphere packing problem have been described and redefined in Hifi and M’Hallah [9], where an instance of these problems can be represented by a set of predetermined items to be packed in one or many larger containers so as to minimize the unused area / space or in some cases to maximize a utility function. Because the packing problem and its variants are generally hard to solve, they are often solved by using approximate methods. To the best of our knowledge, very few papers tackling the three-dimensional sphere packing are available in the literature. One of the paper studying the problem of packing of unequal spheres into a parallelepiped of unfixed height is due to Stoyan et al. [26]. The authors designed a model containing quadratic and linear constraints and, a linear objective function, where it has been heuristically 4

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

solved using a special tree-search procedure that includes a neighborhood search in order to move towards a global optimum. More recently, Scheithauer et al. [21] proposed a generic approach-based algorithm for tackling packing unequal spheres into containers of different types. Hifi and M’Hallah [11] proposed a hybrid method in which beam-search and non-linear programming tools are combined for solving the (two-dimensional) unequal circular packing into the smallest containing circle. M’Hallah et al. [17] proposed a variant of the method discussed in Hifi and M’Hallah [11] for tackling the problem of packing identical spheres into spherical container. In this case, a variable neighborhood search combined with non-linear programming is used for tackling the problem. M’Hallah et al. [18] proposed a straightforward adaptation of the aforementioned method to solve the unit spheres into the smallest container. Hifi and M’Hallah [12] tackled the circular cutting problem (packing circles into a rectangular container) by using a genetic algorithm-based heuristic. The algorithm uses an adapted version of the best local position to check the feasibility of a chromosome, to mutate any infeasible chromosome into a feasible one, and to evaluate the fitness of a feasible chromosome. Huang et al. [16] proposed a greedy algorithm for packing unequal circles into a rectangle, where the hole degree strategy was used for positioning the next circle into the rectangle. More recently, Stetsyuk et al. [23] studied a special packing problem where a set of circles is given and its aims is to pack all these circles into a larger circle of minimum radius. The problem has been tackled as a multiextremal nonlinear programming problem, where it was reduced to unconstrained minimization problem of a nonsmooth function by means of nonsmooth penalty functions. A nice efficient algorithm to search for local extrema and a strategy-based procedure for improving a series of (lower) bounds have been proposed. Note that their proposed procedures uses nonsmooth optimization methods that are based upon Shor’s r-algorithm. Kubach et al. [5, 6] solved two variants of the sphere packing problem, where unequal spheres are to be packed into both cuboid of fixed dimensions and cuboid of variable length. In [5], the authors developed a nice solution procedure, which is based on a simple lookahead strategy which is considered as the core of the algorithm. In [6], the authors proposed a parallel algorithm that can be considered as a straightforward method already considered in [5]. Akeb [1] combines a two-stage look-ahead strategy with Kubach et al.’s [5] method for solving the cuboid of fixed dimensions problem. In [2], the same method has been employed for solving the cuboid of variable length problem Birgin et al. [3] proposed twice differentiable non-linear programming models for both two and three dimensional (circles and spheres) packing into a circular (resp. rectangular) container. In this case, ALGENCAN solver has been used as a powerful tool for approximating the model. Sutou et al. [24] considered the sphere packing problem as a model for tackling the automated radio-surgical treatment planning. This problem was formulated as a nonconvex optimization problem with linear objective function and quadratic constraints. The goal of the problem is to pack a maximum number of spheres that minimize the wastage (unoccupied volume). They proposed an approximate method that is based upon a global optimization approach. Li and Ji [19] used a sphere packing into a cylindrical 3-D system to simulate the granular materials and then a special quasi-dynamics method was proposed in order to generate densely distributed spheres into an enclosed container. Also, Farr [4] studied the random close packing fractions of lognormal distributions in order to simulate granular and 5

ACCEPTED MANUSCRIPT

CR IP T

mesoscopic systems. This problem is related to granular and mesoscopic systems, where various material properties depend on the close packed volume fraction of the constituent particles. He proposed to solve the problem by adapting a special one-dimensional algorithm for predicting the random close packing fraction. Soontrapa and Chen [22] considered the sphere packing process for modeling the carbon atoms. They proposed an adaptive random search technique to achieve the configurations that the total potential energy must be in the global minimum state. Finally, Hifi and Yousef [13, 14, 15] proposed several methods for tackling the unequal sphere packing by combining a series of iterative methods, like, dichotomous search, tree search, hill-climbing and other ones. The authors underlined the importance of some hybridization related to strategies, especially whenever the truncated tree search and local search are combined.

3

A local search-based method for 3DKP

AN US

This section exposes the local search-based method for the sphere packing problem. It starts by describing the method applied for solving 3DKP and later its adaptation for solving 3DODP. The main principle of the local search-based method can be summarized as follows: (i) Starting the search process by an initial solution using a basic greedy procedure (cf. Section 3.1).

M

(ii) Building an improved solution using a series of decision problems (cf. Section 3.2). (iii) Perturbing the search process and re-constructing a new current solution using a basic greedy procedure according to the new order (cf. Section 3.3).

ED

(iv) Steps (ii)-(iii) are repeated until a satisfactory solution is reached.

Building a feasible solution

CE

3.1

PT

In what follows, Section 3.1 exposes the used greedy procedure that yields a feasible solution. Section 3.2 discusses the intensification stage which tries to improve the solution yielded after each iteration. Finally, Section 3.3 details the diversification stage which tries to explore other local optima and to enhance the solutions built.

AC

Generating a solution is equivalent to generate a subset of positions of the spheres in the container P. It can be done by applying the so-called Basic Greedy Procedure (BGP) (cf. Hifi and Yousef [13, 14]) that is based on positioning a series of spheres step by step following the best local position of all candidate spheres. Such a procedure can also be used for building either “a partial (incomplete) solution” or “a complementary” one. Such a strategy can be viewed as an extended version of the hole degree strategy used by Huang et al. [16] for solving the circular packing problem. Let suppose that all spheres are reorganized in decreasing order of their radii. Assume that the bottom-left-back corner of P is positioned at (0, 0, 0), where P is presented by a cuboid bounded by six faces F = {left, top, right, bottom, back, front} (as illustrated in Figure 1). Then, P may be represented in the Euclidean space as follows: • Each placed sphere i ∈ N is centered at the position (xi , yi , zi ). 6

ACCEPTED MANUSCRIPT

• A pair (i, j), i ∈ N and j ∈ N, are characterized by their distance ∆i,j , where it represents a normalized φ-function for two-spheres (as described in Stoyan and Romanova [20]): q ∆i,j = (xi − xj )2 + (yi − yj )2 + (zi − zj )2 − (ri + rj ), ∀(i, j), i ∈ N, j ∈ N, i < j, (12)

AN US

CR IP T

where both i and j can be positioned when ∆i,j ≥ 0 and constraints from (2) to (4) are satisfied.

M

Figure 1: Illustration of (a) the first positioned sphere (1) at position (r1 , r1 , r1 ), (b) the second (candidate) sphere (2) to pack with its eligible positions p12 , p22 and p32 .

ED

Herein, “a partial solution” is equivalent to a subset of positioned spheres in the container P. The first sphere i = 1, where i ∈ N , is packed at the position (ri , ri , ri ) = (r1 , r1 , r1 ) (cf. Figure 1(a)). Let i ∈ N (i ≥ 2,) be the next sphere to be packed into P and PIi be the subset of its eligible positions. Therefore, the following notations are considered:

PT

• Ii : the set of spheres of N already packed in the container P (Figure 1(a) illustrates the first (position) set I1 = {1}).

CE

• Ii : the set of spheres of N which are not yet packed in the container P (Figure 1(b) shows the induced set I1 = {2, 3, . . . , n}).

AC

• PIi : the set of distinct eligible positions for packing the next sphere i + 1. Such a set PIi is built regarding the set Ii (herein, PI1 = {p12 , p22 , p32 }, as illustrated by Figure 1(b)).

An eligible position pki+1 ∈ PIi (k denotes the number of possible positions for positioning the next sphere i + 1) is that touching three different elements e1 , e2 and e3 . An element is either a sphere of N already positioned (representing Ii ) or one of the six faces belonging to F (as illustrated in Figure 1(b), p12 touches the sphere 1 and both left and back faces belonging to F. So, assume that Tpk denotes the set of these three elements (Tp12 = {1, lef t, back}, i+1 Tp22 = {1, lef t, bottom} and Tp32 = {1, back, bottom}). Then, a step of BGP, for i ≥ 2, i ∈ N, proceeds as follows: the distance β(pki+1 , .), pki+1 ∈ PIi , associated to each position for packing the next sphere (namely i + 1), realizes what

7

ACCEPTED MANUSCRIPT

follows: β(pki+1 , j) =  min ∆(pki+1 , j), ∀k.

(13)

j∈Ii ∨j∈F j6∈T k p i+1

Table 1 reports the distance ∆ when the next item is positioned at position p and touches one of the faces of the container P. ∆(p, j), j ∈ F xp − rp yp − rp zp − rp L − xp − rp H − yp − rp W − zp − rp

AN US

j left bottom back right top front

CR IP T

Table 1: The distance between a position p and a face j

ED

M

Hence, BGP starts by positioning the first sphere i = 1 at the bottom-left-back position, i.e., at the position (r1 , r1 , r1 ), while the remaining n − 1 spheres are successively positioned according to β(., .) (cf., Equality (13)). Note that if different eligible positions realize the same value (according to equation (13)), then one can choose the eligible position that realizes the smallest distance to the left side of P. Algorithm 1 describes the framework of the basic greedy procedure, used for reaching a starting feasible solution. BGP (Algorithm 1) can be viewed as a greedy procedure that tries to pack all spheres into the container P. Whenever all the spheres are placed successfully, then an optimal solution is reached for the problem. Otherwise, a subset of spheres is placed with a density of D. Indeed, according to the instances of the literature, when can observe that in some cases the instances used may be easy to solve instead of other ones for which computing good approximations remains difficult. Algorithm 1 . Basic Greedy Procedure (BGP)

PT

Input. An instance of 3DKP.

Output. A feasible solution of density D with a subset Ii ⊆ N of the packed spheres.

CE

Set i := 1. Let Ii = {i}, Ii = N \ Ii and PIi be the set of eligible positions of the next sphere to pack. Update PIi with the eligible positions of each j ∈ Ii . repeat Let pki+1 be the best position for the next sphere according to β(., .) (cf., Equality (13))). Position the (i + 1)th sphere into pki+1 and move the packed sphere to Ii . Update the set of eligible positions PIi . until (PIi = ∅) Exit with a feasible solution of density D and the subset Ii ⊆ N containing the packed spheres.

AC

1: 2: 3: 4: 5: 6: 7: 8: 9:

Indeed, in order to show the hardness of some used instances and how BGP can fill the cuboid of fixed dimensions, two instances taken from the literature (KBG2 with 30 spheres to pack (taken from Kubach et al. [5]) and SYS1KP with 25 spheres to pack (taken from Akeb [1])) are considered for illustrating the phenomenon. Figure 2 illustrates KBG2’s solution provided by BGP whereas Figure 3 shows that of SYS1KP. As one can see, the greedy procedure is able to pack all spheres (cf. Figure 2) realizing a density of 30.0709%

8

Figure 3: Illustration of BGP’s solution for the instance SYS1KP (taken from Akeb [1]) containing 25 spheres (density of the provided solution is equal to 48.4266%).

AN US

Figure 2: Illustration of BGP’s solution for the instance KBG2 (taken from Kubach et al. [5]) containing 30 spheres (density of the provided solution is equal to 30.0709%).

CR IP T

ACCEPTED MANUSCRIPT

(that is an optimal solution) whereas it packs 12 spheres over 25 for the second instance (cf. Figure 3) by realizing a density of 48.4266% (that is an approximate solution).

Improving the quality of solutions

M

3.2

ED

In order to escape from local optima, Wang et al. [30] proposed the so-called quasi-physical quasi-human algorithm. It tries simulate the physical model for tackling the circular packing (packing a set of circles into a containing circle of minimum radius). The approach is composed of two stages according to a set of items (circles) randomly packed into a containing circle:

AC

CE

PT

• The first stage applies the quasi-physical approach that solves a decision problem which checks if a feasible packing of the available circles in the containing circle is possible. Such an approach uses a continuous optimization in order to reduce the overlapping depth related to overlapping between circle(s) and circle(s) and, between circle(s) and container. The quasi-human approach (the second stage) is then employed to make jumps in order to get out of local trapped minima. Because the used process is often trapped in the unfeasible local optima, circles are allowed to relocate new randomly positions according to the packed circles. • The second stage recalls the first stage after reducing the radius of the containing circle (in this case, a simple dichotomous search may be used as a post-procedure for solving this problem). Such process is iterated till obtaining the final container of minimum length.

In our study, we use an adaptation of Wang et al’s approach as a part of the local search-based method for both 3DKP and 3DODP. The proposed adaptation can be viewed as a two-step procedure:

9

ACCEPTED MANUSCRIPT

(i) An insertion procedure is used in order to improve the packing density of starting solution, and (ii) a continuous local optimization procedure is used to repair the unfeasible solution produced during the insertion procedure.

3.2.1

CR IP T

Note that recently, Stoyan et al. [28] have proposed a nice approach, where a specialized insertion strategy has been also used for improving the quality of the solution at hand. Let S be the current feasible solution (not necessary the starting solution reached by BGP). Suppose that N1 ⊆ N be the subset of spheres positioned into the container P and N2 = N \N1 be the remaining non-packed spheres. Then, the following description explains the two-steps procedure used for determining a better feasible solution when it exists. The first step

AN US

The first step tries to determine the largest “unoccupied spaces” according to eligible positions. From the current solution S, containing the set N1 of positioned spheres, one can determine a set of eligible positions for the non-packed spheres. It can be obtained by applying the strategy used by the basic procedure: to determine some positions inside the container, where for a fictive candidate sphere i of radius rf , its eligible position pi+1 is that touching three elements e1 , e2 and e3 . Each element is either a sphere of N1 already positioned (representing the set Ii ) or one of the six faces belonging to F. We then apply the following steps in order to compute the current eligible positions (unoccupied spaces):

M

1. Let rmax be the greatest radius among the unpacked spheres.

ED

2. Let rf = rmax be the radius of the “fictive sphere” and Lp = ∅ be the set of positions.
3. While rf > rmin −  do (a) Compute all eligible positions according to the “fictive sphere”;

The second step

CE

3.2.2

PT

(b) Update Lp and reduce the radius of the current “fictive sphere” (the experimental part has been realized by setting  = 10−2 and decreasing the fictive sphere’s radius as follows: rf = rf × 0.9).

AC

Determining an improved solution (with better density) is equivalent to solve a decision problem. Indeed, given a subset N1 of spheres of radius ri , i ∈ N , and a container of dimensions (L, W, H), is it possible to pack all N1 spheres inside the container without overlap or not. Such a problem can be modeled as a non-linear program that can be solved by using a continuous optimization procedure (as used by Wang et al. [30]). Its goal is to reduce the overlapping depth to zero. When the optimizer stops with overlapping close to zero, the decision problem exits with the value “true”, “false” otherwise. Of course, the true value means that the decision problem is solved, false otherwise. In order to measure the overlapping depth, Table 2 shows the quantity of each overlap generated between items and the faces of a given closed container. Hence, according to both terms of overlapping, following the couple of spheres or a sphere with the container, the quantity of overlap E associated to the current unfeasible

10

ACCEPTED MANUSCRIPT

Between sphere and container

overlapping depth  Oi,x = max  0, ri + |xi | − 21 L Oi,y = max  0, ri + |yi | − 21 H Oi,z = max 0, ri + |zi | − 12 W Oi,j = max {0, ri + rj − δij }

sphere and sphere

solution is given as follows: E=

N −1 X

N X

2 Oi,j +

i=1 j=i+1

N X i=1

CR IP T

Table 2: The overlapping depth between the positioned spheres and between spheres and the container

(Oi,x + Oi,y + Oi,z )

AC

CE

PT

ED

M

AN US

The decision problem (of packing N1 items into a container) admits a local optimum whenever E is reduced to zero. It means that all spheres of N1 can be packed into the container. Otherwise, the procedure fails to provide a feasible solution for the current subset of spheres N1 . Figure 4 illustrates an instance of such a problem, where the container is considered as an “elastic container” and all the spheres as the “smooth elastic solids” which are forced to be included into the container. According to the elasticity mechanism, there exist some conjugated extrusion elastic force whenever some spheres are deformed by other ones or by the container. In this case, a series of complicated distortions would occur due to these elastic forces. During the search process, the overlapping gap (either between the deformed spheres or that realized between spheres and the container) will be reduced to zero. Reducing such deformations to zero induces the convergence of the process to a local optimum (for more details, the reader can refer to Wang et al. [30]).

Figure 4: Illustration of the decision problem and the mechanism simulated for achieving a feasible solution. Algorithm 2 describes the main steps of the Intensification Procedure (noted IP) which is applied for yielding an extended (better) feasible solution for the problem. It can be viewed as a stepwise insertion procedure where, at each step of the procedure, the insertion 11

ACCEPTED MANUSCRIPT

Algorithm 2 . Intensification Procedure (IP) Input. A starting feasible solution S with a subset N1 (resp. N2 ) of packed (resp. non-packed) spheres. Output. An (improved) solution S ? and (updated) subsets N1 (resp. N2 ) of packed (resp. non-packed) spheres. Set N 2 = N2 and i be the first item of N 2 (according to the current order). while (N 2 6= ∅) do Set O = N1 ∪ {i} and solve the decision problem with items of O. if (all spheres of O are packed) then Let S 0 be the new solution. update S ? with S 0 (of density D? ). Set N1 = O. end if Set N 2 = N 2 \ {i}. end while Exit with the best solution S ? and the updated subsets N1 and N2 = N \ N1 .

CR IP T

1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:

AC

CE

PT

ED

M

AN US

respects the feasibility of the provided solution; that is equivalent to solve the decision problem with items of the current set N1 to pack. We recall that IP applies a sequentially insertion according to the available spheres (belonging to the subset N2 ), one by one, into the existing solution. The main loop (Lines from 2 to 10) stops when all the spheres have been examined. Line 3 serves to solve a decision problem which tries to pack all items of the current set O (by applying Wang et al.’s procedure), where a new subset of spheres is considered (that is composed of the subset N1 augmented with a new sphere belonging to the complementary subset N2 − initialized to N 2 for giving a chance to each sphere initially belonging to N2 ). Lines 6 and 7 check the feasibility of the new solution and update that solution with the new positioned spheres in the case where the decision problem stops with a true decision. Line 9 decreases the size of the set of un-packed spheres. Such a process is repeated till examining all spheres initially unpacked. Finally, the procedure exits (line 11) with the best solution S ? (with density D? ) and both updated subsets of packed and unpacked spheres.

Figure 5: BGP’s solution for the instance SYS1KP (taken from Akeb [1]) with density equals to 48.4266%, representing 12 positioned spheres.

Figure 6: The eligible positions generated according to the fictive sphere.

12

Figure 7: Illustration of the solution provided by the second step of insertion procedure: density of the solution is augmented to 51.0489%,representing 14 positioned spheres.

ACCEPTED MANUSCRIPT

In order to illustrate IP’s mechanism, Figures 5 and 7 display two solutions. The first solution (cf. Figure 5) is constructed by using the greedy procedure BGP and the second one (cf. Figure 7) is obtained after running IP. One can observe that despite the good result that can reach BGP, the used process improves the quality of the solution. Indeed, the density of the new solution reached by IP becomes equal to 51.0489% instead of 48.4266%.

3.3

Using a diversification strategy

AN US

CR IP T

Herein, we propose a diversification strategy that consists of removing a subset of spheres from a feasible solution. The removing strategy diversifies the search process by degrading the quality of the solution at hand with the aim of avoiding stagnating in a local optimum. Then, “a partial solution” is obtained and it is completed by applying the basic procedure BGP as a tool for refining the quality of the current “partial solution”, according to the new order associated to the remaining spheres. Such a strategy has been already used with success for solving variants of the knapsack type problems (cf., Hifi [8] and, Hifi and Michrafy [10]). For the knapsack type problems, several binary variables were randomly removed and replaced by an optimal solution corresponding to specified complementary subproblems: the provided optimal solution is reached by applying the Cplex solver as a blackbox solver. Algorithm 3 . Drop and Rebuild Procedure (DRP)

Input. S, a feasible solution with its subset N1 (resp. N2 ) of packed (resp. non-packed) spheres.

ED

Set S ? ←− S; Apply a random destroying strategy (remove α% of packed spheres) from the solution S ? . Place the α% of spheres after the first unpacked sphere. Let S p be the new provided “partial solution”. Complete S p by using BGP and let S 0 be the new solution reached. Return the solution S 0 (with density D0 ) and the updated both subsets N1 and N2 = N \ N1 .

PT

1: 2: 3: 4: 5: 6:

M

Output. S 0 , an (improved) solution and its (updated) subsets N1 (resp. N2 ) of packed (resp. non-packed) spheres.

AC

CE

The main principle of the diversification procedure is described in Algorithm 3. It starts (Line 2) by degrading the quality of the current solution: destroying the current solution (noted S ∗ ), where α% of spheres belonging to S ∗ are randomly removed. Then, a new order is reconsidered (Line 3) by positioning α% of spheres after the first unpacked sphere; that is the first sphere not positioned in the current solution according to the current order. A new (feasible) “partial solution” S p is reached (Line 4) which is completed by applying the basic procedure BGP (cf. Section 3.1). Because a new order is considered, then the procedure tries to build (Line 5) a new feasible solution, noted S 0 . Finally, the procedure returns the resulting solution S 0 and its updating subsets N1 and N2 . Figures 8 and 9 illustrate (i) “the partial solution” realized when the dropping stage (removing 30% of spheres from the packed ones) is applied, and (ii) the final solution obtained when the rebuilding stage is applied. As shown from Figure 8, the first solution degrades the quality of the solution as expected whereas the second solution (cf. Figure 9) seems as a modified solution. If such a solution is compared to that builded by BGP (cf. Figure 5), one can observe that both solutions differ.

13

3.4

Figure 9: Using the second stage of the diversification strategy: the rebuilding stage (a solution realizing a density of 49.3742%).

AN US

Figure 8: Using the first stage of the diversification strategy: the dropping stage (the solution realizes a density equals to 45.021%, where 30% of packed spheres were removed from the current solution).

CR IP T

ACCEPTED MANUSCRIPT

An overview of the local search-based method

AC

CE

PT

ED

M

Algorithm 4 summarizes the main steps of the proposed Local Search-Based Method (noted LSBM). The input of LSBM is an instance of the packing problem composed of a set N of n spheres that would be packed into a container of dimensions (L, W, H). It starts (Line 1) by ordering the spheres of the set N in decreasing order of their radii. At line 2, the first initial feasible solution S 0 is reached by applying the basic procedure BGP (cf. Section 3.1). Lines from 4 to 15 represent the main loop of CM, which is composed of both intensification and diversification stages. The intensification stage (Lines from 5 to 8) applies IP procedure in order to solve a series of decision problems. We recall that the decision problem amounts to find the best density of spheres regarding the (current) specific order of the spheres. Whenever the new solution achieves a better density so far, then the best solution S ? is updated with its subsets N1 (the packed spheres) and N2 (the remaining spheres). The diversification stage (Lines from 10 to 14) applies DRP procedure which combines the destroying operator and the rebuilding strategy. Indeed, the first step removes α% of packed spheres from the current solution S 0 and the second one rebuilds a new diversified solution; that is a solution completing the “partial solution” obtained after destroying the current solution. The complete solution is obtained by using the basic procedure BGP. In the case where the solution is improved (Lines from 12 to 14), i.e., reaching a higher density, the best solution S ? is replaced by Snew with its two subsets. Because the method alternates between these two stages, in this case the diversified solutions undergo the improvement stage (applying IP procedure), i.e., solving a series of decision problems to attempt improving the quality of the solutions. Finally, both stages are repeated until a maximum number of iterations (or fixed runtime) is performed. In order to show how LSBM works, Figure 10 illustrates (from the left-hand to the right-hand) (i) the solution realized when the rebuilding stage is applied, (ii) the fictive positions associated to the remaining (non-positioned) spheres and, (iii) the final solution 14

ACCEPTED MANUSCRIPT

Algorithm 4 : A Local Search-Based Method for the sphere packing problem Input. An instance of packing problem. Output. S ? , the best solution of the problem. 1: Consider that all spheres of N are ordered in decreasing order of their radii. 2: Let S 0 be a starting solution realized by the procedure BGP (cf. Section 3.1). 3: S ? := S 0 // the best solution found so far.

CR IP T

4: while (the stopping criteria is not performed) do 5: Apply IP to S 0 and let Snew be the new reached solution. // Improvement stage 6: if (DSnew > DS ? ) then 7: Snew := S ? and update both subsets N1 and N2 (according to Snew ). 8: end if

PT

ED

M

AN US

9: Snew := S 0 . 10: Apply DRP to S 0 (cf. Section 3.2). // Diversification stage 11: Let Snew be a new solution reached. 12: if (DSnew > DS ? ) then 13: Snew := S ? and update both subsets N1 and N2 (according to S 0 ). 14: end if 15: end while 16: return S ? .

CE

Figure 10: Illustration of LSBMs’ iterations on SYS1KP instance: (i) the rebuilding’s solution (on the left-hand), (ii) the eligible positions generated according to the unpacked spheres (on the middle) and (iii) the final solution achieved when applying the insertion procedure (on the right-hand).

AC

obtained when the insertion procedure IP is applied. As shown from Figure 10, the density of the new solution becomes equal to 52.1885%.

4

Computational results

The purpose of this section is two-fold: (i) to show how to determine a good trade-off between the quality of the obtained solutions when varying the destroying parameter and (ii) to evaluate the effectiveness of the proposed Local Search-Based Method (LSBM) when its achieved results are compared to the best results available in the literature. Moreover, the proposed method was coded in C++, tested on a Pentium 2.4 Mhz (with 2 Gb of RAM) and the other tested methods were also ran on equivalent computers. 15

ACCEPTED MANUSCRIPT

G-Av.

53.489

#d 15 25 24 27 29 35 25.83

Inst. SYS1 SYS2 SYS3 SYS4 SYS5 SYS6

Min DLSBM 52.547 54.172 55.265 54.538 54.322 56.963

G-Av.

54.635

Inst. SYS1 SYS2 SYS3 SYS4 SYS5 SYS6

Min DLSBM 50.476 52.345 54.025 52.294 53.069 54.945

G-Av.

52.859

α = 10 Mean DLSBM 53.2558 55.2099 55.1969 53.2778 54.1031 56.2595 54.551

#d 16.5 24.7 25.7 28.5 29.4 36.2 26.83

α = 30 Mean #d 17 24 27 30 33 36

DLSBM 53.5928 55.8714 56.0701 54.7868 56.1864 57.7119

#d 17.2 24.3 26.8 30.3 33.3 36.4

27.83

55.703

28.05

#d 15 21 24 27 28 35 25.00

α = 50 Mean DLSBM 51.7235 53.8049 54.5342 53.1705 54.0218 56.122 53.896

Max DLSBM 53.901 56.026 56.252 53.631 54.852 56.906

#d 15.2 22.5 24.3 28.8 29.8 35.5

55.261

#d 17 26 25 31 30 35 27.33

Max DLSBM 54.229 56.396 57.191 55.734 57.307 58.205 56.510

#d 18 25 26 31 34 35 28.17

Max DLSBM 52.665 54.508 55.368 54.445 55.711 57.839

26.02

#d 16 25 27 31 32 36

Min DLSBM 52.761 54.353 55.113 53.471 54.664 56.196 54.426

Min DLSBM 53.008 54.099 54.308 53.035 55.081 55.573 54.184

Min DLSBM 50.031 52.557 53.143 52.222 53.437 53.858

#d 17 24 27 31 29 37 27.50

α = 20 Mean DLSBM 53.4298 55.7760 56.1647 53.5515 55.5279 57.0971 55.258

55.089

27.83

52.541

#d 17.1 24.8 26.2 29.5 32 37 27.77

α = 40 Mean #d 17 21 23 30 34 35

DLSBM 53.5395 55.7458 55.8985 54.2323 56.0719 57.3126

#d 17 24.5 25.7 30.1 33 36.8

26.67

55.467

27.85

Max DLSBM 53.901 56.322 56.449 53.611 56.808 57.964 55.843

Max DLSBM 53.653 56.301 56.871 55.856 57.012 58.131

#d 17 24 27 29 33 38 28.00

CR IP T

Min DLSBM 52.623 53.622 53.511 52.678 53.314 55.183

#d 14 22 23 27 29 34

α = 60 Mean DLSBM 50.8633 53.7337 54.1650 52.7863 54.0989 55.5016

AN US

Inst. SYS1 SYS2 SYS3 SYS4 SYS5 SYS6

24.83

53.525

#d 14.4 22.8 24 26.8 29.8 34.9

25.45

56.304

Max DLSBM 51.932 54.316 55.074 53.448 55.601 56.249 54.437

#d 17 25 26 31 34 38

28.50

#d 14 22 24 30 32 36

26.33

Table 3: Behavior of LSBM on instances of 3DKP, when varying the parameter α.

4.1

Behavior of LSBM: 3DKP

AC

CE

PT

ED

M

This part was conducted on two sets of instances related to 3DKP: the first set (noted Set1) contains twelve instances divided into two subsets, where each of them is composed of six instances. The first subset (noted Set1a) contains instances taken from Akeb [1], where the number of spheres to pack varies from 25 to 60 (noted from SYS1KP to SYS6KP). The second subset (noted Set1b) contains instances (noted KBG1, KBG2, KBG3,KBG7, KBG8 and KBG9) taken from Kubach et al. [5], where the number of spheres varies from 30 to 50. These instances proposed and used by Kubach et al. [5] and they were used as benchmarks in Akeb [1] for testing his proposed approach. The second set of instances (Noted Set2) is composed of 96 instances extracted from Kubach et al. [6]. In this case, Set2 was also divided into four subsets (noted set2a, set2b, set2c and set2d) following the number of spheres to pack which varies in the discrete interval {20, 30, 40, 50}. Because the size of instances differs for both sets, i.e., Set1 and Set2, we then fixed the runtime limit, considered as the stopping criterion of the method, to 3600 seconds for instances of Set1 and 7200 seconds for the instances of Set2 (which represent more hardness instances). 4.1.1

Effect of the parameter α

In the preliminary results, two parameters should be taken into account by LSBM (cf., Algorithm 4): the value associated to α and the maximum runtime fixed (which can be considered as a standard runtime limit considered by algorithms of the literature). Herein, LSBM was tested on instances of Set1a by varying the value of α in the discrete interval {10, 20, 30, 40, 50, 60}, where for each instance, ten trials were considered. Over the ten trials, Table 3 exposes for each tested instance the minimum (Min), the average (Mean) and the maximum (Max) density (noted D) reached by LSBM regarding the value associated to α. It also shows the number of spheres packed (noted #d) by LSBM 16

ACCEPTED MANUSCRIPT

4.1.2

Behavior of LSBM on instances of Set1

CR IP T

(minium, average and maximum values, respectively). From Table 3, one can observe that LSBM is able to provide better average solution values for α = 30. We believe that the phenomenon is due to the greedy procedure that uses an intuitive packing which, in some cases, it induces solutions with lesser densities. Indeed, both smallest and greatest values of α may provide a premature convergence of the method and so, the insertion procedure may have a less chance to improve the quality of the solution at hand. Finally, the moderate value of α seems more interesting when combining both the greedy strategy and the insertion one. Hence, α with the value 30 will be adopted for the rest of the paper.

LSBM is first tested on the set Set1 containing 12 instances, represented by two subsets Set1a and Set1b (extracted from Kubach et al. [5] and used as benchmarks in Akeb [1]). All bounds achieved by LSBM are compared to the best bounds reached by TSLA.

n 25 35 40 45 50 60

H 6.9 7.9 6.9 9.9 9.9 9.9

W 5.5 6.5 5.5 8.5 8.5 8.5

L 5 5 5 5 5 5

Search-Based Method Max DLSBM #dLSBM 54.229 18 56.396 25 57.191 26 55.734 31 57.307 34 58.205 35 56.510 28.167

AN US

Inst. SYS1KP SYS2KP SYS3KP SYS4KP SYS5KP SYS6KP G-Av.

Local Mean DLSBM #dLSBM 53.593 17.2 55.871 24.3 56.070 26.8 54.787 30.3 56.186 33.3 57.712 36.4 55.703 28.050

TSLA DT SLA #dT SLA 53.008 15 55.358 25 55.086 23 54.048 29 55.757 32 56.960 36 55.036 26.667

%Improvement Mean Max 1.091 2.252 0.919 1.841 1.755 3.681 1.349 3.025 0.764 2.705 1.303 2.139 1.197 2.607

4.1.2.1

ED

M

Table 4: Performance of LSBM (with α = 30) versus TSLA on instances of Set1a. The value in “boldface” means that the method achieves a better bound and the value in “italic” means that the average bound achieved by the method improves the best solutions of the literature.

Behavior of LSBM on instances of Set1a

AC

CE

PT

Table 4 compares the performance of LSBM to that of TSLA (proposed in Akeb [1]). Columns from 1 to 5 indicate the instance label and its related information. Column 6 (resp. column 7) shows DT SLA (resp. #dT SLA ), the density (resp. number) obtained (resp. packed) by TSLA whereas columns from 8 to 11 show the same results achieved by LSBM, i.e., DLSBM and #dLSBM , respectively. Finally, column 12 (resp. column 13) tallies the maximum (resp. average) percentage improvement (when it happens) yielded by LSBM when compared its results to those reached by TSLA (the percentage improvement −DT SLA is computed as follows: DLSBM × 100). DLSBM The analysis of Table 4 shows what follows: • LSBM outperforms TSLA whenever the runtime limit was fixed to 3600 seconds. Indeed, it is able to improve all the bounds for all instances of Set1a. • When comparing the maximum global average density (cf. the global average value noted G-Av., column 10), one can observe that LSBM realizes 56.510%, which dominates that achieved by TSLA (55.036%). One can also observe that packing more spheres into the container does not automatically lead to a better solution. Indeed, such a phenomenon can be observed on the the instance SYS6KP, where TSLA packed

17

ACCEPTED MANUSCRIPT

36 spheres for realizing a density of 56.960% whereas LSBM positioned 35 spheres (realizing density of 58.205%). • When comparing the average values achieved by LSBM (55.703%), one can also observe that LSBM’s G-Av. (column 8) remains superior than that realized by TSLA (55.036%). As discussed above, for the same instance SYS6KP, LSBM’s average density (57.712%, column 8) remains significantly greater than that of TSLA (56.960%, column 6) and it is able to pack an average of 36.4 spheres.

M

AN US

CR IP T

• Furthermore, as shown from column 13, the maximum percentage improvement (Max) is very impressive for such a problem, since it varies from 1.841% (instance SYS2KP) to 3.681% (instance SYS3KP). Note also that even if one consider the average values related to the ten trials (cf., column 12, Mean), we can see the superiority of LSBM since it is able to realize an average percentage of improvement varying from 0.764 (instance SYS5KP) to 1.755 (instance SYS3KP) with a global average value equals to 1.197.

(b) SYS5KP

(c) SYS6KP

ED

(a) SYS4KP

PT

Figure 11: Illustration of the solutions’ structures of three instances (SYS4KP on the lefthand, SYS5KP on the middle and SYS6KP on the right-hand) of the set Set1a.

CE

Finally, Figure 11 illustrates three final configurations related to the instances SYS4KP (on the left-hand of the figure), SYS5KP (on the middle of the figure) and SYS6KP (on the right-hand of the figure). These structures represent new solutions, where DLSBM (SYS4KP) = 55.856%, DLSBM (SYS5KP) = 57.307% and DLSBM (SYS6KP) = 58.205%. Behavior of LSBM on instances of Set1b

AC

4.1.2.2

Herein, the performance of LSBM is evaluated on instances of the second part of Set1 (noted Set1b) that contains 6 instances (taken from Kubach et al. [5]). Its performance is compared to those realized by the three heuristics: (i) sequential version of B1.16 heuristic with the parameter τ = 0.8 (cf. Kubach et al. [5]), (ii) sequential version of B1.16 heuristic with the parameter τ = 1.0 (cf. Kubach et al. [5]) and, (iii) TSLA (cf. Akeb [1]). We also evaluated the ratios between LSBM and the three considered heuristics, i.e., both sequential versions of B1.16 and TSLA, respectively. First, Table 5 shows the results achieved by both versions of B1.16, TSLA and those reached by the proposed LSBM. Columns from 1 to 5 show the instance’s information. The 18

ACCEPTED MANUSCRIPT

Inst. KBG1 KBG2 KBG3 KBG7 KBG8 KBG9 G-Av.

n 30 30 30 50 50 50

L 11.103 1.99 17.785 13.71 2.207 27.965

H 10 10 10 10 10 10

W 10 10 10 10 10 10

B1.16 Dv1 Dv2 53.657 55.001 30.071 30.071 52.154 52.375 55 55 46.406 46.517 52.54 53.173 48.305 48.690

TSLA DT SLA #dT SLA 55.001 30 30.071 30 52.625 28 55 50 46.342 41 53.662 48 48.784 37.83

LSBM Mean Max DLSBM #dLSBM DLSBM #dLSBM 54.663 29.6 55.001 30 30.071 30 30.071 30 52.220 27.5 53.507 29 55 50 55 50 46.184 35.8 47.416 41 52.968 46.5 53.662 48 48.518 36.57 49.110 38

CR IP T

Table 5: Performance of LSBM versus two versions of B1.16 (with τ = 0.8 and τ = 1.0) and TSLA on instances of Set1b

AN US

results of both versions of B1.16 are displayed in column 6 for the first version (noted Dv1 ) and in column 7 for the second version (noted Dv2 ), respectively. Column 8 tallies those achieved by TSLA (represented by the density DT SLA ) and column 9 shows the number of spheres packed by the same heuristic (noted #dT SLA ). Finally, columns 10 and 11 (resp. columns 12 and 13) display the average densities (resp. the average number of packed spheres) achieved by LSBM for the ten trials. Second and last, in order to show the behavior of LSBM versus both versions of B1.16 and TSLA, we also evaluated the percentage ratios between LSBM and B1.16 and, LSBM −B1.16 −DT SLA and TSLA; both ratios represent DLSBM × 100 and DLSBM × 100, respectively. DLSBM DLSBM Table 6 reports the percentage ratios relating the quality of the solutions achieved by LSBM (column 1 for LSBM versus the first version of B1.16, column 2 for LSBM versus the second version of B1.16 and, column 3 for LSBM versus TSLA). The analysis of the results of both Tables 5 and 6 follows:

LSBM versus B1.16 (τ = 0.8) Mean Max 1.841 2.444 0.000 0.000 0.126 2.529 0.000 0.000 -0.482 2.130 0.808 2.091 0.382 1.532

LSBM versus B1.16 (τ = 1.0) Mean Max -0.618 0.000 0.000 0.000 -0.297 2.116 0.000 0.000 -0.722 1.896 -0.387 0.911 -0.338 0.821

LSBM versus TSLA Mean Max -0.618 0.000 0.000 0.000 -0.776 1.648 0.000 0.000 -0.343 2.265 -1.311 0.000 -0.508 0.652

AC

CE

Inst. KBG1 KBG2 KBG3 KBG7 KBG8 KBG9 G-Av.

PT

ED

M

• First, LSBM outperforms the sequential B1.16 heuristic (with τ = 0.8). Indeed, LSBM realizes a percentage average value of 49.110% (column 12 and the last line of the table noted G-Av.) compared to the average value of 48.3047% (column 12 and the last line of the table noted G-Av.) realized by the sequential heuristic B1.16. Indeed, LSBM dominates B1.16 heuristic (with τ = 0.8) in four cases and matches the solutions for the rest of the instances. In this case, the percentage improvement (%Imp.) varies from 2.091% (KPG9) to 2.529% (KPG3) and both LSBM and B1.16 match the optimal densities for two instances (KBG2 and KBG7).

Table 6: Behavior of the three compared methods on instances of Set1.b

• Second, LSBM has a better behavior that the second version of B1.16 (with τ =1). Indeed, it improves three solutions (instances) among the sixth ones, where LSBM’s percentage improvement varies from 0.9113% (instance KP9) to 2.11561% (instance KPG3). Moreover, both LSBM and B1.16 match the optimal densities for the rest of instances, i.e., KBG1, KBG2 and KBG7. • Third, LSBM performs better than TSLA for the instances of Set1b. Indeed, LSBM provides a better result for two instances (KPG3 and KBG8) and it matches the rest 19

ACCEPTED MANUSCRIPT

of the bounds achieved by TSLA. By excluding the solutions matched by LSBM, the percentage improvement varies from 1.64838% (KPG3) and 2.265% (KPG8). • Furthermore, over the six instances, LSBM is able to reach two new solutions when compared to the best published bounds and matches the rest of the solutions (4). Local Search-Based Method’s cpu min mean max 50.67 589.505 54.42 0 0 0 3.87 522.252 1934.46 0 0 0 57.01 1322.926 113.75 128.34 2176.186 3476.85 39.98 685.14 929.91

CR IP T

Inst. KBG1 KBG2 KBG3 KBG7 KBG8 KBG9 G-Av.

Table 7: Variation of LSBM’s runtime on instances of Set1b

AC

CE

PT

ED

M

AN US

Moreover, as we have already mentioned that the runtime limit was fixed to one hour for all compared methods. For instances of Set1b, we also saved the best (average) runtime for which LSBM provides the best (average) bound, as shown in Table 7. Indeed, over the ten trials, and by excluding both KBG2 and KBG7 (for which the optimality is proven), one can observe that the minimum and maximum global average runtimes (last line, column 2 and column 4, respectively) is equal to 39.98 seconds and 929.91 seconds, respectively. Finally, the global average runtime over the ten trials is equal to 685.14 seconds, which is much less important than the one-hour considered as one of the stopping criteria. We also provide the structure related to two instances of Set1b as illustrated in Figure 12 (KBG3 and KBG8 from the left-hand to the right-hand), where LSBM is able to improve their bounds (we recall that for the rest of instances, all methods are able to provide the optimal solutions, where all spheres are packed). These structures represent new solutions, where DLSBM (KBG3) = 53.507% and DLSBM (KBG8) = 47.416%.

(a) KBG3

(b) KBG8

Figure 12: Illustration of the solutions’ structures of two instances ((a) KBG3 and (b) KBG8) belonging to Set1b.

20

ACCEPTED MANUSCRIPT

4.1.3

Behavior of LSBM on instances of Set2

Df 49.517 50.901 51.148 51.605 50.79275

(τ = 0.8) Dtot ttot 41.798 13 46.519 458 49.439 8259 51.020 7647 47.194 4094.25 (τ = 1.0) Dtot ttot 44.067 1566 48.941 35528 51.565 222400 52.950 203163 49.381 115664.25

AN US

Subset Set2a Set2b Set2c Set2d G-Av.

Parallel B1.16 tf tsucc 16 4 550 0.01 9438 0.01 8273 759 4569.25 190.755 Parallel B1.16 tf tsucc 2208 6 53301 10646 469900 12978 405077 1250 232621.5 6220

M

Df 46.368 49.224 50.553 51.042 49.29675

ED

Subset Set2a Set2b Set2c Set2d G-Av.

CR IP T

In order to study the behavior of the proposed method on more hardness instances, we evaluated its performance on a set (noted Set2) containing 96 instances taken from Kubach et al. [6]. We then compared its obtained results to those reached by Kubach et al.’s parallel heuristic (noted B1.16). Set2a includes the instances noted from KP1 to KP24, Set2b is represented by the instances noted from KP74 to KP96, Set2c includes the instances noted from KP145 to KP240 and, Set2d contains the rest of the instances noted from KP217 to KP240, respectively. Moreover, Kubach et al.’s [6] parallel heuristic has been run without runtime limit whereas LSBM’s runtime limit was fixed to 7200 seconds, where the last bounds achieved by LSBM was considered as the final solution value (density) for the method. Because Kubach et al.’s [6] parallel heuristic uses a parameter τ , which serves to reduce the number of positions to attempt, we then considered both versions: the first one runs with the parameter τ = 0.8 and the second one that considers all positions, i.e., with the parameter τ = 1.0. nbopt 7 4 3 2 4 nbopt 7 10 13 12 10.5

Maximum and average results on ten trials: Local Search-Based Method Df Mean 48.111 49.485 49.762 46.117 48.369

tf

Max 49.632 50.922 51.234 51.615 50.851

Mean 250.75 1330.86 1838.34 3603.63 1755.89

PT

Subset Set2a Set2b Set2c Set2d G-Av.

Max 638.30 3062.74 2542.66 5176.80 2855.12

tsucc Mean Max 15.86 48.09 274.27 595.23 398.34 767.27 159.14 197.46 211.91 402.01

— Mean 42.949 47.917 50.853 48.021 47.435

Dtot Max 44.111 48.954 51.604 52.955 49.406

ttot Mean Max 182.24 466.15 890.61 2034.61 1058.34 1580.99 1881.39 2687.13 1003.15 1692.22

nbopt 7 10 13 12 10.50

CE

Table 8: Behavior of both versions of B1.16 and LSBM on the instances of Set2.

AC

Table 8 shows the results realized by both versions of B1.16 (displayed in the first twoblocks of the table) and those obtained by the proposed method LSBM (the third block of the table). First, for all blocks, column 1 displays the data labels of the four subsets of Set2. Second, Df denotes the average density of each subset of Set2 (Set2a, . . ., Set2d) whenever the method is not able to pack all spheres and tf is the average final runtime that needs the method to achieve the final solutions (Df ). Third, tsucc is the average final runtime that needs each method to pack all spheres (with success), Dtot (resp. ttot ) represents the average density (resp. its corresponding average runtime) for all instances of each subset. Finally, nbopt is the number of instances solved to optimality by the corresponding method. In addition, because LSBM considers ten trials, we then display (in the third block) the average (noted Mean) and the maximum (noted Max) values of Df , tf , tsucc , Dtot , ttot and nbopt , respectively. From Table 8, we can remark what follows. 21

ACCEPTED MANUSCRIPT

• LSBM outperforms both parallel versions of B1.16. Indeed, LSBM’s global average density is equal to 49.406% whereas the first (resp. second) version of the parallel B1.16 is equal to 47.194% (resp. 49.381%).

CR IP T

• LSBM has a better behavior when comparing its obtained results to those reached by the first version of parallel B1.16. Indeed, the percentage gap improvement of LSBM (LSBM−B1.16)×100 (between both average densities) varies from 0.5% (instances of Set2d) to 4.4% (instances of Set2a). The percentage related to filling spheres into the container increases according to the number of spheres to pack, even if the gap related to the improvement doesn’t have the same behavior (for example, see the average densities for the instances of Set2b and Set2c). • LSBM’s results remain better than those reached by the second version of parallel B1.16. Indeed, its gap remains positive, where its average density varies from 41.798% to 51.02% whereas LSBM’s average density increases for realizing the value of 44.111% to 52.955%, respectively.

AN US

• Even if the average runtime of the first version of B1.16 remains interesting (4094.25 seconds), that of LSBM is very impressive (1884.61 seconds). Indeed, both average runtimes show that LSBM is 2.17 times faster than the first parallel version of B1.16.

ED

M

• As pointed above, the second parallel version of B1.16 has a great chance to find solutions with high quality, because all eligible positions of each candidate sphere to pack are considered. Despite its unlimited runtime, LSBM is able to achieve the same number of optimal solutions (at the same time it improves the global average density of the instances of Set2, as discussed above). In fact, LSBM’s average runtime is now 61.37 times faster than that used by the second parallel version of B1.16 for achieving the average densities exposed in the second block of the table (even if the processors used by each method are slightly different).

AC

CE

PT

Second and last, Tables from 9 to 12 report the average results realized by LSBM for the four subsets of instances (Set2a, Set2b, Set2c and Set2d, respectively) over the ten trials. Column 1 (of each table) displays the label of the instance and reports its data information in the three next columns. Column 5 (resp. column 6) reports LSBM’s average density (resp. average runtime) for the ten trials whereas column 7 (resp. column 8) shows the best density (resp. runtime) achieved (resp. needed) by LSBM over the ten trials. From Tables 9 to Table 12, one can observe that LSBM confirms its superiority by achieving several optimal solutions (as shown in column 10 of Table 8: an average of 10.5 optimal solutions realized by LSBM, for which all spheres were packed) almost of the four ones realized by the first parallel version of B1.16 (cf. Tables from 9 to 12, all instances marked with symbol “?”). On the other hand, even LSBM has a better behavior than that the second parallel version of B1.16, it realizes the same number of optimal solutions (marked with the symbol “?”).

4.2

Behavior of LSBM: 3DODP

We recall that both 3DKP and 3DODP are studied. Hence, in order to compare the results achieved by Stoyan et al.’s [27] algorithm, we adapted the proposed method for solving 3DODP. Indeed, a series of LSBM are applied for a predefined interval search [L, L], where L denotes a lower bound and L its upper bound. For this adaptation, the used instances are 22

ACCEPTED MANUSCRIPT

H 30.345 20.23 15.172 10.416 6.944 5.208 13.146 8.764 6.573 7.171 4.781 3.586 1.97 1.97 1.97 11.916 7.944 5.958 5.136 3.424 3.29 1.956 1.956 1.956

W 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10

D 46.4265 46.967 50.4692 44.1237 45.4693 45.5844 46.1502 49.6914 50.438 48.4374 47.1886 47.4227 43.8306 30.573 22.93 48.4367 52.4241 48.7011 51.8485 48.1165 42.934 33.511 22.34 16.755 42.9487

ts 33.53 387.70 1269.20 185.78 184.18 588.84 623.24 829.53 902.98 131.60 345.56 611.21 336.49 0.00 0.00 37.87 1187.47 209.76 3251.79 70.80 0.10 0.01 0.02 0.03 466.15

CR IP T

L 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10

Local Search-Based Method Average Max #d ts D #d 15.6 64.89 47.670 17 17 107.50 47.687 17 18 940.59 51.509 18 15.9 202.81 47.336 17 15.7 76.74 49.249 18 16.4 282.51 47.382 17 16.6 341.61 47.597 17 17.5 188.56 51.630 19 18 449.03 51.175 18 18.2 115.91 50.179 19 17.6 105.62 48.384 18 16.9 209.15 48.684 17 18.6 110.88 45.860 20 20 0.00 30.573 20 20 0.00 22.930 20 17.9 70.94 49.260 18 18.5 218.69 53.994 19 18.8 126.03 49.916 19 18.6 635.36 53.089 19 17.7 126.75 49.011 18 20 0.10 42.934 20 20 0.01 33.511 20 20 0.02 22.340 20 20 0.03 16.755 20 18.06 182.24 44.1106 18.54

AN US

Set2a (n = 20) Inst KP1 KP2 KP3 KP4 KP5 KP6 KP7 KP8 KP9 KP10 KP11 KP12 KP13? KP14? KP15? KP16 KP17 KP18 KP19 KP20 KP21? KP22? KP23? KP24? G-Av.

Table 9: Average and best densities (with their runtimes) reached by LSBM over the ten trials for the instances of Set2a

CE AC

Local Search-Based Method Average Max #d ts D #d 24.3 387.82 48.112 26 26.7 652.83 51.957 28 27.9 1954.72 53.476 28 24.4 1858.68 48.948 28 25.9 682.22 52.037 26 25.9 1003.15 51.693 28 27.1 1664.00 49.58 28 29 1596.04 55.001 30 35.3 3461.10 54.0093 29 28.8 898.44 53.69 29 28.6 306.21 55.001 30 27.3 1381.51 52.928 28 21.3 698.28 44.267 23 30 0.00 40.094 30 30 1.49 30.071 30 26.6 1758.30 49.608 28 29.7 43.38 55.001 30 30 17.70 55.001 30 29.3 363.60 54.997 30 30 414.00 55.003 30 27.9 1582.21 52.981 29 25.4 648.73 49.624 27 30 0.00 35.322 30 30 0.32 26.492 30 27.98 890.61 48.9539 28.54

M

H 35.57 23.714 17.785 14.506 9.671 7.253 25.011 16.674 12.506 9.167 6.111 4.583 2.176 1.99 1.99 22.206 14.804 11.103 9.019 6.012 4.509 1.926 1.926 1.926

W 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10

D 47.3891 50.9428 52.595 47.0779 50.0356 49.1036 48.23256 52.922 52.94614 53.042 53.651 49.3988 43.6295 40.094 30.071 48.3313 54.6358 55.001 54.0305 55.003 50.9798 49.0881 35.322 26.492 47.9173

ED

L 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10

PT

Set2b (n = 30) Inst KP73 KP74 KP75 KP76 KP77 KP78 KP79 KP80? KP81 KP82 KP83? KP84 KP85 KP86? KP87? KP88 KP89? KP90? KP91? KP92? KP93 KP94 KP95? KP96? Av.

ts 605.28 250.51 4413.77 5421.99 1818.54 2026.08 3470.98 3902.78 7110.85 344.71 902.78 5398.43 1023.87 0.00 1.49 4950.59 50.38 17.70 662.84 414.00 5154.99 887.83 0.00 0.32 2034.61

Table 10: Average and best densities (with their runtimes) reached by LSBM over the ten trials for the instances of Set2b extracted from Stoyan et al. [27] and for instances extracted from another paper published

23

ACCEPTED MANUSCRIPT

H 44.958 29.972 22.479 20.267 13.511 10.134 27.612 18.408 13.806 12.816 8.544 6.408 3.316 2.211 1.968 25.748 17.165 12.874 9.655 6.436 4.827 2.433 1.986 1.986

W 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10

D 47.2917 52.7601 53.821 45.3572 51.2145 50.6889 49.2829 54.999 54.999 53.9366 55.001 55.001 52.05369 45.5288 46.333 50.6741 54.831 54.999 54.998 55.004 54.3789 48.7117 44.914 33.686 50.8527

ts 2692.12 4213.43 3397.87 3397.87 1664.17 422.23 3697.03 0.10 19.12 5415.09 13.80 84.28 1264.11 1087.51 280.84 1549.36 1140.20 29.74 9.16 23.01 2958.68 4583.57 0.00 0.55 1580.99

CR IP T

L 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10

Local Search-Based Method Average Max #d ts D #d 31.3 1346.11 47.974 32 34.1 3121.53 53.496 37 35.5 2379.24 54.439 39 29.8 1490.51 46.597 39 35.9 1389.44 53.989 37 32.1 1418.10 53.805 37 35.7 2666.36 49.874 36 40 0.10 54.999 40 40 19.12 54.999 40 37.3 2038.65 55.001 40 40 13.80 55.001 40 40 84.28 55.001 40 35.9 1487.85 53.7899 38 29.7 1092.41 46.752 31 40 280.84 46.333 40 36.5 1834.21 53.001 38 39.7 928.99 55.00 40 40 29.74 54.999 40 40 9.16 54.998 40 40 23.01 55.004 40 39.6 1750.22 55.00 40 35 1996.03 49.854 37 40 0.00 44.914 40 40 0.55 33.686 40 37.00 1058.34 51.6044 38.38

AN US

Set2c (n = 40) Inst KP145 KP146 KP147 KP148 KP149 KP150 KP151 KP152? KP153? KP154? KP155? KP156? KP157 KP158? KP159? KP160? KP161? KP162? KP163? KP164? KP165? KP166 KP167? KP168? Av.

Table 11: Average and best densities (with their runtimes) reached by LSBM over the ten trials for the instances of Set2c

CE AC

Local Search-Based Method Average Max #d ts D #d 40.7 3004.12 47.887 41 48.3 3932.58 53.252 48 48.8 4770.90 54.358 49 43.9 3586.39 48.993 45 47.6 5024.99 53.223 48 46.9 4387.99 52.705 48 46.6 3200.08 51.668 47 50 0.00 55.001 50 50 27.78 55.000 50 50 46.05 55.002 50 50 42.02 55.000 50 49.6 555.16 54.998 50 48.7 3168.16 54.417 49 47 3249.02 51.866 48 41.2 2126.06 47.38 42 48.1 2229.54 52.5 49 50 0.00 55.000 50 50 0.00 55.000 50 50 63.41 55.001 50 50 0.00 55.001 50 50 12.20 54.996 50 49.5 1163.09 55.004 50 45.6 4563.74 51.131 47 50 0.00 46.544 50 48.02 1881.39 52.9553 48.38

M

H 55.93 37.287 27.965 23.291 15.527 11.645 34.766 23.177 17.383 15.043 10.029 7.522 4.414 2.943 2.207 27.42 18.28 13.71 11.337 7.558 5.669 3.334 2.223 1.97

W 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10 5 7.5 10

D 47.3276 52.6881 53.5656 47.9998 52.6161 52.2282 50.2251 55.001 55 55.002 55 54.722 53.6287 51.1253 47.1602 51.9017 55 55 55.001 55.001 54.996 54.7558 49.7962 46.544 52.5536

ED

L 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10

PT

Set2d (n = 50) Inst KP217 KP218 KP219 KP220 KP221 KP222 KP223 KP224? KP225? KP226? KP227? KP228? KP229 KP230 KP231 KP232 KP233? KP234? KP235? KP236? KP237? KP238? KP239 KP240? Av˙

ts 3260.51 5432.78 6589.09 6973.36 6303.38 6997.95 2754.82 0.00 27.78 46.05 42.02 660.27 4899.55 7126.74 1851.69 3404.48 0.00 0.00 63.41 0.00 12.20 1517.80 6527.19 0.00 2687.13

Table 12: Average and best densities (with their runtimes) reached by LSBM over the ten trials for the instances of Set2d by the same authors (Scheithauer et al. [21]). These instances are composed of two sets:

24

ACCEPTED MANUSCRIPT

the first set (noted S1) is composed of six instances, where the number of spheres to pack varies from 25 to 60 (noted from SYS1 to SYS6) and , the second set (noted S2) contains 119 instances, where the number of spheres varies from 20 to 50. Note also that SYS and LSBM perform with a precision  of 10−6 (1 ). 4.2.1

Adaptation of LSBM: set S1

n 25 35 40 45 50 60

KBTG 9.2656 8.9301 8.7178 10.4042 10.9865 11.8399 10.02402

SYS 9.1604 8.8195 8.4108 9.9288 10.5905 11.5708 9.74680

Local Search-Based Method Min Mean Max 9.097764 9.187306 9.250169 8.805289 8.820709 8.838000 8.401346 8.431144 8.500924 10.083717 10.096280 10.109824 10.510374 10.808179 10.968444 11.536260 11.690582 11.899249 9.739125 9.839034 9.927768

AN US

Inst. SYS1 SYS2 SYS3 SYS4 SYS5 SYS6 Av.

CR IP T

Table 13 shows the results obtained by all considered heuristics: KBTG (Kubach et al. [5])), SYS (Stoyan et al. [21, 27]) and the adaptation of LSBM for 3DODP. Column 3 (resp. column 4) of Table 13 tallies KBTGs’ (resp. SYSs’) upper bounds and the last three columns display the best (Min), average (Mean) and wort (Max) results achieved by LSBM. Finally, the last line of the table reports the average solution values achieved by the three methods on all instances of the set S1.

Table 13: Behavior of LSBM and two heuristics of the literature on instances of S1.

ED

M

Note also that the size of each instance varies for instances of set S1, we then tried to consider the equivalent runtimes, as used by both tested methods available in the literature (KBTG and SYS). In this case, LSBM’s runtime limit was fixed to one hour for the instance SYS1, six hours for SYS2 and eight hours for instances from SYS3 to SYS6. We note that the value in “bold-face” means that the best solution value has been obtained by the considered algorithm. The analysis of the results of Table 13 follows.

CE

PT

• On the one hand, one can observa that LSBM outperforms KBTG since it is able to improve all the upper bounds reached by KBTG. In this case, LSBM’s global average value is equal to 9.739125 whereas KBTG’s average global value is equal to 10.02402.

AC

• On the other hand, LSBM remains competitive when comparing its achieved results to those reached by SYS. In this case, LSBM’s global average solution value (9.739125) remains slightly better than that SYS’s average global value (9.74680). One can also observe that LSBM is able to provide five new upper bounds and fails in one occasion to match the best result (values marked in bold-space).

4.2.2

Adaptation of LSBM: set S2

Herein, LSBM’s performance is analyzed on the second set of instances S2. Its obtained results are also compared to those published by Kubach et al. [5]) (KBTG) and Stoyan et al. [21, 27] (SYS). 1

The right precision has been provided by Prof. Y. Stoyan and Prof. G. Yaskov that has been used in Stoyan et al.’s [27] method.

25

ACCEPTED MANUSCRIPT

Set S2 Inst. KBG 1-37 KBG 73-102 KBG 145-177 KBG 217-247 Average

n 20 30 40 50

KBTG 11.1921 14.3243 18.3098 20.4260 16.0630

SYS 11.0087 14.0157 17.8771 19.8412 15.6857

Local Search-Based Method Min Mean Max 10.975066 11.009953 11.045211 13.976988 14.036885 14.092108 17.894299 17.991793 18.132221 19.825432 19.975055 20.072401 15.667946 15.753422 15.835485

%Imp KBTG SYS -2.4445 -0.2687 -2.9824 -0.2345 -2.9344 0.1929 -4.007 -0.0306 -3.0920 -0.0853

CR IP T

Table 14: Performance of LSBM versus KBTG and SYS on instances of S2.

AN US

The results realized by the three methods KBTG, SYS and LSBM are summarized in Table 14 and all detailed results are reported in Tables from 15 to 18. For each table (from Table 15 to Table 18), the first two columns report the instance informations. Column 3 (resp. column 4) displays KBTGs’ (resp. SYSs’) upper bounds taken from Kubach et al. [5] (resp. Stoyan et al. [27]) whereas the last three columns display LSBMs’ upper bounds. Finally, the last line (noted Av.) of each table represents the average solution values realized by the three methods, which were summarized in Table 14. From Table 14 (and Tables from 15 to 18), we can observe what follows: • LSBM versus KBTG: On the one hand, LSBM outperforms KBTG, where it improves all the upper bounds. Indeed, over all tested instances of S2 (cf. Table 14), LSBM’s average global solution value is equal to 15.667946 whereas KBTG realizes a smallest global value of 16.0630. On the other hand, the aforementioned values induced the improvements varying from 2.4445% (for n = 20) to 4.007% (for n = 50).

Conclusion

CE

5

PT

ED

M

• LSBM versus SYS: First, LSBM remains competitive when comparing its achieved results to those reached by SYS. Indeed, LSBM is able to provide 70 out of 119 best upper bounds (representing a percentage of 58.82% of improvement) and fails for the rest of the bounds (representing a percentage of 41.18%). Second and last, even if the approach is not tailored to solving 3DODP, one can observe its superiority in therm of the global average quality (upper bound): 15.667946 for LSBM compared to the SYS’s global average value of 15.6857 (it represents a percentage improvement of 0.0853).

AC

In this paper the three-dimensional sphere packing problem is solved by using a Local Search-Based method that combines three main features: a best-local position procedure, an intensification stage and a diversification stage. The first procedure ensures a starting feasible solution using a basic greedy local strategy. The second stage tries to solve a series of decision problems in order to pack a subset of complementary spheres. The third stage tries to modify the current solution by removing some packed items and replacing them with other spheres. The performance of the proposed method is evaluated on benchmark instances taken from the literature where two versions of the problem were discussed: the three-dimensional knapsack problem and the three-dimensional open dimension problem. For both versions of the problem, the results achieved by the Local Search-Based method were compared to those reached by some methods of the literature. For both versions of the problem, the proposed method remains competitive and succeeded in yielding better 26

ACCEPTED MANUSCRIPT

n 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20

KBTG 35.4136 23.9017 16.2248 12.5449 7.7428 5.9228 15.5238 9.3530 7.0103 7.8308 5.1907 4.2992 13.2412 8.1527 6.3768 5.3784 4.0733 21.5230 13.0966 9.9273 12.1398 7.5517 5.8036 26.2532 17.1909 11.9746 5.2754 3.6257 2.0283 11.1921

SYS 35.0671 23.3937 15.9003 12.3128 7.5094 5.6193 15.2504 9.1665 6.8724 7.7385 5.1217 4.1819 13.2261 8.1040 6.2814 5.3207 3.9641 21.3660 12.9842 9.6635 11.7393 7.3368 5.4127 26.0823 17.0373 11.9082 5.2015 3.6220 1.8681 11.0087

Local Search-Based Method Min Mean Max 35.054555 35.151266 35.216652 23.213585 23.264957 23.335176 15.818059 15.860625 15.900415 12.310173 12.342203 12.399494 7.503538 7.505539 7.510306 5.601681 5.608268 5.616826 15.218954 15.244119 15.263363 9.132202 9.155516 9.177137 6.859323 6.868976 6.889449 7.711843 7.722883 7.734400 5.108434 5.124241 5.134972 4.176256 4.198564 4.207187 13.194877 13.219772 13.245203 8.075537 8.146398 8.261675 6.281391 6.283049 6.283929 5.302099 5.320071 5.324549 3.948935 3.980153 4.045569 21.356880 21.418580 21.496157 12.965426 12.990690 13.020762 9.625079 9.674156 9.681318 11.664111 11.691430 11.744234 7.318732 7.328996 7.330802 5.410875 5.436885 5.441172 25.884028 26.030458 26.132173 16.991775 17.083319 17.174265 11.868965 11.913704 11.961869 5.201388 5.216429 5.232102 3.610196 3.636732 3.672938 1.868004 1.870662 1.877032 10.975066 11.009953 11.045211

%Imp KBTG SYS -1.0242 -0.0358 -2.9643 -0.7759 -2.5714 -0.5199 -1.9068 -0.0213 -3.1887 -0.0781 -5.7325 -0.3145 -2.0031 -0.2066 -2.4178 -0.3756 -2.2010 -0.1906 -1.5425 -0.3457 -1.6104 -0.2597 -2.9439 -0.1351 -0.3511 -0.2366 -0.9555 -0.3525 -1.5189 -0.0001 -1.4391 -0.3508 -3.1493 -0.3840 -0.7778 -0.0427 -1.0117 -0.1448 -3.1399 -0.3992 -4.0782 -0.6446 -3.1832 -0.2469 -7.2581 -0.0337 -1.4263 -0.7660 -1.1719 -0.2679 -0.8900 -0.3306 -1.4229 -0.0022 -0.4294 -0.3270 -8.5812 -0.0052 -2.4445 -0.2687

CR IP T

1 2 3 4 5 6 7 8 9 10 11 12 16 17 18 19 20 25 26 27 28 29 30 31 32 33 34 35 37

AN US

Inst. KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG Av.

Table 15: Performance of LSBM versus the two best heuristics of the literature on instances containing 20 spheres. Local Search-Based Method Min Mean Max 41.099846 41.263541 41.340829 24.632250 24.676425 24.772147 17.773808 17.895364 18.009658 15.990759 16.095403 16.210845 9.837248 9.875865 9.911351 7.368337 7.403275 7.434587 27.980609 28.104894 28.189435 16.689074 16.700316 16.739126 12.483584 12.585667 12.621330 9.099775 9.130884 9.171011 5.984890 6.013917 6.075362 4.733263 4.746001 4.764258 2.510286 2.515243 2.520839 24.682793 24.749500 25.059601 14.456219 14.496628 14.527682 10.624330 10.642847 10.655911 8.818091 8.831826 8.856431 5.906831 5.922563 5.939873 4.639138 4.644469 4.649153 2.099234 2.106508 2.111775 36.307555 36.503263 36.547256 20.993734 21.177488 21.267170 15.580379 15.671204 15.762138 18.779297 18.857900 18.916496 11.106067 11.129604 11.145546 8.398203 8.404652 8.419203 15.491773 15.552465 15.617211 9.518986 9.603429 9.651880 7.227620 7.298768 7.353634 8.495651 8.506628 8.521488 13.976988 14.036885 14.092108

M

SYS 40.8194 24.7246 17.7873 15.9578 9.8169 7.3687 28.2552 16.7146 12.6366 9.1451 5.9927 4.7160 2.4988 25.0700 14.4477 10.6820 8.7404 5.8641 4.6309 2.1068 36.3815 21.0630 15.6891 18.6497 11.0527 8.3332 15.7597 9.7371 7.3754 8.4554 14.0157

ED

KBTG 41.2473 25.2885 18.7148 16.3290 10.3207 7.7341 28.3554 16.9853 12.8157 9.4653 6.1404 4.8995 2.5740 25.0866 14.7129 11.0214 9.0678 6.0593 4.7911 2.1420 36.3819 21.8134 16.0896 19.3571 11.6704 8.8849 15.6612 9.8523 7.4751 8.7927 14.3243

PT

73 74 75 76 77 78 79 80 81 82 83 84 85 88 89 90 91 92 93 94 97 98 99 100 101 102 103 104 105 106

n 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30

%Imp KBTG SYS -0.3588 0.6824 -2.6642 -0.3749 -5.2943 -0.0759 -2.1152 0.2061 -4.9145 0.2068 -4.9640 -0.0049 -1.3395 -0.9814 -1.7750 -0.1530 -2.6604 -1.2257 -4.0169 -0.4981 -2.5984 -0.1305 -3.5121 0.3647 -2.5381 0.4575 -1.6360 -1.5687 -1.7756 0.0589 -3.7374 -0.5428 -2.8318 0.8810 -2.5812 0.7234 -3.2756 0.1776 -2.0372 -0.3604 -0.2048 -0.2037 -3.9043 -0.3299 -3.2683 -0.6978 -3.0768 0.6901 -5.0813 0.4805 -5.7953 0.7740 -1.0937 -1.7295 -3.5016 -2.2914 -3.4241 -2.0447 -3.4965 0.4738 -2.9824 -0.2345

AC

CE

Inst. KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG Av.

Table 16: Performance of LSBM versus the two best heuristics of the literature on instances containing 30 spheres.

bounds for several instances. Acknowledgements. Many tanks to anonymous referees for their helpful comments and pertinent suggestions that contributed to the improvement of the presentation of the paper.

27

ACCEPTED MANUSCRIPT

n 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40

KBTG 52.323 31.2374 23.1462 23.186 14.221 10.662 30.7066 17.9416 13.5376 12.9223 8.3402 6.4583 3.413 2.6084 28.2945 16.9952 12.6447 9.4929 6.4625 4.8768 2.6337 55.2813 33.8326 24.9422 19.4988 12.3795 9.3018 30.3761 18.1909 13.3858 18.30976333

SYS 52.1954 30.3641 21.8031 22.5236 13.5161 10.0959 30.5291 17.7204 13.056 12.7793 8.009 6.2224 3.2183 2.5372 28.1614 16.598 12.2718 9.2759 6.1363 4.7526 2.5698 54.9527 32.6877 24.0768 18.6295 11.6346 8.8331 30.3038 17.7822 13.0775 17.87712

Local Search-Based Method Min Mean Max 52.204195 53.022107 53.806043 30.372953 30.416710 30.446119 21.758152 21.775451 21.804886 22.524545 22.574609 22.601163 13.518773 13.549312 13.635661 10.096200 10.182934 10.280862 30.645055 30.757367 30.810998 17.780510 17.833123 17.890020 13.061335 13.069234 13.092057 12.763720 12.810414 12.885694 7.988800 7.999035 8.009345 6.253833 6.257591 6.260922 3.281004 3.330776 3.340302 2.559330 2.576535 2.585771 28.188112 28.190450 28.192843 16.659357 16.669062 16.709116 12.231450 12.242376 12.260819 9.232740 9.500731 9.696650 6.199129 6.324196 6.363759 4.786237 4.935394 5.762550 2.570710 2.597546 2.606910 55.079761 55.449983 55.884757 32.697014 32.821639 33.700782 23.986052 24.002592 24.010874 18.826771 18.860131 18.898127 11.667000 11.694999 11.750420 8.877861 9.148463 9.311218 30.285715 30.312009 30.344575 17.699522 17.716013 17.782700 13.033126 13.133019 13.240698 17.894299 17.991793 18.132221

%Imp KBTG SYS -0.2276 0.0168 -2.8461 0.0291 -6.3794 -0.2066 -2.9366 0.0042 -5.1945 0.0198 -5.6041 0.0030 -0.2008 0.3784 -0.9060 0.3381 -3.6464 0.0408 -1.2424 -0.1221 -4.3987 -0.2528 -3.2695 0.5026 -4.0230 1.9111 -1.9173 0.8647 -0.3774 0.0948 -2.0159 0.3683 -3.3786 -0.3299 -2.8178 -0.4675 -4.2485 1.0135 -1.8922 0.7028 -2.4503 0.0354 -0.3659 0.2307 -3.4731 0.0285 -3.9863 -0.3783 -3.5695 1.0478 -6.1070 0.2777 -4.7752 0.5042 -0.2984 -0.0597 -2.7762 -0.4671 -2.7060 -0.3405 -2.9344 0.1929

CR IP T

145 146 147 148 149 150 151 152 153 154 155 156 157 158 160 161 162 163 164 165 166 169 170 171 172 173 174 175 176 177

AN US

Inst. KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG Av.

Table 17: Performance of LSBM versus the two best heuristics of the literature on instances containing 40 spheres.

AC

Local Search-Based Method Min Mean Max 64.555790 64.867037 65.011581 36.917431 37.113444 37.174114 26.493228 26.557787 26.584077 24.561734 25.095353 25.486597 15.201655 15.356063 15.434481 11.364147 11.372404 11.399967 37.224077 37.934539 38.184683 21.935569 22.025407 22.129128 16.102862 16.344374 16.412241 14.241450 14.266710 14.299231 9.283430 9.390962 9.545100 6.973666 7.120984 7.185059 4.284771 4.361379 4.382747 3.008577 3.023019 3.070892 2.529171 2.554360 2.578750 28.537304 28.614173 28.724305 17.181054 17.272095 17.399908 12.574792 12.710203 12.783483 10.677965 10.844745 10.970605 7.031318 7.142945 7.195949 5.260015 5.337179 5.358390 3.141408 3.208080 3.255846 2.421601 2.446547 2.451083 61.001394 61.111574 61.163093 35.349801 35.371150 35.393524 25.127625 25.366143 25.576266 25.110056 25.231177 25.314150 15.504461 15.566963 15.600274 11.372817 11.410211 11.447622 39.793777 40.234657 40.658889 19.825432 19.975055 20.072401

M

SYS 64.5559 36.9166 26.4942 24.5354 15.2018 11.3203 37.7602 21.8689 16.1218 14.2417 9.2844 6.9775 4.2864 3.0086 2.527 28.5374 17.1816 12.7099 10.6768 7.0079 5.2443 3.1342 2.4059 60.9059 35.3499 25.1186 25.1033 15.5027 11.4661 39.7893 19.84115

ED

KBTG 65.0059 38.6265 27.9506 25.6006 15.9979 11.957 38.1222 22.5376 16.8623 14.8725 9.8161 7.4086 4.5033 3.1262 2.6422 28.6731 17.4944 13.0998 10.9874 7.2682 5.5444 3.3505 2.5055 61.0556 36.2016 27.0162 26.0798 16.3494 12.1584 39.9661 20.42599667

PT

217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 241 242 243 244 245 246 247

n 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50

CE

Inst. KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG KBG Av.

%Imp KBTG SYS -0.6972 -0.0002 -4.6294 0.0023 -5.5009 -0.0037 -4.2296 0.1072 -5.2379 -0.0010 -5.2169 0.3858 -2.4127 -1.4403 -2.7445 0.3039 -4.7162 -0.1176 -4.4311 -0.0018 -5.7379 -0.0105 -6.2368 -0.0550 -5.1001 -0.0380 -3.9096 -0.0008 -4.4690 0.0859 -0.4759 -0.0003 -1.8238 -0.0032 -4.1751 -1.0744 -2.8979 0.0109 -3.3690 0.3331 -5.4065 0.2988 -6.6560 0.2295 -3.4646 0.6484 -0.0889 0.1565 -2.4096 -0.0003 -7.5159 0.0359 -3.8620 0.0269 -5.4497 0.0114 -6.9075 -0.8202 -0.4330 0.0113 -4.0068 -0.0306

Table 18: Performance of LSBM versus the two best heuristics of the literature on instances containing 50 spheres.

References [1] H. Akeb. A look-ahead-based heuristic for packing spheres into a bin: the knapsack case. Procedia Computer Science, Proceedings of the International Conference on Communications, management, and Information technology (ICCMIT 2015), Vol. 65, pp. 28

ACCEPTED MANUSCRIPT

652-661, 2015. [2] H. Akeb. A two-stage look-ahead heuristic for packing spheres into a three-dimensional bin of minimum length, in Recent Advances in Computational Optimization (S. Fidanova, Ed., Springer International Publishing), Vol. 610, pp. 127–144, 2016.

CR IP T

[3] E. G. Birgin and F. N. C. Sobral. Minimizing the object dimensions in circle and sphere packing problems, Computers and Operations Research, Vol. 35, No. 7, pp. 2357–2375, 2008. [4] R. S. Farr. Random close packing fractions of lognormal distributions of hard spheres, Powder Technology, Vol. 245, pp. 28–34, 2013. [5] T. Kubach, A. Bortfeldt, T. Tilli and H. Gehring. Greedy algorithms for packing unequal sphere into a cuboidal strip or a cuboid. Asia-Pacific Journal of Operational Research, Vol. 28, pp. 739-753, 2011.

AN US

[6] T. Kubach, A. Bortfeldt, T. Tilli and H. Gehring. Parallel greedy algorithms for packing unequal spheres into a cuboidal strip or a cuboid. (Working Paper No 440). Hagen: Department of Management Science, University of Magdeburg, (Diskussionsbeitrag der Fakulty Wirtschaftswissenschaft der Fern University in Hagen), No 440, 2009. [7] J. A. George, J. M. George and B. W. Lamar. Packing different-sized circles into a rectangular container. European Journal of Operational Research, Vol. 84, No 3, pp. 693–712, 1995.

M

[8] M. Hifi. An iterative rounding search-based algorithm for the disjunctively constrained knapsack problem. Engineering Optimization, Vol. 46, No 8, pp. 1109–1122, 2014.

ED

[9] M. Hifi and R. M’Hallah. A literature review on circle and sphere packing problems: models and methodologies. Advances in Operations Research, Article ID 150624, 22 p, 2009 (doi.org/10.1155/2009/150624).

PT

[10] M. Hifi and M. Michrafy. A reactive local search-based algorithm for the disjunctively constrained knapsack problem. Journal of the Operational Research Society. Vol. 57, No 6, pp. 718–726, 2006.

CE

[11] M. Hifi and R. M’Hallah. Beam search and non-linear programming tools for the circular packing problem. International Journal of Mathematics in Operational Research, Vol. 1, No 4, pp. 476-503, 2009.

AC

[12] M. Hifi and R. M’Hallah. Approximate algorithms for constrained circular cutting problems. Computers & Operations Research, Vol. 31, No 5, pp. 675–694, 2004. [13] M. Hifi and L. Yousef. Width beam and hill-climbing strategies for the three- dimensional sphere packing problem, In Proceedings of IEEE, Federated Conference on Computer Science and Information Systems (FedCSIS), ACSIS, Vol. 2, pp.421–428, 2014 (doi: 10.15439/2014F284). [14] M. Hifi and L. Yousef. A dichotomous search-based heuristic for the three-dimensional packing problem, Cogent Engineering, (Taylor & Francis Edition), Vol. 1, pp.1–15, 2015.

29

ACCEPTED MANUSCRIPT

[15] M. Hifi and L. Yousef. Handling lower bound and hill-climbing strategies for sphere packing problems, In S. Fidanova (Ed.), Springer, Recent Advances in Computational Optimization Studies in Computational Intelligence, Vol. 610, pp. 145–164, 2016. [16] W-Q. Huang, Y. Li, H. Akeb C-M. Li. Greedy algorithms for packing unequal circles into a rectangular container. Journal of the Operational Research Society, Vol. 56, pp. 539–548, 2005.

CR IP T

[17] R. M’Hallah, A. Alkandari, and N. Mladenovic. Packing unit spheres into the smallest sphere using VNS and NLP. Computers and Operations Research, Vol. 40, No. 2, pp. 603–615, 2013. [18] R. M’Hallah and A. Alkandari. Packing unit spheres into a cube using VNS. Electronic Notes in Discrete Mathematics, Vol. 39, No. 1, 201–208, 2012.

AN US

[19] Y. Li and W. Ji. Stability and convergence analysis of a dynamics-based collective method for random sphere packing. Journal of Computational Physics, Vol. 250, pp. 373–387, 2013. [20] Y. Stoyan and T. Romanova. Mathematical models of placement optimisation: twoand three-dimensional problems and applications. In G. Fasano and J. Pint´er (Eds), Modeling and Optimization in Space Engineering, Springer Optimization and Its Applications, Vol. 73, pp. 363–388, 2013.

M

[21] G. Scheithauer, Y. Stoyan and G. Yaskov. Packing non-equal spheres into containers of different shapes. Working Paper, MA TH-NM-07-2013, Technische Universit¨at Dresden Herausgeber: Der Rektor, 2013.

ED

[22] K. Soontrapa and Y. Chen. Mono-sized sphere packing algorithm development using optimized Monte Carlo technique. Advanced Powder Technology, Vol. 24, pp. 955–961, 2013.

PT

[23] P. I. Stetsyuk, T. E. Romanova and G. Scheithauer. On the global minimum in a balanced circular packing problem. Optimization Letters, Vol. 10, No 6, pp. 1347– 1360, 2016.

CE

[24] A. Sutou and Y. Dai. Global optimization approach to unequal sphere packing problems in 3D. Journal of Optimization Theory and Applications, Vol. 114, pp. 671–694, 2002.

AC

[25] Y. G. Stoyan and A. Chugay. Packing cylinders and rectangular parallelepipeds with distances between them into a given region. European Journal of Operational Research, Vol. 197, No 2, pp. 446–455, 2008. [26] Y. Stoyan, G. Yaskov, and G. Scheithauer. Packing of various radii solid spheres into a parallelepiped. Central European Journal of Operational Research, Vol. 11, pp. 389– 407, 2003. [27] Y. G. Stoyan, G. Scheithauer and G. N. Yaskov. Packing unequal spheres into various containers. Cybernetics and Systems Analysis, Vol. 52, No 3, pp. 419–426, 2016. [28] Stoyan and G. Yaskov. Packing equal circles into a circle with circular prohibited areas. International Journal of Computer Mathematics, Vol. 89, No 10, pp. 1355–1369, 2012. 30

ACCEPTED MANUSCRIPT

[29] J. Wang. Packing of unequal spheres and automated radiosurgical treatment planning, Journal of Combinatorial Optimization, Vol. 3, No. 4, pp. 453–463, 1999. [30] H. Wang, W. Huang, Q. Zhang and D. Xu. An improved algorithm for the packing of unequal circles within a larger containing circle. European Journal of Operational Research. Vol. 141, pp. 440–453, 2002.

AC

CE

PT

ED

M

AN US

CR IP T

[31] G. W¨ ascher, H. Haussner and H. Schumann. An improved typology of cutting and packing problems. European Journal of Operational Research, Vol. 183, pp. 1109– 1130, 2007.

31