Quantifying Robotic Swarm Coverage

TL;DR

This paper introduces a new, continuous error metric for evaluating coverage in swarm robotics, along with benchmarks and theoretical analysis, to improve performance assessment of control algorithms.

Contribution

It presents a novel error metric for coverage evaluation, along with two benchmarks and theoretical insights, enhancing performance analysis in swarm robotics.

Findings

The error metric is continuously sensitive to swarm distribution changes.

Two benchmarks enable comparison of observed coverage to theoretical and random distributions.

The error metric obeys a central limit theorem, facilitating statistical analysis.

Abstract

In the field of swarm robotics, the design and implementation of spatial density control laws has received much attention, with less emphasis being placed on performance evaluation. This work fills that gap by introducing an error metric that provides a quantitative measure of coverage for use with any control scheme. The proposed error metric is continuously sensitive to changes in the swarm distribution, unlike commonly used discretization methods. We analyze the theoretical and computational properties of the error metric and propose two benchmarks to which error metric values can be compared. The first uses the realizable extrema of the error metric to compute the relative error of an observed swarm distribution. We also show that the error metric extrema can be used to help choose the swarm size and effective radius of each robot required to achieve a desired level of coverage. The second benchmark compares the observed distribution of error metric values to the probability density function of the error metric when robot positions are randomly sampled from the target distribution. We demonstrate the utility of this benchmark in assessing the performance of stochastic control algorithms. We prove that the error metric obeys a central limit theorem, develop a streamlined method for performing computations, and place the standard statistical tests used here on a firm theoretical footing. We provide rigorous theoretical development, computational methodologies, numerical examples, and MATLAB code for both benchmarks.

Tables

Table 1. Table 1: Summary of mean and standard deviation of error metric values for the example in Subsection 4.4.1 .

	Mean	Standard deviation
positions sampled from ${\rho }_{\text{ring}}$	0.4933	0.02484
steady state error values	0.5026	0.02586

Table 2. Table 2: Summary of statistical tests for the example in Subsection 4.4.1 .

F-statistic	1.0831
t-statistic	8.5888

Equations

$$G(z)=\frac{1}{2\pi }\exp \left(-\frac{|z{|}^2}{2}\right),$$

$$K^{\delta }(z)=\frac{1}{{\delta }^d}K\left(\frac{z}{\delta }\right).$$

$${\rho }_N^{\delta }(z;x_1,...,x_N)=\frac{{\sum }_{i=1}^N{G}^{\delta }(z-x_i)}{{\sum }_{i=1}^N{\int }_{\Omega }{G}^{\delta }(z-x_i)dz}.$$

$${\rho }_N^{\delta }(z)\approx \frac{1}{N}\sum _{i=1}^N{G}^{\delta }(z-x_i).$$

$$\underset{\delta \to 0}{\lim }{\rho }_N^{\delta }(z;x_1,...,x_N)=\frac{1}{N}\sum _{i=1}^N{\delta }_{x_i}(z)$$

$$e_N^{\delta }(x_1,...,x_N)={\int }_{\Omega }\left|{\rho }_N^{\delta }(z)-\rho (z)\right|dz.$$

$${\int }_{\Omega }\left|f-g\right|dz=2\underset{B\subset \Omega }{\sup }\left|{\int }_{B}fdz-{\int }_{B}gdz\right|.$$

$${\int }_{\Omega }fdz\le |\Omega |{\left({\int }_{\Omega }{f}^{2}dz\right)}^{1/2},$$

$${\int }_{\Omega }\left|\frac{1}{M}\sum _{j=1}^M{\rho }_N^{\delta }(z,t_j)-\rho (z)\right|dz$$

$${\hat{e}}_N^{\delta }={\int }_{{\Omega }^{-}}\left|{\rho }_N^{\delta }(z)-\rho (z)\right|dz,$$

$$\mu =\sum _{i=1}^M\left|{\int }_{{\Omega }_i}\rho (z)dz-\frac{N_i}{N}\right|.$$

$$\text{minimize }{e}_N^{\delta }(x_1,...,x_N),$$

$$\text{subject to }{x}_i\in \Omega \text{ for }i\in \{1,2,\dots ,N\}.$$

$$e_{\text{rel}}=\frac{e_{\text{observed}}-e^{-}}{e^{+}-e^{-}}.$$

$$\text{minimize }{e}_N^{\delta }(x_1,...,x_N),$$

$$\text{subject to }{x}_i\in \Omega \text{ for }i\in \{1,2,\dots ,N\},\delta >0.$$

$${\rho }_N^{\delta }(z)=\frac{1}{N}\sum _{i=1}^N{K}^{\delta }\left(z-X_i\right),$$

$${\rho }_N^{\delta }=K^{\delta }\ast \left(\frac{1}{N}\sum _{i=1}^N{\delta }_{X_i}\right),$$

$${\delta }^{\ast }(N)=C{N}^{-1/6},$$

$$\int zK(z)dz=0,\int z^{2}K(z)dz<\infty ,$$

$$\delta (N)=O(N^{-1/6})\text{ and }\underset{N\to \infty }{\lim }\delta (N)N^{1/4}=\infty ,$$

$$e_N^{\delta (N)}\approx \mathcal{N}\left(\frac{e(N)}{N^{1/2}},\frac{\sigma (N{)}^2\delta (N{)}^2}{N}\right),$$

$$F_{e_N^{\delta }}(z)={\int }_{{\Omega }^N}{1}_{\{x\mid e_N^{\delta }(x)\le z\}}\prod _{i=1}^N\rho (x_i)dx.$$

Full text

1 INTRODUCTION

Much of the research in swarm robotics has focused on determining control laws that elicit a desired group behavior from a swarm [6], and less attention has been placed on methods for evaluating the performance of these controllers quantitatively. Both [6] and [9] point out the lack of developed performance metrics for assessing and comparing swarm behavior, and [6] notes that existing performance metrics are often too specific to the task being studied to be useful in comparing performance across controllers.

This extended version of the authors’ paper [1] studies an error metric that evaluates one common desired swarm behavior: distributing the swarm according to a prescribed spatial density. Here, we delve deeper into the analysis of the properties of the error metric. Further, we include new examples, including a continuous target distribution to complement the original piecewise constant distribution, in order to better illustrate the utility of our methods.

