# Voronoi-based Multi-Robot Formations for 3D Source Seeking via Cooperative Gradient Estimation

Lara Briñón-Arranz, Martin Abou Hamad, Alessandro Renzaglia

► To cite this version:

Lara Briñón-Arranz, Martin Abou Hamad, Alessandro Renzaglia. Voronoi-based Multi-Robot Formations for 3D Source Seeking via Cooperative Gradient Estimation. ICARCV 2024 - 18th International Conference on Control, Automation, Robotics and Vision, Dec 2024, Dubai, United Arab Emirates. pp.1-6. hal-04762116

HAL Id: hal-04762116

https://inria.hal.science/hal-04762116v1

Submitted on 31 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Voronoi-based Multi-Robot Formations for 3D Source Seeking via Cooperative Gradient Estimation

Lara Briñón-Arranz<fnref n="1" />, Martin Abou Hamad<fnref n="1" /> and Alessandro Renzaglia<fnref n="2" />

**Abstract**— In this paper, we tackle the problem of localizing the source of a three-dimensional signal field with a team of mobile robots able to collect noisy measurements of its strength and share information with each other. The adopted strategy is to cooperatively compute a closed-form estimation of the gradient of the signal field that is then employed to steer the multi-robot system toward the source location. In order to guarantee an accurate and robust gradient estimation, the robots are placed on the surface of a sphere of fixed radius. More specifically, their positions correspond to the generators of a constrained Centroidal Voronoi partition on the spherical surface. We show that, by keeping these specific formations, both crucial geometric properties and a high level of field coverage are simultaneously achieved and that they allow estimating the gradient via simple analytic expressions. We finally provide simulation results to evaluate the performance of the proposed approach, considering both noise-free and noisy measurements. In particular, a comparative analysis shows how its higher robustness against faulty measurements outperforms an alternative state-of-the-art solution.

## I. INTRODUCTION

Cooperative source seeking is a fundamental problem in multi-robot applications because of its relevance in numerous different scenarios such as search and rescue operations, environmental monitoring, and pollution source detection among others [1]. In this problem, a spatially distributed signal field is generated by a source of unknown location. The team of autonomous robots is equipped with sensors able to measure the signal intensity at their locations. The final goal is then to cooperatively localize the source of the signal by exploiting the information gathered by each robot.

Numerous solutions have been proposed over the years to deal with the source-seeking problem with both single and multi-robot systems, but mostly limited to 2D environments. Some works presented gradient-free approaches such as the Speeding-Up and Slowing-Down method inspired by schools of fish seeking darker areas [2], the consensus-like Principal Component Analysis perception algorithm without explicit communication [3], and Particle Swarm Optimization based algorithms [4], [5]. However, the vast majority of the proposed strategies compute a local estimation of the signal gradient based on the collected measurements and use it to drive the robots toward the source location, as illustrated in Fig. 1. Obtaining an accurate estimation from sparse and usually noisy data remains nevertheless a challenging problem to solve. Moreover, only very few works addressed the

problem in 3D environments, which better represent many real-world scenarios. Possible approaches include consensus-based algorithms [6] and cooperative Kalman filters [7]. However, these strategies involve high communication and computational complexities, which are crucial aspects in the presence of small limited robotic platforms, and do not explicitly quantify the quality of the estimation algorithm, nor analyze the impact of measurement noise. Furthermore, in [6], the 3D formation can only be defined for specific numbers of robots, such as $N = \{6, 12, 20, 32, ...\}$.

To ensure formal guarantees on the obtained estimation of the gradient, an alternative interesting solution is to define multi-robot formations or configurations respecting specifically defined geometric constraints. Also in this case, the majority of works consider exclusively 2D environments, as in [8], [9]. Within this line, in our previous work [10], we showed that by adopting particular symmetric multi-robot formations, in both two and three dimensions, the estimates of the gradient and the Hessian matrix of an unknown signal can be obtained via simple averages that involve only the products of the measurements collected by the robots and their relative position with respect to the formation center. Despite the advantages brought by these formations, this approach presents a few drawbacks. Firstly, it requires to have an even number of robots. This may raise some problems in several scenarios, for instance in the case of a defective robot which would force the formation to reconfigure with the remaining $N-1$ robots. Furthermore, each robot needs to know its specific role within the formation to compute its position, requiring a centralized architecture. Finally, the proposed 3D formations are composed of two circles lying

<figure><img src="image_2.png" /><figcaption>Fig. 1. Cooperative source seeking in 3D: a formation of 7 robots follows the estimated gradient of a signal field to converge toward its source.</figcaption></figure>

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, 38000 Grenoble, France. Email: firstname.lastname@grenoble-inp.fr</footnote>

<footnote n="2">²Inria, INSA Lyon, CITI, UR3720, 69621 Villeurbanne, France. Email: firstname.lastname@inria.fr</footnote>
on two parallel planes that limit the spatial dispersion of the available sensors. This may produce some redundancy in the gathered data, especially with large teams.

