Optimized packing multidimensional hyperspheres: a unified approach

1.   Introduction

Optimized packing consists in placing a number of geometrical objects in a larger object called a container. The objects have to be arranged subject to placement conditions, i.e. placed without overlapping (non-overlapping condition) and completely inside the container (containment condition). Packing problems appear in different practical applications and one of the most frequently used placement problems is packing hyperspheres of different dimensions.

In biology and medicine applications of spheres packing include spatial organization of chromosomes in cell nucleus ^[1] and neurons ^[2,3], arrangement of ganglion cell receptive fields on retinal surface ^[4], planning radio-surgical treatment of tumors ^[5,6] and retinal laser coagulation ^[7]. Among various engineering applications one can find, e.g., cable bundling problems ^[8,9,10] and topology optimization in additive manufacturing ^[11,12,13], packing fuel elements in nuclear reactors ^[14] and heat exchangers ^[15]. In physics spheres packing arises in studying structure of nanomateralials ^[16], crystals ^[17], concrete ^[18] and granular materials ^[19], as well as in casting techniques ^[20]. Chemistry applications include packing catalysts in chemical reactors ^[21] or columns for gas distillation and absorption ^[22]. Examples of spherical packing in coding theory one can find in ^[23,24].

Hypersphere packing problems (HPP) form a broad area on the boundary between computational geometry and combinatorial optimization. Classical works of F. Toth ^[25], J. Conway, N. Sloane ^[23], T. Hales ^[26] and M. Vyazovska ^[27] study regular lattice hypersphere placement in space (lattice sphere packing). However, in many practical applications ^[28,29] packing hyperspheres in a bounded domain (container) has to be studied. Packing 2D & 3D spheres with balancing conditions were considered, e.g. in ^{[30,31,32,33,34]} using phi-functions ^[35] for modeling placement conditions. Then global/local optimizers combined with multistart or decomposition techniques were used to solve arising optimization problems. Among the principal characteristics of HPP are there space dimension, number of spheres, shape of the container, metric features (congruence, radii) of hyperspheres, correspondence between the sizes of hyperspheres and the container, etc. HPP is NP-hard ^[36] and thus heuristic approaches are widely used to obtain good approximate solutions in a reasonable computational time. To compare efficiency of different algorithmic approaches, open access collections of benchmark instances are used, see e.g. http://www.packomania.com.

There are a large number of publications on 2D & 3D sphere packing problems. However, the multidimensional (with dimension higher than 3) HPP are much less investigated and still of great interest. In this paper we focus on packing high dimensioned hyperspheres into arbitrary shaped containers subject to prohibited packing zones and minimal allowable distance between hyperspheres and/or container’s boundary.

The paper is organized as follows. Section 2 reviews related papers and highlights the main contributions of the paper. Section 3 provides the general mathematical model for packing multidimensional hyperspheres. Variants of the general mathematical model are stated in Section 4. Section 5 presents six solution strategies used to different classes of the original model. Numerical results are given in Section 6, while Section 7 concludes. Definition of the phi-function is presented in Appendix A.

2.   Related works

One of the most general typologies of Cutting & Packing Problems (C & P) was proposed by Wä scher et al. ^[37]. According to this typology, there are the following main types of HPP: ODP, PP (Placement Problem), KP, IIPP (Identical Item Placement Problem).

Many publications study 2D circle packing problems (CPP) arising in chemistry and geography, biology and production planning, logistics and additive manufacturing ^[38]. Hexagonal packing of equal circles was found to be the densest packing among all possible circle packings ^[39,40]. An optimization method for the open dimension circular packing problem (ODP) based on a descent with respect to groups of variables was presented in ^[41]. The method allows finding feasible solutions. This approach uses the observation that in a dense packing a circle touches either two other circles, or another circle and the container frontier, or the container frontier only. This technique is known as the block-coordinate descent method ^[42]. Huang et al. ^[43] proposed two greedy algorithms for packing circles into a rectangle of fixed dimensions. The first technique selects the next circle to be packed according to the maximum-hole degree rule. The second algorithm improves the later by a self-look-ahead search strategy that determines at each iteration the circle to be packed and its position. CPP with different container shapes, such as circles, squares, rectangles, strips and triangles are considered in ^[44]. Using a finite grid to approximate the container, optimized circle packing problem is transformed in ^[45,46] to a large scale linear integer programming problem. In ^[47] the approach is extended to packing the so-called circular-like objects that can be represented as circles in a certain (not necessary Euclidean) metric. CPP with prohibited zones are investigated in ^{[48,49,50,51]}. Prohibited zones in general lead to nonconvexity, multiconnectedness and/or nonconnectedness of the container.

For circular ODP hybrid algorithms combining beam and binary interval search with an open-strip generation procedure and a multi-start separate-beams strategy were proposed in ^[52,53]. Different models and methods for packing circles and spheres were reviewed in Hifi and R’Hallah ^[54]. The benchmark instances and the best known solutions for packing equal and non-equal circles into containers of different shapes are presented at E. Specht’s website ^[55].

Hales ^[26] obtained the upper bound for the density of packing equal spheres ^[40,56]. For packing unequal spheres in a container Sutou and Day ^[6] proposed a global optimization approach using a nonlinear programing formulation with quadratic constraints and a linear objective. Twice-differentiable models for 2D and 3D packing problems including packing different-sized spheres are presented in ^[57]. Kubach et al. ^[58] adapted the parallel greedy algorithms proposed in ^[52] for the 3D case. A hybrid algorithm for packing unequal circles and spheres into a larger circular (spherical) container is proposed in ^[59]. Stoyan et al. ^[60] proposed a modification of the jump algorithm ^[61] developed for CPP. This approach based on homothetic transformations was used for packing unequal spheres in various containers of minimum sizes (including the spherical container of minimum radius). A method for packing unequal spheres by combining the best-local position procedure with intensification and diversification stages is proposed in ^[29]. In ^[62] the problem of densest packing a given number of equal spheres into multiconnected containers is considered. The algorithm based on the optical-geometric approach and billiard simulation combination is proposed and implemented.

