PACKING CIRCLES WITHIN CIRCULAR CONTAINERS: A NEW HEURISTIC ALGORITHM FOR THE BALANCE CONSTRAINTS CASE

ABSTRACT

In this work we propose a heuristic algorithm for the layout optimization for disks installed in a rotating circular container. This is a unequal circle packing problem with additional balance constraints. It proved to be an NP-hard problem, which justifies heuristics methods for its resolution in larger instances. The main feature of our heuristic is based on the selection of the next circle to be placed inside the container according to the position of the system's center of mass. Our approach has been tested on a series of instances up to 55 circles and compared with the literature. Computational results show good performance in terms of solution quality and computational time for the proposed algorithm.

Keywords:
packing problem; layout optimization problem; nonidentical circles; heuristic algorithm

1 INTRODUCTION

We study how to install unequal disks in a rotating circular container, which is an adaptation of the model for the two-dimensional (2D) unequal circle packing problem with balance behavioral constraints. This problem arises in some engineering applications: development of satellites and rockets, multiple spindle box, rotating structure and so on. The low cost and high performance of the equipment require the best internal configuration among different geometric devices.

This problem is known as layout optimization problem (LOP), and consists in placing a set of circles in a circular container of minimum envelopment radius without overlap and with minimum imbalance. Each circle is characterized by its radius and mass. There, the original three-dimensional (3D) case (the equipment must rotate around its own axis) is simplified: different two-dimensional circles (see, Figure 1(c)) represent three-dimensional cylindrical objects to be placed inside the circular container.

Figure 1
Circular devices inside a rotating circular container and a feasible solution

Figure 1 illustrates the physical problem. Figure 1(a) shows a rotating cylindrical container. The symbol ω and the arrow illustrate the rotation around the axis of the equipment, ω is the angular velocity. In another viewpoint, Figure 1(b) shows the interior of the equipment where distinct circular devices need to be placed. In this example, six cylinders are placed, in which the radii, masses and heights are not necessarily equal.

Research on packing circles into a circular container has been documented and used to obtain good solutions. Heuristic, metaheuristic and hybrid methods are used in most of them. There are only a few publications discussing the disk problems with balance constraints.

LOP is a combinatorial problem and has been proved to be NP-hard (^{Lenstra & Rinnooy Kan, 1979}). This problem was first proposed by (^{Teng et al., 1994}), where a mathematical model and a series of intuitive algorithms combining the method of constructing the initial objects topo-models with the model-changing iteration method are described, and the validity of the proposed algorithms is verified by numerical examples. (^{Tang & Teng, 1999}) presented a modified genetic algorithm called decimal coded adaptive genetic algorithm to solve the LOP. (^{Xu et al., 2007}) developed a version of genetic algorithm called order-based positioning technique, which finds the best ordering for placing the circles in the container, and compare it with two existing nature-inspired methods. (^{Qian et al., 2001}) extended the work (^{Tang & Teng, 1999}) by introducing a genetic algorithm based on human-computer intervention, in which a human expert examines the best solution obtained through the loops of many generations and designs new solutions.

Methods based on particle swarm optimization (PSO) have been applied to the LOP. (^{Li et al., 2004}) developed a PSO method based-mutation operator. This approach can escape from the local minima, maintaining the characteristic of fast speed of convergence. (^{Zhou et al., 2005}) proposed a hybrid approach based on constraint handling strategy suit for PSO, where improvement is made by using direct search to increase the local search ability of the algorithm. (^{Xiao et al., 2007}) presented two nature-inspired approaches based on gradient search, the first hybrid with simulated annealing (SA) method and the second hybrid with PSO method. (^{Lei, 2009}) presented an adaptive PSO with a better search performance, which employs multi-adaptive strategies to plan large-scale space global search and refined local search to obtain global optimum.

(^{Huang & Chen, 2006}) proposed an improved version of the quasi-physical quasi-human algorithm proposed by (^{Wang et al., 2002}) for solving the disk packing problem with equilibrium constraints. An efficient strategy of accelerating the search process is introduced in the steepest descends method to shorten the execution time. In (^{Liu & Li, 2010}) the LOP is converted into an unconstrained optimization problem which is solved by the basin filling algorithm presented by them, together with the improved energy landscape paving method, the gradient method based on local search and the heuristic configuration update mechanism. (^{Liu et al., 2010}) presented a simulated annealing heuristic for solving the LOP by incorporating the neighborhood search mechanism and the adaptive gradient method. The neighborhood search mechanism avoids the disadvantage of blind search in the simulated annealing algorithm, and the adaptive gradient method is used to speed up the search for the best solution. (^{Liu et al., 2011}) developed a tabu search algorithm for solving the LOP. The algorithm begins with a random initial configuration and applies the gradient method with an adaptive step length to search for the minimum energy configuration. (^{He et al., 2013}) proposed a hybrid approach based on coarse-to-fine quasi-physical optimization method, where improvement is made by adapting the quasi-physical descent and the tabu search procedures. The algorithm approach takes into account the diversity of the search space to facilitate the global search, and it also does fine search to find the corresponding best solution in a promising local area. (^{Liu et al., 2015}) presented a heuristic based on energy landscape paving. The LOP is converted into the unconstrained optimization problem by using quasi-physical strategy and penalty function method. Subsequently, the heuristic approach combines a new updating mechanism of the histogram function in an improved energy landscape paving, and a local search for solving the LOP.

In this paper, we propose a new heuristic to solve the LOP. The basic idea of our approach, called center-of-mass-based placing technique (CMPT), is to place each circle according to the current position of the center of mass of the system.

Results for a selected set of instances are found in (^{Huang & Chen, 2006}), (^{Xiao et al., 2007}), (^{Lei, 2009}), (^{Liu & Li, 2010}), (^{Liu et al., 2010}, ²⁰¹¹, ²⁰¹⁵) and (^{He et al., 2013}). To validate our approach, we compare the results of our heuristic with these instances. Computational results show good performance in terms of solution quality and computational time.