The main contribution of this paper is to present a new multi-robot formation strategy able to overcome the limitations of the current state of the art and robustly estimate the gradient of a 3D scalar field to localize its source. In particular, we propose and analyze multi-robot formations that correspond to Centroidal Voronoi partitions of spherical surfaces. The advantages of this novel approach are:

* Suitable multi-robot formations can be defined for any number of robots $N$, with $N \ge 4$.

* The formations can be achieved in a distributed way, easing reconfiguration after possible robots' failures.

* Robots are spatially uniformly distributed resulting in a robust gradient estimation of the signal, that can be obtained via simple algebraic computations and whose error remains upper bounded.

## II. PROBLEM FORMULATION & PRELIMINARIES

The main objective of the cooperative source-seeking problem is to guide a team of robots to navigate through a signal field and localize its source. In this scenario, each robot represents a mobile sensor able to collect measurements of the signal strength emitted by the source. More formally, the signal distribution $\sigma(\mathbf{r}) : \mathbb{R}^3 \to \mathbb{R}^+$ is a function representing the scalar field at location $\mathbf{r}$ achieving its maximum at position $\mathbf{r}^*$ where the source is located and smoothly decreasing to zero far from the source. Denoting the signal gradient at a location $\mathbf{r}$ as $\nabla\sigma(\mathbf{r})$, we consider the following assumption on the signal:

**Assumption 1:** The function $\sigma: \mathbb{R}^3 \to \mathbb{R}^+$ is two times continuously differentiable, i.e., $\sigma \in \mathcal{C}^2$, and all its partial derivatives up to order two are globally bounded. Moreover, $\nabla\sigma(\mathbf{r}^*) = 0$ and $\nabla\sigma(\mathbf{r}) \ne 0, \forall \mathbf{r} \ne \mathbf{r}^*$. This implies that there exists a scalar $L$ such that

$$|\varphi(\mathbf{r}, \mathbf{c})| := |\sigma(\mathbf{r}) - \sigma(\mathbf{c}) - \nabla\sigma(\mathbf{c})^T(\mathbf{r} - \mathbf{c})| \le L \|\mathbf{r} - \mathbf{c}\|^2$$

where $|\cdot|$ denotes the absolute value, $\|\cdot\|$ denotes the Euclidean norm, and $\varphi(\mathbf{r}, \mathbf{c})$ corresponds to the first-order remainder of the Taylor expansion about the point $\mathbf{c}$.

The gradient $\nabla\sigma$ can then be used to drive the sensor network towards the signal source via gradient-ascent (GA) methods by steering the center of the formation as follows:

$$\mathbf{c}(k+1) = \mathbf{c}(k) + \varepsilon \nabla\sigma(\mathbf{c}(k))$$

where $\mathbf{c}(k) \in \mathbb{R}^3$ represents the formation center at iteration $k$ and $\varepsilon > 0$ is the step size of the GA algorithm. However, the gradient is usually not directly available and the objective is thus to design a multi-robot strategy able to cooperatively compute a gradient estimate, denoted by $\nabla\varphi(\mathbf{r})$, based on the collected measurements of the signal strength.

Consider a system composed of $N$ robots whose positions in the inertial global frame are denoted by $\mathbf{r}_i \in \mathbb{R}^3$ with $i = 1, \dots, N$. The robots are equipped with the appropriate sensors to measure the signal field of interest. Let $\sigma(\mathbf{r}_i)$ denote the measurement of the signal strength collected by

robot $i$. The robots can be controlled to stabilize their states, position, and velocity, to a desired reference. As presented in a large number of works dealing with robot control, the non-linear dynamics of the robots can be linearized and/or simplified in order to apply several commonly used control techniques (see [11] for an exhaustive survey) and ensure the tracking of a desired position. In this paper, the desired position $\mathbf{r}_i^{\text{ref}}$ for each robot $i$ in the team is provided by the multi-robot formation generator. The study of the low-level control of the robots is out of the scope of this work, consequently, we assume that the position controller of each robot guarantees a good tracking of this reference position, such that $\mathbf{r}_i \to \mathbf{r}_i^{\text{ref}}$. The robots are then controlled to maintain the desired formation while moving following the formation center. To consider realistic velocity constraints for the robots, we impose that the formation center moves with a limited velocity, i.e., $\|\mathbf{c}(k+1) - \mathbf{c}(k)\| \le \gamma$ where $\gamma > 0$ is a control parameter to be designed taking into account the physical limitations of the robots.

The solution proposed in this paper exploits some of the findings published in our previous work on symmetric formations for gradient estimation [10]. For the sake of clarity, we here briefly recall the necessary results.

Consider an even number of robots $N=2n, n \in \mathbb{N}$ forming a symmetric configuration composed of two parallel circular formations whose centers are aligned with the z-axis. The center point $\mathbf{c}$ is located between the two circles at distance $D \sin\theta_F$ from each one. One half of the robots is uniformly distributed in the upper circular formation, the other half is uniformly distributed in the lower circular formation and thus, their relative vectors with respect to the center point $\mathbf{c}$ are also evenly spaced. The robots' positions can be expressed in spherical coordinates as follows:

