A Literature Review on Circle and Sphere Packing Problems: Models and Methodologies

Abstract

This paper reviews the most relevant literature on efficient models and methods for packing circular objects/items into Euclidean plane regions where the objects/items and regions are either two- or three-dimensional. These packing problems are NP hard optimization problems with a wide variety of applications. They have been tackled using various approaches-based algorithms ranging from computer-aided optimality proofs, to branch-and-bound procedures, to constructive approaches, to multi-start nonconvex minimization, to billiard simulation, to multiphase heuristics, and metaheuristics.

1. Introduction

Cutting and Packing problems are scientifically challenging problems with a wide spectrum of applications. They are very interesting NP-hard combinatorial optimization problems; that is, no procedure is able to exactly solve them in deterministic polynomial time. They are encountered in a variety of real-world applications including production and packing for the textile, apparel, naval, automobile, aerospace, and food industries. They are bottleneck problems in computer aided design where design plans are to be generated for industrial plants, electronic modules, nuclear and thermal plants, and so forth [1].

Cutting and packing problems consist of packing a set of geometric objects/items of fixed dimensions and shape into a region Ω of predetermined shape while accounting for the design and technological considerations of the problem [1]. The packing identifies the arrangement and positions of the geometric objects that determine the dimensions of the containing shape and reach the extremum of a specific objective function [1]. When the geometric objects are 𝜑 objects, the packing problem reduces to a mathematical model whose constraints reflect the conditions of arrangement of 𝜑 objects within the given region, the mutual nonintersection conditions for 𝜑 objects, and any technological constraints that cannot be reduced to purely geometric constraints [1] (as guillotine cuts for rectangular packing). One approach to solve this complex problem is to arrange the 𝜑 objects according to some prescribed order in Ω, and to search for the exact local extremum [1]. However, the search for exact local extrema is time consuming without any guarantee of a sufficiently good convergence to optimum [1].

This paper focuses on problems dealing with circular 𝜑 objects, being two- or three-dimensional. These problems have a wide spectrum of applications ranging from natural sciences, to engineering design, to every day life. They include the coverage of a geographical area with cell transmitters, storage of cylindrical drums into containers or stocking them into an open area, packaging bottles or cans into the smallest box, planting trees in a given region as to maximize the forest density and the distance between the trees, and so forth [2].

Packing circular objects is a challenge to discrete and computational geometry [2]. With a large number of circular objects to pack, the optimal solution is very difficult to find. An optimal solution may be rotated, reflected or the circular objects reordered; hence, the number of equivalent optimal solutions blows up as the number of circular objects increases [2–4]. In addition, one or more of the circular objects may be moved slightly without affecting the optimal solution. In fact, there exists a continuum of optimal solutions [2–4]. Last, there is the issue of computational accuracy and precision.

Obviously, packing circular objects gives rise to optimization problems, but their classification into continuous or discrete problems is fuzzy. Indeed, the positions of the circular objects are continuous whereas the structure of an optimal pattern has a discrete nature [5]. The recent trend is to design solution techniques that tackle (to some extent) these two aspects simultaneously.

Many variants of packing circular objects in the plane have been formulated as nonconvex optimization problems; with the hope that these formulations when solved using available nonlinear programming (NLP) solvers yield high quality approximate solutions [6]. However, most NLP solvers fail to identify global optima. In fact, problems with circular objects are very hard optimization problems. They have a large number of variables and local minima. Thus, they require to be tackled with algorithms which mix local searches with heuristic procedures in order to widely explore the search space [7, 8]. The recent trend has been to adopt approximation methods (heuristics) which combine methods of global (heuristic) search and methods of local exhaustive (exact) search of local minima or their approximations [9].

This paper provides an extensive survey of recent literature on packing circular objects. Literature can be classified into two categories: papers seeking optimal solutions or demonstrating the optimality of some patterns, and papers attempting to improve previously published results [10]. Independently of their classification, these papers search for the best solutions regardless of the margin of improvement or the computational expense [10]. Section 2 considers the case where the region Ω is a regular polygon whereas Section 3 focuses on the case the region Ω is circular. Finally, Section 4 is a summary.

2. Polygonal Region Ω

In this section, we distinguish between the two- and the three-dimensional case.

2.1. Two-Dimensional Case

Packing circles in a two-dimensional geometrical form such as a unit square or a unit-side triangle [11] is the best known type of extremal planar geometry problems [12]. Herein, the cases where the region Ω is a square, a rectangle, and a polygon are discussed.

2.1.1. Packing Circles into a Square

One of the problems that attracted most attention consists in packing 𝑛 identical circles within a unit square with the objective of maximizing the minimum distance between the centers of the circles; that is, of identifying the maximum radius 𝑟 of 𝑛 identical circles that fit into a square whose side length is one (i.e., fixed and known). The problem, which is a geometrical one, can be viewed as a continuous global optimization one. It can be stated as finding the optimal level of the decision variables 𝑟 and (𝑥_𝑖,𝑦_𝑖),𝑖∈𝐼={1,2,…,𝑛}, thus Maximize𝑟Subjectto𝑟≤𝑥_𝑖≤1−𝑟𝑖∈𝐼𝑟≤𝑦_𝑖≤1−𝑟𝑖∈𝐼(𝑥_𝑖−𝑥_𝑗)²+(𝑦_𝑖−𝑦_𝑗)²≥2𝑟1≤𝑖<𝑗≤𝑛, (2.1) where (𝑥_𝑖,𝑦_𝑖),𝑖∈𝐼, denote the coordinates of the center of circle 𝑖, and (𝑥_𝑖−𝑥_𝑗)²+(𝑦_𝑖−𝑦_𝑗)² is the Euclidean distance separating the centers of circles 𝑖 and 𝑗,1≤𝑖<𝑗≤𝑛.

This packing problem is equivalent to the problem of scattering 𝑛 points in a unit square such that the minimal distance 𝑑 between any pair of points is maximized. This scattering problem can be stated as a difference of convex functions program [13] or as an all-quadratic optimization problem [14] or as a continuous nonlinear inequality-constrained global optimization one [15]. The first two formulations are interesting to mathematicians since they are hard optimization problems; however, they are of little practical value [16]. The last formulation, which uses the geometrical aspect of the problem, is more effective. It can be stated as Maximize𝑑Subjectto0≤𝑥^′_𝑖≤1𝑖∈𝐼0≤𝑦^′_𝑖≤1𝑖∈𝐼𝑑≤(𝑥^′_𝑖−𝑥^′_𝑗)²+(𝑦^′_𝑖−𝑦^′_𝑗)²1≤𝑖<𝑗≤𝑛, (2.2) where (𝑥^′_𝑖,𝑦^′_𝑖),𝑖∈𝐼, denote the coordinates of point 𝑖. The optimal solution of the original problem is related to that of the scattering problem as follows: 𝑟=𝑑2(,𝑥𝑑+1)_𝑖=𝑥^′_𝑖(1−2𝑟)+𝑦𝑟,_𝑖=𝑦^′_𝑖(1−2𝑟)+𝑟. (2.3) Maranas et al. [15] solve this scattering problem using MINOS and the GAMS modeling language. They employ a multistart strategy to enhance the probability of finding global solutions for problems with up to 30 circles.

