Shape Optimization for Covering Problems

In his 1915 pioneering paper Nev15, Neville describes a game played at fairs where a large disk is painted on a cloth and five smaller, identical circular disks of thin metal are available. An award is offered to the person who is able to completely cover the large disk with the small disks. Neville then proceeds to show that the problem can be modeled by a system of nonlinear equations, and discusses a numerical method for an approximation of the solution. He also provides a figure for the covering of a large disk with five small disks that looks identical to the solution obtained with our algorithm shown in Figure 1. Many other works followed that dealt with the problem of covering a disk with smaller disks of minimum radius or a convex body with smaller homothetic copies. In general, planar geometry techniques are used in these works to, given a fixed and small number of disks, find the optimal radius or a bound for the optimal radius. See, for example, Kro93 and the references therein.

Figure 1.

Covering of a disk of unit diameter with five identical disks of minimum radius approximately 0.3023, obtained with the method introduced in BLMS21.

Nowadays, this type of problem is called a covering problem. The covering of the whole $n$-dimensional space with minimally overlapping identical balls has often been investigated in parallel to the problem of packing nonoverlapping spheres with the highest possible density CS99. It is well-known that the optimal covering of the plane is achieved by the hexagonal lattice, which realizes the thinnest covering, i.e., the plane covering with the least possible density of disk intersection; see Figure 2.

Figure 2.

In the left subfigure, the centers belong to a square lattice, while in the right subfigure they belong to a hexagonal lattice. The hexagonal lattice is visibly more efficient as it leads to less overlapping between the disks.

 

The covering of a bounded set with overlapping identical balls minimizing the number of balls with a fixed radius, or minimizing their radius with a fixed number of balls, as in Figure 3 in two dimensions, also represents a challenging question with a wide variety of practical applications, ranging from the configuration of a gamma ray machine radiotherapy equipment unit LMZ09 to the placement of base stations DDNS06. Compared to the problem of covering the whole space, in which case the solution is a lattice, covering a bounded set naturally yields a less regular solution, as the covering depends on the shape of the covered set. Numerical investigations hence play an important role in studying such problems. The covering of specific shapes such as squares, rectangles, disks, triangles, and polygons with a fixed number $m$ of small disks has naturally been the focus of various papers on this topic; see for instance HM97.

Figure 3.

Region $A$, in green, to be covered by a union of disks $\Omega (\mathbfit{x},r)$ (left) and $\Omega (\mathbfit{x},r) \cap A$, in gray (right).

 

Shape optimization approach

Even though various numerical methods have been introduced to solve the covering problem, the natural approach of considering the shape of the union of balls as the optimization variable has been generally overlooked, except in a few specific cases, see for instance HLE03. In this framework, the tools of shape optimization and shape calculus are employed to investigate the sensitivity with respect to variations of the union of balls, as the balls’ centers or radii are perturbed.

Shape optimization is the study of optimization problems where the variable is a geometric object, such as a subset of $\mathds{R}^n$ or a manifold; see SZ92. The shape sensitivity analysis is usually performed using strong regularity conditions on the geometry, in order to parameterize the perturbation of the geometry to compute derivatives. For instance, one often works with sets of class $C^k$ with $k\geq 2$, i.e., sets whose boundary can be locally represented by a function of class $C^k$. Still, many relevant shape optimization problems depend on mildly nonsmooth shapes such as curvilinear polygons, which means that their boundary is a union of smooth curves and it can have vertices. The covering of a set $A$ may be naturally formulated as a nonsmooth shape optimization problem, since $A$ may be nonsmooth, and the union of balls covering $A$ can be seen, except for degenerate cases, as a curvilinear polygon, as shown in Figure 4.