$$\begin{align}
\mathbf{r}_i &= \mathbf{c} + D[\sin\theta_i \cos\phi_i, \sin\theta_i \sin\phi_i, \cos\theta_i]^T \nonumber \\
\theta_i &= \begin{cases} \theta_F, & \text{if } i = 2k-1 \\ \pi - \theta_F, & \text{if } i = 2k, \end{cases} \quad k = 1, \dots, n \tag{1}
\end{align}$$

where $\phi_i = 2\pi i/N$ is the azimuthal angle, $D$ the radial distance to the center $\mathbf{c}$, $\theta_i$ is the polar angle and $\theta_F$ is defined such that $\sin\theta_F = \sqrt{2/3}$.

This Symmetric Cylindrical formation has been proposed to satisfy three interesting and convenient properties. Firstly, all the robots are placed at the same distance from the center, i.e., $\|\mathbf{r}_i - \mathbf{c}\| = D, \forall i$. Additionally, as proven in [10], when $N = 2n, n \ge 2$, the relative position vectors of the robots with respect to the center satisfy the following properties:

$$\sum_{i}^{N}(\mathbf{r}_{i}-\mathbf{c})=0 \quad \text{and} \quad \sum_{i}^{N}(\mathbf{r}_{i}-\mathbf{c})(\mathbf{r}_{i}-\mathbf{c})^{T}=\frac{N D^{2}}{3} \mathbf{I} \quad (2)$$

where $\mathbf{I} \in \mathbb{R}^{3 \times 3}$ represents the identity matrix. Considering the Symmetric Cylindrical formation defined by (1), the gradient of a 3D signal distribution $\sigma(\mathbf{r})$ can be estimated as presented in the following theorem:

**Theorem 1: (from [10])** Assume that $\sigma(\mathbf{r}) : \mathbb{R}^3 \to \mathbb{R}^+$ satisfies Assumption 1. Considering a team of $N = 2n$ robots with $n \ge 2$ forming a configuration given by (1) and defining
$$ \hat{\nabla}\sigma(\mathbf{c}) = \frac{3}{ND^2}\sum_{i=1}^{N}\sigma(\mathbf{r}_i)(\mathbf{r}_i - \mathbf{c}) \quad (3) $$

then it holds

$$ \|\hat{\nabla}\sigma(\mathbf{c}) - \nabla\sigma(\mathbf{c})\| \le 3LD. \quad (4) $$

The Symmetric Cylindrical formation proposed in [10] presents however some limitations: 1) an even number of robots is needed to satisfy the symmetric geometrical properties; 2) a central station is required to assign an identifier $i$ to each robot allowing defining their positions; 3) all the robots are located on two parallel circular formations, resulting in a limited distribution over the 3D space.

Inspired by this previous result, our objective is to propose a new formation that overcomes such limitations, i.e., allowing for any number of robots, avoiding the need for a central station, and achieving more dispersed measurements resulting in a higher coverage of the signal.

### III. VORONOI-BASED MULTI-ROBOT FORMATION

To achieve these objectives, we propose a new multi-robot formation strategy where the robots generate a Centroidal Voronoi Tessellation (CVT) of a spherical surface. As we show in this section, this solution allows the system to reach a high spatial dispersion while keeping important geometrical properties. Furthermore, each robot can now identify its position only based on the location of its teammates, i.e., without necessitating a centralized architecture.

Voronoi tessellation is a fundamental concept in Local optimization theory [12] often used to design multi-robot strategies, especially for optimal coverage [13], [14], [15]. Given a bounded set $\Omega \in \mathbb{R}^n$ and $N$ points $\{\mathbf{p}_i\}_{i=1}^N$, we define the Voronoi tessellation of $\Omega$ generated by $\{\mathbf{p}_i\}_{i=1}^N$ as $\{V_i\}_{i=1}^N$: $V_i = \{\mathbf{u} \in \Omega: ||\mathbf{u} - \mathbf{p}_i|| < ||\mathbf{u} - \mathbf{p}_j|| \forall j \neq i\}$, $i = 1, \dots, N$. The obtained tessellation is then called Centroidal Voronoi Tessellation (CVT) if $\mathbf{p}_i = \mathbf{p}_i^*, \forall i = 1, \dots, N$, where $\mathbf{p}_i^*$ is the centroid of the $i$-th region.

An important property of the CVT is that it represents a solution of the optimization problem:

$$ \min_{\mathbf{p} \in \Omega} \sum_{i} J_i(\mathbf{p}) \quad \text{where} \quad J_i(\mathbf{p}) = \int_{V_i} ||\mathbf{u} - \mathbf{p}||^2 d\mathbf{u} \quad (5) $$

proving that the Voronoi-based partition allows a good dispersion of the robots over the environment optimizing its coverage. In our context, this ensures having the different signal field measurements taken from well-spread locations, which is crucial to obtain a reliable gradient estimation.

A commonly adopted solution to obtain a CVT starting from an arbitrary set of generators is the Lloyd's algorithm [16], whose main steps are: calculate the Voronoi partitions and the relative centers of mass corresponding to the current generators configuration; move each generator to the center of mass of its region; repeat this procedure until the convergence or the maximum number of iterations are reached. Although obtaining Voronoi regions, and consequently CVTs, is quite straightforward in 2D, this is not the case for general surfaces in 3D [17]. To obtain a solution for our case, we firstly discretize our environment

