# Stein Coverage: a Variational Inference Approach to Distribution-matching
Multisensor Deployment

Donipolo Ghimire and Solmaz S. Kia, _Senior Member, IEEE_ The authors are with
the Department of Mechanical and Aerospace Engineering, University of
California, Irvine, CA 92697<EMAIL_ADDRESS>This work was
supported by NSF award ECCS 1653838.

###### Abstract

This paper examines the spatial coverage optimization problem for multiple
sensors in a known convex environment, where the coverage service of each
sensor is heterogeneous and anisotropic. We introduce the Stein Coverage
algorithm, a distribution-matching coverage approach that aims to place
sensors at positions and orientations such that their collective coverage
distribution is as close as possible to the event distribution. To select the
most important representative points from the coverage event distribution,
Stein Coverage utilizes the Stein Variational Gradient Descent (SVGD), a
deterministic sampling method from the variational inference literature. An
innovation in our work is the introduction of a repulsive force between the
samples in the SVGD algorithm to spread the samples and avoid footprint
overlap for the deployed sensors. After pinpointing the points of interest for
deployment, Stein Coverage solves the multisensor assignment problem using a
bipartite optimal matching process. Simulations demonstrate the advantages of
the Stein Coverage method compared to conventional Voronoi partitioning
multisensor deployment methods.

###### Index Terms:

Optimization, multisensor deployment, coverage, variational inference, optimal
transport.

##  I Introduction

In recent years, the decrease in the cost of embedded microsensing and
wireless communication technologies has led to increased use of wireless
sensor networks for tasks such as surveillance, data collection, situational
awareness, improved wireless coverage, and environmental monitoring. However,
the operational cost is still prohibitive in deploying an unlimited number of
sensors to achieve full coverage. The coverage problem involves strategically
placing a limited number of sensors across the area of interest in such a way
that the spatial distribution of the sensors can capture all the essential
environmental features. This paper examines a multisensor deployment problem
for coverage over an event distribution in a known convex environment, where
the coverage service of each sensor is heterogeneous and anisotropic, see Fig.
1. The aim is to have the final configuration of the sensors achieve a
topological distribution that is similar to the spatial distribution of the
events induced by the targets on the ground.

Figure 1: An example of a multisensor deployment with a limited number of
sensors available. The sensors are heterogenous represented by the different
sizes and colors. The aim is to place the sensors in positions that will yield
the highest benefit, that is, providing service to more densely populated
areas.

Literature Review: This research focuses on a problem in the field of sensor
deployment for coverage [1, 2, 3, 4, 5]. In coverage problems with a limited
number of sensors, the essence of the proposed solutions is to divide the area
into subregions and assign a sensor to each one, with the goal of maximizing a
certain coverage measure. This problem is related to many other problems, such
as various resource allocation problems and the facility location problem, and
is considered NP-Hard [6, 1].

Among the proposed multisensor coverage strategies, geometric approaches like
Voronoi diagrams, which use expected sensing cost as the Euclidean distance
between sensor position and event location [7, 4] uniquely combines both
deployment and allocation in a distributed manner using gradient descent. This
approach decouples the problem of optimal sensor allocation from the problem
of sensor configuration, as each sensor only needs to consider its own Voronoi
cell and its local neighbors. Beyond partitioning, the Voronoi-based approach
can consider two additional components, such as the sensing model of an agent
and the map over each point in the coverage region, known as the event density
function. A Voronoi-based approach can be used to model the heterogeneity of
sensors using the multiplicative weighted version of Voronoi diagrams(power
diagrams) [8, 9, 10]. On the other hand, the density functions can be used to
bias the importance of the areas for coverage. Although these Voronoi-based
approaches have nice properties such as convexity and dual triangulability,
power diagrams may have empty cells associated with, and the locally optimal
solutions provided by these methods strongly depend on proper initialization
[11, 12]. Another downside of Voronoi partitioning methods is that the cells
created by partitioning extend beyond the sensors’ footprint, implicitly
considering a coverage larger than the immediate sensing range. Therefore, the
sensors do not necessarily get assigned to areas of importance that fall in
the actual sensor footprint of the sensors.

In recent years, an alternative deployment approach has been proposed in the
literature that aims to match the collective quality of service (QoS) of the
sensors with the density distribution of the coverage area [1, 13, 14, 15,
16]. This approach considers a more realistic model for the QoS of the sensors
as a spatial distribution. For example [1] proposed a soft and hard sensor
allocation framework using a statistical coverage quality measure based on
distribution matching using f-divergence between the event detection
probability and spatial sensing model of sensors. Alternatively, [15] and [16]
proposed to cluster the area of interest using techniques such as Gaussian
mixture model (GMM) and K-means clustering, respectively, to identify the
potential deployment points based on the spatial distribution of the coverage
subjects. Then, cast the sensor deployment problem as an optimal assignment
problem.