In many applications of swarm robotics, the swarm must spread across a domain according to a target distribution in order to achieve its goal. Some examples are in surveillance and area coverage [8, 20, 23, 31], achieving a heterogeneous target distribution [16, 4, 14, 34, 17, 10, 41], and aggregation and pattern formation [36, 37, 30, 38, 28, 29]. Despite the importance of assessing performance, some studies such as [34, 38, 28, 29, 41] rely only on qualitative methods such as visual comparison. Others present performance metrics that are too specific to be used outside of the specific application, such as measuring cluster size in [36], distance to a pre-computed target location in [31, 10], and area coverage by tracking the path of each agent in [8]. Ergodicity has also been used as a notion of coverage error [26, 2], but it quantifies only the time-average statistics of the robot trajectories rather than the effective coverage of the swarm at an instant in time. In [30] an $L^{2}$ norm of the difference between the target and achieved swarm densities is considered, but the notion of achieved swarm density is particular to the controllers under study. These existing works do not provide methodology for carrying out computations with the notions they introduce. Moreover, they do not give guidance on how to interpret the values produced by the error metrics they define. Our work seeks to fill these gaps.

We develop and analyze an error metric that quantifies how well a swarm achieves a prescribed spatial distribution. Our method is independent of the controller used to generate the swarm distribution, and thus has the potential to be used in a diverse range of robotics applications. In [25, 40, 18], error metrics similar to the one presented here are used, but their properties are not discussed in sufficient detail for them to be widely adopted. In particular, although the error metric that we study always takes values somewhere between 0 and 2, these values are, in general, not achievable for an arbitrary desired distribution and a fixed number of robots. Therefore, one needs a better understanding of the finer properties of the error metric in order to judge whether its values (and hence the performance of the underlying controller) are “good” or not. We address this by studying two benchmarks,

the extrema of the error metric, and 2. 2.

the probability density function (PDF) of the error metric when robot positions are sampled from the target distribution,

which were first proposed in [25]. Using tools from nonlinear programming for (1) and rigorous probability results for (2), we put each of these benchmarks on a firm foundation. In addition, we provide MATLAB code for performing all calculations at https://git.io/v5ytf. Thus, by using the methods developed here, one can assess the performance of a given controller for any target distribution by comparing the error metric value of robot configurations produced by the controller against benchmarks (1) and (2).

Our paper is organized as follows. Our main definition, its basic properties, and a comparison to common discretization methods is presented in Section 2. Then, Section 3 and Section 4 are devoted to studying (1) and (2), respectively. We suggest future work in Section 5 and conclude in Section 6.

2 QUANTIFYING COVERAGE

One difficulty in quantifying swarm coverage is that the target density within the domain is often prescribed as a continuous (or piecewise continuous) function, yet the swarm is actually composed of a finite number of robots at discrete locations. The approach we take here is to use the robot positions to construct a continuous function that represents coverage (e.g. [24, 42, 2, 18, 29]). It is also possible to use a combination of the two methods, as in [11].

The method we present and analyze is inspired by vortex blob numerical methods for the Euler equation and the aggregation equation (see [12] and the references therein). There are also similarities between our approach and the numerical methods known as smoothed particle hydrodynamics (SPH), which have also been applied in robotics [28, 29, 41]. One main idea behind vortex methods and SPH is to use particles to approximate a continuous function that represents the fluid. We apply this idea but with the opposite aim: namely, we represent discrete points (the robots’ positions) with a continuous function. A similar strategy was alluded to in [18] and used in [25, 40] to measure the effectiveness of a certain robotic control law, but to our knowledge, our work here and in [1] is the first to develop any such method in a form sufficiently general for common use.

This section is devoted to our definition of the error metric and to its basic properties and computational considerations. We also provide a contrast between our method and another common way to measure error, namely, by discretizing the domain (e.g. [4, 14]) in Subsection 2.3. Finally, in Subsection 2.4, we present the setup for the two main examples that we will use throughout the paper.

2.1 Preliminaries

We are given a bounded region111We present our definitions for any number of dimensions $d\geq 1$ to demonstrate their generality. However, in the latter sections of the paper, we restrict ourselves to $d=2$, a common setting in ground-based applications. $\Omega\subset\mathbb{R}^{d}$, a desired robot distribution $\rho:\Omega\rightarrow(0,\infty)$ satisfying $\int_{\Omega}\rho(z)\,dz=1$, and $N$ robot positions $x_{1},\dots,x_{N}\in\Omega$. To compare the discrete points $x_{1},\dots,x_{N}$ to the function $\rho$, we place a “blob” of some shape and size at each point $x_{i}$. The blob or robot blob function can be any function $K:\mathbb{R}^{d}\rightarrow\mathbb{R}$ that is non-negative on $\Omega$ and satisfies $\int_{\mathbb{R}^{d}}K(z)dz=1$. In other contexts the blob function is referred to as a “kernel” [39, 28, 41]. Although more general blob shapes may be useful to represent properties of a robot’s sensing and/or manipulation capabilities, it is often natural to use a $K$ that is radially symmetric, decreases along radial directions, and enjoys certain integrability properties. We discuss these issues further in Subsection 4.1.2.

One choice of $K$ could be a scaled indicator function, for instance, a function of constant value within a disc of radius $1$ and [math] elsewhere. This is an appropriate choice when a robot is considered to either perform its task at a point or not, and where there is no notion of the degree of its effectiveness. For the remainder of this paper, however, we usually take $K$ to be the Gaussian

[TABLE]

which is useful when the robot is most effective at its task locally and to a lesser degree some distance away. A list of other common choices of $K$ and different ways of constructing multivariate blobs from single variable functions can be found in [39, Chapter 2, Chapter 4].

To formulate our definition, we need one more parameter, a positive number $\delta$, which we call blob radius. We then define $K^{\delta}$ as,

[TABLE]

We point out that this rescaling preserves the $L^{1}$ norm of $K$, so that $\int_{\mathbb{R}^{d}}K^{\delta}(z)\,dz=1$ holds for all $\delta>0$.

The shape of $K$ and the blob radius $\delta$ have two physical interpretations as:

• •

the area in which an individual robot performs its task, or

• •

inherent uncertainty in the robot’s position.

Either of these interpretations (or a combination of the two) can be invoked to make a meaningful choice of these parameters.

To define the swarm blob function $\rho_{N}^{\delta}$, we place a blob $G^{\delta}$ at each robot position $x_{i}$, sum over $i$ and renormalize, yielding,

