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 NPhard 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 twodimensional (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 threedimensional (3D) case (the equipment must rotate around its own axis) is simplified: different twodimensional circles (see, Figure 1(c)) represent threedimensional cylindrical objects to be placed inside the circular container.
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 NPhard (Lenstra & Rinnooy Kan, 19795 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. 275291.). This problem was first proposed by (Teng et al., 199413 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): 12721280.), where a mathematical model and a series of intuitive algorithms combining the method of constructing the initial objects topomodels with the modelchanging iteration method are described, and the validity of the proposed algorithms is verified by numerical examples. (Tang & Teng, 199912 TANG F & TENG HF. 1999. A modified genetic algorithm and its application to layout optimization. Journal of Software (in Chinese), 10: 10961102, apud in XIAO RB, XU YC & AMOS M (2007).) presented a modified genetic algorithm called decimal coded adaptive genetic algorithm to solve the LOP. (Xu et al., 200716 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. 39383943.) developed a version of genetic algorithm called orderbased positioning technique, which finds the best ordering for placing the circles in the container, and compare it with two existing natureinspired methods. (Qian et al., 200111 QIAN Z, TENG H & SUN Z. 2001. Humancomputer interactive genetic algorithm and its application to constrained layout optimization. Chinese Journal of Computers, 24(5): 553560, apud in XIAO RB, XU YC & AMOS M (2007).) extended the work (Tang & Teng, 199912 TANG F & TENG HF. 1999. A modified genetic algorithm and its application to layout optimization. Journal of Software (in Chinese), 10: 10961102, apud in XIAO RB, XU YC & AMOS M (2007).) by introducing a genetic algorithm based on humancomputer 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., 20046 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): 897903, apud in XIAO RB, XU YC & AMOS M (2007).) developed a PSO method basedmutation operator. This approach can escape from the local minima, maintaining the characteristic of fast speed of convergence. (Zhou et al., 200517 ZHOU C, GAO L & GAO H. 2005. Particle swarm optimization based algorithm for constrained layout optimization. Control and Decision, 20(1): 3640.) 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., 200715 XIAO RB, XU YC & AMOS M. 2007. Two hybrid compaction algorithms for the layout optimization problem. Biosystems, 90(2): 560567.) presented two natureinspired approaches based on gradient search, the first hybrid with simulated annealing (SA) method and the second hybrid with PSO method. (Lei, 20094 LEI K. 2009. Constrained layout optimization based on adaptive particle swarm optimizer. SpringerVerlag Berlin Heidelberg, 1: 434442.) presented an adaptive PSO with a better search performance, which employs multiadaptive strategies to plan largescale space global search and refined local search to obtain global optimum.
(Huang & Chen, 20063 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): 338344.) proposed an improved version of the quasiphysical quasihuman algorithm proposed by (Wang et al., 200214 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): 440453.) 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, 20109 LIU J, LI G, CHEN D, LIU W & WANG Y. 2010. Twodimensional equilibrium constraint layout using simulated annealing. Computers & Industrial Engineering, 59(4): 530536.) 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., 20109 LIU J, LI G, CHEN D, LIU W & WANG Y. 2010. Twodimensional equilibrium constraint layout using simulated annealing. Computers & Industrial Engineering, 59(4): 530536.) 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., 201110 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): 15721584.) 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., 20132 HE K, MO D, YE T & HUANG W. 2013. A coarsetofine quasiphysical optimization method for solving the circle packing problem with equilibrium constraints. Computers & Industrial Engineering, 66(4): 10491060.) proposed a hybrid approach based on coarsetofine quasiphysical optimization method, where improvement is made by adapting the quasiphysical 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., 20157 LIU J, JIANG Y, LI G, XUE Y, LIU Z & ZHANG Z. 2015. Heuristicbased energy landscape paving for the circular packing problem with performance constraints of equilibrium. Physica A: Statistical Mechanics and its Applications, 431: 166174.) presented a heuristic based on energy landscape paving. The LOP is converted into the unconstrained optimization problem by using quasiphysical 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 centerofmassbased 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, 20063 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): 338344.), (Xiao et al., 200715 XIAO RB, XU YC & AMOS M. 2007. Two hybrid compaction algorithms for the layout optimization problem. Biosystems, 90(2): 560567.), (Lei, 20094 LEI K. 2009. Constrained layout optimization based on adaptive particle swarm optimizer. SpringerVerlag Berlin Heidelberg, 1: 434442.), (Liu & Li, 20109 LIU J, LI G, CHEN D, LIU W & WANG Y. 2010. Twodimensional equilibrium constraint layout using simulated annealing. Computers & Industrial Engineering, 59(4): 530536.), (Liu et al., 20109 LIU J, LI G, CHEN D, LIU W & WANG Y. 2010. Twodimensional equilibrium constraint layout using simulated annealing. Computers & Industrial Engineering, 59(4): 530536., ^{2011}10 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): 15721584., ^{2015}7 LIU J, JIANG Y, LI G, XUE Y, LIU Z & ZHANG Z. 2015. Heuristicbased energy landscape paving for the circular packing problem with performance constraints of equilibrium. Physica A: Statistical Mechanics and its Applications, 431: 166174.) and (He et al., 20132 HE K, MO D, YE T & HUANG W. 2013. A coarsetofine quasiphysical optimization method for solving the circle packing problem with equilibrium constraints. Computers & Industrial Engineering, 66(4): 10491060.). 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_{1} (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 _{1}, y _{1}, x _{2}, y _{2}, ..., x_{n} , y_{n} )^{T} that satisfies the following mathematical formulation.
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 _{3} 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).
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 _{1} i _{2}...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 _{1}, q _{1}, p _{2}, q _{2} ∈ {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 _{1} i _{2}...i_{p} _{ 1} i_{p}i_{p} _{+ 1}...i_{t} _{ 1} i_{t} , we say that the circles i _{1} i _{2}...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 = 1324 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, 19681 COXETER HSM. 1968. The problem of Apollonius. The American Mathematical Monthly, 75(1): 515.).
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
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.
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 CENTEROFMASSBASED 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.
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 = 1324. Among many optional positions we can choose X(α(2)), for example, at the xaxis 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.
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).
Next we describe the two cases of inclusion in the following 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.
Postoptimization
A postoptimization 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 postoptimization 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 _{1} i _{2}...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.
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. (Postoptimization) 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 postoptimization 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 postoptimization process. Recall that, before postoptimization, 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 postoptimization 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 postoptimization is bounded by 2 × (n  1) × (C  2).
Therefore the complexity of placing n circles during the Main routine is bounded by O(n ^{2}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 ^{3}).
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 natureinspired approaches based on simulated annealing and particle swarm optimization (Xiao et al., 200715 XIAO RB, XU YC & AMOS M. 2007. Two hybrid compaction algorithms for the layout optimization problem. Biosystems, 90(2): 560567.; Lei, 20094 LEI K. 2009. Constrained layout optimization based on adaptive particle swarm optimizer. SpringerVerlag Berlin Heidelberg, 1: 434442.), a hybrid approach based on simulated annealing, neighborhood search mechanism and the adaptive gradientmethod (Liu et al., 20109 LIU J, LI G, CHEN D, LIU W & WANG Y. 2010. Twodimensional equilibrium constraint layout using simulated annealing. Computers & Industrial Engineering, 59(4): 530536.), a hybrid tabu search algorithm and gradient method (Liu et al., 201110 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): 15721584.), a series of heuristics based on energy landscape paving, gradient method and local search (Liu & Li, 20108 LIU J & LI G. 2010. Basin filling algorithm for the circular packing problem with equilibrium behavioral constraints. Science China Information Sciences, 53(5): 885895.; Liu et al. 20109 LIU J, LI G, CHEN D, LIU W & WANG Y. 2010. Twodimensional equilibrium constraint layout using simulated annealing. Computers & Industrial Engineering, 59(4): 530536., ^{2015}7 LIU J, JIANG Y, LI G, XUE Y, LIU Z & ZHANG Z. 2015. Heuristicbased energy landscape paving for the circular packing problem with performance constraints of equilibrium. Physica A: Statistical Mechanics and its Applications, 431: 166174.), and a series of algorithms based on quasiphysical approaches, gradient method and local search (Huang & Chen, 20063 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): 338344.; He et al., 20132 HE K, MO D, YE T & HUANG W. 2013. A coarsetofine quasiphysical optimization method for solving the circle packing problem with equilibrium constraints. Computers & Industrial Engineering, 66(4): 10491060.; Liu et al., 20157 LIU J, JIANG Y, LI G, XUE Y, LIU Z & ZHANG Z. 2015. Heuristicbased energy landscape paving for the circular packing problem with performance constraints of equilibrium. Physica A: Statistical Mechanics and its Applications, 431: 166174.).
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, 20063 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): 338344.), (Liu & Li, 20109 LIU J, LI G, CHEN D, LIU W & WANG Y. 2010. Twodimensional equilibrium constraint layout using simulated annealing. Computers & Industrial Engineering, 59(4): 530536.) and (Xiao et al., 200715 XIAO RB, XU YC & AMOS M. 2007. Two hybrid compaction algorithms for the layout optimization problem. Biosystems, 90(2): 560567.).
Data of each instance. First set of instances (Liu & Li, 20109 LIU J, LI G, CHEN D, LIU W & WANG Y. 2010. Twodimensional equilibrium constraint layout using simulated annealing. Computers & Industrial Engineering, 59(4): 530536.).
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 _{1}), the imbalance obtained (second objective function f _{2}), and the running time t (in seconds).
Table 2 shows that our approach proved to be competitive. We obtained the best value for f _{1} 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 _{2} = 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 _{3} 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 _{1}(C _{1} = 4) is iteratively transformed in the Border C _{2}(C _{2} = 12), and finally the latter is iteratively transformed in the Border C _{3}. 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 centerofmassbased 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 FAEPEXUNICAMP (grant 285/15), the second author is grateful to FAPESP, the third author thanks CNPq.
REFERENCES