<figure><img src="image_3.png" /><figcaption>Fig. 2. Constrained CVT on a spherical surface with 30 generator points.</figcaption></figure>

by creating a regular mesh of the spherical surface, following the procedure presented in [18]<fnref n="1" />, to easily compute the Voronoi regions at each iteration. Then, a non-trivial step of implementing the Lloyd's algorithm for a surface in 3D is the computation of the constrained mass centroid $\mathbf{p}_i^c$ of each region $i$. The classical mass centroids $\mathbf{p}_i^*$ are indeed always inside the sphere and a projection on the surface is needed before moving each generator to its new position. However, for a spherical surface, it can be shown that $\mathbf{p}_i^c$ is simply represented by the intersection of the radial line passing through $\mathbf{p}_i^*$ with the spherical surface [17]. Fig. 2 shows an instance of constrained CVT with 30 generators.

Finally, it is worth noticing that, as shown in [13], the Lloyd's algorithm can be implemented in an asynchronous distributed fashion. In this same work, the authors also presented an algorithm to adjust the sensing radius of every robot to ensure that each of them is able to communicate with all its Voronoi neighbors. Although the distribution of the complete approach, i.e., including the gradient estimation, is out of the scope of this paper, this is always a desirable property in presence of a multi-robot system. In our scenario, this would mean that the robots may find and maintain their positions in the formation in a distributed way only based on the information shared with their teammates.

Following the concepts presented previously, we propose a CVT-based Spherical formation to uniformly distribute a group of robots in a 3D formation. To build such a formation, the robots implement a constrained Lloyd's algorithm on the surface of a sphere. Let us define $\mathbf{p}_i^c$ for $i = 1, \dots, N$ the mass centroids corresponding to a constrained CVT on a sphere centered at the origin and of unitary radius. We express the position $\mathbf{r}_i$ of robot $i$ in the CVT-based Spherical formation with respect to the global frame as:

$$ \mathbf{r}_i = \mathbf{c} + D\mathbf{p}_i^c, \quad \forall i = 1, \dots, N \quad (6) $$

where $\mathbf{c}$ and $D$ are the center and the radius of the formation.

All the robots in the CVT-based Spherical formation are placed at the same distance from the center, i.e., $||\mathbf{r}_i - \mathbf{c}|| = D, \forall i$. Their positions, however, are not defined to satisfy by construction the properties in (2). Nevertheless, we can

<footnote n="1">¹Note that this approach does not work for any number of points and is not particularly suitable for small numbers of them.</footnote>
study these properties numerically by considering the robots
positions obtained via the constrained Lloyd's algorithm. Let
define the formation parameters $\bar{\mathbf{r}} \in \mathbb{R}^3$ and $\mathbf{M} \in \mathbb{R}^{3 \times 3}$ as

$$
\bar{\mathbf{r}} := \sum_{i=1}^{N} (\mathbf{r}_i - \mathbf{c}) \quad \text{and} \quad \mathbf{M} := \sum_{i=1}^{N} (\mathbf{r}_i - \mathbf{c})(\mathbf{r}_i - \mathbf{c})^T.
$$

Via numerical analysis, it can be verified that the summation
of the relative position vectors, $\bar{\mathbf{r}}$, remains close to zero<fnref n="2" /> and
that the values of matrix $\mathbf{M}$ are close to those of matrix
$\frac{ND^2}{3} \mathbf{I}$. These properties will be used to study the bound of
the gradient estimation error presented in the next section.

IV. 3D COOPERATIVE SOURCE-SEEKING

In this section, we present and analyze the cooperative
gradient estimation of a 3D signal obtained with Voronoï-
based multi-robot formations. Although the result on gradient
estimation presented in Theorem 1 could be still used with a
CVT-based formation, the estimation error would be greater
due to the geometrical formation properties that are not the
same as for the Symmetric Cylindrical formation. To address
this issue, we propose a new expression for the gradient
estimate by adding two correction terms to obtain a more
accurate estimation while keeping the elegant idea of using
only the summation of the relative position vectors pondered
with the signal measurements and simple computations.

Consider the CVT-based Spherical formation of robots
given by (6). The following theorem is proposed:

**Theorem 2:** Assume that $\sigma(\mathbf{r}) : \mathbb{R}^3 \rightarrow \mathbb{R}^+$ satisfies Assumption 1. Consider a team of $N$ robots forming a CVT-based Spherical formation with radius $D$ centered at $\mathbf{c}$ given by (6) and defining the estimated gradient as follows

$$
\hat{\nabla} \sigma(\mathbf{c}) = \mathbf{M}^{-1} \left( \sum_{i=1}^{N} \sigma(\mathbf{r}_i)(\mathbf{r}_i - \mathbf{c}) - \sigma(\mathbf{c}) \sum_{i=1}^{N} (\mathbf{r}_i - \mathbf{c}) \right) \quad (7)
$$

then it holds