[TABLE]

For brevity, we usually write $\rho_{N}^{\delta}(z)$ to mean $\rho_{N}^{\delta}(z;x_{1},...,x_{N})$. This swarm blob function gives a continuous representation of how the discrete robots are distributed. Note that each integral in the denominator of (2) approaches $1$ if $\delta$ is small or all robots are far from the boundary, so that we have,

[TABLE]

Moreover, $\rho_{N}^{\delta}(z;x_{1},...,x_{N})$ is a smoothed version of $\displaystyle\sum_{i=1}^{N}\delta_{x_{i}}(z)$, which represents the physical robot positions (here, $\delta_{x_{i}}(z)$ denotes the Dirac mass at $x_{i}$). Indeed,

[TABLE]

in the sense of distributions [7]. This point of view gives yet another interpretation of $\rho_{N}^{\delta}$.

We now introduce our notion of error:

Definition 2.1**.**

The error metric $e_{N}^{\delta}$ is defined as:

[TABLE]

This definition, including the definitions of $K^{\delta}$ and $\rho_{N}^{\delta}$, was introduced in this context by the authors in [1]. The relation (3) appeared there as well.

2.1.1 Remarks and Basic Properties

Our error is defined as the $L^{1}$ norm of the difference between the swarm blob function $\rho^{\delta}_{N}$ and the desired robot distribution $\rho$. One could use another $L^{p}$ norm; however, $p=1$ is a standard choice in applications that involve particle transportation and coverage such as [40, 18]. Moreover, the $L^{1}$ norm has a key property: for any two integrable functions $f$ and $g$,

[TABLE]

The other $L^{p}$ norms do not enjoy this property [15, Chapter 1]. Consequently, by measuring $L^{1}$ norm on $\Omega$, we are also bounding the error we make on any particular subset, and, moreover, knowing the error on “many” subsets gives an estimate of the total error. This means that by using the $L^{1}$ norm we capture the idea that discretizing the domain provides a measure of error, but avoid the pitfalls of discretization methods described in Subsection 2.3.

Studies in optimal control of swarms often use the $L^{2}$ norm due to the favorable inner product structure [40]. We point out that the $L^{1}$ norm is bounded from above by the $L^{2}$ norm: indeed, according to the Cauchy-Schwarz inequality, for any function $f$ we have,

[TABLE]

where $|\Omega|$ denotes the area of the bounded region $\Omega$. Thus, if an optimal control strategy controls the $L^{2}$ norm, then it will also control the error metric we present here.

Last, we note that for any $\Omega$, $\rho$, $\delta$, $N$, and $(x_{1},...,x_{N})$, we have $0\leq e_{N}^{\delta}\leq 2$. This was established in Proposition 2.1 of [1]. The theoretical minimum of $e^{\delta}_{N}$ can only be approached for a general target distribution when $\delta$ is small and $N$ is large, or in the trivial special case when the target distribution is exactly the sum of $N$ Gaussians of the given $\delta$, motivating the need to develop benchmarks (1) and (2).

2.1.2 Variants of the Error Metric

The notion of error in Definition 2.1 is suitable for tasks that require good instantaneous coverage. For tasks that involve tracking average coverage over some period of time (and in which the robot positions are functions of time $t$), an alternative “cumulative” version of the error metric is

[TABLE]

for time points $j=1,\dots,M$. This was defined in the authors’ [1], and is similar to the metric used in [40]. It may also be likened to the approach of [2], which expresses the average over time using an integral rather than finite sum and measures the difference from the target distribution using the Kullback-Leibler divergence and Bhattacharyya distance instead of an $L^{1}$ norm. Although the cumulative error metric (6) is, in general, distinct from the instantaneous version of (5), note that the extrema and PDF of this cumulative version can be calculated in the same way as the extrema and PDF of the instantaneous error metric with $MN$ robots. Therefore, in subsequent sections we restrict our attention to the extrema and PDF of the instantaneous formulation without loss of generality.

In addition, [40] considers a one-sided notion of error, in which a scarcity of robots is penalized but an excess is not, that is,

[TABLE]

where $\Omega^{-}:=\{z|\rho_{N}^{\delta}(z)\leq\rho(z)\}$. The definition of $\hat{e}_{N}^{\delta}$, which appears in [1], is particularly useful in conjunction with the choice of $K^{\delta}$ as a scaled indicator function, as $\hat{e}_{N}^{\delta}$ becomes a direct measure of the deficiency in coverage of a robotic swarm. For instance, given a swarm of surveillance robots, each with observational radius $\delta$, $\hat{e}_{N}^{\delta}$ is the percentage of the domain not observed by the swarm.222The notion of “coverage” in [8] might be interpreted as $\hat{e}_{N}^{\delta}$ with $\delta$ as the width of the robot. There, only the time to complete coverage ($t$ such that $\hat{e}_{N}^{\delta}(z;x_{1}(t),...,x_{n}(t))=0$) was considered.

Remarkably, $\hat{e}_{N}^{\delta}$ and $e_{N}^{\delta}$ are related by $e_{N}^{\delta}=2\hat{e}_{N}^{\delta}$. This was established by the authors in Proposition 2.2 of [1]. This relationship implies that $e_{N}^{\delta}$ enjoys the interpretation of being a measure of deficiency, and also allows the techniques introduced here to be directly applied to $\hat{e}_{N}^{\delta}$.

2.2 Calculating $e_{N}^{\delta}$

In practice, the integral in (5) can rarely be carried out analytically, primarily because the integral needs to be separated into regions for which the quantity $\rho_{N}^{\delta}(z)-\rho(z)$ is positive and regions for which it is negative, the boundaries between which are usually difficult to express in closed form. Hence we must rely on numerical integration, or quadrature. Here we study the magnitude $E_{m}$ of the difference between the true value of the error metric $e_{N}^{\delta}$ and an approximation $\tilde{e}_{N}^{\delta}$ as the number of grid points per axis $m$ increases for three elementary quadrature rules: the rectangle rule, the trapezoidal rule, and Simpson’s rule. As quadrature convergence rate proofs typically depend on smoothness of the integrand [13], we did not expect all of these rules to achieve their nominal rates of convergence. Indeed, in Figure 1 we see that the trapezoidal rule converges as $O(m^{-1.58})$ instead of $O(m^{-2})$ and Simpson’s rule converges as $O(m^{-1})$ instead of $O(m^{-4})$. As our applications do not require extremely accurate estimates of the error metric, other error metric values reported throughout the paper have been calculated using the rectangle rule with a moderate number of grid points. For applications that demand increased accuracy, we recommend an adaptive scheme such as that employed by MATLAB’s integral2 function.

