HYBRID OPTIMIZATION FOR THE 3D COVERING PROBLEM WITH SPHERES: INTEGRATING HEURISTICS, PATTERN MINING AND BRANCH AND CUT

ABSTRACT

This paper addresses the challenge of covering three-dimensional solids with spheres of varying radii, allowing for partial overlaps, which is particularly relevant for applications such as gamma knife radiotherapy. In this work a novel approach is introduced, approximating the covering problem in the sense that it focuses on an objective function that maximizes the total volume of the used spheres, but it does not consider the volume of the overlapped regions. This approximation is further refined by discretizing the target volume and utilizing a finite set of potential sphere centers. Exact and hybrid methods are developed to solve this approximation. The exact method employs either a Branch-and-Bound or Branch-and-Cut algorithm, while the hybrid method integrates heuristic solutions and data mining techniques to identify promising sphere configurations, which are then refined using either Branch-and-Bound or Branch-and-Cut on a reduced search space. Data Envelopment Analysis (DEA) is used to determine optimal overlap parameters. Computational experiments demonstrate that the hybrid method outperforms the exact method in both coverage and computational efficiency, highlighting the efficacy of integrating heuristic, data mining, and exact techniques for solving complex optimization problems.

Keywords:
covering problem; data mining; branch and cut

1 INTRODUCTION

Covering problems are often related to determining the position of certain combinatorial objects, for example, squares or cubes, to cover the area or volume of a container. To achieve this objective, a level of intersection between the units of the covering objects is allowed. This class of problems is naturally formulated as a minimization problem. A classical application of covering a two-dimensional region with circles arises in the field of telecommunications ^{18); (15}.

The approach used in this work to tackle the problem of covering a tridimensional region with spheres of different radii combines elements of covering and packing ^{18); (22}, two problems from the literature that consider a container and combinatorial objects, such as spheres and cubes. Therefore, a brief description of the packing problem is given in order to point out features that are important to mention while defining the problem subject of this paper.

Packing problems consist of positioning the maximum number of objects within a container without overlapping ²². A classical application of the packing problem is to stack the maximum number of oranges in a box ⁴. Packing problems are still being studied, as can be seen in ¹³, where a systems-based approach to problem formulation is proposed based on elements, containers, and binary relations among them, such as precedence, dominance and element-to-container correspondence. In ²⁵ the authors present a hybrid quantum-classical framework to address realistic three-dimensional bin packing problems, incorporating practical constraints such as package and bin dimensions, weight constraints, category affinities and item ordering preferences.

The problem of covering a tridimensional container with spheres of different radii consists of positioning combinatorial objects inside a container to cover the maximum of its volume, using the few combinatorial objects as possible, while allowing partial overlap. An important application of this problem is the Gamma knife^{7); (8}, a specific treatment of brain tumors in which equipment is placed on the patient’s head and radiation shots are issued. The region of the cerebral tumor affected by a shot can mathematically be approximated by a spherical object. As the equipment has different sizes, its area of effect can be modeled as spheres of 2, 4, 7, or 9 mm radii. The process of planning where to position the spheres and choosing their radii is very hard, so the objective of studying this problem is to automate the planning process of Gamma Knife treatment, aiming to find a good covering of the tumor in an acceptable time. The mathematical problem itself is continuous, since there are infinite points in the region to be covered by the spheres. The Gamma Knife treatment also considers the dose of shots¹⁴ and this issue is not treated in this work.

Considering the difficulty of solving the covering problem, which is usually modeled as a nonlinear and nonconvex mathematical programming problem, several approaches have been proposed in the literature to obtain tractable approximations. One common strategy consists of discretizing the volume of the container to be covered, as explored in ¹⁴ and further developed in ¹⁹. Alternative formulations have also been investigated, such as adopting a packing perspective to address the covering problem. In ²⁴, the authors proposed heuristic methods to generate feasible solutions, followed by a Branch-and-Bound procedure. In the work ⁶, the authors introduced a nonlinear model focused on minimizing the radiation dose, which was solved using heuristic initialization combined with inexact penalization and least squares techniques. In ¹⁴ it is proposed a nonlinear nonconvex mixed-integer formulation along with several reformulations, addressing objectives such as minimizing the number of spheres and limiting radiation exposure to healthy tissue; these models were solved using Variable Neighborhood Search and commercial solvers. Graph-theoretical approaches were explored in ¹⁹, where maximum weight cliques were used to generate cuts within a Branch-and-Cut framework. Other nonlinear formulations were studied using penalization methods, stochastic search heuristics, and smoothing techniques ^{16); (27); (26}. Metaheuristic approaches have also been considered, including Particle Swarm Optimization for radiation therapy planning ²⁰. More recently, in ⁸ the authors proposed new mathematical formulations for covering solids with spheres of different radii, as well as hybrid solution methods combining exact and heuristic techniques.

The contributions given by this work to the Problem of Covering Solids with Spheres of Different Radii are the expansion of the linear integer formulation proposed in ⁸ and a Branch-and-Cut method. Furthermore, considering the success of hybrid methods in solving combinatorial optimization problems ^{10); (21); (8); (5); (23}, a hybrid method is proposed that combines black-box heuristic, data mining, and branch-and-bound / branch-and-cut.

The remainder of this paper is organized as follows. In Section 2, the mathematical formulation is presented. Section 3 presents the valid inequalities used in the Branch-and-Cut. In Section 4, the proposed hybrid method is explained. Section 5 presents the computational results. Finally, Section 6 presents the final remarks and future work.

2 MATHEMATICAL PROGRAMMING MODEL

Let T ⊂ ℝ³ be a compact set. The covering problem consists of positioning the smallest possible number of spheres (possibly of different radii) to cover T . A covering means that each point p ∈ T must also belong to some sphere. This problem is to be referred to as the Covering Solids with Spheres of Different Radii Problem (CSSDRP). Let B(x, r) be the sphere with center x ∈ ℝ³ and radius r ∈ ℝ₊. Mathematically, we can define the problem as follows:

(CSSDRP) Given a compact set T ⊂ ℝ³ and a finite set of different radii R ⊂ ℝ₊, the problem can be described as finding a set $\left\{B(x,r)\left|x\in {\mathrm{ℝ}}^3,r\in R\right.\right\}$ of spheres of minimum cardinality that covers T .