A package for 3-D Molecular Dynamics Simulations was developed by Bigrin et al ^[63,64,65]. The authors consider molecular simulations as a packing problem where the distance between atoms of different molecules has to be greater than some specified tolerance. The software allows packing millions of atoms, grouped in arbitrarily complex molecules, inside a variety of three-dimensional regions, including intersections of spheres, ellipses, cylinders, planes, or boxes in reasonable time. A review of modeling and solution techniques for packing spheres in various containers is presented in ^[66,67].

Comparing to 2D & 3D case, packing high dimensioned spheres is much less investigated. Random packing hyperspheres of higher dimensions using Monte Carlo method for molecular dynamics simulations is considered in ^[68]. To reach higher packing fractions a compression algorithm ^[69] or a particle scaling algorithm ^[70] are used. A random sequential addition algorithm for packing hyperspheres is proposed in ^[71]. Skoge et al. ^[72] study disordered jammed hard-sphere packings in 4D, 5D and 6D. They use a collision-driven packing generation algorithm ^[73] and obtain the estimates for the packing densities of the maximally random jammed states. The algorithm realizes homothetic transformations of hyperspheres with the common homothetic coefficient for all hyperspheres. Granocentric model for polydisperse sphere packings of high dimension is introduced in ^[74]. The homothetic transformations method for packing equal hyperspheres into a hypersphere of fixed radius is considered in ^[75], each hypersphere being not shared with another. The jump algorithm ^[61] was adopted for packing unequal hyperspheres into a hypersphere of minimum radius in ^[76].

To summarize, various approaches are used for HPP. Among them are modeling of sphere interaction by molecular dynamics and discrete elements methods; lattice and random packings; sequential addition and probabilistic methods; metaheuristic approaches (genetic and simulated annealing techniques, ant-colony and greedy algorithms); linear and nonlinear programming (continuous and integer); branch and bound algorithms for integer problems; hybrid approaches combining heuristics and mathematical programming methods, etc.

In this paper a unified methodology for packing hyperspheres into bounded containers of arbitrary shapes is presented. Additional restrictions, e.g., prohibited zones or distant conditions are also taken into account. Exact mathematical models and corresponding mathematical programming problems are formulated. Solution techniques are proposed and results of numerical experiments for the collections of benchmark problem instances are presented.

The main contributions of the paper are:

1) Phi-function ^[36] based modeling tools for packing hyperspheres into containers with prohibited zones bounded by hyperspheres, hypercylinders and hyperplanes.

2) General mathematical model of HPP for different types of objective function (ODP or KP), metric characteristics of hyperspheres (congruence, radii distribution, constraints on the radii values), shapes of the container (hyperrectangle, hypersphere, hypercylinder, d-polytope), restrictions on minimal allowable distances and prohibited zones.

3) Unified methodology to solve HPP based on efficient starting point algorithms and local optimization methods.

4) Computational results for packing multidimensional hyperspheres into different containers with prohibited zones.

3.   Mathematical model

We consider an optimization packing problem of hyperspheres of different dimensions (2D, 3D and $d$ D, $d \geqslant 4$) in the following formulation.

3.1.   Problem formulation

Let $\Omega ({\rm{ \mathsf{ μ} }})$ be a convex container with $k$ variable metric characteristics ${{\rm{ \mathsf{ μ} }} _1}, {{\rm{ \mathsf{ μ} }} _2}, ..., {{\rm{ \mathsf{ μ} }} _k}$ (sizes of the container). Here ${\rm{ \mathsf{ μ} }} = ({{\rm{ \mathsf{ μ} }} _1}, {{\rm{ \mathsf{ μ} }} _2}, ..., {{\rm{ \mathsf{ μ} }} _k})$. The shape of $\Omega ({\rm{ \mathsf{ μ} }})$ can be a hypersphere, a hypercylinder or a hyperrectangle. In addition, ${n_p}$ prohibition zones ${P_l}$, $l \in {I_p} = \{ 1, 2, ..., {n_p}\}$ are allowed to be arranged in the container. Each prohibited zone ${P_l}$ is a composition of hyperspheres, unbounded hypercylinders and/or hyperhalf-spaces.

In what follows the object $C({\rm{ \mathsf{ μ} }}) = \Omega ({\rm{ \mathsf{ μ} }})\backslash \operatorname{int} ({ \cup _{l \in {I_p}}}{P_l})$ is called a placement domain, where $\operatorname{int} ({ \cup _{l \in {I_p}}}{P_l})$ means the interior of the set $({ \cup _{l \in {I_p}}}{P_l})$. The placement domain with prohibited zones is in general a nonconvex set.

A collection of n hyperspheres ${S_i}({u_i}) = \{ X = ({x_1}, {x_2}, ..., {x_d}) \in {{\bf{R}}^d}:{\left\| {X - {u_i}} \right\|^2} \leqslant r_i^2\},$ $i \in {I_n} = \{ 1, 2, ..., n\}$ is given, where ${r_i}$ is radius of ${S_i}({u_i})$ and ${u_i} = ({x_{i1}}, {x_{i2}}, ..., {x_{id}})$ denotes a vector of variable centers of ${S_i}({u_i})$ for $i \in {I_n}$.

Conditions of packing hyperspheres ${S_i}({u_i})$, $i \in {I_n}$ into the domain $C({\rm{ \mathsf{ μ} }})$ are formulated as follows:

The non-overlapping conditions (2) can be extended regarding minimum allowable distances ${{\rm{ \mathsf{ ρ} }} _{ij}} > 0$ between the hyperspheres:

where

${\rm{ \mathsf{ ρ} }} (a, b)$ is the Euclidean distance between points $a$ and $b$.

If ${I_p} \ne \emptyset,$ then restrictions on prohibited zones should be taken into account:

Packing Problem of Hyperspheres (HPP). Pack hyperspheres from the set ${S_i}({u_i}), i \in {I_n}$ into the placement domain $C({\rm{ \mathsf{ μ} }})$ providing packing conditions (1)-(4) to optimize objective: maximize the packing factor or minimize the container $\Omega ({\rm{ \mathsf{ μ} }})$ volume.

Here the packing factor is defined as the following fraction: the volume of all hyperspheres divided by the volume of the placement domain.

3.2.   Tools of mathematical modeling

To formalize the packing conditions (1)-(4) the phi-function technique ^[36] is used.

The condition (1) can be described by means of phi-functions of the hypersphere ${S_i}({u_i})$ and the object $C{^*}({\rm{ \mathsf{ μ} }}) = {{\bf{R}}^{d + k}}\backslash \operatorname{int} C({\rm{ \mathsf{ μ} }}),$ for $i \in {I_n}.$

Let us define phi-functions of the objects ${S_i}({u_i}) \subset {{\bf{R}}^d}$ and $C{^*}({\rm{ \mathsf{ μ} }})$ for $d \geqslant 2.$

If $C({\rm{ \mathsf{ μ} }}) = {\rm{\{ }}u{\rm{ = }}({x_1}, {x_2}, ..., {x_d}) \in {{\bf{R}}^d}{\rm{: }}0 \leqslant {x_k} \leqslant {h_k}, {\rm{ }}k = 1, 2, ..., d\}$ is a hyperrectangle, then

if $C({\rm{ \mathsf{ μ} }}) = \{ {\rm{ }}u \in {{\bf{R}}^d}:{\rm{ }}{\left\| u \right\|^2} \leqslant r_0^2\}$ is a hypersphere, then

if $C({\rm{ \mathsf{ μ} }}) = \{ {\rm{ }}u \in {{\bf{R}}^d}{\rm{: }}x_1^2 + x_2^2 \leqslant r_0^2{\rm{, }}0 \leqslant {x_k} \leqslant {h_k}{\rm{, }}k{\rm{ = }}3, 4{\rm{, }}...{\rm{, }}d\}$ is a right circular hypercylinder, then

if $C({\rm{ \mathsf{ μ} }}){\rm{ = }}\{ {\rm{ }}u \in {{\bf{R}}^d}{\rm{: }}{A_{1m}}{x_1} + {A_{2m}}{x_2} +... + {A_{dm}}{x_d} + {B_m} \geqslant 0, {\rm{ }}m{\rm{ = }}1, 2, ..., M\}$ is a convex $d$ -polytope, then

The conditions (2) can be stated using phi-functions of the hyperspheres ${S_i}({u_i})$ and ${S_j}({u_j})$

The adjusted phi-functions of the hyperspheres ${S_i}({u_i})$ and ${S_j}({u_j})$ can be used to formalize the distance constraints (3) in the form

The conditions (4) are described by means of phi-functions of ${S_i}({u_i})$ and ${P_l}$ for $i \in {I_n}$ $l \in {I_p}$ ^[77].

3.3.   General mathematical model

Let $u = ({u_1}, {u_2}, ..., {u_n})$ be a vector of placement parameters of the hyperspheres $S_{i}\left(u_{i}\right), \quad i \in I_{n}$ and $t = ({t_1}, {t_2}, ..., {t_n}), {t_i} \in B = \{ 0, 1\}$ be a vector of binary variables for ${S_i}\left({{u_i}} \right), \quad i \in {I_n}$ where

Here ${{\rm{}}_i}({u_i}, {\rm{ \mathsf{ μ} }})$ is a phi-function of the objects ${S_i}({u_i})$ and ${C^*}({\rm{ \mathsf{ μ} }}).$

We denote a vector of variables by ${\rm{ \mathsf{ ω} }} = (u, {\rm{ \mathsf{ μ} }}, t).$

A general mathematical model of HPP can be presented in the following form:

where ${{\rm{ \mathsf{ κ} }}} ({\rm{ \mathsf{ ω} }})$ is the objective function; ${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Phi } _{ij}}({u_i}, {u_j})$ is an adjusted phi-function of the hyperspheres ${S_i}({u_i})$ and ${S_j}({u_j}),$ for $i < j \in {I_n};$ ${\Phi _i}({u_i}, {\rm{ \mathsf{ μ} }}) = \min \{ \Phi _{i0}^c({u_i}, {\rm{ \mathsf{ μ} }}), \Phi _{il}^c({u_i}), {\rm{ }}l \in {I_p}\}$ is a phi-function of objects ${S_i}({u_i})$ and $C{^*}({\rm{ \mathsf{ μ} }}),$ for $i \in {I_n};$ $\Phi _{i0}^c({u_i}, {\rm{ \mathsf{ μ} }})$ is a phi-function of ${S_i}({u_i})$ and ${\Omega ^*}({\rm{ \mathsf{ μ} }}) = {{\bf{R}}^{d + k}}\backslash \operatorname{int} \Omega ({\rm{ \mathsf{ μ} }}),$ $i \in {I_n};$ $\Phi _{il}^c({u_i})$ is a phi-function of ${S_i}({u_i})$ and ${P_l}, i \in {I_n}, l \in {I_p};$ ${g_q}({\rm{ \mathsf{ ω} }}) \geqslant 0, {\rm{ }}q \in {I_q}$ are restrictions on the placement parameters and metric characteristics of the container $\Omega ({\rm{ \mathsf{ μ} }}), {g_q}({\rm{ \mathsf{ ω} }})$, ${\rm{ }}q \in {I_q}$ are continuously differentiable functions.

Let us indicate some basic features of the problem (5):

1) The feasible region $W$ is, in general, a disconnected set with multiply connected components.

2) If the phi-functions describing $W$ contain maximum operators, it can be presented as $W = \bigcup\nolimits_{t = 1}^\eta {{W_t}},$ where the subregions ${W_t}, {\rm{ }}t = 1, 2, ..., \eta$ are described by systems with continuously differentiable functions.