The shape optimization viewpoint on the covering problems opens up new perspectives as the tools of shape calculus become available, which allows us to numerically handle the case of a large number $m$ of balls. We describe now the approach that we have developed in BLMS21 and BLMS22. We focus here on a description in two dimensions for the sake of simplicity, but we emphasize that the theoretical part of the shape optimization approach is relatively independent of the dimension. Let $A$ be an open bounded subset of $\mathds{R}^2$ and $\Omega (\mathbfit{x},r) = \bigcup _{i=1}^m B(x_i,r)$, where $\mathbfit{x}\coloneq \{x_i\}_{i=1}^m$ and $B(x_i,r)$ is an open disk with center $x_i\in \mathds{R}^2$ and radius $r$. We consider the problem of covering $\overline{A}$ using a fixed number $m$ of closed disks $\overline{B(x_i, r)}$ with minimal radius $r$, where $\overline{A}$ denotes the closure of $A$. In other words, we are looking for a vector $(\mathbfit{x}, r) \in \mathds{R}^{2m+1}$ such that $\overline{A} \subset \overline{\Omega (\mathbfit{x}, r)}$ with minimal $r$. The optimization problem can be formulated as

$$\begin{equation} \operatorname *{Minimize}_{(\mathbfit{x},r) \in \mathds{R}^{2m+1}} \; r \; \text{ subject to } \; G(\mathbfit{x},r) = 0, {\tag{1}} \end{equation}$$

where

$$\begin{align} G(\mathbfit{x},r) &\coloneq \operatorname {Area}(A\setminus \Omega (\mathbfit{x},r)) {\tag{2}} \end{align}$$

and $\operatorname {Area}(A\setminus \Omega (\mathbfit{x},r))$ denotes the two-dimensional measure of $A\setminus \Omega (\mathbfit{x},r)$. It can be shown that $G(\mathbfit{x},r)=0$ is equivalent to $\overline{A} \subset \overline{\Omega (\mathbfit{x}, r)}$, thus a solution to problem 1 yields a solution to the covering problem. Covering problems are usually formulated with closed sets, but for the shape optimization framework it is convenient to work with the open sets $A$ and $\Omega (\mathbfit{x},r)$.

The function $G$ can be interpreted as the composition of a so-called shape functional $A\setminus \Omega \mapsto \operatorname {Area}(A\setminus \Omega )$ with a function $(\mathbfit{x},r) \mapsto A\setminus \Omega (\mathbfit{x},r)$. Under some geometric conditions detailed in BLMS21 and BLMS22, the derivative of such a function can be computed using techniques of shape calculus and in particular via the concept of shape derivative SZ92. Using these techniques, we proved in BLMS21 that the partial derivatives of $G$ with respect to the centers $x_i$ and the radius $r$ are given by

$$\begin{align} \partial _{x_i} G(\mathbfit{x},r) & = - \int _{\partial B(x_i,r)\cap \partial \Omega (\mathbfit{x},r)\cap A} \nu (z) \, dz,{\tag{3}}\\ \partial _{r} G(\mathbfit{x},r) & =\int _{\partial \Omega (\mathbfit{x},r)\cap A} \, dz, {\tag{4}} \end{align}$$

where $\nu$ is the outward unit normal vector to $\Omega (\mathbfit{x},r)$. The partial derivative $\partial _{x_i} G(\mathbfit{x},r)$ is a vector with two components since $\nu$ is vector-valued.

Note that the formulas for the partial derivatives of $G$ are boundary integrals involving the unit normal vector $\nu$ to the boundary of $\Omega (\mathbfit{x},r)$. This illustrates a fundamental property in shape optimization known as the structure theorem, which essentially states that the first derivative with respect to the shape of a shape functional only depends on the boundary perturbations in the normal direction. This is a well-known property of the shape derivative which is known since the pioneering work of Hadamard Had68.

In order to prove these results, the main task is to build a transformation $T_t$ between the union of disks $\Omega (\mathbfit{x},r)$ and its perturbation $\Omega (\mathbfit{x}+t\delta \mathbfit{x},r)$, where $\delta \mathbfit{x}$ is a perturbation of the disks’ centers. The transformation needs to be bi-Lipschitz, which means that $T_t$ and its inverse are both Lipschitz functions, in order to perform a change of variables in the integral on $\Omega (\mathbfit{x}+t\delta \mathbfit{x},r)$ that appears in $G(\mathbfit{x}+t\delta \mathbfit{x},r)$. The small perturbation $\Omega (\mathbfit{x}+t\delta \mathbfit{x},r)$ and the corresponding mapping $T_t$ are illustrated in Figure 4. In the case of a perturbation of the radii, one considers the set $\Omega (\mathbfit{x},r+t\delta r)$ and the procedure is similar.