The paper is organized as follows. Section 2 presents a formal definition of the unequal circle packing problem with balance constraints, and some definitions are established. Section 3 describes our heuristic. In Section 4, we present and analyze the experimental results, and Section 5 concludes the paper.

2 PROBLEM FORMULATION

We consider the following layout optimization for the disks installed in a rotating circular container: given a set of circles (not necessarily equal), find the minimal radius of a circular container in which all circles can be packed without overlap, and the shift of the dynamic equilibrium of the system should be minimized. The decision problem is stated as follows.

Consider a circular container of radius r, a set of n circles i of radii r_i and mass m_i , i ∈ N = {1, ..., n}. Let (x, y)^T be the coordinates of the container center, and (x_i , y_i )^T the center coordinates of the circle i. Let f₁ (z) = r be the first objective function, and

the second objective function, which measures the shift in the dynamic equilibrium of the system caused by the rotation of the container. Without loss of generality, we can consider ω = 1. The problem is to determine if there exists a (2n + 3)-dimensional vector z = (r, x, y, x ₁, y ₁, x ₂, y ₂, ..., x_n , y_n )^T that satisfies the following mathematical formulation.

(1)

(2)

where λ, β ∈ (0, 1) are a pair of preset weights, λ + β = 1. Constraint (1) states that circle i placed inside the container should not extend outside the container, while constraints (2) require that two any circles placed inside the container do not overlap each other.

Figure 1(c) illustrates a typical feasible solution to the LOP. The circles are numbered from 1 to 7, r ₃ is the radius of the circle 3, there is no overlap between the circles and the seven circles are completely placed into the larger circle of radius r (radius of the container).

2.1 Definitions

We develop a constructive heuristic guided by a simple strategy. A suboptimal solution is reached after gradually placing a circle at a time inside the container. Each circle is placed in an Euclidean coordinates system on the following evaluation criteria: select as the new position of the circle according to the current center of mass of the system without overlapping with the circles placed earlier; attempt to fill the wasted spaces after placing this circle; and in the end, select as the new coordinates of the container center that completely eliminate the dynamic imbalance ofthe system.

To perform the above criteria we need some notations and definitions. We denote by X(i) = (x_i , y_i )^T the center coordinates of the circle i, by

the Euclidean distance between the center coordinates of the circles i and j, and by Γ(i, j) = {(1 - λ)X(i) + λX(j): 0 ≤ λ ≤ 1} the set of points on the line segment whose endpoints are X(i) and X(j). Figure 2(a) illustrates the set Γ(10, 12).

Figure 2
Two Contact Cyclic Orders and a partial Layout.

Definition 1. (Contact Pair) If d(i, j) = r_i + r_j , we say that {i, j} is a Contact Pair of circles.

Definition 2. (Layout) A partial Layout, denoted by L, is a partial pattern (layout) formed by a subset of the m ≥ 2 of the circle centers, which have already been placed inside the container without overlap. Assume in addition that the container itself is in L. If m = n, then L is a complete layout (or solution).

Figure 2(b) illustrates a partial Layout formed by 16 circles placed inside the container. Among others, {5, 9} and {8, 15} are Contact Pairs.

Definition 3. (Placed Cyclic Order) Let C = i ₁ i ₂...i_{t - 1} i_t , i_p ∈ N, p = 1, ..., t, be a cyclic order of circles, which have already been placed inside the container without overlap. In addition, the intersection of any two sets Γ(i_{p 1}, i_{q 1}) and Γ(i_{p 2}, i_{q 2}) have at most one endpoint incommon, p ₁, q ₁, p ₂, q ₂ ∈ {1, ..., t} . We say that C is a Placed Cyclic Order.

Definition 4. (Contact Cyclic Order) Let C be a Placed Cyclic Order. If the circles are two by two Contact Pairs in C, we say that C is a Contact Cyclic Order.

Given a C = i ₁ i ₂...i_{p - 1} i_pi_{p + 1}...i_{t - 1} i_t , we say that the circles i ₁ i ₂...i_{p - 1} are in counterclockwise order in relation to circle i_p and the circles i_{p + 1}...i_{t - 1} i_t are in clockwise order in relation to circle i_p .

Definition 5. (Main Area) Let C be a Placed Cyclic Order. We say that the area bounded by the union of the line segments Γ(i_p , i_q ), where i_p , i_q are in C, is the Main Area of C, which is denoted by A(C).