3) The number of variables of the problem (5) is $O(n)$ and the number of inequalities is $O({n^2}).$

4.   Variants of the general mathematical model

Depending on the type of the objective function ${{\rm{ \mathsf{ κ} }}} ({\rm{ \mathsf{ ω} }})$ the packing problem (5) can be represented as: 1) the packing problem with variable metric characteristics of the container (ODP ^[37]); 2) the problem formulated as a knapsack problem (KP ^[37]); 3) the problem of packing identical hyperspheres (IIPP) being a partial case of KP.

As a container $\Omega$ is considered a hypersphere (in particular, a circle or a sphere), a hyperrectangle (including a rectangle or a cuboid), a hypercylinder (in particular, a right circular cylinder), a convex $d$ -polytope (in particular, a convex polygon or a convex polyhedron). The prohibited zones may be given as: circles, convex polygons, objects bounded by circular arcs and line segments for $d = 2;$ union of spheres, right circular cylinders, cuboids, convex right polygonal prisms for $d = 3;$ hyperspheres, hypercylinders and/or hyperhalf-spaces for $d \geqslant 4.$

4.1.   Open dimension problem (ODP)

Let there be hyperspheres ${S_i}({u_i})$ with radii ${r_i}, i \in {I_n}$ and a placement domain $C({\rm{ \mathsf{ μ} }}) = \Omega ({\rm{ \mathsf{ μ} }})\backslash \operatorname{int} ({ \cup _{l \in {I_p}}}{P_l}).$ The hyperspheres ${S_i}({u_i}), i \in {I_n},$ have to be packed in $C({\rm{ \mathsf{ μ} }})$ providing restrictions on the minimum allowable distances ${{\rm{ \mathsf{ ρ} }} _{ij}}$ between ${S_i}({u_i})$ and ${S_j}({u_j}), i < j \in {I_n},$ so that the sizes of the container will be minimized.

The mathematical model (5) for ODP takes the form

where $u = ({u_1}, {u_2}, ..., {u_n}), {\rm{ \mathsf{ μ} }} = ({{\rm{ \mathsf{ μ} }} _1}, {{\rm{ \mathsf{ μ} }} _2}, ..., {{\rm{ \mathsf{ μ} }} _k}), {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Phi } _{ij}}({u_i}, {u_j})$ is an adjusted phi-function of the hyperspheres ${S_i}({u_i})$ and ${S_j}({u_j}),$ for $i < j \in {I_n};$ ${\Phi _i}({u_i}, {\rm{ \mathsf{ μ} }})$ is a phi-function of ${S_i}({u_i})$ and $C{^*}({\rm{ \mathsf{ μ} }}), {g_q}(u, {\rm{ \mathsf{ μ} }}) \geqslant 0, {\rm{ }}q \in {I_q}$ are restrictions on the placement parameters of the hyperspheres and metric characteristics of the container $\Omega ({\rm{ \mathsf{ μ} }}), {g_q}(u, {\rm{ \mathsf{ μ} }})$ is at least twice continuously differentiable function.

Let us consider the features of ODP:

1) The number of variables of the problem (6) is $nd + k$ and the number of inequalities is $n + C_n^2 + Q.$

2) The Jacobian and Hessian matrices describing the constraints of the problem (6) are highly sparse.

3) If the container bounded by linear and inversely convex functions, then the minimum of the linear objective function is found at the extreme points of ${W_{ODP}}.$ Each extreme point is a solution of the system of $nd + k$ equations specifying the boundary of ${W_{ODP}}.$

4.2.   Knapsack problem (KP)

Let the container $\Omega$ be given by its fixed metric characteristics. Hyperspheres from the set ${S_i}({u_i}), i \in {I_n}$ should be packed into the placement domain $C$ providing the restrictions on the minimum allowable distances ${{\rm{ \mathsf{ ρ} }} _{ij}}$ between ${S_i}({u_i})$ and ${S_j}({u_j}), i < j \in {I_n},$ so that the packing factor will be maximized.

The mathematical model (5) for KP takes the form

where $u = ({u_1}, {u_2}, ..., {u_n}), t = ({t_1}, {t_2}, ..., {t_n}), {t_i} \in B = \{ 0, 1\}, i \in {I_n};$

${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Phi } _{ij}}({u_i}, {u_j})$ is an adjusted phi-function of the hyperspheres ${S_i}({u_i})$ and ${S_j}({u_j}),$ for $i < j \in {I_n}.$

Let us consider the features of KP:

1) The problem (7) is a mixed integer nonlinear programming (MINLP), the variables ${t_i}, {\rm{ }}i \in {I_n}$ are binary and the variables ${u_i}, {\rm{ }}i \in {I_n}$ are continuous.

2) The objective function $\Psi (u, t)$ is piecewise constant due to the presence of factors ${t_i}$ and is proportional to the sum of volumes of hyperspheres packed into $C$ (${t_i} = 1$).

3) If all hyperspheres from the set ${S_i}({u_i}), i \in {I_n}$ can be packed into $C,$ then $\Psi (u, t)$ attains the global maximum.

4) Local and global maxima are, in general, non-strict.

5) To each value $t \in {{\bf{B}}^n}$ there correspond a set of hyperspheres for which a part of values are ${t_i} = 1$ and the rest are ${t_i} = 0, i \in {I_n},$ the set ${W^t} = \{ u \in {{\bf{R}}^{dn}}:(u, t) \in {W_{KP}}\}$ defining various feasible packings of the hyperspheres regarding ${t_i} = 1, i \in {I_n}.$

Obviously, solving the problem (7) needs to handle all ${2^n}$ elements $t \in {{\bf{B}}^n}$ and define a point ${u^t} \in {W^t}$ for each set ${W^t}.$

4.3.   Identical item packing problem (IIPP)