In BLMS22 we have also computed the second-order partial derivatives $\partial ^2_{x_i,x_j} G(\mathbfit{x},r)$, $\partial ^2_{x_i,r} G(\mathbfit{x},r)$, $\partial ^2_{r,r} G(\mathbfit{x},r)$. Their expressions are more involved compared to 3,4 and we refer to BLMS22 for their detailed description, but their computation is based on similar shape calculus techniques as for the calculation of 3,4.

Figure 4.

A small perturbation of the disks’ centers modifies the shape of the union $\Omega (\mathbfit{x},r)$. The corresponding shape perturbation $\Omega (\mathbfit{x}+t\delta \mathbfit{x},r)$ of $\Omega (\mathbfit{x},r)$ can be parameterized with the help of a bi-Lipschitz transformation $T_t$.

Numerical methods and illustrations

In BLMS21, non-polygonal sets $A$ were considered. We even supposed that the set $A$ was defined by an oracle that, given a point in the plane, answers whether the point is inside $A$ or not. This level of generality forced us to compute $G$ and its gradient approximately. (At that time we were not yet considering second derivatives.) For the calculation of $G$, which can be written as an area integral, we simply partitioned the plane into small squares and counted the area of the squares whose center was at the same time in $A$ and in the union of the disks; see BLMS21, Alg. 4.1. For the calculation of the partial derivatives of $G$ in 3,4, we discretized $\partial B(x_i,r) \cap \partial \Omega (\mathbfit{x},r) \cap A$ for $i=1$, …, $m$ and $\partial \Omega (\mathbfit{x},r) \cap A$, which basically consists of considering a finite number of points on the boundary of each disk and of verifying, for each point, whether it lies within $A$ and is not in the interior of any disk. With the points satisfying these conditions, we approximated the line integrals using a quadrature rule; see BLMS21, Alg. 4.2. This approach allowed us to consider a large variety of sets $A$, but resulted in time-consuming and relatively imprecise calculations that enabled us to find approximate solutions for problems with up to 25 disks. Figure 5 illustrates some solutions; see BLMS21 for details. The highlight of the pictures is that the sets $A$ are not polygons.

It should also be noted that the discrete calculation of $G$ together with the minimization of the radius of the disks produces solutions in which boundary spikes may remain uncovered if the discretization is not sufficiently fine, as can be seen in Figure 5. This is due to the fact that spikes have a large perimeter and a small area, therefore their contribution to $G$ is much smaller than the contribution of the smoother regions of $A$, and a very fine discretization is necessary to obtain a visually satisfying covering. Numerically, a domain $D$ containing $A$ and $\Omega (\mathbfit{x},r)$ is partitioned into small squares, and the area of $A \cap \Omega (\mathbfit{x},r)$ is approximated by the sum of the areas of the squares whose centers are both in $A$ and in $\Omega (\mathbfit{x},r)$. If $A$ has a thin spike, as in Figure 5, the squares have to be very small for any of them to have their center inside the spike. A numerical study showing how the covering improves when the size of the small squares used in the discretization of $G$ decreases is presented in BLMS21.

Figure 5.

Heart covered by 15 radius-minimizing identical disks with $r \approx 0.3476$ (left) and curved star covered by 9 radius-minimizing identical disks with $r \approx 0.1300$ (right).

 

In BLMS22, using shape calculus techniques, we obtained the formulas of the second-order partial derivatives of $G$. Moreover, by limiting the numerical experiments to sets $A$ that are the union of disjoint convex polygons, we were able to calculate $G$, its gradient, and its Hessian exactly, disregarding the machine precision. The main tool for the calculation of $G$ and its derivatives was the calculation of Voronoi diagrams $\mathcal{V}(\mathbfit{x}) \coloneq \left\{ V\!_i(\mathbfit{x}) \right\}_{i=1}^m$, with the Voronoi cells