Nurmela and Oestergard [17] introduce, for the circle packing into a unit square, a different nonlinear programming formulation based on an energy function ∑_𝑖≠𝑗(𝜆/𝑑²_𝑖𝑗)^𝑚, where 𝑑_𝑖𝑗 represents the Euclidean distance between the centers of circles 𝑖 and 𝑗,𝜆 is a scaling factor, and 𝑚 is a positive integer. The energy function is to be minimized subject to the items fitting inside the square and not overlapping. The problem is transformed into an unconstrained one by an appropriate change of variables, and solved using a multistart hybrid line-search algorithm that uses gradient-type directions at the beginning and Newton-type directions near the solution. The local optimization is started from at least 50 randomly generated solutions. The solution obtained by the optimization algorithm is then used as a starting point to solve a system of nonlinear equations with the aim of improving the accuracy of the solution. Problems with 𝑛≤50 are solved. The approach finds some alternative solutions and improves some results. It recognizes some regular patterns that are presumably optimal for small 𝑛, but become nonoptimal for large values of 𝑛; for instance, the square lattice packing of 𝑛=𝑘² is optimal for 𝑘≤6 but not for 𝑘=7. Graham and Lubachevsky [18] extend the patterns of [17] using a billiard simulation that allows them to identify threshold indices above which it is guaranteed that the identified regular patterns become nonoptimal.

Boll et al. [19] propose a two-phase approach. The first phase is an approximation one which moves each point along the appropriately chosen directions with a step size that is exponentially decreased during the run. The second phase is a refining one, where the result of the first phase is the starting point for a billiard simulation. They improved the best known packing for 𝑛=32,37,48, and 50.

Casado et al. [20] subdivide the unit square into 𝑘×𝑘 subsquares, where 𝑘=0𝑥0𝑓8𝑒𝑒⌈√𝑛⌉0𝑥0𝑓8𝑓9. They construct an initial solution by placing the 𝑛 points randomly at the center of the 𝑛 distinct subsquares. They then perturb each point randomly. They may accept a perturbed point even when it is nonimproving (i.e., backtracking is allowed during the search). Based on these results, Szabo [21] develops some new regular patterns of points.

Locatelli and Raber [10] provide an upper bound on the maximum radius of 𝑛 identical circles that can fit into a unit square. They model the problem as a quadratic optimization problem, and introduce two properties that must be satisfied by at least an optimal solution. They develop a clever and efficient algorithm which starts from a general rectangular branch-and-bound, and exploits the special structure of the problem and the properties fulfilled by some of its solutions. They prove the optimality, within a 10⁻⁵ precision, of best known solutions in literature for 𝑛=10,…,31,33,34,35,38, and 39. They improve existing packings for 𝑛=32 and prove its optimality within the given tolerance and for 𝑛=37, but do not prove its optimality. Their algorithm provides a vivid example of the exponential growth of the computational effort as 𝑛 increases. For instance, trying to improve previously computed solutions requires less than 0.5 seconds for 𝑛≤13, more than 2 minutes for 𝑛≤21, and more than 30 minutes for 𝑛≤26. For 𝑛≤28, the computational time exceeds 27 CPU hours.

Markot and Csendes [22] and Szabo et al. [23] offer an efficient algorithm for eliminating large sets of suboptimal solutions of the scattering problem.

Markot and Csendes [24] use interval arithmetic as an accelerating device of an interval branch-and-bound optimization algorithm. Their computational proof lays the foundation for solving the previously open instances of packing 28, 29, and 30 circles.

Markot [25] presents a numerically reliable verified method using interval arithmetic computations, which can be regarded as a “computer-assisted verification method’’ for proving the structural properties of circle packing problems. He proves that the best-known packing structure for 28, 29, and 30 circles results in an optimal and (apart from symmetry and the occurrence of a free circle) unique packing. That is, in all three cases the guessed optimal structure is the unique optimal geometric solution of the particular problem.

Addis et al. [26] reformulate the problem as a mathematical program for which efficient local search procedures exist and identifying feasible solutions is relatively easy; that is, into a mathematical program with a linear objective function but with reverse convex quadratic constraints. In fact, they replace the Euclidean distance in the non-overlap constraints of the scatter point formulation by their squares; that is, 𝑑²≤(𝑥^′_𝑖−𝑥^′_𝑗)²+(𝑦^′_𝑖−𝑦^′_𝑗)²1≤𝑖<𝑗≤𝑛. (2.4) The authors conjecture that the problem possesses a funnel landscape, a feature that is commonly found in molecular conformation problems. They develop a stochastic search algorithm which consists in a population-based version of Monotonic Basin Hopping (MBH). MBH, like iterated search, moves from one local optimum to another without too much disrupting the structure of the original local optimum. It investigates the neighborhood of a local optimum by checking if its neighbors could lead to improving solutions, and stops if the number of nonimproving neighbors exceeds a certain prefixed threshold. A neighbor is obtained by perturbing the positions of all circles and using the perturbed solution as a starting point of a local search approach that would converge to a local optimum. The population-based version of MBH applies MBH to a set of local optima (a population of parents). Each local optimum (parent) produces a new set of solutions (a population of children). Subsequently, each solution of the new population (each child) is compared—based on a dissimilarity measure—with each element in the population of parents looking for the closest one. If the child is better than its closest parent, it replaces it. After all children have been compared with the population of parents and a new set of parents has been obtained, the procedure restarts. This procedure improves 32 best known solutions in the range 53≤𝑛≤130. In fact, a key to the success of the population based MBH is the dissimilarity measure which maintains the diversity of the population (thus guarantees the exploration of wider parts of the solution space) by maintaining a sufficient dissimilarity gap among the individuals of the population [8]. Different dissimilarity measures, mainly based on pairwise distances between cluster circles, are introduced [8]. Tests show that although there is no single dissimilarity measure which dominates the others, it is possible to identify a group of dissimilarity criteria which guarantees the best performance [8].

Van Dam et al. [27] consider a different but closely related problem to the problem of packing identical circles into a unit square. Some solutions of the packing problem have a strong shadow effect; that is, if projected over the 𝑥 or 𝑦 axis, the solutions have a number of distinct projected points that is much lower than the overall number of points. Therefore, the authors consider solutions with the additional constraint that the projected points are well separated. The authors propose both a branch-and-bound (which returns optimal solutions to the modified problem for 𝑛≤70), and a heuristic approach for 𝑛≤1000 points.

2.1.2. Packing Circles into a Rectangle

One of the problems having a wide industrial application is cylinder palletization. It consists in identifying the maximum number of identical circles (with a known radius 𝑟) that can be packed into a rectangle of fixed known dimensions (𝐿,𝑊) where the circles should be totally inside the rectangular box and should not overlap. This problem is classified as a single knapsack problem according to the typology of Wascher et al. [28].