$$
\| \hat{\nabla} \sigma(\mathbf{c}) - \nabla \sigma(\mathbf{c}) \| \le \| \mathbf{M}^{-1} \| \mathrm{NLD}^3 = \| \frac{ND^2}{3} \mathbf{M}^{-1} \| \mathrm{3LD}. \quad (8)
$$

*Proof:* Using the first-order Taylor expansion of each measurement $\sigma(\mathbf{r}_i)$ about the point $\mathbf{c}$ and recalling that $\|\tilde{\mathbf{r}}_i\| = D$ where $\tilde{\mathbf{r}}_i := \mathbf{r}_i - \mathbf{c} = D\mathbf{p}_i^c$ as defined in (6), then the following equation holds for all $i = 1, \dots, N$:

$$
\sigma(\mathbf{r}_i) - \sigma(\mathbf{c}) = \nabla \sigma(\mathbf{c})^T \tilde{\mathbf{r}}_i + \varphi(\mathbf{r}_i, \mathbf{c}),
$$

where $\varphi(\mathbf{r}_i, \mathbf{c})$ denotes the remainder of the Taylor expansion.
Multiplying the previous equation by $\tilde{\mathbf{r}}_i$ and summing over
$i = 1, \dots, N$, we obtain

$$
\sum_{i=1}^{N} \sigma(\mathbf{r}_i) \tilde{\mathbf{r}}_i - \sigma(\mathbf{c}) \sum_{i=1}^{N} \tilde{\mathbf{r}}_i = \sum_{i=1}^{N} \nabla \sigma(\mathbf{c})^T \tilde{\mathbf{r}}_i \tilde{\mathbf{r}}_i + \sum_{i=1}^{N} \varphi(\mathbf{r}_i, \mathbf{c}) \tilde{\mathbf{r}}_i \\
= \left( \sum_{i=1}^{N} \tilde{\mathbf{r}}_i \tilde{\mathbf{r}}_i^T \right) \nabla \sigma(\mathbf{c}) + \Phi(D, \mathbf{c}),
$$

where $\Phi(D, \mathbf{c}) = \sum_{i=1}^{N} \varphi(\mathbf{r}_i, \mathbf{c}) \tilde{\mathbf{r}}_i$. By using the definition of $\mathbf{M}$, the previous equation can be rewritten as

$$
\sum_{i=1}^{N} \sigma(\mathbf{r}_i) \tilde{\mathbf{r}}_i - \sigma(\mathbf{c}) \sum_{i=1}^{N} \tilde{\mathbf{r}}_i = \mathbf{M} \nabla \sigma(\mathbf{c}) + \Phi(D, \mathbf{c})
$$

and therefore the gradient estimation (7) satisfies

$$
\hat{\nabla}\sigma(\mathbf{c}) = \nabla\sigma(\mathbf{c}) + \mathbf{M}^{-1}\Phi(D, \mathbf{c}).
$$

According to Assumption 1 the term $\mathbf{M}^{-1}\Phi(D, \mathbf{c})$ satisfies

$$
\|\mathbf{M}^{-1}\Phi(D, \mathbf{c})\| \le \|\mathbf{M}^{-1}\| \sum_{i=1}^{N} L \|\mathbf{r}_i - \mathbf{c}\|^3 \le \|\mathbf{M}^{-1}\| \mathrm{NLD}^3. \blacksquare
$$

*Remark 1:* The previous result relies on the signal measurement at the center of the formation $\sigma(\mathbf{c})$. Note that, this value can be easily estimated by a simple interpolation computing the average of the measurements collected by the robots, i.e., $\sigma(\mathbf{c}) \approx \hat{\sigma}(\mathbf{c}) = \frac{1}{N} \sum_{i=1}^{N} \sigma(\mathbf{r}_i)$.

This new result allows estimating the gradient of the signal
at the center of CVT-based Spherical formations and proves
that the estimation error is bounded. The bound of the error
is similar to the bound obtained in Theorem 1 but multiplied
by an additional value that depends on $\mathbf{M}$. As previously
stated, for a CVT-based Spherical formation, this matrix is
close to $\frac{ND^2}{3}\mathbf{I}$. Via numerical analysis<fnref n="3" />, it can be shown that
$\|\frac{ND^2}{3}\mathbf{M}^{-1}\| \approx 1$. Consequently, the bound of the estimation
error for both results is very similar.

A. Estimation error analysis

Consider a CVT-based Spherical formation given by (6).
Let us now consider that the multi-robot formation computes
the gradient estimation by using (3). Following similar
derivations as in the proof of Theorem 2, the gradient
estimate (3) satisfies

$$
\hat{\nabla} \sigma(\mathbf{c}) = \frac{3\sigma(\mathbf{c})}{ND^2} \sum_{i=1}^{N} \tilde{\mathbf{r}}_i + \frac{3}{ND^2} \mathbf{M} \nabla \sigma(\mathbf{c}) + \frac{3}{ND^2} \Phi(D, \mathbf{c}).
$$

Let us assume that $\mathbf{M} = \frac{ND^2}{3}\mathbf{I} + \mathbf{M}_\alpha$ where $\mathbf{M}_\alpha$ is a small square matrix<fnref n="4" />. The estimation error can be expressed as