$$\begin{equation*} V\!_i(\mathbfit{x}) \coloneq \{ x \in \mathds{R}^2 \; : \; \| x - x_i \|_2 \leq \min _{1\leq j\leq m}\| x - x_j \|_2 \} \end{equation*}$$

for $i=1$, …, $m$, where the generators $x_i$ are the disks’ centers. Voronoi cells $V\!_i$ are disjoint (possibly unbounded) polyhedra which, intersected with the convex polygons composing $A$, become disjoint convex polygons. In turn, the intersection of each convex polygon with its corresponding disk results in a convex plane figure, whose sides are determined by segments and arcs of a circle, as in Figure 6. The area of each of these plane figures is easy to calculate with basic concepts of plane geometry. This partition of $\Omega (\mathbfit{x},r) \cap A$ as well as access to all vertices, edges, and arcs that compose the partition is all we need to calculate exactly not only $G$ but also its first- and second-order derivatives; see BLMS22, Algs. 4.1–4.3. These tools, by enabling the exact evaluation of $G$ and its first and second derivatives, have allowed us to calculate very precise solutions to problems with up to 100 disks in a fraction of the time required by the method introduced in BLMS21, which uses costly approximations of $G$ and its gradient, and does not use the second derivative of $G$. Figure 7 illustrates some solutions; see BLMS22 for details. The highlight of the pictures is that the corners of the regions are “well covered” by the disks. The description of the tools needed for the exact computation of $G$ and its derivatives makes more or less clear that the choice to consider sets $A$ given by the union of disjoint convex polygons was an arbitrary choice that simplifies the implementation of the calculations. Other types of sets $A$ could be considered, as long as one is willing to implement the required tools of plane geometry.

Figure 6.

Voronoi diagram generated using the centers of the disks as generating points that allows, by partitioning $\Omega (\mathbfit{x},r) \cap A$, to calculate its area.

Figure 7.

Minkowski island fractal (left) and eight-pointed star (right) covered by 100 minimizing radius identical disks with $r \approx 0.2753$ and $r \approx 0.2791$, respectively.

 

So far we have not mentioned how the optimization problems were solved. Problem 1 is a nonlinear programming problem with a single hard-to-compute constraint. The familiar tool for solving problems of this type is the augmented Lagrangian method; see BM14. In particular, we used Algencan ABMS07 which is the implementation of an augmented Lagrangian method with safeguards. Very roughly speaking, an augmented Lagrangian method solves a sequence of subproblems in which the violation of shifted constraints is penalized. In the specific case of the augmented Lagrangian method implemented by Algencan, bound constraints are not penalized and remain in the subproblems. However, since problem 1 has no bound constraints, the subproblems that Algencan solves are unconstrained. In Algencan, leaving aside other issues such as availability of a linear system solver and subproblem size, when the subproblems are unconstrained and second derivatives of the functions defining the problem are available, the subproblems are solved with a globally convergent line search Newton’s method. For this reason we can say that the work developed in BLMS22 is an application of a shape-Newton method in a genuinely nonsmooth setting, a notable fact since Newton’s method is rarely used in shape optimization, even in smooth settings.

It should also be emphasized that we are actually seeking a global solution of problem 1. There are no practical deterministic global optimization methods that are capable of dealing with problem 1. Thus, the stochastic options remain and, among them, the most natural is to use a multistart approach. This is precisely what we did in BLMS21 and BLMS22, using randomly generated starting points. In BGL, based on hexagonal lattices, we developed a way to generate better than random starting points. With that, we managed to improve all the solutions reported in BLMS22 using fewer initial points and, consequently, with lower computational cost.

Asymptotic analysis of the optimal radius

Following the numerical approximation of solutions of covering problems of bounded set, a theoretical question that naturally arises is the asymptotic behavior and bounds on the optimal radius, solution of problem 1, as the number of disks grows to infinity. The possibility of running numerical experiments with large $m$ allows us to verify the sharpness of asymptotic estimates for the optimal radius and to formulate new hypotheses for the bounds. This asymptotic analysis creates a bridge between the problems of covering the whole plane and covering a bounded set, as the optimal arrangement for a given $m$ converges in some appropriate sense toward the solution for the covering of the whole plane, which is the hexagonal lattice