An intuitive way of modeling the (CSSDRP) would require one constraint for each point p ∈ T , ensuring that p also belongs to some of the spheres (and hence results in a covering). Considering that T has infinite points, such a model would have an infinite set of constraints. Following this idea, the authors of ¹⁴ proposed a nonlinear and nonconvex mathematical programming model with an infinite (and uncountable) set of variables and constraints. To avoid this difficulty, an approximate reformulation was proposed in the same work, consisting of replacing T by a finite set of points in T .

In ¹⁹, a slightly different approach was used. Each sphere has a fixed radius, instead of being a variable whose value would be determined by the model, as in ¹⁴. For this reason, the following assumption is necessary: for each r ∈ R, an upper bound in the number of spheres of radius r in the optimal solution of (CSSDRP) is known and is denoted by N _r . Then, we have $n=\sum _{r\in R}{N}_r$ available spheres to be selected to compose the minimum cardinality covering of T . The optimal location of each sphere is also to be determined. Then, for all i in {1,..., n}, consider the following variables: binary variables y _i , assuming value 1 if the i-th sphere takes part in the covering and value 0 otherwise; and real variables x ⁱ ∈ ℝ³, representing the coordinates of the center of the i-th sphere. Let r _i ∈ R be the radius of the i-th sphere, for all i in {1,..., n}. The nonlinear nonconvex model proposed in ¹⁹, is written as in (1)-(4):

$$\begin{gathered} \max \sum _{i=1}^n{c}_i{y}_i \\ \text{s. t.}\left|\right|x^i-x^j|{|}^2\ge (r_i+r_j-{\alpha }_{ij}{)}^2(y_i+y_j-1) \end{gathered}$$ (1)

${\alpha }_{ij}\ge 0,\forall 1\le i<j\le n,$ (2)

$x^i\in T,\forall i\in \{1,\dots ,n\},$ (3)

$y_i\in \{0,1\},\forall i\in \{1,\dots ,n\}.$ (4)

The model (1)-(4) is parameterized, which depends on the values of the parameters c _i ≥ 0, for all i ∈ {1,..., n}. Parameters c are associated with the cost of each sphere in the solution. The parameters α control the maximum amount of overlap allowed between pairs of spheres (see Figure 1). That is, the spheres must be located respecting

$\left|\right|x^i-x^j\left|\right|\ge r_i+r_j-{\alpha }_{ij}.$ (5)

Figure 1
The parameter α _ij limiting the amount of overlap between spheres i and j.

The family of constraints ² ensures condition ⁵ only when both spheres are selected in the solution (when both variables y _i and y _j assume value 1). When all α _ij = 0, no overlap is allowed and (1)-(4) is a model for the packing problem.

It is important to note that the optimal solution of (1)-(4) will not necessarily represent the optimal solution of (CSSDRP). As an example, consider T is a sphere of radius 8 cm. Assume there are two available spheres for the covering, one with radius 2 cm and the other with radius 8 cm. Consider α ₁₂ = 1 cm. The solution for (CSSDRP) is to use only the radius 8 sphere in the center of T . However, the optimal solution given by the proposed model will use both spheres, maximizing the objective function. However, there are choices for the values of the parameters for which the optimal solution of model (1)-(4) will correspond to the optimal solution of (CSSDRP), as stated in Theorem 1 ¹⁹, as follows.

Theorem 1: There are sets of parameters α_ij ≥ 0, for 1 ≤ i < j ≤ n, and c _i , for i ∈ {1,..., n}, for which the optimal solution of (1)-(4) is also the optimal solution of (CSSDRP).

Proof: Let (x^∗, y^∗) be an optimal solution of (CSSDRP). For each pair of spheres in the covering, calculate the distance D_ij between their centers, given by

$D_{ij}=\Vert (x^i{)}^{*}-(x^j{)}^{*}\Vert .$

If D _ij > r _i + r _j , let a _ij = 0. Otherwise, let a _ij = r _i + r _j − D _ij . If sphere k is not used in the optimal covering, i.e. $y_k^{*}=0$, let the parameters α related to sphere k assume any value.

Let, for all i ∈ N,