2.3 The Pitfalls of Discretization

Before concluding this section, we analyze a measure of error that involves discretizing the domain and show that the values produced by this method are strongly dependent on a choice of discretization. In particular, this error approaches its theoretical minimum when the discretization is too coarse and its theoretical maximum when the discretization is too fine, regardless of robot positions.

Discretizing the domain means dividing $\Omega$ into $M$ disjoint regions $\Omega_{i}\subset\Omega$ such that $\bigcup_{i=1}^{M}\Omega_{i}=\Omega$. Within each region, the desired proportion of robots is the integral of the target density function within the region $\int_{\Omega_{i}}\rho(z)dz$. Using $N_{i}$ to denote the observed number of robots in $\Omega_{i}$, we can define an error metric as

[TABLE]

This exact definition appeared in the authors’ [1], and is based on discrete error metrics used in practice, for instance in [4, 14]. It is easy to check that $0\leq\mu\leq 2$ always holds. One advantage of this approach is that $\mu$ is very easy to compute, but there are two major drawbacks.

First, the choice for domain discretization is not unique, and this choice can dramatically affect the value of $\mu$. Indeed, if $M=1$ then $\mu=0$. This follows directly from the definition of $\mu$ and appears in [1] as Proposition 2.3. On the other hand, as the discretization becomes finer, the error approaches its maximum value 2, regardless of the robot positions. More precisely:

Proposition 2.2**.**

Suppose the robot positions are distinct333This is reasonable in practice as two physical robots cannot occupy the same point in space. In addition, the proof can be modified to produce the same result even if the robot positions coincide. and the regions $\Omega_{i}$ are sufficiently small such that, for each $i$, $\Omega_{i}$ contains at most one robot and $\int_{\Omega_{i}}\rho(z)\,dz\leq 1/N$ holds. Then $\mu\rightarrow 2$ as $|\Omega_{i}|\rightarrow 0$.

This proposition and its proof appear in [1] as Proposition 2.4.

Note that the shape of each region is also a choice that will affect the calculated value of $\mu$. Although our approach also requires the choice of some size and shape (namely, $\delta$ and $K$), these parameters have much more immediate physical interpretations, making appropriate choices easier to make.

The second, and perhaps more significant, drawback, is that by discretizing the domain we also discretize the range of values that the error metric can assume. This means that we have simultaneously desensitized the error metric to changes in robot distribution within each region. That is, so long as the number of robots $N_{i}$ within each region $\Omega_{i}$ does not change, the distribution of robots within any and all $\Omega_{i}$ may be changed arbitrarily without affecting the value of $\mu$. On the other hand, the error metric $e_{N}^{\delta}$ is continuously sensitive to differences in distribution.

2.4 Setup for Main Examples

We will perform computations with the following two desired distributions. We have purposefully chosen one that is piecewise continuous, $\rho_{\text{ring}}$, and one that is smooth, $\rho_{\text{ripple}}$.

Definition 2.3**.**

The ring distribution $\rho_{\emph{ring}}$ is defined on the Cartesian plane with coordinates $z=(z_{1}$, $z_{2})$ as follows. Let inner radius $r_{1}=11.4$in, outer radius $r_{2}=20.6$in, width $w=48$in, height $h=70$in, domain $\Omega=\{z:z_{1}\in[0,w],z_{2}\in[0,h]\}$, and region $\Gamma=\{z:r^{2}_{1}<(z_{1}-\frac{w}{2})^{2}+(z_{2}-\frac{h}{2})^{2}<r^{2}_{2}\}$. Define the non-normalized $\rho^{\prime}_{\emph{ring}}(z)$ to be $36\text{ if }z\in\Omega\cap\Gamma$ and $1\text{ if }z\in\Omega\setminus\Gamma$. Then let $C=\int_{\Omega}\rho^{\prime}_{\emph{ring}}(z)\,dz$ and $\rho_{\emph{ring}}(z)=C^{-1}\rho^{\prime}_{\emph{ring}}(z)$ for $z\in\Omega$.

Definition 2.4**.**

The ripple distribution $\rho_{\emph{ripple}}$ is defined on the Cartesian plane with coordinates $z=(z_{1}$, $z_{2})$ as follows. Let width $w=48$in, height $h=70$in, domain $\Omega=\{z:z_{1}\in[0,w],z_{2}\in[0,h]\}$, non-normalized $\rho^{\prime}_{\emph{ripple}}(z)=2+\sin[3\pi(z_{1}^{2}+z_{2}^{2})^{\frac{1}{2}}]+2\frac{z_{1}^{2}}{w^{2}}+\frac{z_{2}^{3}}{h^{3}}$, and normalization factor $C=\int_{\Omega}{\rho^{\prime}_{\emph{ripple}}(z)}\,dz$. Then $\rho_{\emph{ripple}}(z)=C^{-1}\rho^{\prime}_{\emph{ripple}}(z)$ for $z\in\Omega$.

3 ERROR METRIC EXTREMA

In the rest of the paper, we provide tools for determining whether or not the values of $e_{N}^{\delta}$ produced by a controller in a given situation are “good”. As mentioned in Section 2.1.1, it is not possible to achieve $e_{N}^{\delta}=0$ for every combination of target distribution $\rho$, number of robots $N$, and blob radius $\delta$ . Therefore, we would like to compare the achieved value of $e_{N}^{\delta}$ against its realizable extrema given $\rho$, $N$, and $\delta$. Unfortunately, $e_{N}^{\delta}$ is a nonlinear and moreover non-convex function of the robot positions $(x_{1},...,x_{N})$, and thus its extrema may elude analytical expression. Thus, we approach this problem with nonlinear programming.

3.1 Extrema Bounds via Nonlinear Programming

Let $x=(x_{1},...,x_{N})$ represent a vector of $N$ robot coordinates. The optimization problem, first introduced in the authors’ [1], is

[TABLE]

Note that the same problem structure can be used to find the maximum of the error metric by minimizing $-e_{N}^{\delta}$. Given $\rho$, $N$, and $\delta$, we approach these problems using a standard nonlinear programming solver, MATLAB’s fmincon.