$$\begin{equation*} \mathcal{L}_r \coloneq \{k v_r +\ell w_r\ |\ (k,\ell )\in \mathbb{Z}^2\} \end{equation*}$$

with $v_r \coloneq \frac{r}{2}(3,\sqrt {3})$, $w_r \coloneq \frac{r}{2}(3,-\sqrt {3})$. This convergence toward a lattice can be observed numerically using a large number of discs, as in the results in Figure 7. One observes that the disks located farthest from the boundary tend to align in a hexagonal lattice pattern, whereas the disks intersecting the boundary display a more intricate and less predictable behavior.

Kershner Ker39 pioneered the topic in 1939, providing an asymptotic result on the smallest number of disks of fixed radius that are necessary to cover an arbitrary region of the plane, a result that was improved ten years later by Verblunksy Ver49.

We have investigated a similar question recently in BGL for the covering of a general class of sets. Using honeycombs, defined as unions of $m$ regular hexagons whose centers belong to a regular hexagonal lattice, we have shown that the solution $r^*(m)$ to the minimization problem 1 satisfies

$$\begin{equation} r^*(m) = \left[ \frac{2\operatorname {Area}(A)}{3\sqrt {3}m} \right]^{1/2} + R(m) {\tag{5}} \end{equation}$$

where

$$\begin{equation} \underline{R}(m) \leq R(m)\leq \overline{R}(m), {\tag{6}} \end{equation}$$

with the following asymptotic expansions, as $m$ goes to infinity,

$$\begin{align*} \underline{R}(m) &= - \frac{ 2\operatorname {Per}(\partial A)}{3\sqrt {3}m} +O\left(\frac{1}{m^{3/2}}\right),\\ \overline{R}(m) &= \frac{2\operatorname {Per}(\partial A)}{3\sqrt {3} m} +O\left(\frac{1}{m^{3/2}}\right), \end{align*}$$

where $\operatorname {Per}(\partial A)$ denotes the perimeter of the boundary of $A$. We have then computed numerical solutions of covering problems with large $m$ in order to verify the sharpness of these asymptotic expansions. Numerically, we have observed that the asymptotic order $m^{-1}$ of $\overline{R}(m)$ is sharp, but that the constant $2\operatorname {Per}(\partial A) / ( 3\sqrt {3} )$ could probably be improved by a factor roughly equal to $0.2$.

We also observed in all our numerical experiments that the lower bound $\underline{R}(m)$ is always positive, which suggests that our asymptotic estimate for the lower bound is not sharp. To improve the lower bound estimate, a finer study of the behavior of the disks in the neighborhood of the boundary of the optimal covering set $\Omega (\mathbfit{x}^*(m),r^*(m))$ will have to be performed, as this behavior is the primary driver of the asymptotic expansions of the bounds $\underline{R}(m)$ and $\overline{R}(m)$.

Minimizing eigenvalues with respect to a union of disks

So far we have discussed several features of the covering of a bounded set with a union of $m$ disks. This problem does not involve the solution of a partial differential equation. Many shape optimization problems actually involve the solution of a partial differential equation, which usually models a physical phenomenon. For instance, the optimization of eigenvalues of differential operators with respect to geometrical features is a topic of high interest in pure and applied mathematics but also in engineering and natural sciences, such as in structural mechanics for the control of vibration frequency, in mathematical biology, acoustics or electromagnetism.

In particular, the optimization of Laplacian eigenvalues is a popular topic in mathematics as these problems are often simple and elegant to formulate, but are also challenging and require deep mathematical tools from a large spectrum of disciplines such as partial differential equations, spectral theory, and differential geometry. The celebrated Rayleigh–Faber–Krahn inequality, conjectured by Lord Rayleigh in the 19th century and proved several decades later by Faber and Krahn, states that the ball minimizes the first Dirichlet eigenvalue under a volume constraint. Since then, many shape optimization problems of this nature have been considered, such as the minimization of the $k$th eigenvalue of the Dirichlet Laplacian for $k>2$, or the minimization of eigenvalues with other types of partial differential equations and boundary conditions; see Hen06.