The deployment algorithm proposed in this paper falls into the category of
distribution-matching multisensor deployment methods. We present a novel
approach to determining the most suitable locations for sensor deployment by
utilizing a variational inference technique. Variational inference allows
drawing ‘super samples’ from the density distribution of the coverage
subjects. These samples in the statistical sense are the closest
representatives of the density distribution, thus capturing the most important
points in the coverage area. In particular, we use the Stein Variational
Gradient Descent (SVGD) method [17] to draw the ‘super samples’. Compared with
other sampling methods such as Monte Carlo methods, SVGD is known to achieve
good approximation even with a very small number of samples (particles). SVGD
is a deterministic gradient-based sampling algorithm for approximate
inference. Given a probability density function $p(x)$, SVGD finds a set of
samples (particles) to approximate $p(x)$ by a simple iterative update that
resembles a steepest descent gradient method; for more see Section III.
However, our approach is not the immediate application of the SVGD algorithm.
The sampling methods are oblivious to the coverage footprint of the sensors
and can result in finding the candidate deployment points that will result in
significant overlap in the coverage provided by the sensors as shown in Fig.
2. The overlapping multisensor coverage is not only a waste of coverage
resources but also can result in undesirable interference between the sensors’
services as in wireless coverage problems. To consider the effective coverage
footprint of the sensors, we include a repulsive force between the samples in
the SVGD algorithm to spread the samples to avoid overlap for the deployed
sensors. The result is a pruned region of space, where the set of points
determines the regions of interest and helps to create a utility map of the
environment with the event distribution. This approach significantly reduces
our search space for multisensor deployment. Once we generate the low-cost and
constraint-satisfying utility map, we frame our problem as an assignment
problem and solve the problem using the Hungarian Algorithm. We refer to our
proposed deployment algorithm as Stein Coverage. We demonstrate the
performance of our algorithm for heterogeneous sensors with both isotropic and
anisotropic footprints.

_Notation_ : We let ${\mathbb{R}}$, $\mathbb{Z}$, and $\mathbb{N}$ denote the
set of reals, integers, and natural numbers, respectively. The positive and
nonnegative subspaces of these number spaces are indicated by indexes $>0$ or
$\geq 0$, respectively. For any closed and compact set
$\mathcal{C}\subset^{2}$, $\partial\mathcal{C}$ is its boundary.

##  II Problem Definition

We are considering a time-invariant event distribution function that is known
a priori, $p(x):\mathcal{W}\rightarrow\mathbb{R}_{\geq 0}$, which represents
the probability that some event takes place in a two-dimensional space
$\mathcal{W}\subset^{2}$. The probability density represents a finite set of
$m$ targets that are densely scattered over the planar space $\mathcal{W}$.
These targets can be real such as humans or animals, information sources,
pollution spills, forest fires, etc. Suppose
$\mathcal{A}=\\{1,\cdots,m\\}\subset\mathbb{N}$, represents the set of $m$
sensors that can be identical (homogeneous) or non-identical (heterogeneous).
Every sensor $i\in\mathcal{A}$ provides a coverage service, e.g. wireless
coverage, data harvesting, or event detection. The quality of service of each
sensor is described by a spatial coverage distribution
$s_{i}(x|x_{i},\theta_{i})$ where $x_{i}$ is the position and $\theta_{i}$ is
the sensor orientation. The support for this distribution can be
$\mathbb{R}^{2}$, but we assume that there is an effective coverage footprint
$\mathcal{C}_{i}(x|x_{i},\theta_{i})\subset\mathbb{R}^{2}$ around any point
$x_{i}$ at which the sensor is placed. $\mathcal{C}_{i}(x|x_{i},\theta_{i})$
and $s_{i}(x|x_{i},\theta_{i})$ can be considered as rigid bodies.

###### Assumption 1

There exists $\mathsf{R}\in\mathbb{R}_{>0}$ such that
$\|x-x_{i}\|\leq\mathsf{R}$ for any
$x\in\partial\mathcal{C}_{i}(x|x_{i},\theta_{i})$ for all $i\in\mathcal{A}$.

In the statistical sense, the ultimate goal for optimal coverage is finding
$(x_{i},\theta_{i})$ configuration for each sensor $i\in\mathcal{A}$ such that
the collective distribution of
$\frac{1}{m}\sum_{i\in\mathcal{A}}s_{i}(x|x_{i},\theta_{i})$ of the team is as
similar as possible to the distribution of $p(x)$. The similarity between the
distributions can be measured using various statistical measures, e.g.,
Kullback-Leibler (KL) divergence [18]. However, this optimization problem is
highly nonlinear and computationally challenging. Therefore, the literature
resorts to solving the problem under various assumptions, such as an isotropic
or uniform coverage distribution for sensors or seeking suboptimal solutions.

When the sensors are homogeneous and their coverage distribution is uniform
over their effective coverage range, one of the well-known frameworks to solve
the multisensor coverage problem is to assign each sensor to a weighted
centroid point of Voronoi partition
$\mathcal{V}(.)=\\{\mathcal{V}_{1},\mathcal{V}_{2},\cdots\mathcal{V}_{N}\\}$
or power diagrams
$\mathcal{P}(.)=\\{\mathcal{P}_{1},\mathcal{P}_{2},\cdots,\mathcal{P}_{N}\\}$
obtained by using a continuous-time version of Lloyd’s algorithm. This method
drives the weighted centroid points to observe events distributed according to
$p(x)$ by minimizing a locational optimization cost function based on a
Euclidean metric. This continuous optimization can be very costly when the
domain is very large and the decision variables are high dimensional. There is
also no guarantee of convergence to the global minimum of the coverage merit
function [1, 8].

