An Optimized Covering Spheroids by Spheres

: Covering spheroids (ellipsoids of revolution) by di ﬀ erent spheres is studied. The research is motivated by packing non-spherical particles arising in natural sciences, e.g., in powder technologies. The concept of an ε -cover is introduced as an outer multi-spherical approximation of the spheroid with the proximity ε . A fast heuristic algorithm is proposed to construct an optimized ε -cover giving a reasonable balance between the value of the proximity parameter ε and the number of spheres used. Computational results are provided to demonstrate the e ﬃ ciency of the approach.


Introduction
Covering a certain region by simple shapes has various applications. Our interest in covering problems is motivated by packing particles arising in science and engineering applications. Packing problems are widely used in modeling liquid and glass structures [1,2], representing granular materials [3], packing beds, and cermet [4], as well as in many other applications (see, e.g., [5][6][7][8][9][10]). Many algorithms have been proposed for packing spherical particles (see, e.g., [11][12][13][14][15][16] and the references therein). However, to have a better and more adequate representation of particulate microstructure, packing non-spherical particles must be considered. The approaches proposed to handle non-spherical shapes can be divided in three large groups [17]. Techniques in the first group use analytical equations for the shapes of the particles, e.g., ellipsoids, and the main modeling problem is to state analytically non-overlapping and containment conditions [18][19][20][21][22][23]. The second approach is based on tessellating the container/particles shapes with a grid and then approximating them by corresponding collections of grid nodes (see, e.g., [15,24] and the references therein). This way, detecting the overlap is reduced to verify if two shapes share the same node. However, to get a reasonable approximation, fine grids must be used resulting in large-scale and memory-consuming problems. In the third approach, the shape of the particle is represented approximately by a collection of spheres having different sizes and positions (see, e.g., [25][26][27]). Then, detecting the overlap is reduced to verifying the overlapping for two spheres from different particle collections. This problem is much simpler than detecting overlapping in the first approach. The third approach can be considered as a compromise between the simplicity of the tessellating techniques and the rigor of the methods based on analytical shapes presentation. The efficiency of the third approach depends on the quality of the multi-spherical approximation of the shape of the particle. The number of the spheres representing the shape of the particle must be small enough to simplify the non-overlapping (pairwise) test. At the same time, the union of the spheres in the collection must represent the original shape of the particle closely enough.
The paper is organized as follows. Section 2 presents the basic definitions and formulations used in the paper as well as the two-stage conceptual approach. Solution techniques are presented in detail in Section 3. Numerical results are presented in Section 4, while Section 5 presents the conclusions. Expressions for the parameters used in Section 3 are derived in the Appendix A.

Basic Constructions
The following definitions are used throughout the paper. Let be a given spheroid (ellipsoid of revolution) with c = b.

Definition 1.
A set Λ ⊂ R 3 is called a cover (outer approximation) for the spheroid E if Λ ⊇ E.
The following extended spheroid, referred to as an ε-spheroid, can be considered as an outer approximation of the original spheroid E with a given error ε ≥ 0.
where ε is an error of the cover.
In this paper, the ε-cover of the spheroid E by spheres S k , k = 1, . . . , M is studied.
From a practical point of view, a "good" ε-cover Λ(ε, M) must use a few circles M and provide a small approximation error ε. However, we may expect that for sufficiently small ε, decreasing ε in Λ(ε, M) results in increasing M, and vice versa. Thus, we cannot minimize ε and M simultaneously. Instead, our problem is formulated as follows: Find Λ(ε, M), i.e., positions, radii and total number M of spheres, providing a reasonable balance between ε and M.
The following heuristic two-stage approach is proposed to construct the optimized ε-cover. Stage 1. For a given ε, find the minimum number M * of the spheres covering the given spheroid E; i.e., find Λ(ε, M * ). Here, the decision variables are the positions and radii of the circles and their total number M.