Let there be a set $\Lambda$ of congruent hyperspheres ${S_i}({u_i}), i \in {I_\lambda } = \{ 1, 2, ..., \lambda \}$ with radius $r$ and a placement domain $C.$ The number of hyperspheres ${{\rm{ \mathsf{ σ} }} ^*} \leqslant \lambda$ from the set $\Lambda$ which are packed in $C$ without mutual overlappings should be maximized.

Given the sphere congruence, the mathematical model of KP (7) takes the following form:

where $u = ({u_1}, {u_2}, ..., {u_\lambda }), t = ({t_1}, {t_2}, ..., {t_\lambda }), {t_i} \in {\bf{B}} = \{ 0, 1\}, i \in {I_\lambda },$

${\Phi _i}({u_i})$ is a phi-function of the hypersphere ${S_i}({u_i})$ and the object $C{^*} = {{\bf{R}}^d}\backslash \operatorname{int} C.$

It should be noted that the maximum value of the objective function is equal to the number of spheres packed into the container and denoted by ${{\rm{ \mathsf{ σ} }} ^*} = \Psi ({u^*}, {t^*}).$

5.   Solution algorithms

The overall solution methodology includes analyses of the problem statement, initial data and restrictions; deriving and analysis of mathematical models for HPP in different dimension; constructing initial feasible packings (starting points or approximate solutions) and methods for local and global optimization. The basic elements of the methodology for solving HPP are shown in Figure 1.

5.1.   Constructing starting feasible packings

The following methods to construct starting feasible points are proposed: the random packing method, the lattice packing method, the greedy algorithm (modification of the block-coordinate descent method) and the residual optimization method.

The random packing method is used to quickly obtain feasible points for $d=4,5, \ldots, 19$ by choosing random values of hyperspheres centers in the Cartesian or hyperspherical coordinate systems and verifying feasibility.

The lattice packing method is used for 2D and 3D packing into a convex container without prohibited zones. It is assumed that the hyperspheres diameters are much smaller than the container’s size. For the hexagonal lattice packing each circle is tangent to 6 other circles and each sphere is tangent to 12 other spheres. Translations of the hexagonal lattices are realized to generate denser packings of circles (spheres) into the container and to construct various starting points which further can result in different local extrema.

The greedy algorithm is a modification of the block-coordinate descent method ^[42] with specific criteria. It can be used to obtain promising starting points in a reasonable time. In the greedy algorithm the problem (5) is reduced to the sequence of problems

The residual optimization method is used either for $d = 2, 3$ and the container with prohibited zones or for $d \geqslant 20.$ Let ${X^0} \notin {W_g} \subset {{\bf{R}}^\tau }$ where ${W_g} = \{ X \in {{\bf{R}}^\tau }:{g_l}(X) \geqslant 0, {\rm{ }}l \in {L_m}\}$ is the feasible region of the packing problem stated as a nonlinear programming problem, ${L_m} = \{ 1, 2, ..., m\}, \tau$ is the number of variables, $m$ is the number of restrictions. Define a set of indices of violated constraints: $L' = \{ l \in {L_m}:{g_l}({X^0}) < 0\}.$ To search for a feasible point the following problem is solved:

A point $({\chi ^0}, {X^0}) \in W'$ with ${\chi ^0} = \min \{ {g_l}({X^0}), {\rm{ }}l \in L'\}$ is chosen as a starting point. If the problem has a global maximum point $(\chi *, X*)$ where $\chi * = 0,$ then $X* \in {W_g}.$ If in the local (or global) maximum point $\chi * < 0,$ then $X* \notin {W_g}.$

5.2.   Local optimization

To find local extrema the interior point solver IPOPT (Interior Point Optimizer) ^[78]) is coupled with decomposition methods, modifications of the feasible direction method, the block-coordinate descent method ^[42] and the descent method based on analysis of Lagrange multipliers.

The combined decomposition and IPOPT method allows reducing the original problem with many nonlinear constraints to a sequence of non-linear programming problems with a significantly smaller number of non-linear constraints. The interior-point method used in IPOPT is a straight-dual interior point method (barrier function method) with proven convergence.

The block-coordinate descent method is used for ODP with more than 10, 000 variables. The variables are divided into several groups. Selecting one group of variables, the other variables are considered fixed. To solve the subproblems IPOPT or modifications of the feasible direction method are used. After solving all subproblems the best result is supposed to be an approximation to a local extremum of ODP. To improve the objective function value the process is repeated several times for different groups of variables.

The descent method based on analysis of Lagrange multipliers is applied to problems with inverse convex and linear constraints and is based on the analysis of Lagrange multipliers. It is used together with the block-coordinate descent method and the active constraints strategy. The variables are updated iteratively as follows:

To select the descent direction, the vector of Lagrange multipliers ${\lambda ^k} = (\lambda _1^k, \; \lambda _2^k, \; ..., \; \lambda _m^k)$ is obtained from the following system of linear equations

where $A({u^k})$ is the Jacobian matrix.

The constraint corresponding to the minimum negative multiplier is relaxed. One moves along the feasible region frontier decreasing the objective function until at least one of the non-active constraints is violated. The process continues until all Lagrange multipliers become nonnegative thus fulfilling the necessary and sufficient condition of the extremum.

5.3.   Global optimization

The problem HPP is NP-hard and thus the exhaustive search of all local extrema is not efficient. In this work methods of searching on a subset of the feasible set and cutting off its unpromising subsets are used.

The decremental neighborhood method (DNS ^[79]) is based on statistical properties of the objective on permutations of the spheres or the solution tree. Modifications of the method aim to organize direct search of the spheres sequences or the tree branches. The method is used to search for solutions of HPP and to construct promising starting points.