Our goal is to deploy the sensors in a distribution-matching fashion at points
in the coverage area where their collective service distribution best matches
the distribution of the coverage subjects. To do this, we aim to obtain a set
of discrete samples from the most informative regions using a functional
gradient method while accounting for overlapping coverage between sensors.
SVGD, Kernel herding, and Markov Chain Monte Carlo (MCMC) are some techniques
that can be employed to facilitate the important sample selection process. Our
method uses the Kernelized Stein discrepancy, a new discrepancy statistic for
measuring differences between two probability distributions, that enables an
inference problem, which looks for

$\displaystyle\vspace{-2.5 in}\min_{q}\,\mathcal{KL}(q(x)\,\|\,p(x))\quad
s.t.\quad q(g(x)\geq 0)=1$ (1)

where $q$ is the empirical distribution of a set of points generated from an
arbitrary (Gaussian, Uniform, Beta, etc) distribution, and $g(x)$ is the soft
constraint imposed on the sensor locations. Our method is a spatially
regulated SVGD that introduces a metric that encourages the generated sensor
locations to maintain a desired minimum distance between each other. This
ensures that the sensors are strategically placed and collectively cover the
dense and important regions of the environment. Similarly, for heterogeneous
sensors, these final sensor configurations resolve an assignment problem,
where the assignment cost for each sensor is determined by the dissimilarity
of the coverage utility and the distribution of the targets near the potential
assignment location.

##  III A brief overview of Stein Variational Gradient Descent Preliminary
material

In this section, for the convenience of the reader and also to introduce the
notation and definitions that we use in developing our algorithm, we give a
brief overview of the Stein variational gradient descent (SVGD) method
following [17] and [19]. SVGD is a deterministic sampling algorithm that
iteratively transports a set of samples (particles) to approximate a given
distribution $p(x)$ on an open set $X\subset\mathbb{R}^{d}$, based on a
gradient-based update that guarantees optimally decreasing the KL divergence
within a function space [19]. Note that we will be using samples and particles
interchangeably. To achieve this, we initialize the samples with some simple
distribution $q(x)$, and update them via map

$T(x)=x+\epsilon\,\phi(x),$ (2)

where $\epsilon$ is a small step size, and $\phi(x)$ is a perturbation
direction, or velocity field, which should be chosen to maximally decrease the
KL divergence of the sample distribution with the target distribution $p(x)$;
this is framed by [17] as solving a functional optimization problem. Given
some structural properties of this optimization problem, which allows drawing
connection to ideas in the Stein’s method [20] used for proving limit theorems
or probabilistic bounds in theoretical statistics, Liu and Wang [17] showed
that the functional optimization problem that gives us $\phi(x)$ reduces to