Stage 2.
For M * obtained at Stage 1, find a cover of the given spheroid E with the minimal value ε * ≤ ε; i.e., find Λ(ε * , M * ). Here, the decision variables are the positions and radii of the circles and the proximity parameter ε.
In what follows, constructing the cover Λ(ε * , M * ) is referred to as an Optimized Spherical Covering (OSC). Detailed descriptions of problems arising at Stages 1,2 and their solutions are presented in the next section.

Solution Algorithm
Let an ellipse E with semi-axes a and b be the projection of the spheroid E on the plane XOY. It is assumed that the centers v k = (x k , y k , z k ) of the spheres S k of variable radii r k belong to the axis OX, i.e., v k = (x k , 0, 0).
The parametric description of the ellipse E presented in [28] is used: Note that the point x(s, t) = cos t a (a 2 − b 2 ) belongs to the axis OX for s = 0. Denote by Λ(ε, M) the cover that is the projection of Λ(ε, M) on the plane XOY. By the definition, the cover Λ(ε, M) is an outer approximation of the ellipse E, i.e., Λ(ε, M) = M ∪ k=1 C k (x k , 0). Suppose that the odd number of the circles is used for the cover, and the first circle has its radius b + ε and is centered at the point (0, 0, 0), as shown at Figure 1a. Then, due to the symmetry, the ellipse E can be covered by a collection of M = 2m + 1 circles. Here, the right-hand side of the E can be covered by the circles C k (x k , 0), k = 0, 1, . . . , 2m with x k+1 > x k , x k ≥ 0, and similarly the left-hand side.
For the case of an even number of circles, the ellipse E be covered by a collection of 2m circles C k (x k , 0), k = 1, . . . , 2m. Let the first circle be centered at the point (x 0 , 0, 0), as shown at Figure 1b. To find the parameters x 0 , r 0 of the first circle for 0 ≤ t ≤ π 2 , the following optimization problem is solved: For the given t 0 corresponding parameters x 0 , r 0 are defined in the form of Using the parametric description of the ellipse [28], the parameter x 0 can be obtained analytically as follows (we would like to thank the anonymous referee for pointing out the possibility of getting the explicit solution): The details are presented in the Appendix A.
Then, due to the symmetry, the right-hand side of the ellipse E can be covered by a collection of m circles C k (x k , 0), k = 1, . . . , m, for x k+1 > x k , x k > 0 and similarly for the left-hand side.
In what follows, the solution algorithm for the OSC problem is presented only for the odd number of spheres, since the difference between the odd and even cases is only in the parameters of the first sphere.
Then, due to the symmetry, the right-hand side of the ellipse E can be covered by a collection In what follows, the solution algorithm for the OSC problem is presented only for the odd number of spheres, since the difference between the odd and even cases is only in the parameters of the first sphere.

Solution Algorithm for Stage 1
Let the value of 0  be given. At Stage 1, the objective is to find the minimum number of circles Step 2. Get the intersection point of the ellipse E and the circle k C solving the following system: Using the parametric description of the ellipse [28], the point is presented analytically by the following formula (see the appendix for details): Step 3. Set 1 kk .
Step 4. Solve the following nonlinear programming problem:

Solution Algorithm for Stage 1
Let the value of ε > 0 be given. At Stage 1, the objective is to find the minimum number of circles C k (x k , 0), k = 0, 1, . . . , 2m covering the ellipse E centered at (0, 0) for the given value of ε.
Step 2. Get the intersection point p k = (p kx , p ky ) of the ellipse E and the circle C k solving the following system: where r k is the radius of the circle C k centered at the point (x k , 0). Using the parametric description of the ellipse [28], the point p k = (p kx , p ky ) is presented analytically by the following formula (see the Appendix A for details): Step 3. Set k = k + 1.
Step 4. Solve the following nonlinear programming problem: and get its optimal solution (r * k , t * k ). Here, the feasible set V is defined by the following system of inequalities: where inequality (5) assures that the point p k−1 is inside the circle C k of radius r k centered at (x k , 0); inequality (6) guarantees that the circle C k of radius r k centered at (x k , 0) is inside the ellipse E ε with semi-axes a ε = a + ε, b ε = b + ε; inequality (7) reflects the monotonous decrease of the corresponding angle t k ; and the last inequality describes the natural constraint for the radius of the circle C k . Find . Using the parametric description of the ellipse [28], the parameters x k and r k are presented analytically in the following form (see the Appendix A for details): Step 5. If the point (a, 0) is inside the circle C k of radius r k centered at (x k , 0), then the solution is obtained and the algorithm stops. Otherwise, go to Step 2. Figures 2 and 3 illustrate an iterative procedure for the algorithm described above.
( , 0) Using the parametric description of the ellipse [28], the parameters k x and k r are presented analytically in the following form (see the appendix for details):