The branch and bound algorithm is adjusted to solve KP. To construct the solution tree, indexes of the binary variables are sorted taken into account hyperspheres having the equal radii. Each level of the tree corresponds to the radius of the certain subset of equal hyperspheres. The number of nodes of each level corresponds to the number of hyperspheres having the same radius. Each terminal node of the solution tree corresponds to a certain subset of hyperspheres from the given set. Pruning rules are based on the lower and upper bounds of the objective function value. The lower bound corresponds to the best currently found value of the objective. The starting value is set to 0 and dynamically updated while the solution tree constructing. The value corresponding to the best known packing factor is chosen as an upper bound. A node is considered as unpromising if it does not meet the corresponding lower and upper bounds.

The homothetic transformations method is based on allowing hyperspheres radii be variable. The sphere ${S_i}({u_i})$ with radius ${r_i}$ is denoted as ${S_i}({u_i}, {r_i})$ Problem IIPP (8) can be reduced to a sequence of problems of nonlinear programming with linear objective functions ^[75]:

where ${X^{\rm{ \mathsf{ σ} }} } = ({v^{\rm{ \mathsf{ σ} }} }, {u^{\rm{ \mathsf{ σ} }} }), {v^1} = {r_1}, {u^1} = {u_1}, {v^2} = ({r_1}, {r_2}), {u^2} = ({u_1}, {u_2}), {v^{\rm{ \mathsf{ σ} }} } = ({r_1}, {r_2}, ..., {r_{\rm{ \mathsf{ σ} }} }), {u^{\rm{ \mathsf{ σ} }} } = ({u_1}, {u_2}, ..., {u_{\rm{ \mathsf{ σ} }} }), {\rm{ \mathsf{ σ} }} = 3, 4, ..., {{\rm{ \mathsf{ σ} }} ^*}, {\varsigma _{\rm{ \mathsf{ σ} }} }({v^{\rm{ \mathsf{ σ} }} }) = \sum\nolimits_{i = 1}^{\rm{ \mathsf{ σ} }} {{r_i}}, {{\rm{ \mathsf{ σ} }} ^*}$ is an optimized value of the objective function of the problem (8),

${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Phi } _{ij}}({r_i}{\rm{, }}{r_j}{\rm{, }}{u_i}{\rm{, }}{u_j}{\rm{)}}$ is an adjusted phi-function of the hyperspheres ${S_i}({u_i}, {r_i})$ and ${S_j}({u_j}, {r_j})$ with variable radii ${r_i}$ and ${r_j}, i < j \in {I_{\rm{ \mathsf{ σ} }} };$ ${\Phi _i}({u_i}, {r_i})$ is a phi-function of the hypersphere ${S_i}({u_i}, {r_i})$ with variable radius ${r_i}$ and the set $C{^*}.$

Due to linearity of the objective function ${\varsigma _{\rm{ \mathsf{ σ} }} }({X^{\rm{ \mathsf{ σ} }} })$ local maxima of the problems (9) are attained at extreme points of the feasible region ${W_{\rm{ \mathsf{ σ} }} }.$

The homothetic transformations are also used to pack unequal hyperspheres (ODP) and form the core of the method of directed sorting of local extrema. Given a local minimum point of the problem (6), an auxiliary problem is formulated considering the hyperspheres radii as variables to define, while metric characteristics of the container are fixed ^[76]:

where $\Gamma (r) = \sum\nolimits_{i = 1}^n {r_i^n}$,

${r_{\min }} = \min \{ r_i^0, i \in {I_n}\}, {r_{\max }} = \max \{ r_i^0, i \in {I_n}\}, r_i^0, i \in {I_n}$ are the starting values of radii of the hyperspheres.

The objective function (10) is formed to maximize the volume (area) of the container. Solution of the problem (10) provides a starting point for the problem (6). The corresponding objective function value is at least as good as the previous local minimum. The process continues until there is no an improvement of the objective function.

The method is effective if the hyperspheres radii are uniformly distributed in the range between their minimum and the maximum values. The greater is the heterogeneity of the radii, the less useful is the method.

The multistart method is used together with all the methods of local and global optimization except for the branch and bound algorithm. The best local extremum is taken as an approximation to the global extremum.

To transform KP to ODP a container with variable homothetic ratio can be introduced. Thus, the methods developed for ODP can be applied for KP as well.

5.4.   Basic strategies

After analyzing all the mathematical models proposed, six basic strategies for solving HPP are highlighted: three for ODP (5) (Figure 2) and three for KP (IIPP) (6) (Figure 3). These strategies depend on: a) objective function type, b) problem dimension, c) metric characteristics of hyperspheres (congruence, radii distribution and values), d) container’s shape; e) prohibited zones in the container and/or minimal allowable distances.

Strategy 1. ODP-SA-DNS-permutations-IPOPT is based on DNS for sphere permutations and is used to pack a small number of spheres (10-60) with unequal radii significantly different one from another. Consistent statistical optimization is applied to select promising sphere permutations. Using randomly generated sequences and the greedy algorithm, a set of extreme points of the feasible region is generated. The best extreme point and the corresponding permutation are chosen as the center of the neighborhood on the permutation set. The process continues decreasing neighborhood radius. A set of points with minimum values of the objective function are taken as starting points, and a set of local minima of the problem are calculated by the decomposition methods and IPOPT. The best local minimum is selected as an approximation to the global minimum of the problem.

Strategy 2. ODP-SA-DNS-tree-IPOPT is another modification of DNS for packing 10-150 equal circles. A special solution tree that defines the extreme points of the feasible region is constructed by the greedy algorithm. Vertices of the tree are associated with extreme points (that can be obtained by the method of eliminating unknowns) and sequences of numbers. The Euclidean metric is introduced on the discrete set formed by the sequences. For directed search of extreme points the decremental neighborhood method is applied. The extreme points are, as in the previous modification, taken as starting points. With the decomposition method and IPOPT corresponding local minima are obtained. The best local minimum is considered as an approximation to the global minimum.