A limitation of all general nonlinear programming algorithms is that successful termination produces only a local minimum, which is not guaranteed to be the global minimum. Since there is no obvious formulation of this problem for which a global solution is guaranteed, we use the local minimum as an upper bound for the global minimum of the error metric. Heuristics such as multi-start (running the optimization many times from several initial guesses and taking the minimum of the local minima) can be used to make this bound tighter. We use $e^{-}$ to denote the best local minimum found and $e^{+}$ to denote the best local maximum found. The bounds $e^{-}$ and $e^{+}$ serve as benchmarks against which we can compare a value of the error metric achieved by a given control law. This is reasonable, because if a configuration of robots with a lower value of the error metric exists but eludes numerical optimization, it is probably not a fair standard against which to compare the performance of a general controller.

Examples of extremal swarm blob functions found using this approach are shown in Figures 2, 3, and 4.

3.1.1 Relative Error

We introduce the notion of relative error, a quantitative way of assessing the performance of a robot distribution controller. This notion first appeared in the authors’ [1].

Definition 3.1**.**

Let $e_{\emph{observed}}$ denote the error value of a robot configuration. The relative error is defined as,

[TABLE]

In order to apply this definition, we need to specify how to find $e_{\text{observed}}$. If the robot positions $x_{1},...,x_{N}$ produced by a given controller are constant, then $e_{\text{observed}}$ can simply be taken as $e_{N}^{\delta}(x_{1},...,x_{N})$. In general, however, the positions $x_{1},...,x_{N}$ may change over time, and it is natural to define $e_{\text{observed}}$ based on the steady state value of $e_{N}^{\delta}$. We suggest defining the steady state settling time $t_{s}$ to be the time at which the error has settled to within $2\%$ of its asymptotic value. Then, we propose taking $e_{\text{observed}}$ to be the third-quartile value observed for $t>t_{s}$, which we denote $e_{\text{Q3}}$.

We suggest that if $e_{\text{rel}}$ is less than 10%, the performance of the controller is quite close to the best possible, and if this ratio is 30% or higher, the performance of the controller is rather poor.

We summarize the procedure described in this section:

• •

Step 1: Given $\rho$, $N$ and $\delta$, compute $e^{-}$ and $e^{+}$.

• •

Step 2: Compute $e_{\text{observed}}$ for the desired swarm controller. That is,

• –

Step 2.a: Calculate $e_{N}^{\delta}(t)$ from robot trajectories $x_{1}(t),...,x_{N}(t)$ produced by the controller.

• –

Step 2.b: Find the steady state settling time $t_{s}$.

• –

Step 2.c: Measure the third quartile value of the error metric $e_{Q_{3}}$ for $t>t_{s}$. This is $e_{\text{observed}}$.

• •

Step 3: Evaluate $e_{\text{rel}}$ according to Definition 3.1.

We emphasize that $e^{-}$ and $e^{+}$ are independent of the particular controller; therefore, Step 1 has to be completed only once when comparing different controllers or different outcomes of running one controller.

3.1.2 Example

We apply this method to assess the performance of the controller in [25], which guides a swarm of $N=200$ robots with $\delta=2$in (the physical radius of the robots) to achieve the ring distribution from Definition 2.3. These calculations were originally performed in the authors’ [1].

**Step 1: Compute $e^{-}$ and $e^{+}$. ** To determine an upper bound on the global minimum of the error metric, we computed 50 local minima of the error metric starting with random initial guesses, then took the lowest of these to be $e^{-}=0.28205$. An equivalent procedure bounds the global maximum as $e^{+}=1.9867$, produced when all robot positions coincide near a corner of the domain. The corresponding swarm blob functions are depicted in Figures 2 and 4. Note that the minimum of the error metric is significantly higher than zero for this finite number of robots of nonzero radius, emphasizing the importance of performing this benchmark calculation rather than using zero as a reference.

Step 2: Find $e_{\text{observed}}$. Under the stochastic control law of [25], the behavior of the error metric over time appears to be a noisy decaying exponential. Therefore, we fit to the data shown in Figure 7 of [25] a function of the form $f(t)=\alpha+\beta\exp(-\frac{t}{\tau})$ by finding error asymptote $\alpha$, error range $\beta$, and time constant $\tau$ that minimize the sum of squares of residuals between $f(t)$ and the data. By convention, the steady state settling time is taken to be $t_{\text{s}}=4\tau$, which can be interpreted as the time at which the error has settled to within $2\%$ of its asymptotic value [3]. The third quartile value of the error metric for $t>t_{\text{s}}$ is $e_{Q3}=0.5157$.

Step 3: Find $e_{\text{rel}}$. Using these values for $e_{\text{Q3}}$, $e^{-}$, and $e^{+}$, we calculate $e_{\text{rel}}$ according to Equation 9 as $13.71\%$.

Remark 3.2**.**

Although the sentiment of the $e_{\text{rel}}$ benchmark is found in [25], we have formalized the calculation and made three important improvements to make it suitable for general use. First, we have placed the calculation of $e^{-}$ and $e^{+}$ into a framework appropriate for nonlinear programming, so that the calculation is objective and repeatable. On the other hand, in [25], values analogous to $e^{-}$ and $e^{+}$ were found by “manual placement” of robots. Second, we suggest a general way of evaluating error in time-dependent controllers. In particular, although [25] refers to steady state, a definition is not suggested there. Adopting the 2% settling time convention allows for an unambiguous calculation of $e_{Q_{3}}$ and other steady state error metric statistics. It also provides a metric for assessing the speed with which the control law effects the desired distribution. Finally, we suggest using the third quartile value $e_{Q_{3}}$ in the calculation, in contrast to the minimum observed value of the error metric used in [25]. Since $e_{Q_{3}}$ better represents the distribution of error metric values achieved by a controller, our notion of relative error is more representative of the controller’s overall performance.

These changes account for the difference between our calculated value of $e_{\text{rel}}=13.71\%$ and the report in [25] that the error is “$7.2$% of the range between the minimum error value… and maximum error value”. Our substantially higher value of $e_{\text{rel}}$ indicates that the performance of this controller is not very close the best possible. We emphasize this to motivate the need for our second benchmark in Section 4, which may be more appropriate for a stochastic controller like the one presented in [25].

3.2 Error Metric for Optimal Swarm Design