Solution Algorithm for the Stage 2
At this stage, the minimum value of  is obtained for a given number 21 m  of covering circles.

Solution Algorithm for the Stage 2
At this stage, the minimum value of ε is obtained for a given number 2m + 1 of covering circles. Let E ε be the ellipse centered at (0, 0) with semi-axes a ε = a + ε, b ε = b + ε (see Figure 4).
The following nonlinear programing problem is used to optimize ε: ε.

Computational Results
Computational experiments were run on an AMD Athlon 64 X2 5200+ computer. Local optimization was performed by the IPOPT solver [29], which is available at an open access noncommercial software depository (https://projects.coin-or.org/Ipopt).
To illustrate the algorithm performance, examples for ellipses and spheroids are considered. For each problem instance, the computational time required to construct the approximated cover is less than a second.

Instance 1.
The ellipse E with semi-axes a =1.3, b =1 is given.
(a) For ε = 0.3, the solution with one covering circle is presented in Figure 5a: (b) For ε = 0.2, the solution with an even number of covering circles is shown in Figure 5b:  Figure 5c:

Computational Results
Computational experiments were run on an AMD Athlon 64 X2 5200+ computer. Local optimization was performed by the IPOPT solver [29], which is available at an open access noncommercial software depository (https://projects.coin-or.org/Ipopt).
To illustrate the algorithm performance, examples for ellipses and spheroids are considered. For each problem instance, the computational time required to construct the approximated cover is less than a second.
, the solution with one covering circle is presented in Figure 5a: , the solution with an even number of covering circles is shown in Figure 5b:
(a) The optimized cover with an even number of circles is presented in Figure 6a:

Conclusions
The problem of optimized multi-spherical covering for spheroids is introduced. This problem is motivated by the packing non-spherical particles arising in naturals sciences and engineering. The simple heuristic approach is proposed to construct an optimized covering, providing a reasonable balance between the number of spheres and the error of approximation. Computational experiments indicate that the proposed approach constructs good feasible coverings very fast: in less than a second. The multi-spherical approximations obtained in the paper provide a basis for fast optimized packing spheroids in different containers, using algorithms proposed, e.g., in [30][31][32]. The two-stage scheme presented in Section 2 can be applied for the general optimized multi-spherical covering. An extension of the covering approach to the case of more complex objects [33][34][35][36] is an interesting direction for the future research. Some results on this topic are on the way.

Acknowledgement:
The authors would like to thank anonymous referees for careful reading the paper and constructive comments. The first two authors were partially supported by the "Program for the State Priority Scientific Research and Technological (Experimental) Development of the Department of Physical and Technical Problems of Energy of the National Academy of Sciences of Ukraine" (#6541230).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
We would like to thank the anonymous referee for pointing out that the explicit expressions can be obtained for k r , k x and ( , )

Conclusions
The problem of optimized multi-spherical covering for spheroids is introduced. This problem is motivated by the packing non-spherical particles arising in naturals sciences and engineering. The simple heuristic approach is proposed to construct an optimized covering, providing a reasonable balance between the number of spheres and the error of approximation. Computational experiments indicate that the proposed approach constructs good feasible coverings very fast: in less than a second. The multi-spherical approximations obtained in the paper provide a basis for fast optimized packing spheroids in different containers, using algorithms proposed, e.g., in [30][31][32]. The two-stage scheme presented in Section 2 can be applied for the general optimized multi-spherical covering. An extension of the covering approach to the case of more complex objects [33][34][35][36] is an interesting direction for the future research. Some results on this topic are on the way.