$$
\hat{\nabla}\sigma(\mathbf{c}) - \nabla\sigma(\mathbf{c}) = \frac{3\sigma(\mathbf{c})}{ND^2}\bar{\mathbf{r}} + \frac{3}{ND^2}\mathbf{M}_{\alpha}\nabla\sigma(\mathbf{c}) + \frac{3}{ND^2}\Phi(D, \mathbf{c})
$$

and the following inequality holds

$$
\| \hat{\nabla} \sigma(\mathbf{c}) - \nabla \sigma(\mathbf{c}) \| \le \frac{3\sigma(\mathbf{c})}{ND^2} \| \bar{\mathbf{r}} \| + \frac{3}{ND^2} \| \mathbf{M}_{\alpha} \| \| \nabla \sigma(\mathbf{c}) \| \\
\qquad + \frac{3}{ND^2} \| \Phi(D, \mathbf{c}) \| \le 3LD + \Phi_{\alpha}
$$

where $\Phi_{\alpha} = \frac{3}{ND^2} (\sigma(\mathbf{c})\|\bar{\mathbf{r}}\| + \|\mathbf{M}_{\alpha}\|\nabla\sigma(\mathbf{c})\|)$ denotes the additional term which depends on the CVT-based formation properties and on the signal distribution.

We can notice that, as expected, $\Phi_{\alpha}$ decreases with both
the number of robots $N$ and the radius $D$. We then recall
that the terms $\sigma(\mathbf{c})$ and $\|\nabla\sigma(\mathbf{c})\|$ are both bounded while
$\|\bar{\mathbf{r}}\|$ and $\|\mathbf{M}_{\alpha}\|$ are close to zero. Moreover, if we focus
on the proximity of the signal source, we can see that the
term $\|\nabla\sigma(\mathbf{c})\|\|\mathbf{M}_{\alpha}\|$ rapidly goes to zero because $\nabla\sigma(\mathbf{c}) \to 0$
when $\mathbf{c} \to \mathbf{r}^*$. However, for strong signals, $\sigma(\mathbf{c})$ can be a

<footnote n="2">²For the formations used in this work, with a sphere of radius $D = 1$, $\bar{\mathbf{r}} \approx 10^{-4} \mathbf{1}$ where $\mathbf{1}$ denotes the unit vector.</footnote>

<footnote n="3">³For the formations used in this work, $||\frac{ND^2}{3}\mathbf{M}^{-1}|| \le 1.04$.</footnote>

<footnote n="4">⁴For the formations used in this work, $||\mathbf{M}_{\alpha}|| \le 0.3$ with a radius $D = 1$.</footnote>
large value and, even if $||\bar{\mathbf{r}}||$ is close to zero, this might lead to a value of $\Phi_{\alpha}$ that remains larger than the bound of the estimation error obtained in Theorem 2.

## B. Noise analysis

To consider more realistic scenarios, we assume that the signal measurements collected by the sensors are corrupted by Gaussian zero-mean white noise, i.e.:

$$y(\mathbf{r}_i) = \sigma(\mathbf{r}_i) + v_i(\mathbf{r}_i), \quad v_i(\mathbf{r}_i) \sim \mathcal{N}(0, \mathbf{v}^2)$$

where $y(\mathbf{r}_i)$ is now the noisy measurement collected by robot $i$ and $\mathbf{v}^2$ represents the variance of the noise. This situation can be the result of the noise produced by the physical limitations of the sensors and their electronic components. The noise can also model possible small-scale spatial variations of the robots' positions due to turbulence or local perturbations. Let us assume that the noise is independent in each robot measurement, i.e., $E[v_i(\mathbf{r}_i)v_j(\mathbf{r}_j)] = 0, i \neq j$, which is a realistic assumption if $||\mathbf{r}_i - \mathbf{r}_j|| > D_{min}$ where $D_{min}$ represents the spatial correlation distance of the disturbance.

Since the noise is additive in the measurement, there will be an additional additive estimation error when computing the gradient estimate defined in (3). This additional error caused by the noisy measurements has zero mean and it has been shown that its expected standard deviation is a monotonically decreasing function of the formation radius and the number of robots (see [10]).

Another advantage of CVT-based formations is that the minimum distance between two robots is larger than for Symmetrical Cylindrical ones and this difference increases by increasing the number of robots. To illustrate this difference, let us study the minimum distance $d_{min} = \min\{||\mathbf{r}_i - \mathbf{r}_j||, \forall i, j, i \neq j\}$ for the two formations. For the Symmetric Cylindrical formation, $d_{min}$ can be analytically obtained by simple trigonometric computations. Considering the formation (1), we have $d_{min} = 2D\cos(\theta_F)\sin(\frac{2\pi}{N})$. For the CVT-based formation, we obtain this minimum distance via numerical computations. Table I presents the values of $d_{min}$ for the two formations with $D = 1$ and different values of $N$.