^{1}COXETER HSM. 1968. The problem of Apollonius. The American Mathematical Monthly, 75(1): 515.

^{2}HE K, MO D, YE T & HUANG W. 2013. A coarsetofine quasiphysical optimization method for solving the circle packing problem with equilibrium constraints. Computers & Industrial Engineering, 66(4): 10491060.

^{3}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): 338344.

^{4}LEI K. 2009. Constrained layout optimization based on adaptive particle swarm optimizer. SpringerVerlag Berlin Heidelberg, 1: 434442.

^{5}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. 275291.

^{6}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): 897903, apud in XIAO RB, XU YC & AMOS M (2007).

^{7}LIU J, JIANG Y, LI G, XUE Y, LIU Z & ZHANG Z. 2015. Heuristicbased energy landscape paving for the circular packing problem with performance constraints of equilibrium. Physica A: Statistical Mechanics and its Applications, 431: 166174.

^{8}LIU J & LI G. 2010. Basin filling algorithm for the circular packing problem with equilibrium behavioral constraints. Science China Information Sciences, 53(5): 885895.

^{9}LIU J, LI G, CHEN D, LIU W & WANG Y. 2010. Twodimensional equilibrium constraint layout using simulated annealing. Computers & Industrial Engineering, 59(4): 530536.

^{10}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): 15721584.

^{11}QIAN Z, TENG H & SUN Z. 2001. Humancomputer interactive genetic algorithm and its application to constrained layout optimization. Chinese Journal of Computers, 24(5): 553560, apud in XIAO RB, XU YC & AMOS M (2007).

^{12}TANG F & TENG HF. 1999. A modified genetic algorithm and its application to layout optimization. Journal of Software (in Chinese), 10: 10961102, apud in XIAO RB, XU YC & AMOS M (2007).

^{13}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): 12721280.

^{14}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): 440453.

^{15}XIAO RB, XU YC & AMOS M. 2007. Two hybrid compaction algorithms for the layout optimization problem. Biosystems, 90(2): 560567.

^{16}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. 39383943.

^{17}ZHOU C, GAO L & GAO H. 2005. Particle swarm optimization based algorithm for constrained layout optimization. Control and Decision, 20(1): 3640.
Publication Dates

Publication in this collection
MayAug 2016
History

Received
20 Feb 2015 
Accepted
12 May 2016