Correia et al. [29] introduce a new upper bound for the optimal number of circles that fit inside a rectangle. In fact, existing bounds for similar problems cannot be adapted to the circle packing case. For instance, the bound which consists of the ratio of the area of the rectangle to the area of a circle is inefficient due to the existence of inevitable unused space among the circles. Other bounds exploit the geometric characteristics of packing rectangles into a rectangle; thus, are unsuitable for this case. Finally, using the area of the square hull of the circles does not yield an upper bound for the optimal value, although it can be used as a lower bound. The developed bound uses a better estimate of the usable area of the rectangle. It subtracts from the overall area of the rectangle an external area corresponding to the unused spaces along the pallet borders and an internal area corresponding to the unused spaces among circles, where the circles have regular arrangements. The analysis of the results obtained for five distinct sets of problems show that the developed bound dominates existing upper bounds, and that its gap from a lower bound, obtained by applying simulated annealing [30] to the cylinder packing problem, is very small.

Birgin et al. [31] tackle this problem by solving a series of decision problems. Each problem investigates the feasibility of packing 𝑘 identical circles of known radius 𝑟 into (𝐿,𝑊). If such a packing is feasible, 𝑘 is incremented by one and the decision problem is solved again. The algorithm stops when the decision problem yields an infeasible packing. The decision problem is formulated as Minimize_𝑖<𝑗max(0,4𝑟²−(𝑥_𝑖−𝑥_𝑗)²−(𝑦_𝑖−𝑦_𝑗)²)²Subjectto𝑟≤𝑥_𝑖≤𝐿−𝑟𝑖∈𝐼^′={1,…,𝑘}𝑟≤𝑦_𝑖≤𝑊−𝑟𝑖∈𝐼^′. (2.5) The decision problem is feasible if the value of the objective function is zero, and infeasible otherwise. In fact, if any pair of circles (𝑖,𝑗) overlap, the Euclidean distance separating their centers will be less than the sum of their radii: (𝑥_𝑖−𝑥_𝑗)²+(𝑦_𝑖−𝑦_𝑗)²<2𝑟, (2.6) which is equivalent to 4𝑟²−(𝑥_𝑖−𝑥_𝑗)²−(𝑦_𝑖−𝑦_𝑗)²>0, (2.7) and implies that max{4𝑟²−(𝑥_𝑖−𝑥_𝑗)²+(𝑦_𝑖−𝑦_𝑗)²,0}>0. (2.8) If the value of the objective function given by (2.5) is zero, then every one of the terms of the summation equals zero. Thus, none of the pairs of circles overlap.

The problem given by (2.5) is a nonlinear piecewise continuous optimization problem whose objective function has continuous first derivatives but discontinuous second derivatives. That is, the objective function changes its analytical definition on the boundary that separates the regions. This is evidently a disadvantage for minimization algorithms based on quadratic models, like Newton's method, which enjoys good convergence properties. The authors overcome the discontinuity issue by means of a continuous regularized Hessian which ensures more stability of the iterations. They solve the resulting problem using GENCAN, an augmented Lagrangian method for the minimization of a smooth function with general constraints. It uses, for the direction chosen at each step inside the faces, a truncated-Newton approach which employs the regularized Hessian approximations. This means that the search vector is an approximate minimizer of the quadratic approximation of the function in the current face. Conjugate gradients are used to find this direction. GENCAN with the regularized Hessian turns out to be much more efficient than the same method using the discontinuous Hessian. The method is restarted with 50 000 initial randomly generated solutions in search for a global minimum.