Let $\Omega \subset \mathds{R}^2$ be a bounded open set and introduce the function spaces

$$\begin{equation*} H^1(\Omega ) \coloneq \{v\in L^2(\Omega )\ |\ \nabla v\in L^2(\Omega ) \}, \end{equation*}$$

$$\begin{equation*} H^1_0(\Omega ) \coloneq \{v\in H^1(\Omega )\ |\ v = 0\text{ on }\partial \Omega \}, \end{equation*}$$ where $L^2(\Omega )$ denotes the space of functions whose square is integrable on $\Omega$. The first Dirichlet Laplacian eigenvalue is defined as

$$\begin{equation*} \lambda (\Omega ) \coloneq \min _{u\in H^1_0(\Omega )\setminus \{0\}} \frac{\displaystyle \int _{\Omega } |\nabla u|^2}{\displaystyle \int _{\Omega } u^2}. \end{equation*}$$

The corresponding eigenfunction $u$ satisfies

$$\begin{align} -\Delta u &= \lambda (\Omega ) u \text{ in }\Omega ,{\tag{7}}\\ u &= 0 \text{ on }\partial \Omega , {\tag{8}} \end{align}$$

and we impose the normalization condition $\|u\|_{L^2(\Omega )}=1$. The name “Dirichlet” refers to the boundary condition 8. We consider the following eigenvalue minimization problem:

$$\begin{align} \begin{split} \text{Minimize } & \lambda (\mathbfit{x}) \\ \text{ with respect to } & \mathbfit{x}\coloneq \{x_i\}_{i=1}^m \in \mathds{R}^{2m}, \end{split} {\tag{9}} \end{align}$$

where $\lambda (\mathbfit{x}) \coloneq \lambda (\Omega (\mathbfit{x}))$ denotes the solution to 7,8 with $\Omega = \Omega (\mathbfit{x}) \coloneq \cup _{i=1}^m B(x_i,r)$, with $m$ and $r$ fixed, $x_i\in \mathds{R}^2$ for all $i\in \{1,\dots ,m\}$.

The minimizers of problem 9 produce an interesting geometrical configuration of the centers $\{x_i\}_{i=1}^m$. First, when $\mathbfit{x}^\star \coloneq \{x_i^\star \}_{i=1}^m$ is a solution of problem 9, $\Omega (\mathbfit{x}^\star )$ must be connected due to a monotonicity property of Dirichlet eigenvalues. Second, $\Omega (\mathbfit{x}^\star )$ achieves an equilibrium between two competing tendencies. On the one hand, $\Omega (\mathbfit{x}^\star )$ strives to maximize its area, as the first Dirichlet eigenvalue tends to decrease as the area of the domain increases. On the other hand, $\Omega (\mathbfit{x}^\star )$ seeks to minimize the angles (measured from the interior of $\Omega (\mathbfit{x}^\star )$) at circle intersections, and to prevent small gaps from appearing, as these create strong singularities in the eigenfunction $u$ in the neighborhood of the circle intersections, which increase the eigenvalue. Here, “singularities” roughly means that $\nabla u$ is unbounded in these neighborhoods. Furthermore, considering that the minimizer of $\lambda$ with respect to a free-form shape under a volume constraint is a disk due to the Rayleigh–Faber–Krahn inequality, one expects the small disks $B(x_i^\star ,r)$ to agglomerate and approximate a large disk as $m\to \infty$, and to minimize their overlapping while avoiding gaps in $\Omega (\mathbfit{x}^\star )$, i.e., $\Omega (\mathbfit{x}^\star )$ should be simply connected. Note that the minimizing of overlapping is also characteristic of solutions of covering problems BLMS21BLMS22CS99, as discussed in the previous sections.

In BFHL23 we have developed an algorithm to find approximate solutions of problem 9. The approach is similar to the method described in the previous sections for the covering problem. Here we can also compute the derivative of the eigenvalue $\lambda (\mathbfit{x})$ with respect to $\mathbfit{x}$ using techniques of shape calculus. We have shown that the partial derivative of the eigenvalue is given by