<table><caption>TABLE I
MINIMUM DISTANCE $d_{min}$ WITH RESPECT TO $N$.</caption><thead><tr><th>N</th><th>6</th><th>8</th><th>10</th><th>20</th></tr></thead><tbody><tr><td>Symmetric Formation d<sub>min</sub></td><td>1.00</td><td>0.81</td><td>0.68</td><td>0.37</td></tr><tr><td>CVT-based Formation d<sub>min</sub></td><td>1.41</td><td>1.14</td><td>1.06</td><td>0.76</td></tr></tbody></table>

This property is important because it strengthens the assumption about the independence of noisy measurements and it proves a better uniform distribution of the robots around the center of the formation. As shown in the next section, this leads to an increased robustness to asymmetric noisy measurements, for example, when one robot sensor has a level of noise greater than the rest of the team.

## V. SIMULATION RESULTS AND ANALYSIS

In this section, we present numerical simulations to illustrate the performance of the proposed cooperative estimation algorithms applied to the source-seeking problem in 3D. The

<figure><img src="image_4.png" /><figcaption>Fig. 3. Evolution of the three components of the gradient estimate $\hat{\nabla}\sigma(\mathbf{c})$ for a team of $N = 7$ robots with CVT-based Spherical formations compared to the real computed gradient of $\sigma_1(\mathbf{r})$.</figcaption></figure>

signals considered in the simulations are Gaussian functions whose surface sets are ellipsoids given by $\sigma_1(\mathbf{r}) = A_1e^{-\mathbf{r}^T\mathbf{S}_1\mathbf{r}}$ where $A_1$ determines the intensity of the signal strength and $\mathbf{S}_1 \in \mathbb{R}^{3 \times 3}$ is a symmetric matrix. The maximum corresponding to the source is located at $\mathbf{r}^* = [0, 0, 0]^T$.

In our analysis, we study two different scenarios: first the case of a strong signal without any noise and a second case with a weaker signal where measurements are affected by noise. The signal functions for these two scenarios are: $\sigma_1(\mathbf{r})$ with $A_1 = 100$ and $\mathbf{S}_1 = 10^{-4} \begin{bmatrix} 1 & 1 & 1 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \end{bmatrix}$ and $\sigma_2(\mathbf{r})$ with $A_2 = 1$ and $\mathbf{S}_2 = 10^{-4} \begin{bmatrix} 20 & 5 & 2 \\ 5 & 20 & 2 \\ 2 & 2 & 100 \end{bmatrix}$ respectively. For all the results, simulation parameters are set to $\varepsilon = 1$ and $\gamma = 0.1$.

Fig. 1 illustrates a CVT-based formation with $D = 4$ and an odd number of robots ($N = 7$) performing a GA algorithm where the gradient of $\sigma_1$ is estimated by (7). The red line represents the trajectory of the formation center, the blue circles represent the initial robots' positions and the green ones an intermediate state at iteration $k = 500$. The gradient estimation for this scenario is shown in Fig. 3. The exact computation of the gradient $\nabla\sigma(\mathbf{c})$ is compared with two gradient estimates $\hat{\nabla}\sigma(\mathbf{c})$ computed by (3) and by (7) respectively with the same CVT-based formation. The data show how the correction term presented in (7) provides an estimation closer to the exact gradient.

We then analyze the effect of the radius $D$ on the gradient estimation error. Fig. 4 shows the evolution of the norm of this error, i.e., $||\hat{\nabla}\sigma(\mathbf{c}) - \nabla\sigma(\mathbf{c})||$, where $\hat{\nabla}\sigma(\mathbf{c})$ is computed by (7), for three different values of $D$. This study supports the analytical results given in (8): as expected the norm of the error increases monotonically with the radius. In particular, we can see that, although $D = 4$ remains significantly closer to $D = 1$ than $D = 7$, its final error at convergence is larger and the minimum error is obtained with $D = 1$.

Let us now consider the second scenario of a weaker signal with robots' measurements affected by noise. To show the robustness of our solution to noisy measurements, we compare the obtained results with the Symmetric Cylindrical formation presented in [10]. Note that comparisons with other alternative solutions are difficult to obtain due to the very limited attention that the current literature has dedicated
<figure><img src="image_5.png" /><figcaption>Fig. 4. Evolution of the norm of the gradient estimation error, i.e., ||∇σ(c) - Vσ(c)||, with a CVT-based Spherical formation for different values of the formation radius D.</figcaption></figure>

to the case of 3D fields with noisy signal measurements (see discussion in Section I). Fig. 5 shows a statistical analysis of the evolution of the distance ||**c** − **r**\*|| for a series of 100 trials to compare the influence of noisy measurements on the source-seeking convergence for a Symmetric Cylindrical formation and a CVT-based Spherical formation of N = 8 robots. The top figure shows the results when v = 0.1 for all the robots in the formation. In this situation, the behavior of the two formations is almost identical. The figure at the bottom shows the results where one robot obtains faulty measurements, modeled by a higher level of noise v = 0.5, while for all the others v = 0.1. This asymmetry in the noise level has an important influence on the gradient estimation when computed by the Symmetric Cylindrical formation that strongly relies on its geometrical properties to ensure a good gradient estimate. Consequently, the convergence toward the source is affected resulting in a greater final error and a much larger standard deviation along its entire trajectory. On the other hand, the CVT-based formation is less affected and presents a behavior not far from the equal noise case. This represents a significant improvement in robustness against unexpected faulty measurements that can be of paramount importance in real applications.