Strategy 3. ODP-random-JA-IPOPT is used to pack 10-300 hyperspheres with unequal radii varying gradually from the smallest to the largest (ODP). The radii are supposed to be variable. Feasible starting points are constructed by the random packing and the residual optimization method. The corresponding local minima are calculated by the decomposition method using IPOPT. Using the auxiliary problems, hyperspheres radii are redistributed filling up the container volume (area). The directed search of local extrema is applied yielding an approximation to the global minimum of ODP. This approach allows to make use of discrete and continuous nature of HPP.

Strategy 4. KP-random-tree-IPOPT is suitable to pack 10-30 unequal hyperspheres from a given set of hyperspheres into a fixed size container. A tree of solutions is constructed using the exhaustive search of all choices of the hyperspheres from the set (values of binary variables). A set of pruning rules reducing the number of tracing vertices is proposed. Based on the truncated tree subsets of spheres from the given set are formed. The subset of hyperspheres which can be packed into the container with the maximum sum of hypersphere volumes defines an approximation to the global maximum of KP. To verify if the hyperspheres are packed completely into the container auxiliary nonlinear programming problems are solved.

Strategy 5. IIPP-random / lattice-IPOPT is used to pack 10-5000 equal hyperspheres. It is based on the homothetic transformations method. First, a starting number of hyperspheres is chosen. To verify whether the hyperspheres are packed into the container an auxiliary problem is solved. Starting points are constructed either by the random packing, or the lattice packing, or the residual optimization depending on characteristics of the hyperspheres and the shape of the container. If a feasible packing is obtained, the number of hyperspheres increases by 1 and a new auxiliary problem is solved. If all starting points failed, then the number of hyperspheres obtained on the previous step is taken as an approximation to the global maximum. The strategy can be also applied to pack over 5000 hyperspheres. In this case, the block-coordinate descent method allows to calculate an approximation to a local maximum of the auxiliary problem. The duration of the process is limited by the computing resources available. For global optimization the block-coordinate descent method together with the multistart method is used.

Strategy 6. IIPP-random-SA-Lagrange is suitable for packing more than 5, 000 equal spheres. In Strategy 6 optimization is performed by the block-coordinate descent method. Feasible starting points for groups of variables are generated randomly. Descent directions are obtained by analyzing Lagrange multipliers. Iterations are terminated if the algorithm fails to generate a starting point for the current sphere.

6.   Computational results

Numerical experiments were implemented for the known collections of benchmark instances presented in ^[28,55,80]. Also, new instances were proposed for various container shapes and for packing hyperspheres into a hypersphere ($d \geqslant 4$). To solve linear and nonlinear programming problems the local solvers HOPDM ^[81], BPMPD ^[82], IPOPT ^[78] were used. To test the developed methods for IIPP ($d \leqslant 7$) the global solver GAMS / BARON ^[83] was applied. We used Intel Core I5 750 processor (2.5 GHz), 6 Gb RAM.

The proposed approach for packing hyperspheres was tested on instances available at the specialized websites and/or published in the corresponding papers. For the problems HPP (ODP, KP, IIPP) the dimension varied from $d = 2$ (circles) to $d = 24$ (hyperspheres).

The benchmark instances used for packing unequal circles into a minimal length rectangle (ODP) are available at the E. Specht website ^[55]. Strategy 3 was used to solve the problem. Two groups of instances were tested. The first consists of 25 instances with 20-300 circles. The comparison with the best-known results is shown in Table 1. Here the first column provides instance name, the second gives the number of circles, the third and the fourth indicate the minimal length obtained by the proposed approach and the best-known minimal length, correspondingly. Two last columns give the rectangle width and the packing factor.

The second group contains 128 instances with 25, 50, 75, 100 circles thus giving in total 153 instances in both groups. In 81 of 153 instances the best-known results were improved. Based on the numerical experiments we may conclude that the proposed technique is the most efficient for instances with 50-300 circles with radii evenly distributed in a given range. Figure 4 shows packing 300 circles in the rectangle of minimum length (Strategy 3). More examples are available in ^[55].

The other set of benchmark instances corresponds to packing equal circles in a circle of minimal radius given in ^[55] in the form of ODP. To perform experiments ODP was reduced to a sequence of IIPP solved in accordance with Strategy 5. The best-known results ^[55] were improved for some instances having from 1077 to 5000 circles.

Table 2 presents corresponding results, while Figure 5 shows packing of 5000 circles in the optimized circular container.

Thus, Strategy 5 using lattice packings as starting points is effective for more than 1000 circles where the lattice packing structure is perturbated marginally.

Packing spheres into a sphere for the benchmark instances with ${r_i} = i, i = 1, 2, ..., n, 20 \leqslant n \leqslant 35$ was introduced in ^[28]. The results are presented in Table 3. Here the first two columns provide the instance name and the number of spheres to pack. Two next columns define the radius $r$ obtained in ^[28] and the objective value of $r*$ found for the spherical container by Strategy 3. At the end of the Table 3 new results for the instances corresponding to $40 \leqslant n \leqslant 100$ (not reported in ^[28]) are presented.

Packing for 100 spheres (Example ZHXF100) is shown in Figure 6.

More results of packing spheres in containers of more complex geometries are illustrated in Figure 8. The results confirm viability of Strategy 3 for different container shapes. Data for these instances are provided in ^[80].

Packings for 100 spheres into a cylinder of minimum height and in an annular cylinder of minimum outer radius (ODP, Strategy 3) are shown in Figure 7.

Some numerical results on packing hyperspheres for $d \geqslant 4$ without prohibited zones one can find in ^[75,76,80].

New examples of packing unequal hyperspheres into a hypersphere of minimum radius with prohibited zones (ODP, Strategy 3) are provided below.

Example 1. $n = 23, d = 3, 4, 8, 16, 24, \{ {r_i}, i = 1, ..., 23\} =$ {1.2, 1.4, 1.7, 1.9, 1.2, 1.4, 1.7, 1.9, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0} are radii of ${S_i}({u_i}).$ The prohibited zone composed by two hyperspheres of radius $r = 3$ centered at $(2, 3, 2, 1, ..., 1, 1)$ and $(3, 2, 2, 1, ..., 1, 1)$ is given. The computational accuracy is 10⁻⁶. Table 4 presents the rounded values of the objective function. CPU is limited by 5 minutes.