$$\begin{align} \partial _i \lambda (\mathbfit{x}) = -\int _{\mathcal{A}_i} |\nabla u|^2 \nu , {\tag{10}} \end{align}$$

where $\mathcal{A}_i \coloneq \partial B(x_i,r) \cap \partial \Omega (\mathbfit{x})$ for $i\in \{1,\dots ,m\}$ and $\nu$ is the outward unit normal vector to $\Omega (\mathbfit{x},r)$. Notice the similarity between 10 and formula 3 which was obtained in the case of the covering problem. This similarity is due to the structure theorem of shape optimization that we have discussed above. There is nevertheless an important difference between 10 and 3, as $\partial _i \lambda (\mathbfit{x})$ depends on the gradient of the eigenfunction $\nabla u$. Indeed, $\nabla u$ is unbounded in the neighborhood of the circle intersections due to the nonsmoothness and the concavity of the shape at these points, and this makes the accurate numerical approximation of $\partial _i \lambda (\mathbfit{x})$ challenging.

Figure 8 shows the results obtained by our algorithm for $r=1$ and $m\leq 10$. We can clearly see how the shape of $\Omega (\mathbfit{x}^\star )$ converges to a disk as the number $m$ of disks grows. We also observe symmetries and regular patterns formed by the disks’ centers at the solutions. However, the solutions sometimes have less regularity and symmetries than expected. For instance, in the case $m=4$, the disks’ centers have two symmetries, but they are the vertices of a rhombus and not of a square. In the cases $m=6,7,8$, the results suggest that the exterior layer of disks should form a regular pentagon, hexagon, and heptagon, respectively. The case $m=5$ is interesting as it only possesses one axis of symmetry, and the pentagon formed by the disks’ centers is not regular. The structure of the solution is also remarkably similar to the solution of the covering problem with five disks shown in Figure 1. Starting from $m\geq 9$, regular patterns are more difficult to observe as the shape of $\Omega (\mathbfit{x}^\star )$ becomes more complex due to the appearance of two disks in the inner layer.

As in the covering problem, the asymptotic behavior of $\Omega (\mathbfit{x}^*(m))$ is an interesting theoretical question. As discussed in the previous sections, for the covering problem, the disks’ centers converge in some sense toward the hexagonal lattice. For the eigenvalue minimization problem, numerical results also suggest that the optimal placement of the disks’ centers probably converges, in some appropriate sense, toward a subset of the hexagonal lattice. This conjecture is also supported by a simple argument: asymptotically, the hexagonal lattice configuration is a trade-off between maximizing the total area of the union of balls, and keeping $\Omega (\mathbfit{x}^*(m))$ simply connected, i.e., avoiding any gaps in $\Omega (\mathbfit{x}^*(m))$.

Figure 8.

Numerical approximations of minimizers of the first Dirichlet eigenvalue for $m = 2,3,4,5,6,7,8,9,10$.

Conclusions and future research

Shape calculus and optimization are a powerful set of techniques for the sensitivity analysis of functions depending on the geometry. There exists an extensive literature in the smooth setting, but shape calculus still requires an active development in the nonsmooth case. Nonsmooth shape optimization has a variety of relevant applications such as the modeling of evolving nonsmooth sets and the optimization of complex geometries such as the union and intersection of moving components, curvilinear polygons, tessellations, generalized Voronoi diagrams and minimization diagrams BLM23. Applied to covering problems, it provided a new perspective on the problem and allowed us to design efficient numerical methods. Here we have presented results with a union of balls of identical radius, but the shape optimization approach is versatile and union/intersection of sets of various shapes can be treated in a similar way, also in dimension greater than two.

Of particular interest is nonsmooth shape optimization involving partial differential equations, as irregular geometries appear naturally in applications. Both theoretical and numerical challenges arise in this context, one of the main issues being the singularities that appear in the corners of the domain, which need to be carefully studied and handled numerically. The eigenvalue problem presented in this notice represents a first step in this direction, and other nonsmooth shape optimization problems involving partial differential equations will be investigated in the near future.