A variant of this problem is to consider nonidentical circles. George et al. [32] propose a nonlinear mixed integer program for packing nonidentical circles (of known radii) into a rectangle (𝐿,𝑊), of fixed known dimensions. The model regards (𝑥_𝑖,𝑦_𝑖),𝑖∈𝐼, as continuous variables and uses binary variables 𝛿_𝑖,𝑖∈𝐼, to signal the inclusion/exclusion of a circle from the pattern. The model can be stated as Maximize^𝑛_𝑖=1𝑐_𝑖𝛿_𝑖Subjectto𝛿_𝑖𝑟_𝑖≤𝑥_𝑖≤𝛿_𝑖(𝐿−𝑟_𝑖𝛿)𝑖∈𝐼_𝑖𝑟_𝑖≤𝑦_𝑖≤𝛿_𝑖(𝑊−𝑟_𝑖)𝑖∈𝐼(𝑥_𝑖−𝑥_𝑗)²+(𝑦_𝑖−𝑦_𝑗)²≥𝛿_𝑖𝛿_𝑗(𝑟_𝑖+𝑟_𝑗𝑥)1≤𝑖<𝑗≤𝑛,_𝑖≥0,𝑦_𝑖≥0,𝛿_𝑖∈{0,1}𝑖∈𝐼, (2.9) where 𝑐_𝑖,𝑖∈𝐼, is the profit obtained when including circle 𝑖 of radius 𝑟_𝑖 in rectangle (𝐿,𝑊). This model involves 𝑛 binary decision variables, 2𝑛 continuous decision variables, and 2𝑛+(𝑛(𝑛−1)/2) functional constraints reinforcing the fact that a packed circle is totally enclosed in the rectangle, and that no pair of packed circles overlap. This model is very difficult to solve using any off-the-shelf nonlinear mixed integer mathematical solver. Indeed, there are usually many different packing configurations which give the same objective function value. These feasible solutions differ only in the pattern in which the circles are packed. In addition, the feasible region is discontinuous; that is, an improving feasible solution occurs only when an additional whole circle can be placed in the rectangle. Moreover, the problem exhibits many local optima, and in many cases a globally optimal point is relatively rare compared to the number of second best solutions. In other words, the objective function consists of a small number of spikes in a large, relatively flat surface. Evidently, this model could not solve even small-sized problems. Therefore, the authors solved the problem using constructive heuristics which produce either unstable or stable configurations. A configuration is stable if every circle is adjacent to two sides of the containing rectangle, or to one side of the rectangle and to another larger circle, or to two circles, and is unstable otherwise. These heuristics apply a set of common sense rules: (a) first fit decreasing packing, (b) corner packing, (c) edge packing, (d) packing identical circles in a cluster, (e) generating stable configurations, (f) random positions, (g) spin-out, and (h) shake down. The heuristics are coupled with a genetic algorithm which chooses the position of every circle. The first fit rule (i.e., rule (a) is reinforced for all solutions; that is, the circles are sorted in a nonincreasing order of their diameters for all solutions. The authors show experimentally that the quasi-random rule (i.e., rule (f) coupled with the genetic algorithm yield the best solutions. Hifi et al. [33] adapt simulated annealing to this variant of the problem. They start from a superoptimal non feasible solution and restore feasibility by moving some circles out of the strip.

Another problem of interest is to find the dimensions of the rectangle of minimal area that can contain 𝑛 nonoverlapping congruent identical circles (of known fixed radius), where the height to width ratio of the rectangle is variable. In fact, when the aspect ratio of the enclosing rectangle is fixed, the densest packing minimizes both the area and the perimeter of the rectangle. However, by allowing a variable aspect ratio, the two optima may differ for the same number of circles 𝑛.

Lubachevsky and Graham [34] use a computational technique that consists of two independent algorithms: a restricted search and a compactor simulation. The restricted search operates on the assumption that the desired minimum is achieved on a set of configurations which is much smaller than the set of all possible configurations. The set is restricted to include only hexagonal patterns, square-grid patterns, their hybrids, and the patterns obtained by removing some circles from these patterns. A configuration is defined by 6 parameters: the number 𝑤 of circles in the longest row, the number of rows arranged in a hexagonal alternating pattern, those rows—among the hexagonally arranged rows—that have 𝑤−1 circles each, the number of rows that are stacked in the square-grid pattern, those among the square-grid rows that consist of 𝑤−1 circles each, and the number of “monovacancies’’ or holes. These six integers are subject to a set of constraints which induce the pattern and its enclosing rectangle. The restricted search procedure does not necessarily produce the global optimum. The compactor simulation generates a random starting feasible but sparse configuration with the 𝑛 circles lying inside a (large) rectangle. Then, the simulation imitates a compactor with each side of the rectangle pressing against the circles, so that the circles are forced towards each other until they jam. Possible circle-circle or circle-boundary conflicts are resolved using a simulation of hard collisions so that no overlap or boundary-penetrating circles occur during the process. The simulation is repeated many times, with different starting circle configurations, and the packing resulting in the smallest perimeter of the enclosing rectangle is retained as a candidate for the optimal packing. The compactor simulation makes no assumption regarding the resulting packing pattern; that is, the circles are free to choose any final configuration as long as it is jammed. However, this results in the simulation needing a large number of initial configurations prior to reaching a good candidate packing (of the same quality as that obtained by the restricted search). For example, it may take a fraction of a second to find the minimum perimeter packing of 15 circles by the restricted search but days with thousands of attempts to produce the same answer by simulation. Using these two techniques, the authors conjecture the minima for 𝑛≤5000.

Lubachevsky and Graham [35] note that in many of the identified best configurations, the circles form the usual regular square-grid or hexagonal patterns or their hybrids. However, for most values of 𝑛≤5000, for example, for 𝑛=7,13,17,21,22,26,31,37,38,41,43,…,4997,4998,4999,5000, they prove that the optimum cannot possibly be achieved by such regular arrangements. In general, the deviations of the optimal packing from the regular patterns are predictable small localized modifications to regular patterns. However, some of the identified best configurations show unexpected substantial extended irregularities. These correspond, for the range 𝑛≤62, to the instances with 𝑛=𝑘(𝑘+1)+1,𝑘=3,4,5,6,7, that is, for 𝑛=13,21,31,43, and 57. The height-to-width ratio of minimum perimeter rectangles containing 𝑛 congruent circles tends to 1 as 𝑛→∞.

Another problem of interest is packing circles into a rectangle of fixed length 𝐿 and minimum width 𝑊 without exceeding the dimensions of the rectangle or having overlapping circles; that is, the circle two-dimensional open dimension problem [28].

Stoyan and Yaskov [36] propose a mathematical model for the nonidentical circles variant of the problem Minimize𝑊Subjectto𝑟_𝑖≤𝑥_𝑖≤𝐿−𝑟_𝑖𝑟𝑖∈𝐼_𝑖≤𝑦_𝑖≤𝑊−𝑟_𝑖𝑖∈𝐼(𝑥_𝑖−𝑥_𝑗)²+(𝑦_𝑖−𝑦_𝑗)²≤(𝑟_𝑖+𝑟_𝑗)²1≤𝑖<𝑗≤𝑛𝑊≥0,𝑥_𝑖≥0,𝑦_𝑖≥0,𝑖∈𝐼. (2.10) The authors discuss the peculiarities as well as the characteristics of the feasible region of the resulting multidimensional multiextremal problem. They then solve it using a combination of branch-and-bound and the reduced gradient method. They search for the global extremum of the linear objective function by investigating all extreme points, which are sorted via a constructed search tree. Each extreme point is uniquely determined by an end-node of the search tree where the node corresponds to a quadratic system of 2𝑛+1 equations. They move from one extreme point to another having a smaller objective function value via a computed step size. They illustrate the application of their approach via a set of examples, and indicate that they find a global minimum if 𝑛≤8. However, for 𝑛>8, their approach can obtain, in general, only some approximation of the global minimum.

For the same problem and model formulation, Stoyan and Yaskov [9] offer an original method of jumping from one local minimum to another in search for an approximation of a global minimum. For a given local minimum, the method chooses, by means of Lagrangian multipliers, two nonidentical circles and interchanges their positions. That is, the radii of these two circles can be taken as variables; thus, the method is not suitable for packing identical circles. The method increases the dimension of the problem but permits the transition from one local minimum to another such that the transition decreases the objective function value. In fact, the solution method is based on simultaneously increasing the problem's dimension and choosing the movement direction from one local minimum to another via a reduced gradient method as well as on the concept of active inequalities and the Newton method. Since some systems of active inequalities can be inconsistent, a way to overcome the inconsistencies is considered. The final solution depends on the choice of the first extreme point. Therefore, obtaining a good approximation of a global minimum requires applying the afore-described approach with different starting extreme points. Computational tests show that the method, which can be viewed as a tunneling method—for escaping from local minimizers—coupled with an efficient local solver and a multistart strategy, is sufficiently effective for 𝑛 up to 35 circles. Indeed, it performs better than the branch-and-bound method of [36].

For the same problem, Hifi and M'Hallah [37] propose a constructive heuristic and a genetic algorithm. Both algorithms search for a good order of the pieces and use an adaptation of the best local position procedure (BLPP) [38, 39] to pack the circles. The adapted BLPP (ABLPP) positions circles in the upper left-most available position but takes advantage of the circularity of the pieces to explore more promising positions. ABLPP is simpler and faster than BLPP since a positioned circle cannot be further translated, and a circle can be positioned in holes generated by already packed ones.

To solve the same problem, Huang et al. [40] first study a decision problem which checks the possibility of packing 𝑛 circles into a rectangle of fixed dimensions. They solve the decision problem using one of two greedy algorithms: B1.0 and B1.5. Algorithm B1.0 selects the next circle to be positioned according to the maximum-hole degree (MHD) rule, which is inspired from human activity in packing. In fact, MHD, positions the current circle into the maximum hole degree corner position among all its feasible positions in Ω, that is, without overlapping any of the already packed circles. A corner position is a feasible position of the current circle in Ω if the current circle is tangent to two already packed circles or tangent to an already packed circle and to the perimeter of Ω. Algorithm B1.5 improves B1.0 with a self-look-ahead search strategy that determines, at each iteration, the circle to be packed and its position. The authors solve the original problem of determining the minimum width of a rectangle of fixed length and containing the 𝑛 circles, by applying a dichotomous search to get rapidly a good enough upper bound for the width. They, then, gradually reduce this upper bound until neither B1.0 nor B1.5 finds a successful configuration. The final upper bound is then taken as the minimal width 𝑊 of the rectangle.

Cui [41] formulates the same problem as a constrained circular cutting problem with the objective of minimizing material waste, and proposes a heuristic solution strategy that generates T-shaped cutting patterns.

Castillo et al. [42] model the same problem using 4𝑛 linear constraints, 𝑛(𝑛−1)/2 nonlinear constraints (as in [36]), and a set of simple bounds. The objective function can be of minimizing either the area of the containing rectangle or one of its dimensions. The authors apply off-the-shelf generic global optimization techniques. To further improve their results, they apply a posteriori strategy that, given a near-optimal initial arrangement, swaps all pairs of adjacent sized circles till no possible improvement is possible.

For the same problem, Akeb and Hifi [43] propose three heuristics. The first one positions the circles ordered in a nonincreasing order of their radii using the best local position rule. The second applies the first heuristic a fixed number of times using a different order of the circles each time. Finally, the third combines beam search with a dichotomous search for the width of the rectangle.

Birgin and Sobral [44] consider the problem of minimizing the dimensions of an object that contains 𝑛 nonoverlapping items where the object can be a two- or three-dimensional rectangle, triangle, square, or circle, and the items are circular. They develop twice differentiable nonlinear programming models that reduce the computational cost of reinforcing the nonoverlap constraints. In fact, the model has the classical constraints requiring each of the 𝑛 circles to be totally embedded in the containing circle, but replaces the nonoverlapping constraints by a single constraint: _𝑖<𝑗max(0,(𝑟_𝑖+𝑟_𝑗)²−(𝑥_𝑖−𝑥_𝑗)²−(𝑦_𝑖−𝑦_𝑗)²)²=0. (2.11) This constraint has two main advantages. Not all 𝑂(𝑛²) terms of the left hand side of the equation need to be evaluated (i.e., the number of pairs that contribute to the sum is 𝑂(𝑛)), and identifying the pairs of items that contribute to the sum can be undertaken in 𝑂(𝑛). The models are solved, starting from several randomly generated initial solutions, using a local solver that is based on an augmented Lagrangian method for smooth general constrained minimization.

2.1.3. Packing Circles in a Compact Polygon Set

Stoyan and Patsuk [45] consider the problem of covering a compact polygonal set by identical circles of minimal radius. They develop a mathematical model for the problem based on Voronoi polygons and investigate its characteristics. They then apply a modified version of the Zoutendijk feasible directions method to search for local minima, and design a special approach to choose the starting points. They illustrate the success of their proposed approach with many computational examples.

2.2. Three-Dimensional Case

Packing three-dimensional objects such as spheres arises in various branches of the industry [46]. For example, random packing of geometric objects has been used as a model to represent the structure of liquids and glassy materials; to study phenomena such as electrical conductivity, fluid flow, stress distribution, and other mechanical properties of granular materials; and to investigate processes such as sedimentation, compaction and sintering [46, 47]. Furthermore, they have applications in medicine for radio-surgical treatment planning [48], in powder metallurgy for three-dimensional laser cutting [49], in arranging and loading containers for transportation [31], in cutting different natural crystals, in layout of computers, buildings, and so forth. Herein, we consider the case where the region Ω is a cube, a parallelepiped, a cylinder, and a three-dimensional polytope.

Gensane [50] describes an adaptation of the billiard simulation for finding the largest radius 𝑟 of 𝑛identical spheres that can fit inside a cube. Billiard simulation is a stochastic method that simulates the idealized movement of billiard balls inside a domain, with the initial centers of the balls and their directions being randomly fixed. The obtained configuration is the result of these probabilistic choices. To improve the convergence of the stochastic algorithm, he introduces systematic perturbations in it. He considers four different versions of the simulation. In Algorithm 1, the hard spheres do not move along straight lines but randomly around themselves, and realize a random walk. In Algorithm 2, the magnitude of the stochastic movements is progressively reduced when the growing spheres yield a jam. Algorithm 3 introduces simultaneous perturbations of all the spheres. Algorithm 4 alternates the three types of perturbations. The author applies the perturbed billiard algorithm and obtains all the optimal and best known patterns for 𝑛 up to 32. He improves the previous best known solutions for all 11≤𝑛≤26 except 𝑛=13,14, and 18. He proves the existence of the displayed configurations for 𝑛=11,12,15,17,20,21,22,26, and 32, by constructing them explicitly. He conjectures that the minimum distance between spheres' centers of the optimal patterns is constant in the range 29≤𝑛≤32. He shows that contrary to the two-dimensional case, his and earlier billiard algorithms are unable to produce all optimal configurations with a good accuracy without using perturbations.

Stoyan and Yaskov [51] deal with the optimization problem of packing identical spheres of radius 𝑟 into a parallelepiped(𝐿,𝑊,𝐻) of minimal height 𝐻. They construct a mathematical model: Minimize𝐻Subjectto𝑟≤𝑥_𝑖≤𝐿−𝑟𝑖∈𝐼𝑟≤𝑦_𝑖≤𝑊−𝑟𝑖∈𝐼𝑟≤𝑧_𝑖≤𝐻−𝑟𝑖∈𝐼(𝑥_𝑖−𝑥_𝑗)²+(𝑦_𝑖−𝑦_𝑗)²+(𝑧_𝑖−𝑧_𝑗)²≤4𝑟²1≤𝑖<𝑗≤𝑛, (2.12) where (𝑥_𝑖,𝑦_𝑖,𝑧_𝑖) are the coordinates of sphere 𝑖. This model has a linear objective function but linear and quadratic constraints. The feasible region is generally disconnected with its connected component being 𝑛(𝑛−1)/2 connected. Its frontier is formed by 6𝑛 linear and 𝑛(𝑛−1)/2 quadratic inversely convex surfaces. Based on the peculiarities of this multi-extremal, NP-hard problem, the authors offer a solution strategy which includes a special search tree construction, a modification of the Zoutendijk method of feasible directions to calculate local minima, and a modification of the decremental neighborhood method to search for an approximation to the global minimum.

Stoyan and Yaskov [52] use a similar approach for packing identical spheres of radius 𝑟 into a right circular cylinder of known radius 𝑟 and minimal height 𝐻. The model they use is similar to that used in [51] except that the constraints regarding the 𝑥_𝑖 and 𝑦_𝑖 coordinates of a sphere 𝑖,𝑖∈𝐼, are replaced by 𝑟≤𝑥_𝑖≤2𝑟−𝑟𝑖∈𝐼,𝑟≤𝑦_𝑖≤2𝑟−𝑟𝑖∈𝐼, (2.13) respectively. The authors obtain the best results to date for 𝑛=498,499, and 500. Their approach is very effective for 𝑛≤500, and can handle instances with 𝑛≤2000 in very short computational times.

Stoyan et al. [53] use techniques similar to those considered in [9] to identify a packing of 𝑛nonidentical spheres, each of radius 𝑟_𝑖,𝑖∈𝐼, into a parallelepiped of fixed length 𝐿 and width 𝑊 but of variable height 𝐻 with the objective of minimizing 𝐻. The model is straightforward: Minimize𝐻Subjectto𝑟_𝑖≤𝑥_𝑖≤𝐿−𝑟_𝑖𝑟𝑖∈𝐼_𝑖≤𝑦_𝑖≤𝑊−𝑟_𝑖𝑟𝑖∈𝐼_𝑖≤𝑧_𝑖≤𝐻−𝑟_𝑖𝑖∈𝐼(𝑥_𝑖−𝑥_𝑗)²+(𝑦_𝑖−𝑦_𝑗)²+(𝑧_𝑖−𝑧_𝑗)²≤(𝑟_𝑖+𝑟_𝑗)²1≤𝑖<𝑗≤𝑛. (2.14) The authors provide numerical results with up to 60 spheres.Yaskov et al. [54] tackle the problem of maximizing the number of identical spheres that can be packed into a cylindrical composed domain. They construct a mathematical model based on the concept of Φ-functions [55], and design a solution algorithm based on a modification of the optimization method by groups of variables.

Wang [48] formulates mathematically the automated radio-surgical treatment planning problem as the packing of spheres into a three-dimensional region with a packing density greater than a given threshold level. He proves that this packing problem is NP-complete and proposes an approximate algorithm to solve it.

Sutou and Dai [56] assimilate the automated radio-surgical treatment planning problem to packing nonidentical spheres in a three-dimensional polytope with the objective of maximizing the sum of the volumes of the spheres packed in the polytope. They formulate the problem as a nonconvex optimization one with quadratic constraints and a linear objective function. On the basis of the special structures associated with this problem, they propose a variety of algorithms which improve markedly the existing branch-and bound algorithm for the general nonconvex quadratic program. They incorporate heuristic algorithms into the branch and bound to strengthen its efficiency. The computational study demonstrates the efficiency of the proposed algorithm for limited problem sizes.

Stoyan and Chugay [46] consider the problem of packing cylinders and parallelepiped shapes into a three-dimensional region so that the height of the occupied part of the region is minimal and the distances between each pair of items, and the distance between each packed item and the frontier of the region must be greater than or equal to given distances. A mathematical model of the problem is built and its characteristics are investigated. Methods for fast construction of starting points, searching for local minima, and a special non-exhaustive search of local minima to obtain good approximations to a global minimum are offered.

3. Circular Region Ω

When the 𝜑 objects and the region Ω are two-dimensional circles, the problem is referred to as the circle packing problem (CPP). CPP has many important applications in manufacturing, logistics, networks, facility layout, and materials science [42]. For example, in the automobile industry, design engineers have to estimate the size of the hole to be drilled on the body of the car and through which they plan to pass the bundle of wires that connect car's sensors to the display board [57]. The hole has to be large enough to allow all wires to pass, but as small as possible to avoid unnecessarily weakening the body [57]. CPP is also encountered in the manufacturing of sprockets for the motorcycle industry [58]. Similarly, it is of interest to the telecommunication/electrical/oil companies and refineries which have to pass bundles of different types of cables, pipes, and insulated pipes through cylindrical shapes over very long distances. The smaller the diameters of the cylinders, the cheaper the cost is. Finally, CPP emerges in material science where it is used to interpret topological relationships encountered when analyzing the normal grain growth in two dimensions [59] and to model certain absorption patterns of molecules [60].

This section distinguishes between the case the circles are identical and the case the circles are nonidentical.

3.1. Packing Identical Circles

Packing identical circles into a unit circle has been viewed as a generic problem with little industrial relevance. Yet, solving it remains challenging. Mladenovic et al. [6] apply a general reformulation descent heuristic (RD) to the problem of identifying the largest radius of identical circles that can be packed into a unit containing circle. RD iterates switching from solving CPP expressed in Cartesian coordinates to solving it expressed in polar coordinates and vice versa until no further improvement is obtained. The model expressed in Cartesian coordinates is Maximize𝑟Subjectto𝑥²_𝑖+𝑦²_𝑖≤(1−𝑟)²𝑖∈𝐼(𝑥_𝑖−𝑥_𝑗)²+(𝑦_𝑖−𝑦_𝑗)²≥4𝑟²(𝑖,𝑗)∈𝐼×𝐼,𝑖<𝑗−𝑟+𝑟_𝑖≤𝑥_𝑖≤𝑟−𝑟_𝑖𝑖∈𝐼−𝑟+𝑟_𝑖≤𝑦_𝑖≤𝑟−𝑟_𝑖𝑖∈𝐼𝑟≥0; (3.1) whereas its formulation in polar coordinates is Maximize𝑟Subjectto𝜌_𝑖𝜌+𝑟≤1𝑖∈𝐼_2𝑖+𝜌²_𝑗−2𝜌_𝑖𝜌_𝑗cos(𝜙_𝑖−𝜙_𝑗)≥4𝑟²𝑖∈𝐼,𝑗∈𝐼,𝑖<𝑗𝑟≥0,𝜌_𝑖≥0,𝜙_𝑖∈[0,2𝜋]𝑖∈𝐼; (3.2) where (𝜌_𝑖,𝜙_𝑖) denote the polar coordinates of circle 𝑖,𝑖∈𝐼. In fact, 𝜌_𝑖 is the distance of (𝑥_𝑖,𝑦_𝑖) from the origin of the polar coordinate system (which coincides herein with the center of the containing circle), and 𝜙_𝑖 is the angle of the sector delimited by the line joining (𝑥_𝑖,𝑦_𝑖) to the origin of the polar coordinates and the horizontal axis. That is, 𝑥_𝑖=𝜌_𝑖cos(𝜙_𝑖) and 𝑦_𝑖=𝜌_𝑖sin(𝜙_𝑖). RD allows nonlinear programming solvers using first-order information to escape stationary points. The experimental results for instances with up to 100 identical circles show that RD is 100 times faster than second-order nonlinear programming methods, and that the best known solution is obtained in 40% of the cases while in other cases the error never exceeds 1%.

Birgin and Sobral [44] solve their new model, starting from several initial solutions, using a local solver that is based on an augmented Lagrangian method for smooth general constrained minimization. Their reported solutions for instances with up to 100 identical circles coincide with the best known solutions (up to a prescribed tolerance).

Pinter [61] discusses the Lipschitz Global Optimizer (LGO) software which integrates global and local scope search methods to handle a very general class of nonlinear optimization models. In particular, he reviews the key features and basic usage of the LGO implementation linked to the General Algebraic Modeling System (GAMS). Then, he illustrates its application to packing 𝑛 identical circles into a unit circle where the objective is to determine the largest radius of the small circles.

3.2. Packing Nonidentical Circles

One variant of CPP consists in packing 𝑛 nonidentical circles without overlap into the smallest containing circle 𝐶 where each circle 𝑖 is characterized by its radius 𝑟_𝑖. The goal is to search for the best packing of the 𝑛 circles into 𝐶, where the best packing minimizes waste. CPP, according to the typology of Wascher et al. [28], is a two-dimensional variant of the Open Dimension Problem since all small items (which are circular) have to be packed and the extension of the large object (which is circular) is not given but has to be minimized. CPP is equivalent to finding the coordinates (𝑥_𝑖,𝑦_𝑖) of every circle 𝑖,𝑖∈𝐼, the radius 𝑟 and coordinates (𝑥,𝑦) of 𝐶, such that no pair of circles (𝑖,𝑗)∈𝐼×𝐼 and 𝑖<𝑗 overlap. Formally, the problem can be stated as finding the optimal level of the decision variables 𝑟,(𝑥,𝑦), and (𝑥_𝑖,𝑦_𝑖),𝑖∈𝐼, that is Minimize𝑟Subjectto(𝑥_𝑖−𝑥)²+(𝑦_𝑖−𝑦)²≤(𝑟−𝑟_𝑖)²𝑖∈𝐼(𝑥_𝑖−𝑥_𝑗)²+(𝑦_𝑖−𝑦_𝑗)²≥(𝑟_𝑖+𝑟_𝑗)²(𝑖,𝑗)∈𝐼×𝐼,𝑗<𝑖𝑟≥𝐿𝐵_𝑟, (3.3) where 𝐿𝐵_𝑟=∑_𝑖∈𝐼𝑟²_𝑖. The first set of constraints reinforce the complete containment of every circle within 𝐶. The second set reinforces the no overlap constraint of any pair of distinct circles. Finally, the last constraint provides a positive lower bound for the radius of the containing circle. It substitutes the nonnegativity constraint whose elimination from the model makes CPP unbounded. CPP has a linear objective function but non linear constraints. Even though its set of feasible solutions has piecewise smooth surface frontiers, CPP is a difficult problem to solve [62]. For all practical purposes, CPP can be solved by setting either (𝑥,𝑦)=(0,0) or any other pair of coordinates (𝑥_𝑖,𝑦_𝑖)=(0,0).

As in the case of packing circles into a square or a rectangle, Birgin and Sobral [44] reformulate the non-overlap constraints of CPP using the single constraint given by (2.11).

Alternatively, the problem can be formulated in polar coordinates: Minimize𝑟Subjectto𝜌_𝑖+𝑟_𝑖𝜌≤𝑟𝑖∈𝐼_2𝑖+𝜌²_𝑗−2𝜌_𝑖𝜌_𝑗cos(𝜙_𝑖−𝜙_𝑗)≥(𝑟_𝑖+𝑟_𝑗)²𝑖∈𝐼,𝑗∈𝐼,𝑖<𝑗,𝑟≥0,𝜌_𝑖≥0,𝑟_𝑖≥0,𝜙_𝑖∈[0,2𝜋]𝑖∈𝐼, (3.4) where (𝜌_𝑖,𝜙_𝑖) denote the polar coordinates of circle 𝑖,𝑖∈𝐼. In fact, 𝜌_𝑖 is the distance of its center from the center (𝑥,𝑦) of 𝐶, and 𝜙_𝑖 is the angle of the sector delimited by the line joining (𝑥_𝑖,𝑦_𝑖) to (𝑥,𝑦) and the horizontal axis. That is, 𝑥_𝑖=𝜌_𝑖cos(𝜙_𝑖)+𝑥 and 𝑦_𝑖=𝜌_𝑖sin(𝜙_𝑖)+𝑦.

CPP is a difficult problem to solve [62]. It poses several challenges. It cannot be tackled effectively by purely analytical approaches [42]. It embeds two extremely difficult problems: a pure continuous optimization problem, and a combinatorial one [63]. In addition, it has an infinite number of alternative optima [62], an exponential number of local optima which are not globally optimal, and an uncountable set of stationary points [6] (i.e., solutions that satisfy the KKT conditions but which do not correspond to local optima [63]). Moreover, CPP’s solution is sensitive to the choice of the initial solution (as is the case in nonlinear optimization) [6]. Finally, the larger the number of items, the larger the number of local minima and stationary points; thus, the simple multistart global optimization strategies need more local minimizations to reach a global minima [44].

CPP has been tackled by three distinct approaches based on: decision problems (where packings are obtained using the maximum hole degree rule) coupled with some bounding mechanism for 𝑟, constructive heuristics, and nonlinear programming.

The first approach fixes a value for the radius 𝑟 of the containing circle 𝐶, and solves—based on the concept of maximum hole degree or of simulation of elastic forces—a decision problem which checks if a feasible packing of the 𝑛 circles in 𝐶 is possible. When such a packing is feasible, the radius is decreased, and the process is repeated until no feasible packing is obtained.

Using this approach, Wang et al. [64] describe a quasi-physical quasi-human algorithm which mimics the physical model in which a number of smooth cylinders are packed inside a container. A quasi-human strategy is then proposed to trigger a jump for a stuck object in order to get out of local minima. The algorithm can be assimilated to an adaptive Tabu search. It randomly generates an initial pattern where every circle has its center inside the containing circle. It measures the infeasibility of a solution, and translates a circle, whose position is infeasible, along both axes. The translation distance is a function of an adaptive step size and the gradient of the objective function of the current pattern. Since the iterative process converges quickly to (in)feasible local optima, circles are allowed to jump in search for a new position within the containing circle.

Wenqi and Yan [65] formulate CPP as a potential energy function by simulating a system of elastic solids. They position all circles randomly inside the containing circle. If this configuration has no overlapping circles, a feasible solution is at hand. Otherwise, the elastic repulsion forces generated by the overlaps drive the overlapping circles to restore their shape and size. The circles move along straight lines, colliding with each other and with the containing circle until the composition of elastic forces is decreased to zero. If the amount of overlap is also decreased to zero, then the process stops with a feasible solution. Otherwise, the process restarts. The authors then give a quasi-physical algorithm that reaches a global minimum of the potential energy function by simulating the movements of all overlapping circles. They improve the quasi-physical algorithm by using two strategies: an early-escape strategy that helps the search escape from local minima, and a good-prospect strategy that generates good initial solutions for the search.

Sugihara et al. [57] apply a shrink and shake strategy. They construct a sufficiently large circle that contains all circles, then reduce its radius by translating circles using a shrink and shake iterative strategy. They further speed their strategy by using the circle Vonoroi diagram to determine the translation distance of any protruding circle—while keeping all circles inside the containing circle.

Zhang and Deng [66] adopt the model of Wang et al. [64], and use a hybrid approach consisting of simulated annealing to explore the neighborhood of the current solution, and tabu search to implement the jumps. When exploring the neighborhood of the current solution, one of the circles whose position is infeasible is translated and the degree of infeasibility of the neighbor is computed. A neighboring solution that reduces the degree of infeasibility becomes the incumbent solution whereas a nonimproving solution is accepted with a given probability, which decreases as the search becomes more selective. As this method often converges to infeasible local optima very quickly, a tabu search is adopted to allow circles causing infeasibility to jump out of their current position and randomly get a new position within the containing circle.

Huang et al. [67] solve CPP using two heuristics A1.0 and A1.5, which are extensions of the approach of [68]. They use the notion of the MHD of a position of a circle to be packed given a fixed value of 𝑟 and the positions of already packed circles. A1.0 chooses for the circle to be packed the feasible position having the highest hole degree. It is run 𝑛(𝑛−1) times; each time starting with a different pair of the 𝑛 circles. A1.5 is a modified version of A1.0. It applies a self look-ahead search for every feasible corner position of the circle to be packed. Given the set of 𝑖 packed circles, A1.5 packs 𝐶_𝑖+1 in every feasible corner position and uses A1.0 to pack all remaining items. If it successfully packs all items, A1.5 stops; else, it chooses for 𝐶_𝑖+1 the feasible corner position yielding the infeasible solution having the highest density and pursues packing the remaining circles. Experimental results, on benchmark instances with up to 100 circles, show that A1.5 has a good performance in terms of solution quality and computational time for packing unequal circles.

Huang and Chen [69] propose an improved version of the algorithm of Wang et al. [64] for solving CPP with equilibrium constraints. An efficient strategy of accelerating the search process is introduced in the gradient method.

Lu and Huang [70] incorporate the concept of maximum cave degree corner occupying actions into an improved PERM. The comparison of the proposed algorithm to existing ones shows that it is less efficient than Zhang and Deng's [66] for several large-scale identical circle instances—as is the case for all heuristics based on the concept of maximum hole degree—but obtains new lower bounds for several benchmark instances of nonidentical circles.

Akeb et al. [71] solve CPP using an adaptive beam search algorithm that combines beam search with the notion of local position distance and dichotomous search. It uses a width first beam search where the decisions at each node of the developed tree are based on the maximum hole degree that uses the local minimum distance. It escapes local minima by employing diversification strategies.

Akeb et al. [72] solve CPP using a binary search to determine 𝑟, and a beam search to check the feasibility of packing the 𝑛 circles into 𝐶 when the radius is 𝑟. A node of level ℓ,ℓ=1,…,𝑛, of the beam search tree corresponds to a partial packing of ℓ circles into 𝐶. The potential of each node of the tree is assessed via a lookahead strategy that, starting with the partial packing of the current node, assigns each unpacked circle to its maximum hole degree position. The beam search stops either when the lookahead strategy identifies a feasible packing or when it has fathomed all nodes.

The second approach is a constructive approach that successively packs items and searches for the smallest radius of the containing circle. For example, Hifi and M'Hallah [3] propose a dynamic, adaptive, local search algorithm which iterates through two complementary steps that are dynamically coupled in search for the smallest circle. At each iteration, the algorithm identifies the set of potential best local positions of a circle 𝐶_𝑖,𝑖∈𝑁, given the positions of the previously packed circles, and determines for each of these positions the coordinates and radius of the smallest containing circle. The best local position minimizes the radius of the current containing circle. That is, every time an additional circle is packed, both the center and the radius of the containing circle are dynamically updated, and the smallest containing circle is known. The procedure considers the circles in a nonincreasing order of their radii.

Hifi and M’Hallah [62] propose a three-phase approximate algorithm. During its first phase, the algorithm successively packs the ordered set of circles. It searches for each circle’s “best’’ position given the positions of the already packed circles where the best position minimizes the radius of the current containing circle. During its second phase, the algorithm tries to reduce the radius of the containing circle by applying (i) an intensified search—based on a reduction search interval—and (ii) a diversified search—based on the application of a number of layout techniques. During its third and last phase, the algorithm introduces a restarting procedure that explores the neighborhood of the current solution in search for a better ordering of the circles.

The third approach tackles CPP using random multistart global optimization run with off-the shelf optimizers. For example, Pinter and Kampas [73] present numerical results obtained using LGO. Castillo et al. [42] apply various off-the-shelf generic global optimization techniques, and compare their performance. They further improve the results of the generic solvers by implementing a posteriori strategy that, given a near-optimal initial arrangement, swaps all pairs of adjacent sized circles till no possible improvement is possible.

Addis et al. [63] mix standard local optimization routines with local moves between minima while reinforcing solution dissimilarity but reducing the solution space. The resulting approach obtains the best known solution for problems with up to 50 circles and 𝑟_𝑖=𝑖,𝑖=1,…,𝑛.

Hifi and M’Hallah [4] show that the combinatorial structure and the continuous optimization aspects of CPP should not be treated individually, but must be considered simultaneously. They base their recommendation on the comparison of the performance of two newly designed algorithms: BS1 and BS2. BS1 is a two-stage approach. The first stage uses a beam search (along with the packing heuristic of [3] and nonlinear optimization) to identify the best ordering of the circles. The second stage considers the circles in the order identified in stage 1, and uses a beam search to find the best position of each circle. BS2 embeds both searches into a single search-tree, where each level of the tree has two sub-levels: one for the circle to be packed, and one for the packing position of the circle. The proposed BS1 and BS2 are constructive in nature but use global optimization techniques too. They search for the best ordering of the circles and the best position of a circle given those already packed as to minimize the radius of the current containing circle. They further improve each local partial solution by applying global optimization techniques. They map most of the search space, guide the search into the most promising neighborhoods, allow escape from local minima, and reinforce solution dissimilarity.

Al-Modahka et al. [74] present an adaptive hybrid algorithm that addresses the combinatorial structure of CPP via a tabu search (TS), and its continuous optimization aspects via a combination of nested partitioning (NP) and nonlinear optimization. The hybrid TS/NP algorithm exploits the advantages of TS to undertake a local search aimed at identifying a good permutation of the circles whereas NP undertakes a global search to identify their respective best positions. The provided results are further modified/improved using some diversification strategies. Numerical examples show that the adaptive hybrid algorithm is both efficient and robust.

4. Conclusion

This paper reviews the recent literature on the NP hard optimization problem of packing circular items/objects into Euclidean regions of the plane. This challenging and difficult problem has many real-world applications, and has received a lot of attention. It has been tackled by constructive approaches coupled with meta-heuristics or local search methods, by simulation mimicking some physical or molecular phenomena, by branch- and bound-type approaches, and by nonconvex optimization. Any of these approaches tackles the continuous and the discrete structures of the problem using good standard local optimization routines, local moves between local minima, and enforcement of dissimilarity between local minima. The success of many of these approaches is due to the advancement of computer science related technology; subsequently, of the capacity of many nonlinear solvers to deal with large sized instances.

Acknowledgments

The authors thank two anonymous referees for their helpful comments. This research was partially supported by Kuwait University Research Grant SS02/06. The second author is grateful for their support.