Figure 9 illustrates projections on $O{x_1}{x_2}{x_3}$ of the results of Example 1 for $d = 3, 4, 8, 16, 24,$ the spheres packed are colored in purple and the prohibited zone is colored in blue.

Example 2. $n = 23, d = 3.$ Radii of spheres are taken from Example 1 and the prohibited zone is composed by three hyperspheres of radius $r = 1$ centered at $(0, 0, 1), (0, 1, 0)$ and $(1, 0, 0).$ The best objective value is ${r^*} = 4.6674.$

Example 3. $n = 23, d = 3.$ Radii of spheres are taken from Example 1 and the prohibited zone is the hypersphere of radius $r = 4$ centered at $(2, 2, 2).$ The best objective value is ${r^*} = 5.2471.$

Figure 10 provides illustrations for Examples 2 and 3. The spheres packed are colored in purple while the prohibited zone is colored in blue.

Computational results have shown that increasing dimension of hyperspheres (for $d > 3$) leads to dramatically increasing the number of hyperspheres’ contacts while decreasing the packing factor. This correlates the research ^[71].

Sometimes a “degeneration” of the problems with the small number of hyperspheres arises when increasing the hypersphere dimension. In the case the objective value is determined by the position of the largest hypersphere. Insignificant changes of the hyperspheres’ positions with the smaller radii do not affect the objective function value.

Increasing dimension of hyperspheres also makes harder finding the starting feasible solutions.

Figure 11 shows the degenerated case for 2D spheres: the red sphere is the prohibited zone, the blue sphere is the largest sphere packed. The objective value (the spherical container radius) is determined by the position of the blue sphere.

Therefore to avoid these situations the use of Strategy 3 (for $d > 3$) is recommended.

7.   Concluding remarks

In this paper smart technologies to solve HPP are proposed. The main factors affecting the computational time are the number of hyperspheres to be packed and their dimension.

Six basic strategies were proposed and tested for different classes of the hyperspheres packing problems. In many practical applications, complex container shapes with prohibited zones arise. To cope with this problem, it is necessary to state analytically corresponding placement conditions. The phi-function modeling approach was used to state non-intersection and containment for hyperspheres in a convex container with prohibited zones.

The other line for future research is developing techniques based on homothetic transformations of hyperspheres and the container. Also, creating new approaches for constructing feasible starting packing is very important in most packing iterative algorithms. An alternative way is studying the combinatorial structure of the packing problem ^[84] to carry out a directed search for local solutions in the configuration space of geometric objects ^[85].

An interesting direction for the future research is approximating complex objects or containers by simpler shapes ^[86], e.g. by spheres ^[87,88]. Most of subproblems arising in the proposed approach are large-scale optimization problems. Using agammaegation and/or decomposition techniques based on the special structure of constraints ^{[89,90,91,92]} may reduce computational time and increase the quality of solutions.

Acknowledgments

This work was partially supported by the Technological Institute of Sonora (ITSON), Mexico through the Research Promotion and Support Program (PROFAPI 2020) and by the Nuevo Leon State University (UANL), Mexico. The authors would like to thank E. Specht for maintaining the open source collection of circle packing instances and J. Gondzio and C. Meszaros for providing access to their solvers, HOPDM and BPMPD. T. Romanova was 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 interests

The authors declare no conflict of interest.

Appendix A

Let ${T_1} \subset {{\bf{R}}^d}$ and ${T_2} \subset {{\bf{R}}^d}$ be two objects. The positions of the objects ${T_1}$ and ${T_2}$ in ${{\bf{R}}^d}$ are defined by the corresponding vectors of placement parameters ${u_1}$ (${u_2}$).

A continuous and everywhere defined function $\Phi ({u_1}, {u_2})$ is called a phi-function for objects ${T_1}$ and ${T_2}$ if

Phi-functions allow distinguishing the following three cases: ${T_1}({u_1}) ~{\rm{and}}~ {T_2}({u_2})$ are intersecting so that ${T_1}({u_1}) ~{\rm{and}}~ {T_2}({u_2})$ have common interior points; ${T_1}({u_1}) ~{\rm{and}}~ {T_2}({u_2})$ do not intersect, i. e. ${T_1}({u_1}) ~{\rm{and}}~ {T_2}({u_2})$ do not have any common points; ${T_1}({u_1}) ~{\rm{and}}~ {T_2}({u_2})$ are tangent, i. e. ${T_1}({u_1}) ~{\rm{and}}~ {T_2}({u_2})$ have only common frontier points.

The inequality $\Phi ({u_1}, {u_2}) \geqslant 0$ describes the non-overlapping constraint, i.e. $\operatorname{int} {T_1}({u_1}) \cap \operatorname{int} {T_2}({u_2}) = \emptyset,$ and the inequality ${\Phi ^*}({u_1}, {u_2}) \geqslant 0$ describes the containment constraint $A({u_A}) \subset$ $B({u_B}),$ i.e. $\operatorname{int} {T_1}({u_1}) \cap \operatorname{int} T_2^*({u_2}) = \emptyset,$ where $T_2^* = {{\bf{R}}^d}\backslash \operatorname{int} T_2^{}.$

Let the minimum allowable distance ${\rm{ \mathsf{ ρ} }} > 0$ between objects ${T_1} ~{\rm{and}}~ {T_2}$ be given.

A continuous and everywhere defined function $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Phi } ({u_1}, {u_2})$ is called an adjusted phi-function for objects ${T_1} ~{\rm{and}}~ {T_2}$ if

The inequality $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\Phi } ({u_1}, {u_2}) \geqslant 0$ describes the distance constraint, i.e. $dist({T_1}({u_1}), {T_2}({u_2})) \geqslant {\rm{ \mathsf{ ρ} }}.$