$\displaystyle\mathbb{D}(q\|p):=\max_{\phi\in\mathcal{H}}\\{\mathbb{E}_{q}[\mathcal{S}_{p}\phi]\leavevmode\nobreak\
\leavevmode\nobreak\ \text{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\
\|\phi\|_{\mathcal{H}}\leq 1\\}$ (3)

where $\mathcal{H}$ is a normed function space chosen to optimize over. Here,
$\mathbb{D}(q\|p)$ is called Stein discrepancy, which provides a discrepancy
measure between $q$ and $p$, with the proven properties that
$\mathbb{D}(q\|p)=0$ if $p=q$ and $\mathbb{D}(q\|p)>0$ if $q\neq p$ given
$\mathcal{H}$ is sufficiently large. Furthermore,

$\mathcal{S}_{p}\phi(x):=\nabla\log p(x)^{\top}\phi(x)+\nabla\cdot\phi(x),$
(4)

where $\nabla\cdot\phi(x)=\sum_{k=1}^{d}\partial_{k}\phi_{k}(x)$, and
$\mathcal{S}_{p}$ is a linear operator that maps a vector-valued function
$\phi(x)$ to a scalar-valued function $\mathcal{S}_{p}\phi(x)$, and
$\mathcal{S}_{p}$ is called the Stein operator in connection with the so-
called Stein’s identity, which shows that
$\mathbb{E}_{q}[\mathcal{S}_{p}\phi]=0$ if $q=p$.

Algorithm 1 Stein Variational Gradient Descent

Input: The score function $\nabla_{x}\log p(x)$.

Goal: A set of particles $\\{x_{i}\\}_{i=1}^{n}$ that approximates $p(x)$.

Initialize a set of particles $\\{x_{i}^{(0)}\\}_{i=1}^{n}$; choose a positive
definite kernel $k(x,x^{\prime})$ and step-size.
For iteration l do

$\displaystyle x_{i}^{(\ell+1)}\leftarrow
x_{i}^{(\ell)}+\phi^{\star}(x_{i}^{(\ell)})\quad\forall\,i\in\\{1,\ldots,n\\}$

where,

$\displaystyle\phi^{\star}(x)=\frac{1}{n}\\!\sum\nolimits_{j=1}^{n}\\!\left[\nabla\log
p(x_{j})k(x_{j},x)\\!+\\!\nabla_{x_{j}}k(x_{j},x)\right].$ (5)

The optimization problem (3) is an infinite-dimensional functional
optimization, which searches for the transformation function $\phi$ in (2)
that maps samples drawn from $q$ to $p$. It is critical to select a space
$\mathcal{H}$ that is both sufficiently rich and also ensures computational
tractability in practice. Literature has shown that the Kernelized Stein
discrepancy (KSD) which selects a Reproducing kernel Hilbert space (RKHS) to
be $\mathcal{H}$, not only meets this expectation but also goes one step
further by yielding a _closed-form_ solution [17].

To demonstrate this closed-form solution, let $\mathcal{H}_{0}$ be a RKHS of
scalar-valued functions with a positive definite kernel $k(x,x^{\prime})$, and
$\mathcal{H}=\mathcal{H}_{0}\times\cdots\times\mathcal{H}_{0}$ the
corresponding $d\times 1$ vector-valued RKHS. Then it can be shown that the
optimal solution of (3), denoted by $\phi^{\star}_{q,p}$, is

$\displaystyle\phi^{\star}_{q,p}(.)\propto\mathbb{E}_{x\sim q}[\nabla\log
p(x)\,k(x,.)+\nabla_{x}k(x,.)]$ (6)

The key advantage of KSD is its computational tractability: it can be
empirically evaluated with samples drawn from $q$ and the gradient $\nabla\log
p$, which is independent of the normalization constant in $p$, see [17]. In
SVGD algorithm, we iteratively update a set of samples using the optimal
transform just derived, starting from a certain initialization. Let
$\\{x^{i}_{\ell}\\}^{n}_{i=1}$ be the samples at the $\ell$-th iteration. In
this case, the exact distributions of $\\{x^{i}_{\ell}\\}^{n}_{i=1}$ are
unknown or difficult to keep track of, but can be best approximated by their
empirical measure
$\hat{q}_{\ell}^{n}(\text{d}x)=\sum_{i}\delta(x-x_{\ell}^{i})\text{d}x/n$.
Therefore, it is natural to think that $\phi^{\star}_{\hat{q}^{n}_{\ell},p}$,
with $q$ in (6) replaced by $\hat{q}^{n}_{\ell}$, provides the best update
direction to move the samples (and equivalently $\hat{q}^{n}_{\ell}$ ) “closer
to" $p(x)$. Implementing this update (5) iteratively, we get the main SVGD
algorithm in Algorithm 1.

##  IV Stein Coverage: a variational inference approach to multisensor
deployment

To deploy the robot team $\mathcal{A}$, our goal is to identify representative
samples $n\geq N$, called points of interest (PoIs), from the spatial
distribution of events $p(x)$. Once the PoIs are identified, we propose a
mechanism that assigns PoIs to the sensors in accordance with their
heterogeneous spatial coverage ability. We call our proposed deployment
algorithm Stein Coverage.

Figure 2: The figure on the left displays the traditional SVGD approach, where
the red dots represent the samples placed on a given Gaussian Mixture Model.
The dotted circles represent the sensor footprints that overlap when the
sensors are installed at the sample points drawn by SVGD. The figure on the
right, however, shows the sample placement using our proposed spatially
regulated SVGD algorithm, which reduces the overlap of the sensor footprints
when placed on the sample points.

###  IV-A Designing PoIs for robot deployment

Let $\mathcal{S}=\\{1,\cdots,n\\}\subset\mathbb{N}$ represent the label set of
PoIs, and $\mathbf{x}:=(x_{1},x_{2},\cdots,x_{N})\subset\mathcal{W}^{n}$ be
the location of PoIs. The SVGD algorithm can be used to draw these PoIs, but
the naive implementation of the SVGD is unaware of the coverage footprint of
the sensors, as illustrated in Fig. 2. This can lead to overlapping coverage,
which is an inefficient use of the limited coverage resources.

Given the effective coverage footprint of the sensors, under Assumption 1, if
any two samples are placed in a way that

$\displaystyle g=\mathbb{E}[\left\|x_{i}-x_{j}\right\|^{2}]-\mathsf{R}^{2}\geq
0,\quad i\neq j\in\mathcal{N},$ (7)

the overlap between the coverage of any sensors deployed at $x_{i}$ and
$x_{j}$ is negligible. We will impose constraint (7), by informed design of
the kernel function $k(x,x^{\prime})$ of the SVGD Algorithm (1).

Intuitively, the update in (5) pushes the samples (particles) toward the high
probability regions of the event probability via the gradient term $\nabla\log
p$, while maintaining a degree of diversity via the second term
$\nabla_{x_{j}}k(x_{j},x)$ to prevent the particles from collapsing together
in the local Maximum a posteriori (MAP) solutions [21]. In other words,
$\nabla_{x_{j}}k(x_{j},x)$ acts as a deterministic repulsive force [21, 22].
For our application, we want to choose a kernel that can capture the
underlying geometry of the target distribution while also satisfying the
pairwise distance constraint (7).

For a kernel to be a viable choice in the Stein Variational Gradient Descent
formulation, $k(x,.)$ should be in the Stein class of $p$ satisfying

$\lim_{\left\|x\right\|\rightarrow\infty}p(x)\,k(x,.)=0.$ (8)

Light-tailed distributions satisfy this condition, i.e. smooth functions
$k(x,.)$ with proper boundary conditions and bounded target distribution
$p(x)$ [17]. The kernel function that we choose and that satisfies (8), is the
Mahalanobis-Gaussian kernel function.

$\displaystyle
k(x,x^{\prime})=\exp(-\frac{1}{h}(x_{i}-x_{j})^{\top}\,W^{-1}\,(x_{i}-x_{j})),$
(9)

where we want a positive definite matrix $W$ to scale the pairwise distance
inside a smooth kernel function.

We formulate the weighting matrix $W$ as a variance induced by the pairwise
distance constraint (7). Notice that

$\displaystyle\sum\nolimits_{i,j}\left\|x_{i}-x_{j}\right\|^{2}=2\,n\sum\nolimits_{i=1}^{n}\left\|x_{i}-\mu\right\|^{2}.$

Now, if the samples are iid, we obtain

$\displaystyle\mathbb{E}(\left\|x_{i}-x_{j}\right\|^{2})=2\,\text{Var}(X_{\ell}).$
(10)

where $\mathbb{E}(\left\|x_{i}-x_{j}\right\|^{2})$ can be represented as the
empirical sum of $n^{2}$ pairwise differences and the empirical variance of
the samples is calculated using
$\text{Var}(X_{\ell})=\frac{1}{n}\sum\nolimits_{i=1}^{n}\left\|x_{i}-\mu\right\|^{2}$.
Here, $X_{\ell}$ represents the samples of the empirical distribution
approximating the target distribution using the update (2), and $\mu$ is the
mean of these samples. In our formulation, we want the empirical variance that
is proportional to the pairwise difference to converge to
$\text{Var}(X_{\ell})\to\text{Var}(X)\geq\mathsf{R}^{2}I$ after $\ell$
iterations. $\text{Var}(X)$ is a $d\times d$ dimensional matrix. Since (5), we
are performing the functional gradient descent in RKHS, we need the weighted
matrix, $W$ to be represented in the RKHS. Thus, we map the empirical variance
using the radial basis function to the RKHS, and define
$W=\exp(-\frac{1}{h}\text{Var(X)})$.

Figure 3: The assignment problem, where each sensor $i\in\mathcal{A}$ with
spatial coverage distribution $s_{i}(x)$ finds the best placement position
${x}_{i}$ using the cost matrix $C_{ij}$.

######  Remark 1 (Algorithm initialization)

Since we want to diversify the PoIs around the event distribution such that
there is no overlap, we use our kernel, (9) to generate a repulsive force that
repels the samples to satisfy the no-overlap condition. However, in doing so
we may not have samples in the local modes (MAP), which represent the high
region of probability. Thus, we select a few particles in our initialization
for which the $\nabla_{x}k(x_{i}^{l},x)$ vanishes and reduces to a typical
gradient ascent for MAP. This can be done when we let the parameter
$h\rightarrow 0$ for (9), the term $\nabla_{x}k(x_{i}^{l},x)$ from (5)
vanishes.  $\Box$

_Selecting PoIs to deploy the robots_ : In summary, to select the PoIs to
deploy the robots, we will use the SVGD Algorithm 1 with the kernel function
(9) we designed. Our algorithm’s initialization is guided by Remark 1. The
selected PoIs are denoted by
$\mathcal{X}^{\star}=\\{x_{1}^{\star},\cdots,x_{n}^{\star}\\}.$

###  IV-B Assigning the PoIs to the robots

For isotropic homogeneous sensors, the number of samples/PoIs extracted by the
SVGD process can be set to be the same as the number of sensors ($n=N$). This
means that the samples/PoIs will be the final deployment locations for the
sensors, and sensors can be assigned to these points randomly.

However, when sensors are heterogeneous and/or have an isotropic or
anisotropic quality of service distribution $s(x|x_{i},\theta_{i})$,
$i\in\mathcal{A}$, we need a mechanism to assign sensors to the PoIs and
decide the deployment orientation. In what follows, we assume that we search
over a finite number of deployment orientations for each sensor, that is,
$\theta_{i}\in\Theta=\\{\bar{\theta}_{1},\cdots,\bar{\theta}_{M}\\}\subset[0,2\pi]$.
We cast our assignment problem as a bipartite matching problem; see Fig. 3.
For every sensor $i\in\mathcal{A}$ and every PoI $j\in\\{1,\cdots,M\\}$, let

$\displaystyle
C_{ij}^{\star}=\min_{\theta_{i}\in\Theta}\\{\mathcal{K}\mathcal{L}(s(x|x_{j},\theta_{i})\|p(x))\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ \leavevmode\nobreak\
x\in\mathcal{C}_{i}(x|x_{i},\theta_{i})\\},$

which gives us the minimum dissimilarity among all possible orientations
$\Theta$ between the service distribution of sensor $i$ and $p(x)$ over the
effective coverage footprint of the sensor at the given PoI $j$. In what
follows, we let $\theta_{i}^{j\star}$ be the orientation of sensor $i$
corresponding to $C_{ij}^{\star}$. $C_{ij}^{\star}$ and corresponding
$\theta_{i}^{j\star}$ can also be computed using the discrete sampling method
proposed in our previous paper [16]. With $C_{ij}^{\star}$’s at hand, we cast
the assignment problem as

$\displaystyle\bm{\mathbf{Z}}^{\star}=\arg\min\sum\nolimits_{j\in\mathcal{S}}\\!\sum\nolimits_{i\in\mathcal{A}}\\!\\!\\!Z_{i,j}C_{i,j}^{\star},$
(11a) $\displaystyle Z_{i,j}\in\\{0,1\\},\quad
i\in\mathcal{A},\leavevmode\nobreak\ j\in\mathcal{S},$ (11b)
$\displaystyle\sum\nolimits_{i\in\mathcal{A}}\\!\\!\\!Z_{i,j}=1,\leavevmode\nobreak\
\forall i\in\mathcal{A},$ (11c)
$\displaystyle\sum\nolimits_{j\in\mathcal{S}}\\!\\!\\!Z_{i,j}\leq
1,\leavevmode\nobreak\ \forall j\in\mathcal{S},.$ (11d)

The optimization problem (11) is in the form of a standard assignment problem
and can be solved using existing algorithms, such as the Hungarian algorithm
[23]. The solution of the assignment problem gives the final configuration of
the sensors $(x_{i},\theta_{i}^{\star})$.

##  V Simulation and Experimental Analysis

This section provides a set of experiments and simulations to study the
performance of the proposed framework for homogeneous and heterogeneous
sensors with both isotropic and anisotropic configurations. The framework
provides a final configuration of the sensor that maximizes the spatial
coverage of the environment. To illustrate the collective behavior, we first
consider a $50\times 50$ square environment. The event distribution function
$p:[0,50]^{2}\rightarrow\mathbb{R}_{>0}$ is given by a GMM described by

$\displaystyle
p(x|\psi)\quad=\quad\sum\nolimits^{K}_{k=1}\pi_{k}\mathcal{N}(x|\mu_{K},\Sigma_{k}).$
(12)

where $K\in_{>1}$ is the number of basis of the GMM and
$\psi=(\pi_{k},\mu_{K},\Sigma_{k})$ are the mixing weights, mean, and
covariance, respectively, for $k\in\\{1,\cdots,K\\}$.

Figure 4: Sensor deployment for homogeneous sensors: (left) proposed Stein
Coverage method, (right) Voronoi partitioning. The plots in the bottom row
show the experimental deployment at the Robotarium [24], a remotely accessible
multi-agent research testbed.

###  V-A Homogeneous isotropic sensors

First, we consider a group of $N=10$ isotropic homogeneous sensors with a
probabilistic sensing model of
$s_{i}(x|x_{i},\theta_{i})=\mathcal{N}(x|\mu_{i},\Sigma_{i})$ where
$\mu_{i}=x_{i}$ and $\Sigma_{i}=\sigma^{2}I$, $\sigma=2.5$. Here, the
orientation of the robot, $\theta_{i}$, does not play a role as the coverage
is isotropic. We let the radius of the effective sensing footprint of each
sensor be $2\sigma=5$, i.e.,
$\mathcal{C}_{i}(x|x_{i},\theta_{i})=\\{x\in\mathcal{W}\,|\,\|x-x_{i}\|^{2}\leq
5\\}$, which results in $\mathsf{R}=5$. We use our proposed Stein Coverage
algorithm to determine the location of the sensors to cover the event density
as shown on the left side of Fig. 4. To compare our proposed method with
conventional Voronoi partitioning, we obtain the desired Voronoi partitions
$\mathcal{V}(.)=\\{\mathcal{V}_{1},\cdots,\mathcal{V}_{n}\\}$, from minimizing

$\displaystyle\sum\nolimits^{n}_{i=1}\int_{V_{i}(x)}\left\|x-x_{i}\right\|^{2}p(x)dx.$

where $p(x)$ is given in (12), and is used to bias the importance of the
location optimization function. The result is shown in the plots on the right-
hand side of Fig. 4. Because the sensors are homogeneous, there is no need to
solve the optimal assignment (11). Comparing the two configurations shows us
that Stein Coverage is able to capture the high-density region of the density
function $\phi$, compared to Voronoi partitioning. This is because SVGD uses a
kernelized (functional) gradient descent method, to move the particles where
the KL divergence is minimized in terms of Stein discrepancy, the update is
better at capturing the mean and variance of Gaussian distributions. However,
Voronoi partitioning based on Euclidean distance might not capture the
underlying structure of the distribution well.

(a) Sensor deployment for heterogeneous sensors: (left) Stein Coverage,
(right) Voronoi partitioning using power diagrams.

(b) Comparison of sensor deployment:(left) proposed Stein Coverage, (right)
Voronoi partitioning using power diagrams, when the number of sensors is less
than the number of clusters.

Figure 5: Sensor deployment for heterogeneous sensor.

###  V-B Heterogeneous isotropic and anisotropic sensors

Next, we present a simulation for a group of $N=10$ heterogeneous sensors.
Here, we consider a scenario where the sensors are isotropic with uniform
coverage everywhere, so we can compare our results with the power diagrams
method, a generalized Voronoi partitioning deployment for heterogeneous
sensors with uniform quality of service distribution.

(a) Sensor deployment for anisotropic heterogeneous sensors. The level sets
represent the multimodal distribution of the target points (blue dots) and the
ellipsoid encircling the level sets represents the configuration of
anisotropic footprint calculated using the proposed algorithm.

(b) Sensor deployment for heterogeneous anisotropic sensors when the number of
sensors is less than the number of clusters.

Figure 6: Heterogeneous and anisotropic sensor deployment

We assume that every sensor’s effective footprint is
$\mathcal{C}_{i}(x|x_{i},\theta_{i})=\\{x\in\mathcal{W}|\|x-x_{i}\|\leq
r_{i}\\}$ for $i\in\mathcal{A}$ with corresponding $r_{i}$ belonging to the
randomly generated ordered set
$\\{2.2,2.3,2.5,2.8,2.8,3.0,3.3,3.4,3.6,3.8\\}\subset\mathbb{R}_{>0}$. We
compare this method with coverage for heterogeneous sensors, using power
diagrams, $\mathcal{P}(.)=\\{\mathcal{P}_{1},\cdots,\mathcal{P}_{n}\\}$
generalized Voronoi partitioning with additive weights$(r_{i}\in$
$\mathbf{r})$ obtained from minimizing

$\displaystyle\sum\nolimits^{n}_{i=1}\int_{\mathcal{P}_{i}}(\left\|x-x_{i}\right\|^{2}-r_{i}^{2})p(x)dx.$

The result of implementing the Stein Coverage and deployment using the power
diagrams is shown in Fig. 5(a). As seen, the Stein Coverage achieves a more
effective placement, with sensors deployed to the important regions and in
accordance to their coverage footprint. The difference between the two methods
is even more evident when the number of sensors is significantly less than
important coverage subject clusters; an example scenario is shown in Fig.
5(b). As seen in Fig. 5(b), a Voronoi-based solution is far less effective in
this case, because of its implicit consideration of coverage beyond the
effective footprint of the sensors (enforcing the footprint as a soft
constraint).

Next, we present a simulation for $N=10$ heterogeneous sensors, where the
footprints are anisotropic with probabilistic sensing models of
$s_{i}(x|x_{i},\theta_{i})=\mathcal{N}(x|\mu_{i},\Sigma_{i})$ where
$\mu_{i}=x_{i}$ and
$\bm{\mathbf{\Sigma_{i}}}=\left[\begin{smallmatrix}\sigma^{2}_{i,11}&\sigma^{2}_{i,12}\\\
\sigma^{2}_{i,21}&\sigma^{2}_{i,22}\end{smallmatrix}\right]$ ,where
$\sigma_{i,11}\in\\{2,2,2,1,1,1,2,1,2,2\\}$,
$\sigma_{i,12}\in\\{1,0.5,1,2,0,0,0,0,1,0\\}$,
$\sigma_{i,21}\in\\{1,1,0.1,0,1,0,0,0,1,0\\}$ and,
$\sigma_{i,22}\in\\{2,2,1,1,2,1,1,2,2,2\\}$. Here, the orientation of the
robot, $\theta_{i}$, is important and should be taken into account when
designing the final sensor deployment configuration. We let the radius of the
effective sensing footprint of each sensor $i\in\\{1,\cdots,N\\}$ be
$2\max\\{\sigma_{i,1},\sigma_{i,1}\\}$, and use $\mathsf{R}=4$.

The density distribution $p(x)$ is a GMM with $K=10$ basis, as shown in Fig.
6(a). In this case, because the sensors are anisotropic, we want to find the
configuration $(x_{i},\theta_{i})$, $i\in\\{1,\cdots,N\\}$, to maximize the
coverage. In implementing Stein Coverage, we use a sampling-based method
proposed in our previous work [16] to compute $C_{ij}^{\star}$ and
$\theta_{i}^{j\star}$. Fig. 6(a) shows the deployment configuration
$(x_{i},\theta_{i})$ when the number of sensors is equal to the number of
basis of the GMM distribution $p(x)$. On the other hand, Fig. 6(b),
illustrates the deployment configuration $(x_{i},\theta_{i})$, when we have a
limited number of sensors ($N=4$) compared to the number of the basis
($K=10$). As seen, Stein Coverage deploys the sensors to most important areas
and determines their orientation to provide the best local distribution-
matching between the sensors and the coverage subjects distribution.

##  VI Conclusion

In this paper, we consider the spatial coverage optimization problem of
multiple sensors and propose a novel framework inspired by Stein Variational
Gradient Descent to configure a group of sensors to optimally cover a given
event distribution in convex environments. We use a new statistical distance
measure to determine potential sensor assignment locations. We characterized a
locally optimal coverage configuration with soft workspace statistics and
compared it with hard workspace statistics like Voronoi partitioning. We
compared these two different approaches using numerical simulations for both
heterogeneous and homogeneous sensors. In numerical simulations, we were able
to show how the coverage quality of the proposed framework is better compared
to the Voronoi partitioning approaches. Finally, we extend our work to deploy
the sensor with an anisotropic footprint using the method we proposed in our
previous paper. Future work will explore the practical extension of our
framework to non-convex environments and 3D coverage for environmental
monitoring, surveillance, and search and rescue.

## References

  * [1] O. Arslan, “Statistical coverage control of mobile sensor networks,” IEEE Transactions on Robotics, vol. 35, no. 4, pp. 889–908, 2019.
  * [2] G. Notomista, S. Mayya, Y. Emam, C. Kroninger, A. Bohannon, S. Hutchinson, and M. Egerstedt, “A resilient and energy-aware task allocation framework for heterogeneous multirobot systems,” IEEE Transactions on Robotics, vol. 38, no. 1, pp. 159–179, 2022.
  * [3] A. Farinha, R. Zufferey, P. Zheng, S. F. Armanini, and M. Kovac, “Unmanned aerial sensor placement for cluttered environments,” IEEE Robotics and Automation Letters, vol. 5, no. 4, pp. 6623–6630, 2020.
  * [4] J. Cortes, 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.
  * [5] J. Aponte-Luis, J. A. Gómez-Galán, F. Gómez-Bravo, M. Sánchez-Raya, J. Alcina-Espigado, and P. M. Teixido-Rovira, “An efficient wireless sensor network for industrial monitoring and control,” Sensors, vol. 18, no. 1, 2018.
  * [6] N. Megiddo and K. J. Supowit, “On the complexity of some common geometric location problems,” SIAM Journal on Computing, vol. 13, no. 1, pp. 182–196, 1984.
  * [7] M. Schwager, D. Rus, and J.-J. Slotine, “Decentralized, adaptive coverage control for networked robots,” The International Journal of Robotics Research, vol. 28, no. 3, pp. 357–375, 2009.
  * [8] O. Arslan and D. E. Koditschek, “Voronoi-based coverage control of heterogeneous disk-shaped robots,” in 2016 IEEE International Conference on Robotics and Automation (ICRA), pp. 4259–4266, 2016.
  * [9] A. Kwok and S. Martinez, “Deployment algorithms for a power-constrained mobile sensor network,” in 2008 IEEE International Conference on Robotics and Automation, pp. 140–145, 2008.
  * [10] M. Santos and M. Egerstedt, “Coverage control for multi-robot teams with heterogeneous sensing capabilities using limited communications,” in 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5313–5319, 2018.
  * [11] A. Krause, A. Singh, and C. Guestrin, “Near-optimal sensor placements in Gaussian processes: Theory, efficient algorithms and empirical studies.,” Journal of Machine Learning Research, vol. 9, no. 2, 2008.
  * [12] A. Likas, N. Vlassis, and J. J. Verbeek, “The global k-means clustering algorithm,” Pattern Recognition, vol. 36, no. 2, pp. 451–461, 2003.  Biometrics.
  * [13] D. Inoue, Y. Ito, and H. Yoshida, “Optimal transport-based coverage control for swarm robot systems: Generalization of the voronoi tessellation-based method,” IEEE Control Systems Letters, vol. 5, no. 4, pp. 1483–1488, 2021.
  * [14] V. Krishnan and S. Martínez, “A multiscale analysis of multi-agent coverage control algorithms,” Automatica, vol. 145, p. 110516, 2022.
  * [15] Y.-f. Chung and S. S. Kia, “A distributed service-matching coverage via heterogeneous agents,” IEEE Robotics and Automation Letters, vol. 7, no. 2, pp. 4400–4407, 2022.
  * [16] D. Ghimire and S. S. Kia, “Optimal multi-sensor deployment via sample-based quality-of-service distribution matching,” in 2023 European Control Conference (ECC), pp. 1–6, 2023.
  * [17] Q. Liu, J. Lee, and M. Jordan, “A kernelized stein discrepancy for goodness-of-fit tests,” in International Conference on Machine Learning, pp. 276–284, PMLR, 2016.
  * [18] S.-i. Amari and A. Cichocki, “Information geometry of divergence functions,” Bulletin of the polish academy of sciences. Technical sciences, vol. 58, no. 1, pp. 183–195, 2010.
  * [19] Q. Liu, “Stein variational gradient descent as gradient flow,” Advances in neural information processing systems, vol. 30, 2017.
  * [20] C. Stein, “Approximate computation of expectations,” in Lecture Notes-Monograph Series, IMS, 1986.
  * [21] Q. Liu and D. Wang, “Stein variational gradient descent: A general purpose bayesian inference algorithm,” Advances in neural information processing systems, vol. 29, 2016.
  * [22] G. Detommaso, T. Cui, Y. Marzouk, A. Spantini, and R. Scheichl, “A stein variational newton method,” in Advances in Neural Information Processing Systems (S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, eds.), vol. 31, Curran Associates, Inc., 2018.
  * [23] H. W. Kuhn, “The hungarian method for the assignment problem,” Naval Research Logistics Quarterly, vol. 2, no. 1-2, pp. 83–97, 1955.
  * [24] S. Wilson, P. Glotfelter, L. Wang, S. Mayya, G. Notomista, M. Mote, and M. Egerstedt, “The robotarium: Globally impactful opportunities, challenges, and lessons learned in remote-access, distributed control of multirobot systems,” IEEE Control Systems Magazine, vol. 40, no. 1, pp. 26–44, 2020.