VI. CONCLUSION

This paper presented a new strategy to design suitable multi-robot formations that allow an accurate and robust estimation of the gradient of a 3D signal field generated by a source of unknown position. This estimation, obtained via simple analytic expressions that ensure a low computation cost, is then leveraged to guide the multi-robot formation through the field to finally localize the source. Results in simulations supported the analytic analysis of the proposed approach. Moreover, they showed how it is able to outperform an alternative approach, providing a higher robustness with respect to faulty measurements.

For future work, our primary intent is to build upon these
results to obtain the full distribution of this solution, includ-
ing limited communication capabilities, while maintaining
theoretical guarantees on the gradient estimation.

REFERENCES

[1] B. Bayat, N. Crasta, A. Crespi, A. M. Pascoal, and A. Ijspeert, "Environmental monitoring using autonomous vehicles: a survey of recent searching techniques," *Current opinion in biotechnology*, 2017.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Mean and standard deviation over 100 trials of the distance between the center and the source location for both Symmetric and CVT-based formations in presence of noisy measurements with v = 0.1 for all robots (top) and with v = 0.5 only for one robot and v = 0.1 for the others (bottom).</figcaption></figure>

[2] W. Wu and F. Zhang, “A speeding-up and slowing-down strategy for distributed source seeking with robustness analysis,” *IEEE Transactions on Control of Network Systems*, vol. 3, no. 3, pp. 231–240, 2015.

[3] S. Al-Abri, W. Wu, and F. Zhang, “A gradient-free three-dimensional source seeking strategy with robustness analysis,” *IEEE Transactions on Automatic Control*, vol. 64, no. 8, pp. 3439–3446, 2018.

[4] R. Zou, V. Kalivarapu, E. Winer, J. Oliver, and S. Bhattacharya, “Particle swarm optimization-based source seeking,” *IEEE Transactions on Automation Science and Engineering*, vol. 12, pp. 865–875, 2015.

[5] J. Zhang, Y. Lu, Y. Wu, C. Wang, D. Zang, A. Abusorrah, and M. Zhou, “PSO-based sparse source location in large-scale environments with a uav swarm,” *IEEE Transactions on Intelligent Transportation Systems*, 2023.

[6] R. Fabbiano, C. Canudas De Wit, and F. Garin, "Source localization by gradient estimation based on poisson integral," *Automatica*, vol. 50, no. 6, pp. 1715–1724, 2014.

[7] W. Wu and F. Zhang, “Cooperative exploration of level surfaces of three dimensional scalar fields,” *Automatica*, vol. 47, no. 9, 2011.

[8] A. Kahn, J. Marzat, H. P. Lahanier, and M. Kieffer, “Cooperative estimation and fleet reconfiguration for multi-agent systems,” *IFAC-PapersOnLine*, vol. 48, no. 5, pp. 11–16, 2015.

[9] A. Renzaglia and L. Briñón-Arranz, "Search and localization of a weak source with a multi-robot formation," *Journal of Intelligent & Robotic Systems*, vol. 97, no. 3-4, pp. 623–634, 2020.

[10] L. Briñón-Arranz, A. Renzaglia, and L. Schenato, “Multirobot symmetric formations for gradient and hessian estimation with application to source seeking,” *IEEE Transactions on Robotics*, vol. 35, no. 3, pp. 782–789, 2019.

[11] T. P. Nascimento and M. Saska, "Position and attitude control of multi-rotor aerial vehicles: A survey," *Annual Reviews in Control*, vol. 48, pp. 129–146, 2019.

[12] A. Okabe, *Spatial tessellations*. Wiley Online Library, 1992.

[13] J. Cortés, S. Martinez, T. Karatas, and F. Bullo, "Coverage control for mobile sensing networks," *IEEE Transactions on Robotics and Automation*, vol. 20, no. 2, pp. 243–255, 2004.

[14] A. Renzaglia, J. Dibangoye, V. Le Doze, and O. Simonin, “A common optimization framework for multi-robot exploration and coverage in 3D environments,” *Journal of Intelligent & Robotic Systems*, vol. 100, pp. 1453–1468, 2020.

[15] F. Pratissoli, B. Capelli, and L. Sabattini, "On coverage control for limited range multi-robot systems," in *IEEE/RSJ International Conference on Intelligent Robots and Systems*, 2022, pp. 9957–9963.

[16] S. Lloyd, “Least squares quantization in PCM,” *IEEE Transactions on Information Theory*, vol. 28, no. 2, 1982.

[17] Q. Du, M. D. Gunzburger, and L. Ju, "Constrained centroidal Voronoi tessellations for surfaces," *SIAM Journal on Scientific Computing*, vol. 24, no. 5, 2003.

[18] M. Deserno, "How to generate equidistributed points on the surface of a sphere," *If Polymerforshung (Ed)*, vol. 99, no. 2, 2004.