So far we have taken $N$ and $\delta$ to be fixed; we have assumed that the robotic agents and size of the swarm have already been chosen. Now, we briefly consider the use of the error metric as an objective function for the design of a swarm. Adding $\delta>0$ as a decision variable to (8) yields the minimization problem,

[TABLE]

Solving (10) for several fixed values of $N$ provides insight into the number of robots and what effective working radius are needed to achieve a given level of coverage for a particular target distribution. Visualizations of such calculations are provided in Figure 5 and 6, and Figure 7 shows the optimal value of $\delta$ as a function of swarm size.

Note that “breakthroughs”, or relatively rapid decreases in the error metric, can occur once a critical number of robots is available; these correspond with a qualitative change in the distribution of robots. For example, at $N=22$ in Figure 5 the robots are arranged in a single ring; beginning with $N=25$ we see the robots start to be arranged in two separate concentric rings of different radii and the error metric begins to drop sharply. On a related note, there are also “lulls” in which increasing the number of robots has little effect on the minimum value of the error metric, such as between $N=44$ and $N=79$. We leave for future work the task of finding a theoretical explanation for these breakthroughs and lulls. In the meantime, computations like these can help a swarm designer determine the best number of robots $N$ and effective radius of each $\delta$ to achieve the required coverage.

4 ERROR METRIC PROBABILITY DENSITY FUNCTION

In the previous section, we described how to find bounds on the minimum and maximum values for error and use these as benchmarks against which to measure the performance of a given control law. However, stochastic control laws may tend to produce robot positions with observed error $e_{\text{observed}}$ well above the minimum $e^{-}$. For these controllers, $e_{\text{rel}}$ may be too stringent of a measure for assessing whether the controller’s performance is “good” or “bad”, and so we propose an alternative.

According to the setup of our problem, the goal of any control law is for the robots to achieve the desired distribution $\rho$. Thus, it is natural to compare the outcome of a given control law to simply picking the robot positions at random from the target distribution $\rho$. We take the robots’ positions $X_{1}$, …, $X_{N}$, to be independent, identically distributed bivariate random vectors in $\Omega\subset\mathbb{R}^{2}$ with probability density function $\rho$. We place a blob of shape $K$ at each of the $X_{i}$ (previously we took $K$ to be the Gaussian $G$), so that the swarm blob function is,

[TABLE]

where $K^{\delta}$ is defined by (1). Equation (11) appears in the authors’ [1]. We point out that the right-hand sides of (11) and (3) agree upon taking $K$ to be the Gaussian $G$ in (11) and the robot locations $x_{i}$ to be the randomly selected $X_{i}$ in (3). Moreover, we have the relationship,

[TABLE]

where $*$ denotes convolution. This equality is the reason (4) holds.

In this context, the error $e_{N}^{\delta}$ is a random variable, the value of which depends on the particular realization of the robot positions $X_{1}$, …, $X_{N}$, and thus has a well-defined probability density function (PDF) and cumulative distribution function (CDF). We denote the PDF and CDF by $f_{e_{N}^{\delta}}$ and $F_{e_{N}^{\delta}}$, respectively. The performance of a stochastic robot distribution controller can be quantitatively assessed by calculating the error values it produces in steady state and comparing their distribution to $f_{e_{N}^{\delta}}$. We introduce this benchmark in Subsection 4.4.

4.1 Kernel Density Estimation

The approach we take in this section is closely linked to the statistical theory of kernel density estimation (KDE). For an overview, see, for example, the book [32]. Broadly, the aim of this area of study is to find the underlying density $\rho$ from the values of identically distributed random variables sampled from this density. The expression (11) is the so-called kernel density estimator of $\rho$, and $\rho_{N}^{\delta}$ is considered as an approximation to $\rho$. In a sense, this is the reverse point of view from the one we take, since we are given $\rho$ and seek to learn about the distribution of the $X_{i}$s. Nevertheless, we are able to apply ideas from this well-developed field of study to our context.

There are many notions of error that are used in the context of KDE. We refer the reader to [39, Chapter 2] for a summary. Our notion $e_{N}^{\delta}$ corresponds to integrated absolute error (IAE). The use of IAE instead of other notions, particularly ones involving the $L^{2}$ norm, was advocated for in [15].

We conclude this subsection by applying two results from KDE theory to ideas we’ve presented so far. Then, in Subsection 4.2 we present rigorous results that show that the error metric has an approximately normal distribution when the $X_{i}$ are sampled from $\rho$. As a corollary we obtain that the limit of this error is zero as $N$ approaches infinity and $\delta$ approaches 0. Subsections 4.3 and 4.4.1 include a numerical demonstration of these results.

The theoretical results presented in the next subsection not only support our numerical findings, but they also allow for faster computation. Indeed, if one did not already know that the error when robots are sampled randomly from $\rho$ has a normal distribution for large $N$, tremendous computation may be needed to get an accurate estimate of this probability density function. On the other hand, since the results we present prove that the error metric has a normal distribution for large $N$, we need only fit a Gaussian function to the results of relatively little computation. In Subsection 4.4.1 we present an example calculation demonstrating the utility of this approach.

4.1.1 Optimal $\delta(N)$

In Subsection 3.2, we included $\delta$ as a parameter in the optimization problem (8) and numerically obtained an estimate for the optimal $\delta$ given a swarm size $N$. This problem has also been studied in the context of KDE. In [33] it was found that the optimal choice is,

[TABLE]

where $C$ is a constant that does not depend on $N$. This is equation (2.3) of [33] in the 2-dimensional case. Here, “optimal” means that this choice of $\delta$ minimizes the expectation of the $L^{1}$ error between $\rho$ and $\rho_{N}^{\delta}$. This result holds under the hypotheses that $K$ satisfies (13) and $\rho$ has bounded second derivatives.

We remark that the optimal $\delta(N)$ minimizing (10) that was found in Subsection 3.2 and shown in Figure 7 scales with $N^{p}$ where $p\approx-0.4\neq-1/6$. Despite the difference in exponent, these computations do not contradict the results of [33] as the latter considers the radius $\delta$ that minimizes the expected error value, and hence takes into account robot positions that are not optimal. On the other hand, the values of $\delta(N)$ computed in Subsection 3.2 are optimal only in concert with the optimal robot positions. The contrast between the two situations suggests that for the design of a swarm as described in Subsection 3.2, the sensing/manipulation radius required to maximize coverage may depend on whether the controller is stochastic or deterministic.