Figure 2(a) illustrates a Contact Cyclic Order i, j} in Figure 2(b), in addition to are two by two Contact Pairs. In , it is illustrated a Contact Cyclic Order C formed by 4 circles. Note that all circles in C (dashed lines) are completely placed on A( formed by 12 circles placed inside the container. Note that all {). This is a feature of our approach, since several Contact Cyclic Order are obtained by circling each other. This approach is an important requirement, since it can yield a more compact layout.

Let N be a subset of circles placed inside the container and || be the cardinality of by the coordinates .. We denote the centroid of ⊂

Definition 6. (Border) Let L be a partial Layout and C be a Contact Cyclic Order. If thecenter of each circle in L belongs to A(C), we say in addition that C is the Border of the partial Layout L.

Figure 2(b) illustrates the Border L (16 circles). Note that all circle centers in L belong to A(C = 1|3|2|4 is not a Border of L, since X(5) ∉ A(C). (12 circles) of the partial Layout ). On the other hand,

We consider two cases of inclusion for placing circles. In the first case, we require that the circle k to be included must touch at least two previously placed circles. After this, in the second case, we require that another circle ℓ to be included occupies the wasted spaces after placing the circle k. This is a reasonable requirement, since it will generally yield a more compact layout than one defined by separate circles.

These two cases of inclusion can be explained by a partial Layout of the LOP example with seven existing circles illustrated in Figure 3. In Figure 3(a), it is shown the first case of inclusion. There are two positions to place the circle 9 (dashed lines) touching the Contact Pair {1, 7}, and two positions to place the circle 6 (dashed lines) touching the Contact Pair {2, 3}. Each position can be obtained by the solutions of the following particular case of the problems of Apolonio (^{Coxeter, 1968}).

(3)

Figure 3
Two cases of inclusion: (a) External Placement and (b) Internal Placement.

We denote by St(k, i_p , i_q ) the coordinates of the solution of the System (3) which does not belong to A(C). Note that the System (3) has two real solutions whenever d(i_p , i_q ) ≤ r_k . In Figure 3(a), by choosing St(9, 7, 1) as the coordinates of the circle 9, we obtain a feasible layout. However, it is not enough to choose St(6, 2, 3) as the coordinates of the circle 6, because the circle 6 overlaps the circles 1 and 5. + + 2

In our approach, we always select the coordinates St(k, i_p , i_q ) ∉ A(C) in order to place the new circle k touching the Contact Pair {i_p , i_q } in the Border C (in Figure 3(a), we have k = 6 and {i_p , i_q } = {2, 3}), but due to the potentially large differences in the radii, it is possible to occur overlap with the circles in the Border C. As it is illustrated in Figure 4(a), we get around this situation by repositioning the circle k to the coordinates of the solution of new System (3) for k, Figure 3(b), we have k = 6, (in and = 5). This first case of inclusion and the possible reposition define the following placement approach. and = 1 and , where now the circle k touches the circles

Figure 4
Example of the CMPT routine.

Definition 7. (External Placement) Let L be a partial Layout and C be the Border of L. An External Placement is the placement of a circle k inside the container, so that there is no overlap, its center does not belong to A(C), and it becomes Contact Pair with at least two circles in C. We denoted an External Placement by p_E (k).

The External Placement is always selected outside A(C), however if there is overlap on C, the repositioning of the new circle k (as explained above) is done in the following routine.

Procedure 1: External Placement routine

Input: a circle k, a Contact Pair {i_p , i_q }, a partial Layout L and the Border C

Output: an External Placement p_E (k), p and q

Step 1. Calculate St(k, i_p , i_q ) by System (3) and p_E (k) ← St(k, i_p , i_q ). If the circle k does not overlap any circles in C stop, otherwise go to Step 2.

Step 2. While there is overlap between the circle k and the circles in C repeat. If the circle k overlaps the circle C, p ← C, q ← , and choose pE(k) as the solution of the System (3) that is furthest from the centroid of the circles in L with respect to the Euclidean distance. furthest with respect to the counterclockwise order of the Border furthest with respect to the clockwise order of the Border , and if the circle k overlaps the circle

First, if the new circle k does not overlap any circles in Border C, the External Placement routine selects pE(k) = St(k, ip, iq). In our approach, this case is the most convenient way to place the next circle. However, if there is overlap, in Step 2 the routine identifies such circles () in order to reposition the circle k further from the centroid of the partial Layout, eventually avoiding any kind of overlap.and

To obtain a more compact layout, after including the circle k, it is checked the possibility of including another circle to occupy the wasted spaces after placing the circle k. We check among the remaining circles outside the container (preferably the largest one) if there is a circle ℓ that can be placed into the container in a centralized position without overlap. Each centralized position is the centroid coordinates of a certain set of circles which includes the circle k and the two circles touching the circle k. Figure 4(b) illustrates the second case of inclusion, which we can investigate the possibility of positioning a circle in the wasted space after placing the circle 9 touching the Contact Pair {1, 7} (centroid of the circles {1, 2, 3, 4, 5, 6}). of the circles {1, 7, 9}), and in the wasted space after placing the circle 6 touching the circles 1 and 5 (centroid

Definition 8. (Internal Placement). Let L be a partial Layout, C be the Border of L, and {i_p , k}, {k, i_q } be Contact Pairs in C, where k is the previous circle included. An Internal Placement is the placement of a circle ℓ inside the container, so that its center belongs to A(C), ℓ does not overlap with any other circle, and the center of ℓ is placed at the centroid coordinates of the subset . We denote an Internal Placement by pI(, where {k, ip, iq} ⊆ ,ℓ), meaning that ℓ is to be placed at the centroid coordinates XC of .

Let L be a partial Layout and C be the Border of L. In our algorithm, each positioning in the first case of inclusion is always done by looking at the Contact Pairs in the Border C. Suppose that the remaining circle k is selected to be placed touching the Contact Pair {ip, iq} in C. The placement of the circle k causes the addition of one element in L and one index in C, and perhaps the removal of some indices from C. This will be represented by the following operation.

O-+(k, ip, ip + s)(C) = C + {k} - {ip + 1, ..., ip + s - 1} = i1i2...ipkip + s...it - 1it,

where 1 ≤ s ≤ ⌊(t - 2)/2⌋.

With this choice for s there are fewer indices between p and p + s than p + s and p.

The operation O-+(k, ip, ip + s)(C) applied to C means that the circle k was placed inside the container touching the circles ip and ip + s without overlap. Then the index k is added to C, the subset of indices {ip + 1, ..., ip + s - 1} between p and p + s is removed from C, and the coordinates X(k) are added to the partial Layout L. Note that if s = 1 there is no removal of indices from C, and the index k is inserted in C between the indices ip and ip + 1.

The possible placement of the circle ℓ after the placement of the circle k only causes the possible addition of the coordinates X(ℓ) to the partial Layout L.

In our approach, we require the imbalance of the system be zero. It seems intuitive that this requirement may result in a good solution. We denote the center of mass of the system by

then one can shift the center of the rotating circular container to the center of mass of the system to have zero imbalance. This shift is made at each outer iteration and at the end of the algorithm, but it may increase the envelopment radius. Thus, if the Layout L represents a complete solution of the LOP, we denote the radius r of the container by

Moreover, the index where R(L) is reached is denoted by kmax ≡ arg(R(L)).

3 CENTER-OF-MASS-BASED PLACING TECHNIQUE (CMPT)

We present a new placing technique which yields compact layouts and quality solutions in an efficient manner. Let α = (α(1), α(2), ..., α(n)) be a permutation of (1, 2, ..., n). We place the circles in the partial Layout L one by one according to the order defined by this permutation. Given a order of inclusion, the first circles α(1), α(2), α(3) and α(4) must be positioned as follows.

Procedure 2: Initial layout routine

Input: the circles α(1), α(2), α(3) and α(4)

Output: an initial Layout L and the initial Border C

Place the circle α(1) at coordinates X(α(1)) = (0, 0). Choose the coordinates X(α(2)), such that α(2) touches α(1). For each circle α(3) and α(4), solve the System (3) and place them at coordinates X(α(3)) = (xα(3), yα(3)) = St(α(3), α(1), α(2)) and X(α(4)) = St(α(4), α(1), α(2)) without overlap. L = {X(α(1)), X(α(2)), X(α(3)), X(α(4))} and C ← α(1)|α(3)|α(2)|α(4).

Figure 2(b) illustrates the initial L = {X(1), X(2), X(3), X(4)} and the initial Border C = 1|3|2|4. Among many optional positions we can choose X(α(2)), for example, at the x-axis with the coordinates (rα(1) + rα(2), 0).

Suppose we have already placed the circles α(1), α(2), ..., α(k - 1), we describe our approach for placing the circle α(k), and after that, we verify the possibility of placing another circle α(ℓ), k < ℓ ≤ n.

When we place the circle α(k) (where k>4, see ^{Procedure 2}), we require that the circle touches at least two previously placed circles (see Figure 3(a) and ^{Procedure 4}). This will generally yield a more compact layout. However, we can increase the compactness of the layout if the wasted spaces after placing the circle α(k) can be occupied by another circle α(ℓ) (see Figure 3(b) and ^{Procedure 4}).

We observe that for each additional circle, the envelopment radius of a layout is generally enlarged. In order to minimize the rate of growth of this radius during inclusions, we must properly choose a new position for circle α(k) which yields a smaller envelopment radius. Our strategy CMPT attempts to reduce the rate of growth of the envelopment radius by including every circle around the coordinates of the center of mass of the system, which is updated during each outer iteration.

This strategy consists of shifting the origin of the Euclidean plane to the current center of mass of the system. Then we require that the circle α(k) touches the circles of a Contact Pair arbitrarily chosen among the elements of the Border C, taking into consideration the quadrants of the Euclidean plane. This approach is performed according to the following routine.

Procedure 3: CMPT routine

Input: a partial Layout L and the Border C

Output: the sets Q1, Q2, Q3 and Q4

Step 1. Calculate the coordinates of the center of mass XCM of the circles in L, and translate the origin of the Euclidean plane to XCM.

Step 2. Include each Contact Pair {ip, ip + 1} of C in the set Qh if the center of ip belongs to the quadrant h of the Euclidean plane, for h = 1, 2, 3, 4.

Given the Border C, the ^{Procedure 3} only separates the Contact Pairs in C according to the quadrants of the Euclidean plane with origin shifted to the current center of mass of the system.

Figure 4 illustrates the ^{Procedure 3}. In Figure 4(a) we observe that the coordinates of the center of mass XCM of the system do not coincide with the coordinates of the origin X(α(1)) = (0, 0). We wish to place the next circle α(k) around the coordinates XCM in order to mitigate the growth of the envelopment radius. By dividing the plane into quadrants, we can obtain a Border C with a more circular shape. We see in Figure 4(b) that if we place each new circle α(k) at a different quadrant of the Euclidean plane (with the origin shifted to XCM), then the layout is moreevenly distributed.

The choice of different quadrants (a Contact Pair in Qh (h = 1, 2, 3, 4)) to position the next circle α(k), and the operation O+- on the Border C lead to a updated Border C more similar to a circular shape. This will generally yield a more compact layout, because the wasted space between the Main Area A(C) and envelopment radius is minimized (see the example in Figure 5).

Figure 5
A suboptimal solution: 45 circles inside the container.

Next we describe the two cases of inclusion in the following routine.

Procedure 4: Inclusion routine

Input: a circle α(k), a permutation α, the sets Qh, h = 1, 2, 3, 4, a Contact Pair {ip, iq} in a set Qh, a partial Layout L and the Border C

Output: a partial Layout L, the Border C and the sets Qh, h = 1, 2, 3, 4

Step 1. Obtain an External Placement pE(α(k)) and the new values for p and q by ^{Procedure 1}. If there are fewer indices in the Border C between q and p than those between p and q, then p, p ← q and q ← ← .

Step 2. C ← O+-(α(k), ip, iq)(C), X(α(k)) ← pE(α(k)), L ← L ∪ {X(α(k)) and Qh ← Qh\{{ip, ip + 1}, ..., {iq - 1, iq}}, h = 1, 2, 3, 4 (note that q = p + s, where 1 ≤ s ≤ ⌊(t - 2)/2⌋).

Step 3. If it is possible to obtain an Internal Placement pI(k), i_p , i_{p + 1}, ..., i_{q - 1}, i_q } and a circle α(ℓ) (preferably the largest) in the permutation α, k < ℓ ≤ n, then X(α(ℓ)) ← p_I (L ← L ∪ {X(α(ℓ))}, and exclude α(ℓ) from α., α(ℓ)) for the set , α(ℓ)), = {α(

The Inclusion routine attempts to place the new circles in a more compact layout. First, it computes an External Placement for the next circle α(k) by ^{Procedure 1} and updates the values for p and q in order to obtain fewer indices between p and q than those between q and p. In Step 2, the Border C is updated by the operation O ⁺ _-, where the indices between p and q are removed from C (by Step 1, those indices correspond to the internal circles from the Border when compared to the indices between q and p), and α(k) is added to C. The circle α(k) is placed inside the container and all Contact Pairs between p and q (including {i_p , i_{p + 1}} and {i_{q - 1}, i_q }) are removed from the sets Q_h , h = 1, 2, 3, 4. This choice guarantees that the operation O ⁺ _- will exclude the internal circles from the Border C. Finally, a search to place another circle α(ℓ) (Internal Placement) is performed.

Post-optimization

A post-optimization is performed after the algorithm builds a complete solution (represented by Layout L), which contemplates improvements via circle repositioning at the Border C of L. This post-optimization process causes changes in C, where an index is removed and then it is repositioned in C by operation O ⁺ _-. The removal of the index from C will be represented by the following operation.

D(i_p )(C) = i ₁ i ₂...i_{p - 1} i_{p + 1}...i_{t - 1} i_t .

The operation D(i_p )(C) applied to C means that the circle i_p is deleted from its position. We delete the current X(i_p ) from the Layout L, and we test if a new position p_E (i_p ) for i_p improves the radius R(L) of the container.

Main routine

We choose to position each circle inside the container according to the following main procedure.

Main routine

Input: a permutation α = (α(1), α(2), ..., α(n))

Output: a Layout L (complete solution)

Step 1. (Initialization) Obtain the initial Layout L and the initial Border C by ^{Procedure 2}, k ← 5.

Step 2. (CMPT) Obtain the sets Q_h , h = 1, 2, 3, 4 by ^{Procedure 3}.

Step 3. (Layout construction) While there are Contact Pairs in any Q_h and circles outside the container, repeat for each h = 1, 2, 3, 4.

If Q_h ≠ ∅, choose an arbitrary {i_p , i_q } ∈ Q_h and include the circle α(k) and the possible circle α(ℓ), k < ℓ ≤ n, by ^{Procedure 4} and k ← k + 1.

Step 4. If there are circles outside the container, return to Step 2. Otherwise, go to Step 5.

Step 5. (Post-optimization) L, C, R(k _max in k ← k _max.), compute and ← ← ←

Step 5.1. (, excluding {ik - 1, ik} and {ik + 1, ik}., ) ← X(ik), delete the current X(ik) and repeat Step 5.2. for each Contact Pair {ip, iq} of

Step 5.2. Obtain an External Placement pE(ik) and the new values for p and q by ^{Procedure 1}, and X(ik) ← pE(ik). If the radius of the container is improved, then D(i_k )(O ⁺ _-(i_k , i_p , i_q )(L ← ← , C ← ), ← and return to Step 5.),

Step 5.3. X(i_k ) ← (, , whose container center is the center of mass of the system.) and finish the routine with the complete solution L ←

Given a permutation α, the ^{Main routine} builds an initial Layout L in Step 1 by placing the first four circles as in ^{Procedure 2}. Next, in Step 2 the main aspect of our approach is performed by ^{Procedure 3} (CMPT routine), where the Euclidean plane is divided into four parts and the subsets Qh (h = 1, 2, 3, 4) of Contact Pairs are obtained. Next, Step 3 is repeated by looking at each subset Qh and while there are circles remaining to be placed. In this step, an arbitrary Contact Pair in Qh is chosen and the two cases of inclusions are performed by ^{Procedure 4}. After we finish placing all circles inside the container, we obtain a complete solution L and its Border C. Then, in Step 5, a post-optimization is performed via circle repositioning at the Border C, which attempts improvements in the envelopment radius. In the end, the center of the container is shifted to center of mass of the system, which achieves zero imbalance.

Order of placement of the circles

As previously described, a permutation α = (α (1), α(2), ..., α(n)) of (1, 2, ..., n) is used as an input in our algorithm to generate a layout by specifying the order in which the circles are placed. Since there exist n! possible permutations for n circles, we need an appropriate technique in order to search in such a large space. Preliminary tests show that the wasted spaces after placing circles are minimized with greater efficiency when the order of addition of the circles favors those of larger radii.

Let α = ( α(1), α(2), ..., α(n)) be a sequence obtained by considering their radii in descending order of the circles, i.e., rα(k) ≥ rα(j), 1 ≤ k < j ≤ n. Choose an integer b, 1 ≤ b ≤ n and subdivide the terms of the sequence α in ℓ = ⌊n/b⌋ blocks. Thus, it is possible to obtain a subsequence bℓ), ..., α(n - 1), α(n) last elements of α. With this procedure, several subsequences to place the different circles may be generated. Actually, there are ((ℓ!)^b )(n - (bℓ))! possibilities, so that 1 ≤ ((ℓ!)^b )(n - (bℓ))! ≤ n!. Thus, when b = n, we only obtain the sequence b = 1, we can generate at most n! distinct subsequences. In our numerical experiments, for each instance of dimension greater than or equal to 10 we chose b = 5 to generate such sequences. = α, and when of α to be used as an input to the algorithm, by permuting the positions of the first α(1), ..., α(ℓ) elements of α, the α(ℓ + 1), ..., α(2ℓ) elements of α, and so on, until we permute the positions of α(

Complexity

The analysis of the real computational time of the ^{Main routine} is difficult, because it does not depend only on the number of circles, but also on the diversity of the circle radii and the number of circles in a current Border C, as well as the implementation. Here, we analyze the upper bound of the complexity of the ^{Main routine}, when it finds a complete solution L with Border C, such that |C| = λn, where 0 < λ ≤ 1, including the post-optimization process. Recall that, before post-optimization, the circles in Border C are two by two Contact Pair.

Given a partial Layout L with m circles already placed inside the container and n -m circles outside. Let |C| be the number of circles in the Border C of the partial Layout L.

The strategy CMPT in ^{Procedure 3} checks the position of |C| circles in the Euclidean plane, which is done in O(|C|).

When we position the circle i_k (where k > 2, see ^{Procedure 2}) touching two circles in Border C, |C| existing circles can define 2 × |C| positions, since two existing circles define two possible positions for the third i_k . To determine an External Placement for i_k , we must check the overlap with A(C) or with any circles in C. This is the same that we check the overlap with each circle in L, that is, m circles (a good implementation can reduce the number of checks). Because we assess 2 × |C| positions when we place the circle i_k , each time checking for overlaps m times, then the complexity to obtain an External Placement is about 2 × m × |C|.

After placing the circle i_k , we must check if there is a circle i _ℓ outside the container to be placed in an Internal Placement, then we must check the overlaps among n - m circles and a subset C, which is done in about |n - m) ≤ |C| × (n - m). in | × (

In the post-optimization process, we select one circle in C and assess at most (|C| - 2) Contact Pairs in C to try to improve of the envelopment radius by checking at most 2 × (|C| - 2) External Placements, thus the complexity of the post-optimization is bounded by 2 × (n - 1) × (|C| - 2).

Therefore the complexity of placing n circles during the ^{Main routine} is bounded by O(n ²|C|). After placing the new circle i_k , the operation O ⁺ _- modifies the Border C. This operation controls the size of C during the iterations. Since |C| ≤ n, the theoretical upper bound is O(n ³).

4 EXPERIMENTAL RESULTS

In this section, we measure the quality and performance of our algorithm on a series of instances up to 55 circles from literature. We tested three sets of instances from the literature. We compare our approach with a series of hybrid nature-inspired approaches based on simulated annealing and particle swarm optimization (^{Xiao et al., 2007}; ^{Lei, 2009}), a hybrid approach based on simulated annealing, neighborhood search mechanism and the adaptive gradientmethod (^{Liu et al., 2010}), a hybrid tabu search algorithm and gradient method (^{Liu et al., 2011}), a series of heuristics based on energy landscape paving, gradient method and local search (^{Liu & Li, 2010}; ^{Liu et al. 2010}, ²⁰¹⁵), and a series of algorithms based on quasi-physical approaches, gradient method and local search (^{Huang & Chen, 2006}; ^{He et al., 2013}; ^{Liu et al., 2015}).

These methods search for the optimal layout by directly evolving the positions of every circle, as well as considering imbalance. We use the benchmark suite of 21 instances of the problem described in Table 1 to test our algorithm. For each instance we present the range for r_i and m_i . A more detailed description of the instances can be found in (^{Huang & Chen, 2006}), (^{Liu & Li, 2010}) and (^{Xiao et al., 2007}).

Table 1
Data of each instance. First set of instances (^{Liu & Li, 2010}).

Size	Radii	Mass		Size	Radii	Mass
5	[20.71, 50]	[20.71, 50]		7	[8.5, 12]	[72.25, 144]
9	[0.17, 0.41]	[0.17, 0.41]		40	[81, 120]	[6, 14]
5	[6, 11]	[36, 121]
			Second set of instances (Huang & Chen, 2006)
Size	Radii	Mass		Size	Radii	Mass
7	[20, 20]	[400, 400]		17	[5, 25]	[25, 625]
12	[23.72, 100]	[562.64, 104]		37	[20, 20]	[400, 400]
15	[1, 15]	[1, 225]		50	[20, 20]	[400, 400]
			Third of instances (Xiao et al., 2007)
Size	Radii	Mass		Size	Radii	Mass
10	[5, 23]	[20, 93]		35	[7, 24]	[10, 99]
15	[6, 24]	[12, 98]		40	[6, 23]	[12, 99]
20	[5, 24]	[11, 94]		45	[6, 24]	[11, 99]
25	[6, 24]	[11, 96]		50	[5, 24]	[10, 99]
30	[6, 24]	[12, 97]		55	[6, 24]	[13, 99]

The routines were implemented in MATLAB language, and executed on a PC with an Intel Core i7, 7.7 GB of RAM and Linux operating system.

Except for the instances 1.1 and 1.3, both of them with 5 circles, we decide to generate 7! = 5040 distinct permutations α as input for the algorithm in each instance, i.e., we fixed in 5040 the number of executions of the ^{Main routine} for each instance and the best solution found was selected. This amount of tests proved adequate for our comparisons.

The results from the first, second and third sets of instances are presented in Tables 2, 3 and 4 respectively, where we compare our approach with those described in each indicated reference. The results are shown for the size of the instances, the best radius of the container obtained (first objective function f ₁), the imbalance obtained (second objective function f ₂), and the running time t (in seconds).

Table 2
Numerical results for the third set of instances

Instance	n	(Huang & Chen, 2006)			(Lei, 2009)			(Liu & Li, 2010)
		f 1	f 2	time(s)	f 1	f 2	time(s)	f 1	f 2	time(s)
1.1	5	120.710	7.63 × 10-4	3.82	-	-	-	120.710678	3.89 × 10-9	0.594
1.2	9	-	-	-	-	-	-	1.000000	9.7 × 10-9	2.891
1.3	5	22.327	3.90 × 10-4	2.91	-	-	-	22.249000	3.12 × 10-9	1.313
1.4	7	31.981	4.11 × 10-3	76.33	31.924	1.40 × 10-5	427	31.854600	1.5 × 10-9	54.828
1.5	40	742.750	5.40 × 10-4	12.27	769.819	3.25 × 10-4	1724	725.043500	2.47 × 10-9	7.187
Instance	n	(Liu et al., 2010)			(Liu et al., 2011)
		f 1	f 2	time(s)	f 1	f 2	time(s)
1.1	5	120.7106782	4.30 × 10-9	0.696	120.710678	3.25 × 10-9	0.350
1.2	9	-	-	-	1.000000	2.38 × 10-9	0.225
1.3	5	22.2463	3.52 × 10-9	0.469	22.246203	2.63 × 10-9	0.506
1.4	7	31.8412	2.03 × 10-10	58.797	31.841134	5.89 × 10-9	27.766
1.5	40	719.6782	2.86 × 10-9	154.594	716.610932	1.77 × 10-9	63.263
Instance	n	(He et al., 2013)			(Liu et al., 2015)			Our algorithm
		f 1	f 2	time(s)	f 1	f 2	time(s)	f 1	f 2	time(s)
1.1	5	120.710678	2.11 × 10-10	0.046	120.710678	3.49 × 10-10	-	120.710618	0	1.925
1.2	9	1.000000	2.00 × 10-12	0.172	1.000000	9.92 × 10-10	-	1.000000	0	2.365
1.3	5	22.246203	3.44 × 10-10	0.15	22.246203	9.22 × 10-10	-	23.624138	0	1.812
1.4	7	31.841134	1.80 × 10-10	98.031	31.841133	1.13 × 10-9	-	31.905109	0	110.912
1.5	40	709.812500	1.85 × 10-10	434.734	715.978266	3.15 × 10-7	-	745.118349	0	829.156

Table 3
Numerical results for the second set of instances

Instance	n	(Huang & Chen, 2006)			(Liu et al., 2010)			(Liu, 2011)
		f 1	f 2	time(s)	f 1	f 2	time(s)	f 1	f 2	time(s)
2.1	7	60.00	3.6 × 10-3	7.78	60.0000	1.59 × 10-9	1.031	60.000000	7.38 × 10-9	0.206
2.2	12	215.470	9.5 × 10-3	44.78	215.470054	2.74 × 10-7	11.562	215.470054	3.77 × 10-7	3.042
2.3	15	39.780	7.6 × 10-3	91.92	39.123400	2.97 × 10-9	130.362	39.065105	4.13 × 10-9	365.797
2.4	17	49.720	5.1 × 10-3	157.92	49.507800	6.83 × 10-8	123.407	49.368194	2.50 × 10-8	597.719
2.5	37	135.176	6.7 × 10-3	18.29	135.175410	1.23 × 10-8	0.203	135.175410	4.91 × 10-8	0.426
2.6	50	159.570	8.0 × 10-3	348.97	158.967800	4.81 × 10-9	453.625	158.967698	3.49 × 10-9	1078.386
Instance	n	(He et al., 2013)			(Liu et al., 2015)			Our algorithm
		f 1	f 2	time(s)	f 1	f 2	time(s)	f 1	f 2	time(s)
2.1	7	60.000000	2.66 × 10-9	0.015	60.000000	5.13 × 10-9	-	60.000000	0	0.118
2.2	12	215.470054	2.92 × 10-9	0.391	215.470054	1.69 × 10-7	-	215.470054	0	42.125
2.3	15	39.037635	7.78 × 10-10	36.343	39.065105	1.12 × 10-8	-	40.090538	0	167.248
2.4	17	49.338750	1.51 × 10-9	10.562	49.368194	1.42 × 10-8	-	50.246873	0	325.080
2.5	37	135.175410	5.84 × 10-9	0.234	135.175410	3.35 × 10-9	-	145.624231	0	0.402
2.6	50	158.963672	5.14 × 10-9	268.078	158.967698	3.65 × 10-9	-	163.206698	0	0.309

Table 4
Numerical results for the third set of instances

Instance	n	(Xiao et al., 2007) (SA)			(Xiao et al., 2007) (PSO)
		f 1	f 2	time(s)	f 1	f 2	time(s)
3.1	10	69.96	0	3237	59.93	0	2898
3.2	15	68.77	0	8320	67.65	0	8659
3.3	20	83.09	0	18431	83.06	0	20035
3.4	25	83.97	0	34032	84.24	0	36815
3.5	30	99.58	0	54565	99.89	0	62360
3.6	35	102.86	0	76760	102.71	0	86537
3.7	40	115.15	0	128112	115.58	0	122390
3.8	45	120.63	0	1674848	119.67	0	153006
3.9	50	125.82	0	198071	126.19	0	199050
3.10	55	138.22	0	198071	138.89	0	244171
Instance	n	(Liu et al., 2015)			Our algorithm
		f 1	f 2	time(s)	f 1	f 2	time(s)
3.1	10	59.92	5.40 × 10-9	-	59.96	0	28
3.2	15	67.39	6.44 × 10-9	-	67.59	0	39
3.3	20	82.99	2.48 × 10-9	-	83.85	0	82
3.4	25	82.98	5.10 × 10-9	-	83.68	0	107
3.5	30	98.97	4.14 × 10-9	-	98.95	0	134
3.6	35	102.32	4.50 × 10-9	-	105.15	0	155
3.7	40	115.00	5.82 × 10-9	-	114.15	0	96
3.8	45	119.07	7.64 × 10-9	-	117.69	0	118
3.9	50	124.98	4.45 × 10-9	-	124.32	0	136
3.10	55	135.13	2.19 × 10-9	-	137.37	0	154

Table 2 shows that our approach proved to be competitive. We obtained the best value for f ₁ on instance 1.1, tied in instance 1.2, while obtained results 6.1%, 0.2% and 4.0% worse than the best result on instances 1.3, 1.4 and 1.5 respectively.

Similar results were obtained in the second set test, as can be seen in Table 3. Our algorithm tied in instance 2.1 and 2.2, while obtained results at most 7.7% worse than the best result on instances 2.3, 2.4, 2.5 and 2.6.

In Table 4 we compare our approach with three other algorithms. The first set of data are from a version of simulated annealing (SA), the second set of data are from the same reference, but one of them is a version of particle swarm optimization (PSO), while the third set of data is a heuristic based on energy landscape paving. Again, our approach proved to be competitive. In relation to the envelopment radius we obtained better results in 5 out of the 10 instances. We only obtained worse results in five cases, but they were on average approximately 0.63% worse than the results from literature for such instances.

Overall, running time obtained by our algorithm can be considered good. Since the center of the rotating circular container is shifted to the center of mass of the system we always have f ₂ = 0, making our solutions more interesting than the others for this first set of instances.

Figure 5 illustrates a typical solution obtained by our algorithm for an instance of 45 circles. Note that the large Border C ₃ have 27 circles, i.e., 60% of the size, and when we carefully read the CMPT routine, we can see that the initial Border C ₁(|C ₁| = 4) is iteratively transformed in the Border C ₂(|C ₂| = 12), and finally the latter is iteratively transformed in the Border C ₃. In this example there were only two inclusions by Internal Placement.

The computational results show that the proposed algorithm is an effective method for solving the circular packing problem with additional balance constraints.

5 CONCLUSIONS

We have presented a new heuristic called center-of-mass-based placing technique for packing unequal circles into a 2D circular container with additional balance constraints. The main feature of our algorithm is the use of the Euclidean plane with origin in the center of mass of the system to select a new circle to be placed inside the container. We evaluate our approach on a series of instances from the literature and compare with existing algorithms. The computational results show that our approach is competitive and outperforms some published methods for solving this problem. We conclude that our approach is simple, but with high performance. Future work will focus on the problem of packing spheres.

ACKNOWLEDGEMENTS

The authors are indebted to the anonymous reviewers for their helpful comments. The first author wishes to thank CAPES and FAEPEX-UNICAMP (grant 285/15), the second author is grateful to FAPESP, the third author thanks CNPq.

REFERENCES

• ¹
  COXETER HSM. 1968. The problem of Apollonius. The American Mathematical Monthly, 75(1): 5-15.
• ²
  HE K, MO D, YE T & HUANG W. 2013. A coarse-to-fine quasi-physical optimization method for solving the circle packing problem with equilibrium constraints. Computers & Industrial Engineering, 66(4): 1049-1060.
• ³
  HUANG W & CHEN M. 2006. Note on: An improved algorithm for the packing of unequal circles within a larger containing circle. Computers & Industrial Engineering, 50(3): 338-344.
• ⁴
  LEI K. 2009. Constrained layout optimization based on adaptive particle swarm optimizer. Springer-Verlag Berlin Heidelberg, 1: 434-442.
• ⁵
  LENSTRA JK & RINNOOY KAN AHG. 1979. Complexity of packing, covering and partitioning problems. In: SCHRIJVER A. (Ed.), Packing and Covering in Combinatorics. Econometric Institute. Mathematisch Centrum, Amsterdam, pp. 275-291.
• ⁶
  LI N, LIU F & SUN D. 2004. A study on the particle swarm optimization with mutation operator constrained layout optimization. Chinese Journal of Computers (in Chinese), 27(7): 897-903, apud in XIAO RB, XU YC & AMOS M (2007).
• ⁷
  LIU J, JIANG Y, LI G, XUE Y, LIU Z & ZHANG Z. 2015. Heuristic-based energy landscape paving for the circular packing problem with performance constraints of equilibrium. Physica A: Statistical Mechanics and its Applications, 431: 166-174.
• ⁸
  LIU J & LI G. 2010. Basin filling algorithm for the circular packing problem with equilibrium behavioral constraints. Science China Information Sciences, 53(5): 885-895.
• ⁹
  LIU J, LI G, CHEN D, LIU W & WANG Y. 2010. Two-dimensional equilibrium constraint layout using simulated annealing. Computers & Industrial Engineering, 59(4): 530-536.
• ¹⁰
  LIU J, LI G & GENG H. 2011. A new heuristic algorithm for the circular packing problem with equilibrium constraints. Science China Information Sciences, 54(8): 1572-1584.
• ¹¹
  QIAN Z, TENG H & SUN Z. 2001. Human-computer interactive genetic algorithm and its application to constrained layout optimization. Chinese Journal of Computers, 24(5): 553-560, apud in XIAO RB, XU YC & AMOS M (2007).
• ¹²
  TANG F & TENG HF. 1999. A modified genetic algorithm and its application to layout optimization. Journal of Software (in Chinese), 10: 1096-1102, apud in XIAO RB, XU YC & AMOS M (2007).
• ¹³
  TENG HF, SUN SL, GE WH & ZHONG WX. 1994. Layout optimization for the dishes installed on rotating table - the packing problem with equilibrium behavioural constraints. Science in China Mathematics (Series A), 37(10): 1272-1280.
• ¹⁴
  WANG H, HUANG W, ZHANG Q & XU D. 2002. An improved algorithm for the packing of unequal circles within a larger containing circle. European Journal of Operational Research, 141(2): 440-453.
• ¹⁵
  XIAO RB, XU YC & AMOS M. 2007. Two hybrid compaction algorithms for the layout optimization problem. Biosystems, 90(2): 560-567.
• ¹⁶
  XU YC, XIAO RB & AMOS M. 2007. A novel genetic algorithm for the layout optimization problem. In: 'Evolutionary Computation, 2007. CEC 2007. IEEE Congress on', IEEE, pp. 3938-3943.
• ¹⁷
  ZHOU C, GAO L & GAO H. 2005. Particle swarm optimization based algorithm for constrained layout optimization. Control and Decision, 20(1): 36-40.

Publication Dates

• Publication in this collection
  May-Aug 2016

History

• Received
  20 Feb 2015
• Accepted
  12 May 2016