$c_i=\left\{\begin{array}{rl}1,& \text{if}{y}_i^{*}=1\\ -1,& \text{if}{y}_i^{*}=0.\end{array}\right.$

Then (x ^∗ , y ^∗) is also an optimal solution for model (1)-(4), using the values chosen for the parameters α and c. In fact, it is easy to check that it is a feasible solution. We are left to show that it has the highest objective function value. We have

$\sum _{i\in N}{c}_i{y}_i=\sum _{i\in N|y_i^{*}=1}{y}_i+\sum _{i\in N|y_i^{*}=0}-y_i.$

This sum will reach its highest value if we take y _i = 1 whenever $y_i^{*}=1$ and take y _i = 0 otherwise. □

However, there is no information about the adequate values of the parameters. Theorem 1 states only their existence. This work focuses on solving the model (1)-(4), regardless of these parameter values. It is a non-linear and non-convex model and its solution is hard to obtain by current commercial solvers. Possible values for α _ij are studied using Data Envelopment Analysis.

In ¹⁹, model (1)-(4) is revisited from a graph-theory perspective, and this new formulation will be used for later comparisons. The fundamental concepts of this approach will be introduced.

2.1 Discretization Procedure

The discretization procedure takes as input the compact set T ⊂ ℝ³ and a real parameter δ > 0. This procedure aims at eliminating the nonconvexity of Constraints (2) present in the model described in (1) - (4). The points in the intersection of T with a grid of size δ will be used as the discretization of T . This procedure is represented in Figure 2.

Figure 2
In the left, the set T and the δ -spaced grid of gray points. In the right, the discretization of T , corresponding to the black dots.

2.2 Mathematical Formulation

In this section, the new mathematical formulation is introduced, based on the formulation of (9).

The discretization procedure generates a finite set of p points. Once the coordinates of these points are known, the distance between each pair of points can be calculated in advance. Proceeding in this way, the left side of constraints (2) becomes a constant. Let P = {1,..., p}. Define

$D_{ij}=\left|\right|p_i-p_j\left|\right|,i,j\in P.$ (6)

According to the above ideas, the following reformulation is proposed.

$\max \sum _{i\in P}\sum _{r\in R}{c}_i^r{y}_i^r$ (7)

$\text{s. t.}{D}_{ij}\ge (r+s-{\alpha }_{rs})(y_i^r+y_j^s-1),\forall i,j\in P,\forall r,s\in R,$ (8)

$\sum _{r\in R}{y}_i^r\le 1,\forall i\in P,$ (9)

$y_i^r\in \{0,1\},\forall i\in P,\forall r\in R.$ (10)

where $y_i^r$ is a binary variable that indicates whether or not a sphere centered on i ∈ P and of radius r ∈ R is in the solution and $c_i^r$ is the cost of such a sphere. The family of inequalities (8) limits the amount of overlap among each pair of spheres. The restrictions (9) ensure that for each point i of the discretization, only one radius is selected, or none. The last set of constraints, the set (10), describes the domain of variables y.

Since D _ij and (r + s − α _rs ) are positive constants, for all i, j ∈ P and for all r, s ∈ R, the restrictions (8) can be strengthened, tightening the linear relaxation bound. Isolating $y_i^r$ and $y_j^s$ and using the fact that these variables are binary, constraints (8) can be written as

$$\begin{gathered} y_i^r+y_j^s\le ⌊\frac{D_{ij}}{(r+s-{\alpha }_{rs})}⌋+1, \\ \forall i,j\in P,\forall r,s\in R|i\ne j. \end{gathered}$$ (11)

Constraints (11) become redundant when D _ij ≥ (r + s − α _rs ), once, in this case, the right side of the inequalities is at least two. Meanwhile, the constraints (11) are meaningful when D _ij < (r + s − α _rs ), causing the first term on the right side of the inequalities to equal zero. Thus, the constraint (11) can be rewritten as

$$\begin{gathered} y_i^r+y_j^s\le 1,\forall i,j\in P,\forall r,s\in R \\ \text{such that}{D}_{ij}<(r+s-{\alpha }_{rs})\wedge i\ne j. \end{gathered}$$ (12)

Once the modifications were made, the model becomes the following:

$$\begin{gathered} \max \sum _{i\in P}\sum _{r\in R}{c}_i^r{y}_i^r \\ \text{s.t.}{y}_i^r+y_j^s\le 1,\forall i,j\in P,\forall r,s\in R|D_{ij}<(r+s-{\alpha }_{rs})\wedge i\ne j \end{gathered}$$ (13)

$y_i^r\in \{0,1\},\forall i\in P,\forall r\in R.$ (14)

It is important to point out that in Constraint (12), i ≠ j, otherwise this family of constraints will refer to the same sphere. This model can still be strengthened by adding constraints (15).

$$\begin{gathered} \sum _{r^{\text{'}}\in R|r^{\text{'}}\ge r}{y}_i^{r^{\text{'}}}+\sum _{j\in P\setminus \left\{i\right\} \\ s\in R|s\ge r \\ {D}_{ij}<(r+s-\alpha \text{'})/2 \\ {\alpha }_{rs}\ge \alpha \text{'}}{y}_j^s\le 1\forall i\in P,\forall r\in R, \end{gathered}$$ (15)

where, $mi{n}_{s\ge r}{\alpha }_{rs}\le \alpha \text{'}\le ma{x}_{s\ge r}{\alpha }_{rs}$, forall s ∈ R | s ≥ r.

Analyzing the first part of the Constraints (15), suppose there is a sphere centered at i ∈ P with radius r ∈ R that cannot be at the solution together with a sphere centered at j ∈ P with radius s ∈ R because they violate Constraints (8), which means these two spheres exceed the allowed intersection. Now, suppose the two centers i and j ∈ P are the same as well as the radius s ∈ R, but sphere i has radius r’ ≥ r. Clearly, if the new sphere centered at i has radius r’ ≥ r, its intersection with sphere centered at j will be equal or larger the former one, again violating Constraints (8). This result can be considered for all spheres with r’ ≥ r radius and is represented by the first summation in Constraints (15).

Now, assume that there is one sphere centered at i ∈ P with radius r ∈ R. Suppose, for example, that j and j’, are taken as centers of spheres with radius, respectively, s, s’ ≥ r at distances D _ij < (r + s − α _rs )/2 and D _ij ′ < (r + s’ − α _rs ′)/2. These two sphere centers are shown to be infeasible as a solution, unless α _rs ≤ α’ is considered for all r, s ∈ R, which is the case in this paper since α _rs = β . This occurs because even if the centers j and j’ are diametrically opposite, their distance D _jj ′ would be strictly less than (s + s’ − α _ss ′ ), violating the constraints (8). Thus, it occurs to all points that are at a distance D _ij < (r + s − α _rs )/2 around the sphere centered at i ∈ P with radius r ∈ R and this result is stated in the second summation in Constraints (15).

The formulation that follows aggregates the Constraints (15) and is referred to as Enhanced Formulation from now on throughout this text.

$$\begin{gathered} \max \sum _{i\in P}\sum _{r\in R}{c}_i^r{y}_i^r \\ \text{s.t.}\left(10\right),\left(12\right),\left(15\right) \end{gathered}$$

3 VALID INEQUALITIES

This section presents inequalities that can be applied at each Branch-and-Bound iteration, resulting in a Branch-and-Cut. Considering that several constraints of the integer programming formulation are a case of clique cuts ²⁸, following them as a base for the Branch-and-Cut seems rather natural. Hereafter, the construction of the prohibition graph is introduced to explain how the inequalities arose and how their joint work with Branch-and-Bound turned it into a Branch-and-Cut thus accelerating the process of determining the centers of the spheres.

Let the grid of points resulting from the discretization procedure described in Section 2.1 correspond to a graph whose vertices are divided into four different vertices, one for each possible radius, and initially are totally disconnected. An edge is added among two vertices when their radii do not satisfy the distance constraint D _ij ≥ (r + s − α _rs ), which means that, related to these radii, $y_i^r$ and $y_j^s$ are not feasible together at the solution of the enhanced formulation. Enlarging the search for variables that are not feasible when considered together at the solution and taking into account the weights of the vertices, if a clique is identified on the graph, only one of its vertices may be at the solution, and clique cuts model that situation. Due to the NP-hardness to find a maximum weight clique, a heuristic procedure developed by James McCaffrey ¹⁷ was used.

Clique cuts interact with Branch-and-Bound as follows. At each Branch-and-bound iteration, each viable variable $y_i^r$ has its values resulting from the relaxed or the integer solution, which will represent the weights on the prohibition graph. At this point, a heuristic procedure is executed to find a maximum weight clique. At the moment one clique with a weight greater than one is found, a clique cut is inserted into the model, and all variables that belong to the cut have their weights updated to 0. This procedure is repeated until no clique is found anymore, and another Branch-and-Bound iteration is executed.

4 HYBRID METHOD

This section presents the hybrid method developed in this work, which combines heuristics, data mining, and the modified enhanced formulation (presented in Section 4.3). The objective of this hybrid approach is to find high-quality solutions within a reasonable computational time frame. Initially, a brief overview of the method is provided, followed by a detailed presentation of the algorithm and an explanation of each component’s role and their interaction.

The heuristic component contributes significantly by rapidly generating several high-quality solutions. Subsequently, data mining techniques are applied to these solutions to extract meaningful patterns, specifically subsets of sphere centers. These identified patterns are then incorporated into the modified enhanced formulation to guide the solution process. The Algorithm 1 below details the steps of the proposed hybrid method.

Algorithm 1:
Hybrid Method

In Algorithm 1, ES is a set of high-quality solutions and PS is the maximal pattern found by the mining Algorithm FPMax*, both of which are explained in more detail in Section 4.2.

4.1 Heuristic - LocalSolver

The heuristic component of the hybrid method used the LocalSolver framework (www. localsolver.com), a powerful tool in mathematical programming known for its unique hybrid approach. LocalSolver stands out from traditional solvers by integrating a variety of optimization techniques. These include local search techniques, constraint propagation, inference techniques, linear mixed-integer programming, and non-linear programming. This amalgamation of methods enables LocalSolver to effectively tackle complex problems, encompassing combinatorial and continuous optimization challenges.

Unlike traditional solvers, LocalSolver is adept at handling large-scale problems with millions of variables, a task often beyond the scope of conventional linear mixed-integer solvers. This capability was a key factor in its selection for this project.

Implementation details were carried out using LocalSolver’s C++ API, tailored to the complexity and scale of our problem. The problem was modeled in LocalSolver based on the formulation presented in Section 2.2. This specific formulation was chosen for its classical structures, which are well-suited to the capabilities of LocalSolver. These structures allow LocalSolver to efficiently explore and exploit the solution space, enhancing overall solution quality and computational performance.

4.2 Data Mining

According to ¹², data mining refers to the process of extracting a small set of valuable information from large volumes of data, typically represented in the form of patterns that capture relevant structural characteristics of the dataset ¹. In the context of this work, the role of data mining is not to perform a comprehensive pattern analysis, but rather to support the hybrid optimization framework by identifying structural information that can be directly exploited by the exact solution phase.

The Miner component is designed to identify sets of sphere centers that consistently co-occur in high-quality solutions generated by the heuristic phase. To this end, each high-quality solution is added to an Elite Set, which constitutes the transaction database for the mining process. Each solution is represented as a set of pairs {i, r}, where i denotes the center of a sphere and r its radius. Under this representation, the mining task consists of identifying groups of centers that appear together with high frequency across the Elite Set.

This task is performed using the FPmax* algorithm (11), which extracts maximal frequent item-sets. The worst-case computational complexity of FPmax* is exponential in the number of items, as is the case for frequent pattern mining algorithms in general. However, in the proposed hybrid method, the mining process is applied to a relatively small Elite Set composed only of high-quality solutions, which significantly limits the size of the transaction database. As a consequence, the computational cost of the mining step is negligible when compared to the subsequent exact optimization phases.

Alternative pattern-mining techniques, such as Apriori-based algorithms, FP-Growth, or methods for extracting closed frequent itemsets, could also be employed to identify recurrent structures in the Elite Set. However, these approaches often generate multiple overlapping or redundant patterns, which would require additional filtering or selection mechanisms before integration into the exact optimization model. In contrast, FPmax* directly returns maximal patterns, providing a compact representation of large sets of mutually compatible sphere centers.

This property makes FPmax* particularly suitable for the proposed hybrid framework, since the mined patterns can be incorporated directly into the modified enhanced formulation as structural constraints, without introducing additional decision layers or post-processing steps. The objective is not to claim optimality or superiority of FPmax* over other mining algorithms in general, but rather to adopt a method whose output is naturally aligned with the requirements of the subsequent exact optimization step.

4.3 Modified Enhanced Formulation

The necessary modifications that were made in the Enhanced Formulation are presented in this section. These modifications were made in the Enhanced Formulation to make it consider the maximal frequent sets, found by the Miner, as input data. The new model presented is the result of the modifications.

$$\begin{gathered} \max \sum _{i\in P}\sum _{r\in R}{c}_i^r{y}_i^r \\ \text{s.t.}\left(10\right),\left(12\right),\left(15\right) \\ \sum _{r\in R}{y}_i^r=1\forall i\in \text{PS} \end{gathered}$$ (16)

where the set of constraints (16) ensures that all centers of spheres i that are in the maximal pattern must be used by the Modified Enhanced Formulation.

5 COMPUTATIONAL EXPERIMENTS

This section delineates computational experiments, but before they are performed, some considerations must be made. As referenced in Theorem 1, optimal values exist for certain parameters, although they remain unidentified. To address this, the parameter $c_i^r$ is selected as r ³, representing the volume of spheres, a decision rooted in the geometric considerations of the coverage problem. This assumes uniform spherical coverage, a point that warrants consideration when interpreting the results. The set of radii values R = {2, 4, 7, 9} is used, chosen for its potential applicability in varied real-world scenarios.

The remainder of the section is organized into two subsections. The first subsection uses Data Envelopment Analysis (DEA) to examine various values of the overlapping parameter. DEA is a nonparametric method in operations research, which is crucial to estimate efficiency frontiers for different decision-making units (DMUs). Two important references on DEA are ^{2); (3}. The first includes a comprehensive review and comparative discussion of the fundamental DEA models. The second reference provides a comprehensive review and discussion of basic DEA models, extensions of these foundational methods presented in this issue, and a set of DEA applications in banking, engineering, health care, and services. In this work, the DMUs are represented by different values of the overlapping parameter, and the DEA is utilized to identify the most efficient parameter values. The second subsection is dedicated to evaluating the four proposed methods using the overlapping parameters determined from the DEA. The methods are tested in 15 instances, each named after the dimensions of a parallelepiped, which forms the set T . The results of these computational experiments will be presented through a combination of tables and graphical representations, which will facilitate a comprehensive analysis of the efficacy of the proposed methods.

5.1 Data Envelopment Analysis

In Section 2, the parameters of the problem were defined by setting α _rs = β for all r, s ∈ R. This choice is motivated by the role of α _rs in controlling the amount of overlap allowed between pairs of spheres. The overlap is governed by the constraint

$\Vert x_i-x_j\Vert \ge r+s-\beta ,$

which imposes a lower bound on the distance between the centers of two spheres of radii r and s. When β = 0, this constraint reduces to ∥x _i − x _j ∥ ≥ r + s, corresponding to the classical packing condition in which no overlap is allowed. As β increases, the minimum required distance between sphere centers decreases, allowing progressively larger intersections.

In this work, the set of admissible radii is fixed as R = {2, 4, 7, 9}. Hence, for any pair of spheres, the minimum possible value of r + s is equal to 4. When β = 2, the distance constraint allows the centers of two spheres to be separated by as little as half of this minimum value, allowing deep overlaps and, in extreme cases, configurations in which one sphere is almost entirely contained within another. Such configurations correspond to an excessive and unrealistic overlap from both a geometric and application-oriented perspective, particularly in the context of Gamma Knife treatment planning.

For this reason, β = 0 represents the case of no-overlap, while large values of β lead to degenerate solutions with limited practical relevance. Consequently, the values of β were restricted to the interval [0.1, 1.9] in the computational experiments, where a meaningful balance between coverage and overlap control can be achieved.

To determine the most suitable values of β for computational tests, DEA methodology was used, considering execution time as input and percentage coverage as output for each instance at a given β . DEA identified β = 0.2 and β = 1.1 as efficient values for the five instances, as illustrated in Figure 3.

Figure 3
Number of instances for which each beta was efficient.

The selected values β = 0.2 and β = 1.1, derived from the DEA analysis, were later used in the computational tests for all instances, which included exact and hybrid methods.

5.2 Computational Results

In this section, the computational results are presented. All algorithms were implemented in C++ using CPLEX 12.6 and LocalSolver 6.0, and executed on an Intel® XeonCore® CPU X5675 @ 3.07GHz with 48 GB of RAM. The maximum execution time was set to one hour for all methods.

The computational benchmark comprises 15 randomly generated parallelepiped instances with side lengths ranging from 5 mm to 30 mm. This instance family follows the same benchmark adopted in ⁸, which is the most closely related work on covering solids with spheres of different radii. This dimensional range corresponds to clinically plausible target volumes in the Gamma Knife planning context. Considering significantly larger dimensions would shift the problem to a different computational regime due to the rapid growth of the discretization grid and the induced search space, and would no longer be representative of the intended application.

Regarding the computational setup, it is important to emphasize that the proposed hybrid framework is not fully single-threaded. Different threading configurations were adopted for each component according to their intrinsic characteristics and to ensure a fair experimental evaluation. The exact components, namely the Branch-and-Bound (B&B) and Branch-and-Cut (B&C) procedures implemented in CPLEX, were executed using multithreading, which reflects their standard mode of operation in practice. In contrast, the heuristic component based on Local-Solver was deliberately executed using a single thread, since its performance is highly sensitive to parallel execution settings and because this choice improves reproducibility while avoiding solver-specific parallelization effects. Therefore, the proposed approach should be understood as a mixed-thread framework, in which the heuristic phase is single-threaded and the subsequent exact improvement phase fully exploits multithreading capabilities.

Since the heuristic component of the hybrid method is stochastic, each instance was solved multiple times under identical computational limits. Results for the hybrid approaches are therefore reported using aggregated indicators, including average solution values and average computational times. When variability was observed, coefficients of variation are reported to characterize performance stability.

In practice, the solution values produced by the hybrid methods exhibited negligible or no variability across runs within the reported precision. As a consequence, inferential statistical tests are not informative in this setting, since they rely on the presence of meaningful variance to assess performance differences. Instead, robustness is supported by the empirical stability of the solutions, by consistent dominance patterns across instances, and by the aggregation of results over heterogeneous instances and overlap settings.

In addition to solution values and execution times, we also report the discretization and mesh sizes (number of grid points) associated with each instance. These values provide a proxy for instance difficulty and highlight the heterogeneity of the benchmark, supporting robustness conclusions across markedly different search-space sizes.

Using the values of the overlap parameter β defined in Section 5.1, the Branch-and-Bound, Branch-and-Cut, and hybrid methods were executed, and the results are summarized in Tables 1 and 2. The two tables share the same structure and differ only in the value of the overlap parameter considered: β = 0.2 in Table 1 and β = 1.1 in Table 2. Each table reports, for each instance, the solution value (Sol), computational time (T), and duality gap (Gap) for the exact methods, as well as the average solution (Avg Sol), average computational time (Avg T), and coefficient of variation of the computational time (CV _T ) for the hybrid methods. The symbol “-” indicates that no feasible solution was found within the imposed time limit. Finally, a summary row reports the average solution value and average execution time for each method, providing an overall comparison of their performance.

Table 1
Computational results for exact and hybrid methods for β = 0.2.

Instances	Exact Method						Hybrid Method
	B & B			B & C			B & B			B & C
	Sol	T	Gap	Sol	T	Gap	Avg Sol	Avg T	CV T	Avg Sol	Avg T	CV T
5-20-8	144	3600	0.5	168	101.93	0	168	124.90	0.04	168	123.72	0.04
6-12-10	160	234.72	0	160	116.14	0	160	133.02	0.02	160	131.47	0.03
7-10-19	232	3600	0.5556	248	3600	0.0887	248	112.92	0.01	248	420.30	0.01
8-7-17	176	3600	0.4545	176	703.1	0	176	105.15	0.04	176	111.13	0.03
8-19-16	336	3600	1.6039	-	3600	-	464	1056.73	0.02	424	3600	0.00
9-16-23	-	3600	-	-	3600	-	584	3600	0.00	544	3600	0.00
10-16-18	-	3600	-	-	3600	-	512	3600	0.00	496	3600	0.00
11-19-6	240	3600	0.5444	256	1362.53	0	256	132.92	0.03	256	341.41	0.02
13-7-10	176	1272.56	0	176	483.11	0	176	106.07	0.01	176	108.56	0.04
14-12-10	216	3600	1.5185	-	3600	-	336	229.02	0.01	336	1281.28	0.02
15-20-14	-	3600	-	-	3600	-	735	3600	0.00	735	3600	0.00
18-15-18	-	3600	-	-	3600	-	855	3600	0.00	855	3600	0.00
19-10-13	-	3600	-	-	3600	-	448	3600	0.00	432	3600	0.00
19-12-5	144	3600	1.4259	216	987.02	0	216	112.62	0.02	216	127.85	0.01
21-9-14	-	3600	-	-	3600	-	464	3600	0.00	456	3600	0.00
Average	202.67	3220.48	-	200	2410.25	-	386.53	1580.89		378.53	1856.381

Table 2
Computational results for exact and hybrid methods for β = 1.1.

Instances	Exact Method						Hybrid Method
	B & B			B & C			B & B			B & C
	Sol	T	Gap	Sol	T	Gap	Avg Sol	Avg T	CV T	Avg Sol	Avg T	CV T
5-20-8	336	0.361	0	336	0.352	0	336	103.76	0.03	336	103.34	0.05
6-12-10	256	1926.08	0	256	377.19	0	256	110.75	0.06	256	218.78	0.05
7-10-19	352	3600	0.7223	408	3600	0.3013	456	3600	0.00	408	3600	0.00
8-7-17	312	3600	0.2316	320	3600	0.1447	336	439.31	0.02	328	3600	0.00
8-19-16	-	3600	-	-	3600	-	824	3600	0.00	680	3600	0.00
9-16-23	-	3600	-	-	3600	-	992	3600	0.00	992	3600	0.00
10-16-18	-	3600	-	-	3600	-	848	3600	0.00	808	3600	0.00
11-19-6	408	3600	0.3217	448	3600	0.0814	456	627.19	0.04	448	3600	0.00
13-7-10	296	3600	0.2560	304	3600	0.1627	312	2312.16	0.02	304	3600	0.00
14-12-10	-	3600	-	-	3600	-	584	3600	0.00	480	3600	0.00
15-20-14	-	3600	-	-	3600	-	1168	3600	0.00	1168	3600	0.00
18-15-18	-	3600	-	-	3600	-	1303	3600	0.00	1303	3600	0.00
19-10-13	-	3600	-	-	3600	-	696	3600	0.00	696	3600	0.00
19-12-5	424	3600	0.04789	424	2768.27	0	424	198.87	0.03	424	3600	0.00
21-9-14	-	3600	-	-	3600	-	808	3600	0.00	752	3600	0.00
Average	340.57	3248.43		356.57	3089.72		653.27	2412.80		625.53	3141.47

The data presented in Tables 1 and 2 offer a comparative analysis of the Branch-and-Cut method against the conventional Branch-and-Bound approach. It is observed that, across the majority of the computational experiments, the Branch-and-Cut method matched or surpassed the Branch-and-Bound method in terms of solution quality. Moreover, it demonstrated a higher efficiency in proving optimality for a greater number of instances, thereby highlighting its effectiveness. However, an interesting divergence is observed when considering β = 0.2. In this scenario, the Branch-and-Bound method successfully found feasible solutions for specific instances, namely 8 − 19 − 16 and 14 − 12 − 10, where the Branch-and-Cut did not.

In terms of hybrid methodologies, the hybrid framework with Branch-and-Bound showed a marked superiority over the exact approaches, particularly in aspects of computational time and solution quality. Specifically, when β = 0.2 was applied, this hybrid method matched the optimal solutions for nine instances. A similar trend was observed for β = 1.1, although with a reduced count of six instances that met this criterion. These observations underscore the nuanced performance variations between these algorithms under different parameters and problem instances, thereby informing future algorithms’ selection and development strategies in this domain.

Another important remark is that each component of the hybrid methods was important to guarantee its efficiency. Localsolver was able to find high-quality solutions, while Modified Enhanced Formulation was responsible for improving them. On average, the Modified Enhanced Formulation was able to improve the solutions obtained by LocalSolver, on average, by 11.23%.

Tables 3 and 4 detail the percentages of coverage for each method, in relation to the parameters β = 0.2 and β = 1.1, respectively. This percentage of coverage is derived from a mesh of points within the solid, set at a 0.1 unit distance from each other. The coverage is calculated as the ratio of the number of mesh points covered by at least one sphere to the total number of mesh points. The last column in both tables denotes the total number of mesh points for each instance.

Table 3
Covering results for exact and hybrid methods - β = 0.2.

Instances	Exact Method		Hybrid Method		Grid Points
	B & B	B & C	B & B	B & C
5-20-8	0.5277620	0.6211825	0.6211765	0.6211765	826200
6-12-10	0.6441184	0.6441036	0.6441036	0.6441036	745481
7-10-19	0.5947647	0.5857232	0.6163326	0.5929086	1362490
8-7-17	0.6354542	0.6111970	0.6499109	0.6499109	983421
8-19-16	0.4890330	-	0.6755290	0.6261309	2477790
9-16-23	-	-	0.6257332	0.5933639	3369730
10-16-18	-	-	0.6303687	0.6030057	2943241
11-19-6	0.6165971	0.6679415	0.6685050	0.6696562	1286490
13-7-10	0.6544287	0.6640710	0.6695277	0.6695277	939401
14-12-10	0.4485890	-	0.6405461	0.6405542	1723161
15-20-14	-	-	0.6523061	0.6523061	4258200
18-15-18	-	-	0.6594723	0.6594723	4946911
19-10-13	-	-	0.6325822	0.6229008	2513890
19-12-5	0.4054823	0.6062824	0.6022610	0.6022610	1172490
21-9-14	-	-	0.6139009	0.6175394	2694510

Table 4
Covering results for exact and hybrid methods - β = 1.1.

Instances	Exact Method		Hybrid Method		Grid Points
	B & B	B & C	B & B	B & C
5-20-8	0.9562116	0.9562116	0.9562116	0.9562116	826200
6-12-10	0.9040995	0.9001249	0.9066093	0.9025421	745481
7-10-19	0.7859346	0.8502051	0.9016308	0.8459959	1362490
8-7-17	0.8694476	0.8792674	0.9077028	0.8944359	983421
8-19-16	-	-	0.9083675	0.8613841	2477790
9-16-23	-	-	0.8984699	0.8984699	3369730
10-16-18	-	-	0.8842558	0.8625743	2943241
11-19-6	0.8646006	0.9163981	0.9223748	0.9129057	1286490
13-7-10	0.8657485	0.8797234	0.8968779	0.8848032	939401
14-12-10	-	-	0.9140057	0.8585832	1723161
15-20-14	-	-	0.8814816	0.8814816	4258200
18-15-18	-	-	0.8824430	0.8824430	4946911
19-10-13	-	-	0.8682870	0.8682870	2513890
19-12-5	0.9285683	0.9256446	0.9217725	0.9217725	1172490
21-9-14	-	-	0.8971037	0.8711350	2694510

For β = 0.2, as indicated in Table 3, the average coverages for the Branch-and-Bound and Branch-and-Cut methods are 0.56 and 0.63, respectively, with corresponding variances of 3.86 × 10⁻² and 9.56 × 10⁻⁴. In the context of hybrid methods, the coverage averages are 0.64 for Branch-and-Bound and 0.63 for Branch-and-Cut, with variances of 4.97 × 10⁻⁴ and 6.93 × 10⁻⁴, respectively.

Similarly, for β = 1.1 in Table 4, the coverage averages for the exact methods using Branch-and-Bound and Branch-and-Cut are 0.88 and 0.9, with variances of 3 × 10⁻³ and 1.2 × 10⁻³, respectively. The hybrid methods exhibit coverage averages of 0.9 for Branch-and-Bound and 0.88 for Branch-and-Cut, with variances of 4.48 × 10⁻⁴ and 8.21 × 10⁻⁴, respectively.

It is worth remarking that model (1)-(4) maximizes coverage only under the parameter conditions shown in Theorem 1, thereby aligning the optimal solution with that of the (CSSDRP). However, in this study, the parameter c was selected as r ³, representing the volume of the sphere and favoring larger spheres in the solution. Thus, while the results in Tables 3 and 4 affirm the applicability of the model to solve (CSSDRP), they are not utilized for comparative analysis of the methods.

Considering the aggregate factors of solution cost, computational time, and variance, it is concluded that the hybrid method employing Branch-and-Bound outperforms both the exact methods and the hybrid method using Branch-and-Cut.

6 CONCLUSION

In this paper, the problem of covering a solid with spheres of different radii was studied. The proposition of a family of new constraints led to a new formulation of the problem. Furthermore, regarding the solution of the problem, to improve the performance of Branch-and-Bound, the use of clique cuts, which consequently led to a Branch-and-Cut algorithm was proposed. Both Branch-and-Bound and Branch-and-Cut were combined with LocalSolver and a Frequent Itemset Mining technique, which led to the proposed hybrid methods aimed at achieving high-quality solutions with lower computational time.

A new family of constraints was derived that led to an improved formulation. In addition, the study proposed the integration of clique cuts into the Branch-and-Bound method, evolving it into an efficient Branch-and-Cut algorithm. These exact methods were further amalgamated with LocalSolver and Frequent Itemset Mining techniques, forming hybrid methodologies aimed at achieving superior solutions within reduced computational times.

As solution methods need a value for the β parameter, which is related to the overlap of spheres, a study using the DEA methodology was carried out to determine which values would produce better coverage results. The efficient betas obtained from this study were used in all computational experiments.

The tests concerning the exact methods, namely Branch-and-Bound and Branch-and-Cut were conducted, and the proposed Branch-and-Cut obtained better results proving the optimality of 6 instances out of 15 instances, against only 2 optimal proofs given by Branch-and-Bound, considering β = 0.2. Regarding β = 1.1, the Branch-and-Cut proved 3 optimal solutions versus 2 proved by Branch-and-Bound. It is worth mentioning that when the Branch-and-Cut did not solve the problem, it found better gaps than the Branch-and-Bound in most cases.

From a computational point of view, Branch-and-Cut used 74.8% of the computational time needed by Branch-and-Bound for β = 0.2 and for β = 1.1 this value is of 95.1%. These results show the superiority of the Branch-and-Cut approach in addition to computational time.

Another set of tests was carried out using Branch-and-Bound and Branch-and-Cut combined with LocalSolver and a Frequent Itemset Mining technique. In this hybrid scenario, the version using Branch-and-Bound delivered better solutions than the version using Branch-and-Cut. In addition, the hybrid method with Branch and Bound used 85.2% of the computational time needed by the version with Branch and Cut for β = 0.2 and for β = 1.1 used 76.8%.

In future work, one may invest in different formulations to deal with the overlapping and exploit the same strategy developed in this work. Also, it is possible to consider the development of other hybrid methods, since this work showed that hybrid methods are efficient in solving the Problem of Covering Solids with Spheres of Different Radii. Another line o investigation may follow the idea of exploring the percentage of the mined pattern that belongs to optimal solutions, which may lead to rethinking the data mining process.

Data Availability

Do not apply.

References

• ¹ BARBALHO H, ROSSETI I, MARTINS SL & PLASTINO A. 2013. A Hybrid Data Mining GRASP with Path-relinking. Computers & Operations Research., 40(12): 3159-3173.
• ² CHARNES AA, SEIFORD L, COOPER W & LEWIN A. 1995. Data Envelopment Analysis: Theory, Methodology and Application. Kluwer Academic Publishers, USA.
• ³ COOPER W, SEIFORD L & ZHU J. 2011. Handbook on Data Envelopment Analysis. International Series in Operations Research & Management Science, Springer US.
• ⁴ CORNFORTH D. 2002. Evolution in the Orange Box - A New Approach to the Sphere-Packing Problem in CMAC-Based Neural Networks. pp. 333-343.
• ⁵ DE SOUSA MM, GONZÁLEZ PH, OCHI LS & DE LIMA MARTINS S. 2021. A hybrid iterated local search heuristic for the traveling salesperson problem with hotel selection. Computers & Operations Research, 129: 105229.
• ⁶ FERRIS MC, JINHO L & SHEPARD DM. 2002. An optimization approach for radiosurgery treatment planning. SIAM Journal on Optimization, 13(3): 921.
• ⁷ GANZ JC. 1997. Gamma knife surgery. Springer-Verlag Wien, Austria.
• ⁸ GONZÁLEZ PH, MACAMBIRA AFU, PINTO RV, SIMONETTI L, MACULAN N & MICHELON P. 2020. New proposals for modelling and solving the problem of covering solids using spheres of different radii. RAIRO-Operations Research, 54(3): 873-882.
• ⁹ GONZÁLEZ PH, PINTO RV, SIMONETTI L, MACULAN N, MACAMBIRA AFUS & MICHELON P. 2016. Uma nova abordagem para o problema de recobrimento de sólidos com esferas de diferentes raios. Anais do XLVIII SBPO, pp. 4179-4189.
• ¹⁰ GONZÁLEZ PH, SIMONETTI L, MICHELON P, MARTINHON C & SANTOS E. 2016. A variable fixing heuristic with Local Branching for the fixed charge uncapacitated network design problem with user-optimal flow. Computers & Operations Research , 76: 134-146.
• ¹¹ GRAHNE G & ZHU J. 2004. Reducing the main memory consumptions of FPmax* and FPclose. In: Proc. Workshop Frequent Item Set Mining Implementations (FIMI 2004, Brighton, UK), Aachen, Germany, Citeseer.
• ¹² HAN J, PEI J & KAMBER M. 2011. Data Mining: Concepts and Techniques. The Morgan Kaufmann Series in Data Management Systems, Elsevier Science.
• ¹³ LEVIN M. 2018. Bin Packing Problems (Promising Models and Examples). Journal of Communications Technology and Electronics, 63: 655-666.
• ¹⁴ LIBERTI L, MACULAN N & ZHANG Y. 2009. Optimal configuration of gamma ray machine radiosurgery units: the sphere covering subproblem. Optimization Letters, 3(1): 109-121.
• ¹⁵ LONGANI P & KANTABUTRA S. 2014. Maximum coverage of radius-one sources in the wireless communication. In: 2014 International Computer Science and Engineering Conference (ICSEC), pp. 29-32.
• ¹⁶ LUBKE DC. 2014. Cobertura de Corpos por Esferas Utilizando Suavização Hiperbólica. Master’s thesis, COPPE/UFRJ, Rio de Janeiro.
• ¹⁷ MICROSOFT. 2011, Web Page. https://msdn.microsoft.com/en-us/magazine/hh547104.aspx
  » https://msdn.microsoft.com/en-us/magazine/hh547104.aspx
• ¹⁸ NANDI M. 2023. Coverage problem in two and three dimension. p. 11.
• ¹⁹ PINTO RV. 2015. O Problema de Recobrimento de Sólidos por Esferas de Diferentes Raios. Ph.D. thesis, COPPE/UFRJ, Rio de Janeiro.
• ²⁰ RIOFRIO D, LUAN S, ZHOU J & MA L. 2015. Particle Swarm Optimization for Radiation Therapy Planning. In: Proceedings of the 6th ACM Conference on Bioinformatics, Computational Biology and Health Informatics, Association for Computing Machinery, New York, NY, USA, pp. 250-257.
• ²¹ SANTOS E, OCHI LS, SIMONETTI L & GONZÁLEZ PH. 2016. A hybrid heuristic based on iterated local search for multivehicle inventory routing problem. Electronic Notes in Discrete Mathematics, 52: 197-204.
• ²² SCHEITHAUER G & TERNO J. 1993. Modeling of packing problems. Optimization, 28: 63-84.
• ²³ SOUTO G, MORAIS I, MAURI GR, RIBEIRO GM & GONZÁLEZ PH. 2021. A hybrid matheuristic for the Two-Stage Capacitated Facility Location problem. Expert Systems with Applications, 185: 115501.
• ²⁴ SUTOU A AND DAI Y. 2002. Global Optimization Approach to Unequal Sphere Packing Problems in 3D. Journal of Optimization Theory and Applications, 114(3): 671-694.
• ²⁵ V ROMERO S, OSABA E, VILLAR-RODRIGUEZ E, OREGI I & BAN Y. 2023. Hybrid approach for solving real-world bin packing problem instances using quantum annealers. Scientific Reports, 13(1).
• ²⁶ VENCESLAU M, VENCESLAU H, PINTO R, DIAS G & MACULAN N. 2021. A new penalty/stochastic approach to an application of the covering problem: the Gamma Knife treatment. RAIRO - Operations Research, 55: S787-S810.
• ²⁷ VENCESLAU MBK. 2015. O Problema de Recobrimento Mínimo de um Corpo em Três Dimensões por Esferas de Diferentes Raios. Ph.D. thesis, COPPE/UFRJ, Rio de Janeiro.
• ²⁸ WOLSEY LA. 2020. Integer programming. John Wiley & Sons.

• Funding
  This research was partially supported by Federal University of Rio de Janeiro (23079.227622/2023-70), Fundação COPPETEC, Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq (304680/2023-0, 304043/2025-7) and Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro - FAPERJ (E-26/201.341/2022).

Edited by

• Editor responsible for the review
  Editor-in-Chief: Annibal Parracho Sant’Anna.

Publication Dates

• Publication in this collection
  06 Mar 2026
• Date of issue
  2026

History

• Received
  01 Oct 2025
• Accepted
  09 Jan 2026