4.1.2 Choice of $K$

In the context of KDE, there is a notion of how efficient a given kernel $K$ is. This notion is formulated in terms of $L^{2}$ error between $\rho$ and $\rho_{N}^{\delta}$, and therefore the precise details do not directly apply to our situation. However, we point out that it is known [39, Section 2.7] that, so long as the kernel $K$ is non-negative, has integral 1, and satisfies

[TABLE]

then the particular choice of $K$ has only a marginal effect on efficiency. We leave to future work the extension of this concept to our context.

4.2 Theoretical Central Limit Theorem

It turns out that, under appropriate hypotheses, the $L^{1}$ error between $\rho$ and $\rho^{\delta}_{N}$ has a normal distribution with mean and variance that approach zero as $N$ approaches infinity. In other words, a central limit theorem holds for the error. The first such result was obtained in [15] in the one-dimensional case. Previously, similar results, such as those in [5], were available only for $L^{2}$ notions of error, which are easier to analyze.

For any result of central limit theorem type to hold, $\delta$ and $N$ have to be compatible. Thus, for the remainder of this subsection $\delta$ will depend on $N$, and we display this as $\delta(N)$. We have,

Theorem 4.1**.**

Suppose $\rho$ is twice continuously differentiable, and $K$ is zero outside of some bounded region and radially symmetric. Then, for $\delta(N)$ satisfying

[TABLE]

we have

[TABLE]

where $\sigma^{2}(N)$ and $e(N)$ are deterministic quantities that are bounded uniformly in $N$.444Here $\mathcal{N}(\mu,\sigma^{2})$ denotes the normal random variable of mean $\mu$ and variance $\sigma^{2}$, and we use the notation $\approx$ to mean that the difference of the quantity on the left-hand side and on the right-hand side converges to zero in the sense of distributions as $N\rightarrow\infty$.

From Theorem 4.1 it is easy to deduce:

Corollary 4.2**.**

Under the hypotheses of Theorem 4.1, the error $e_{N}^{\delta(N)}$ converges in distribution to zero as $N\rightarrow\infty$.

Corollary 4.2, Theorem 4.1, and the proof of Theorem 4.1 (which follows from [22, Theorem, page 1935]) appear in the authors’ [1].

Remark 4.3**.**

There are a few ways in which practical situations may not align perfectly with the assumptions of Theorem 4.1. However, we posit that in all of these cases, the differences are numerically insignificant. Moreover, we believe that a stronger version of Theorem 4.1 may hold, one which applies in our situation. We now briefly summarize these three discrepancies and indicate how to resolve them.

First, we defined our density $\rho_{N}^{\delta}$ by (2), but in this section we use a version with denominator $N$. However, as explained above, the two expressions approach each other for small $\delta$, and this is the situation we are interested in here. We believe that the result of Theorem 4.1 should hold with (2) instead of (11).

Second, in one of the two examples we consider here, the desired density $\rho$ is only piecewise continuous but not twice differentiable. We point out that an arbitrary density $\rho$ may be approximated to arbitrary precision by a smoothed out version, for example by convolution with a mollifier (a standard reference is [7, Section 4.4]). Moreover, the main results of [19] give a version of Theorem 4.1 that require that $\rho$ is only Lebesgue measurable (piecewise continuous functions satisfy this hypothesis). However, the results of [19] only apply in 1 dimension. Generalizing the results of [19] to several spatial dimensions is an interesting open problem in KDE, but is beyond the scope of our work here.

Third, in our computations we use the kernel $G$, which is not compactly supported, for the sake of simplicity. Similarly, this kernel can be approximated, with arbitrary accuracy, by a compactly supported version. Making these changes to the kernel or target density would not affect the conclusions of numerical results.

4.3 Numerical Approximation of the Error Metric PDF

In this subsection we describe how to numerically find $f_{e_{N}^{\delta}}$ and $F_{e_{N}^{\delta}}$. According to Theorem 4.1, for sufficiently large $N$, one could simply use random sampling to estimate the mean and standard deviation, then take these as the parameters of the normal PDF and CDF. However, for moderate $N$, we choose to begin by estimating the entire CDF and confirming that it is approximately normal. To this end, we apply:

Proposition 4.4**.**

We have,

[TABLE]

Here $\boldsymbol{1}$ denotes the indicator function. Proposition 4.4 and its proof appear in the authors’ [1].

Notice that, since each of the $x_{i}$ is itself a 2-dimensional vector as the $X_{i}$ are random points in the plane, the integral defining the cumulative distribution function of the error metric is of dimension $2N$. Finding analytical representations for the CDF quickly becomes infeasible for large swarms. Therefore, we approximate (15) using Monte Carlo integration, which is well-suited for high-dimensional integration [35]555Quasi-Monte Carlo techniques, which use a low-discrepancy sequence rather than truly random evaluation points, promise somewhat faster convergence but require considerably greater effort to implement. The difficulty is in generating a low-discrepancy sequence from the desired distribution, which is possible using the Hlawka-Mück method, but computationally expensive [21].. Next, we fit a Gauss error function ($\operatorname{erf}(\cdot)$, the integral of a Gaussian $G$) to the data. If the fitted curve matches the data well, we differentiate to obtain the PDF. We remark that we express the integral in (15) in terms of an indicator function in order to express the quantity of interest in a way that is easily approximated with Monte Carlo integration.

4.3.1 Computation

We apply Proposition 4.4 to numerically find the PDF $f_{e^{\delta}_{N}}$ and the CDF $F_{e^{\delta}_{N}}$ for $\rho=\rho_{\text{ring}}$ (see Definition 2.3), $N=200$, and $\delta=2\text{in}$.

We approximate $F_{e_{N}^{\delta}}$ using $M=1000$ Monte Carlo evaluation points; this is shown by a solid gray line in Figure 8. The numerical approximation appears to closely match a Gauss error function as theory predicts. Therefore an analytical $\operatorname{erf}(\cdot)$ curve, represented by the dashed line, is fit to the data using MATLAB’s least squares curve fitting routine lsqcurvefit. To obtain $f_{e_{N}^{\delta}}$, the analytical curve fit for $F_{e_{N}^{\delta}}$ is differentiated, and the result is also shown in Figure 8.

In addition, we calculate the mean and standard deviation of $e_{N}^{\delta}$ to be, respectively, 0.4933 and 0.02484.

4.4 Benchmark for Stochastic Controllers

In this subsection we describe how to use the tools developed so far to assess the performance of a given stochastic controller for a given desired distribution $\rho$, number of robots $N$, and blob radius $\delta$. We then demonstrate the method via an example.

We propose using two standard statistical tests, the two-sample t- and F- tests [27], to assess whether the performance of the control law is comparable to sampling robot positions from the target distribution. More precisely, these tests provide criteria for the rejection of the null hypothesis that the mean and variance of the steady state values of the error metric produced by a given control law are indistinguishable from mean and variance of $f_{e_{N}^{\delta}}$. This is summarized in the following procedure, which assumes that $f_{e_{N}^{\delta}}$ is approximately normal. This assumption is reasonable due to Theorem 4.1.

• •

Step 1: Numerically find the mean and variance of $f_{e^{\delta}_{N}}$ using either of the methods described in Subsection 4.3.

• •

Step 2: Find the mean and variance of the steady state error metric values produced by the controller.

• •

Step 3: Perform two-sample t- and F-tests to compare the means and variances, respectively.

• •

Step 4: Assess the results.

• –

Step 4.a: If the tests fail to refute the null hypotheses, the performance of the controller is consistent with sampling robot positions from the target distribution.

• –

Step 4.b: If the two-sample t-test suggests that the means are significantly different, find the 95% confidence interval of their difference to assess the magnitude of the controller’s performance surplus or deficiency.

Just as for the benchmark $e_{\text{rel}}$ defined in Section 3, the first step is independent of choice of controller. Thus, if the goal is to compare several controllers, the computations of Step 1 need to be performed only once.

4.4.1 Example

We apply this method to assess the performance of the controller in [25], where the desired distribution is $\rho=\rho_{\text{ring}}$ (see Definition 2.3), there are $N=200$ robots, and the radius is $\delta=2\text{in}$. These calculations were originally performed in the authors’ [1].

Step 1, the computation of the mean and variance of $e_{N}^{\delta}$ when robot positions are sampled at random from $\rho_{\text{ring}}$, was completed in Subsection 4.3.1. For Step 2, we use the data presented as Figure 7 of [25] to calculate that the distribution of steady state error metric values produced by their controller has a mean of $0.5026$ with a standard deviation of $0.02586$. These values are summarized in Table 1.

Next we perform the two-sample F-test and two-sample t-test. We take the null hypothesis to be that the distribution of these error metric values is the same as $f_{e_{N}^{\delta}}$. The results are summarized in Table 2.

The F-statistic of 1.0831 fails to refute the null hypothesis, indicating that there is no significant difference in the standard deviations. On the other hand, a two-sample t-test rejects the null hypothesis with a t-statistic of 8.5888, indicating that the steady state error is not distributed with the same population mean as $f_{e_{N}^{\delta}}$. Nonetheless, the 95% confidence interval for the true difference between population means is computed to be merely $(0.00717,0.01141)$. This shows that the mean steady state error achieved by this controller is unlikely to exceed that of $f_{e_{N}^{\delta}}$ by more than 2.31%. Therefore, we find the performance of the controller in [25] to be acceptable given its stochastic nature, as the error metric values it realizes are only slightly different from those produced by sampling robot positions from the target distribution.

Remark 4.5**.**

As with $e_{\text{rel}}$ of Section 3, the sentiment of this benchmark is preceded by [25]. However, here we have presented a concrete, repeatable, and objective approach that makes this benchmark suitable for general use. In particular, we have introduced into this context the use of appropriate statistical tests for comparing two approximately normal distributions. On the other hand, [25] relies on visual inspection, noting, “the error values [from simulation] mostly lie between …the 25th and 75th percentile error values when robot configurations are randomly sampled from the target distribution”. The fact that $f_{e_{N}^{\delta}}$ is approximately Gaussian (according to Theorem 4.1) not only allows for the use of these statistical tests, but also provides a way to greatly speed up computation. Indeed, the calculations in [25]666According to the caption of Figure 2 of [25], the figure was generated as a histogram from 100,000 Monte Carlo samples. took two orders of magnitude more computation than we perform in Subsection 4.3.1.

5 FUTURE WORK

Several open problems remain regarding the computational and theoretical aspects of our proposed performance assessment methods. For instance, the benchmark calculations scale quickly with the size of the swarm, and therefore they become unfeasible for a standard personal computer with swarms on the order of 1000 robots. To extend our proposed methods to larger swarms, computational techniques should be improved. Possibilities include the derivation of an analytical representation of the error metric PDF in terms of the target distribution, more sophisticated optimization algorithms to compute the extrema, and alternate formulations of the optimization problem that reduce the number of variables or introduce guarantees on the global optimality of a resulting solution. As for theoretical extensions, analysis of the shape and size of the robot blob function, as well as a rigorous understanding of the optimal radius, are intriguing topics for further study. Specifically, if a certain control law were characterized as “good” with one choice of blob function, would the same conclusion be reached with a different blob function? These questions were raised by the authors’ [1], and here we provide possible routes to answering them, as well as connections to known results in the theory (see Subsection 4.1).

6 CONCLUSION

This work deals with the performance assessment of spatial density control of robotic swarms. We introduce a sensitive error metric and two corresponding benchmarks, one based on the realizable extrema of the error metric and one based on its probability density function, which provide methods for evaluating the performance of control algorithms. The parameters of the error metric correspond to physical properties such as robot shape and size, and can therefore be tailored to the particular system under study. The benchmark calculations then provide reference points to which the performance of different control schemes can be compared. The error metric and these benchmarks were first carefully defined in the authors’ [1]; here, we have made these definitions more transparent, and, in Subsections 3.1 and 4.4, we have added step-by-step methods for performing the relevant computations. This will allow swarm designers and practitioners to more easily use these benchmarks to quantitatively decide which control scheme best fits their application. As demonstrated by the examples, including the new example featured in Subsections 3.1 and 3.2, these evaluation methods work well for a variety of target distributions and can be applied to judge the performance of both deterministic and stochastic control algorithms.

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the support of NSF grant CMMI-1435709, NSF grant DMS-1502253, the Dept. of Mathematics at UCLA, and the Dept. of Mathematics at Harvey Mudd College. OT acknowledges support from the Charles Simonyi Endowment at the Institute for Advanced Study. The authors also thank Hao Li (UCLA) for his insightful contributions regarding the connection between this error metric and kernel density estimation.