Scalar Field Estimation with Mobile Sensor Networks

TL;DR

This paper presents a novel method for estimating scalar fields using mobile sensors, leveraging adaptive control and Lyapunov techniques to ensure stability and convergence of the estimation process.

Contribution

It introduces a new estimation algorithm for scalar fields with mobile sensors, incorporating stability analysis and motion planning for guaranteed convergence.

Findings

Proposed an adaptive estimation algorithm with stability guarantees.

Proved convergence of parameter estimates to true values.

Demonstrated effectiveness through theoretical analysis.

Abstract

In this paper, we consider the problem of estimating a scalar field using a network of mobile sensors which can measure the value of the field at their instantaneous location. The scalar field to be estimated is assumed to be represented by positive definite radial basis kernels and we use techniques from adaptive control and Lyapunov analysis to prove the stability of the proposed estimation algorithm. The convergence of the estimated parameter values to the true values is guaranteed by planning the motion of the mobile sensors to satisfy persistence-like conditions.

Tables

Table 1. TABLE I : Parameters of the simulated scalar field

$c_{i,x}$	$0.20$	$0.35$	$0.60$	$0.85$	$0.70$	$0.75$	$0.15$	$0.35$
$c_{i,y}$	$0.25$	$0.26$	$0.18$	$0.30$	$0.75$	$0.90$	$0.75$	$0.60$
$a^i$	$2.0$	$1.0$	$1.5$	$1.8$	$1.2$	$1.6$	$2.5$	$1.1$

Table 2. TABLE II : Comparison of algorithms for p = 100 𝑝 100 p=100 parameters.

${\sigma }_i=0.04$	$T$ (sec)	${\Vert e\Vert }_2$	${\sigma }_i=0.05$	$T$ (sec)	${\Vert e\Vert }_2$
Algorithm S$1$	$3.1$	$0.045$	Algorithm S$1$	$3.9$	$0.012$
Algorithm S$2$	$3.1$	$0.054$	Algorithm S$2$	$3.7$	$0.053$
Algorithm S$3$	$3.1$	$0.048$	Algorithm S$3$	$3.7$	$0.028$

Table 3. TABLE III : Comparison of algorithms for p = 196 𝑝 196 p=196 parameters.

${\sigma }_i=0.03$	$T$ (sec)	${\Vert e\Vert }_2$	${\sigma }_i=0.04$	$T$ (sec)	${\Vert e\Vert }_2$
Algorithm S$1$	$6.6$	$0.031$	Algorithm S$1$	$8.9$	$0.008$
Algorithm S$2$	$6.6$	$0.059$	Algorithm S$2$	$8.8$	0.073
Algorithm S$3$	$6.6$	$0.053$	Algorithm S$3$	$8.8$	$0.039$

Equations

$$\varphi (q)$$

$$=\sum _{i=1}^p{\mathcal{K}}^i(q)a^i$$

$$\hat{f}(q)$$

$${\mathcal{K}}^i(q)=\phi (\Vert c_i-q\Vert )=\exp \left\{-\frac{\Vert c_i-q{\Vert }^2}{{\sigma }_i^2}\right\}$$

$$S:={\int }_{\mathcal{Q}}\mathcal{K}(q)\mathcal{K}(q{)}^{\top }dq$$

$${\int }_{\mathcal{Q}}|\mathcal{K}(q{)}^{\top }v{|}^{2}dq\ge 0$$

$$\left(\begin{array}[]{cccc}\mathcal{K}^{1}(c_{1})&\mathcal{K}^{1}(c_{2})&\dots&\mathcal{K}^{1}(c_{p})\\ \mathcal{K}^{2}(c_{1})&\mathcal{K}^{2}(c_{2})&\dots&\mathcal{K}^{2}(c_{p})\\ \dots&\ddots&\dots&\vdots\\ \mathcal{K}^{p}(c_{1})&\mathcal{K}^{p}(c_{2})&\dots&\mathcal{K}^{p}(c_{p})\\ \end{array}\right)$$

$${\int }_Q|\mathcal{K}(q{)}^{\top }v{|}^{2}dq>0\text{mbox}foranyv\ne 0$$

$${\mathcal{S}}_T:={\int }_0^T\mathcal{K}(x(t))\mathcal{K}(x(t){)}^{\top }dt$$

$$\dot{\Lambda }$$

$$\dot{\lambda }$$

$$\dot{\hat{a}}=-\Gamma \left(\Lambda \hat{a}-\lambda \right),$$

$$S(T):={\int }_0^T\mathcal{K}(\tau )\mathcal{K}(\tau {)}^{\top }d\tau$$

$$S(t)={\int }_0^t\mathcal{K}(\tau )\mathcal{K}(\tau {)}^{\top }d\tau$$

$$V=\frac{1}{2}{\tilde{a}}^{\top }{\Gamma }^{-1}\tilde{a}$$

$$\dot{V}={\tilde{a}}^{\top }{\Gamma }^{-1}\dot{\hat{a}}$$

$$\dot{V}$$

$${\mathcal{A}}_j^{\epsilon }:=\left\{q\in \mathcal{Q}:|{\mathcal{X}}^j(q)|-\sum _{i=1,i\ne j}^p|{\mathcal{X}}^i(q)|>\epsilon \right\}.$$

$${\bar{\mathcal{S}}}_T={\int }_0^T\mathcal{X}(x(t))\mathcal{X}(x(t){)}^{\top }dt.$$

$${\int }_0^T{u}^{\top }\mathcal{X}(x(t))\mathcal{X}(x(t){)}^{\top }udt<\delta$$

$${\int }_0^T{u}^{\top }$$

$$\ge {\int }_{t_{i1}}^{t_{i2}}|{\mathcal{X}}^{\top }u{|}^{2}dt={\int }_{t_{i1}}^{t_{i2}}|\sum _{j=1}^p{\mathcal{X}}^j{u}^j{|}^{2}dt$$

$$\ge {\int }_{t_{i1}}^{t_{i2}}(|{\mathcal{X}}^i{u}^i|-|\sum _{j=1,j\ne i}^p{\mathcal{X}}^j{u}^j|{)}^{2}dt$$

$$\ge {\int }_{t_{i1}}^{t_{i2}}(|{\mathcal{X}}^i{u}^i|-\sum _{j=1,j\ne i}^p|{\mathcal{X}}^j{u}^j|{)}^{2}dt$$

$$\ge {\int }_{t_{i1}}^{t_{i2}}((|{\mathcal{X}}^i|-\sum _{j=1,j\ne i}^p|{\mathcal{X}}^j|)|u^i|{)}^{2}dt$$

$$\ge {\int }_{t_{i1}}^{t_{i2}}{\epsilon }^2|u^i{|}^{2}dt=(t_{i2}-t_{i1}){\epsilon }^2|u^i{|}^2.$$

$$\Vert \mathcal{K}(x)-\mathcal{K}(c_j){\Vert }_{\infty }<\frac{(1-ϵ)}{2(p-1)\Vert {\mathrm{K}}^{-1}{\Vert }_{\infty }},$$

$${\mathcal{X}}^i(x)-{\mathcal{X}}^i(c_j)={\left[{\mathrm{K}}^{-1}(\mathcal{K}(x)-\mathcal{K}(c_j))\right]}^i$$

$$\left[\begin{array}[]{c}y_{1}\\ y_{2}\end{array}\right]=B_{j}\left(\mathcal{X}(x)-\mathcal{X}(c_{j})\right)$$

$$B_{j}=\left[\begin{array}[]{ccccccc}0&\dots&0&1&0&\dots&0\\ 1&\dots&1&0&1&\dots&1\end{array}\right]$$

$$\Vert y{\Vert }_{\infty }$$

$$\le \Vert B{\Vert }_{\infty }\Vert {\mathrm{K}}^{-1}{\Vert }_{\infty }\Vert \mathcal{K}(x)-\mathcal{K}(c_j){\Vert }_{\infty }$$

$$\Vert \mathcal{K}(x)-\mathcal{K}(c_j){\Vert }_{\infty }<\frac{(1-ϵ)}{2(p-1)\Vert {\mathrm{K}}^{-1}{\Vert }_{\infty }}$$

$$\mathcal{K}(q)=\left[\begin{array}[]{c}\mathcal{K}^{(1)}(q)\\ \mathcal{K}^{(2)}(q)\\ \vdots\\ \mathcal{K}^{(N)}(q)\end{array}\right],\qquad a=\left[\begin{array}[]{c}a^{(1)}\\ a^{(2)}\\ \vdots\\ a^{(N)}\end{array}\right]$$

Full text

I Introduction

Multi-robot systems consists of network of robots which cooperate to perform tasks such as consensus, formation control etc. [1, 2, 3, 4]. Mobile sensor networks consists of network of robots mounted with sensors deployed to perform some distributed sensing task such as monitoring, coverage etc [5]. In this paper, we consider the problem of estimation of an unknown scalar field using mobile sensor networks. There have been many works related to scalar field estimation in literature. Several works have studied field estimation using wireless sensor networks. See for example [6, 7]. In [8] the scalar field is assumed to be modelled using a partial differential equation and finite element methods are used for estimating the field. In [9, 10, 11, 12, 13, 14] the field is modelled as spatial random process and estimated using samples from the sensor nodes. In [15] field reconstruction is posed as an optimization problem constrained by linear dynamics and a gradient-based method is used to solve the problem. In [16], the scalar field is assumed to be linearly parameterized in terms of Gaussian basis functions and the measurements from the sensors are fused together to form an estimate for the scalar field.

In most of these cases, the sensors are assumed to be fixed and distributed over the region of interest. Usually a large number of sensors are required to be installed for achieving enough spatial resolution. Using mobile sensor networks can be highly advantageous since they can move around the region of interest and collect measurements adaptively, the number of sensors required is greatly reduced. In [17, 18], scalar field estimation is done with mobile sensor network by fusing sensor measurements using consensus filters. In [19], information about a scalar field is obtained by exploring the level surfaces of the field using a mobile sensor network. In [20], a static sensor network is used along with a mobile robot to estimate a scalar field by combining the robot measurements with the sensor network measurements and planning the robot trajectory to minimize some reconstruction error. However the method we propose in the current work is motivated by the coverage control problem [5, 21, 22, 23, 24].

In the coverage problem, we are interested in controlling the robots so that the robots attain an optimal configuration or a near optimal configuration with respect to a scalar field. In [5], this is achieved by minimizing a cost function which gives a measure of how good the coverage is. In [21], the authors extended the coverage algorithm for the case where the scalar field is unknown. The scalar field is assumed to be linearly parameterized with unknown constant parameters. In order to achieve the coverage goal, the robot needs to adapt the unknown parameters so that the estimated scalar field is close to the actual field. The exact estimation (asymptotically) of the density function parameters require a time integral quantity to be positive definite, which is a sufficient richness condition for the robot trajectories. See [21] for more details. In general, the robot trajectories need not meet this condition since the trajectories of the robots are decided based on the gradient of the coverage cost function, not on estimating the density function parameters. However, it is crucial to estimate the true values of those parameters since the estimation of the unknown scalar field is often the primary objective for a robotic sensor network and it may lead more efficient deployment of robots. For example, in case of radiation spill, if we have a good estimate of the radiation concentration, we may directly deploy agents to regions of high concentration.

Thus in this work, we look at a slightly different problem closely related to and motivated by the coverage problem discussed above. Our primary aim in this paper is to accurately estimate the scalar field not the coverage. The unknown scalar field is approximated using positive definite radial basis functions and we use a similar adaptive approach as that in [21] for parameter estimation.

In Section II, we discuss the problem statement in detail. In section III, we consider the single mobile sensor case, followed by the mobile sensor network case in Section IV. In Section V we discuss the case where the centres of the radial basis functions are not known exactly, but only to within an $\epsilon$-accuracy. We present some simulations to verify the results in section VI. We conclude the paper with Section VII.

II Preliminaries and Problem Statement

We denote the set of positive real numbers by $\mathbb{R}_{+}$. The components of a vector $v$ are denoted using superscripts $v^{i}$. Subscripts on vector quantities refer to the agent or mobile sensor the quantity is associated to. For example, $v_{i}$ refers to a quantity associated with agent $i$.

We consider a compact region $\mathcal{Q}\subset\mathbb{R}^{n}$ with $N$ mobile sensors. The position of the sensors is denoted by $x_{i};\,\,\,i=1,2,\dots,N$. There also exists a continuous scalar field $\phi:\mathcal{Q}\to\mathbb{R}_{+}$ over $\mathcal{Q}$ which is unknown. The objective is to estimate the unknown scalar field using $N$ mobile sensors assuming the sensors can measure the value of the scalar field at their respective positions. We assume that the unknown scalar field can be represented by positive definite radial basis functions (RBF). In other words, we assume the density function can be parameterized as

[TABLE]

where $a\in\mathbb{R}^{p}$ is a constant vector, and $\mathcal{K}(q)=\left[\mathcal{K}^{1}(q)\,\,\,\mathcal{K}^{2}(q)\,\,\,\dots\,\,\,\mathcal{K}^{p}(q)\right]^{\scriptscriptstyle\top}$ with $\mathcal{K}^{i}:\mathcal{Q}\to\mathbb{R}_{+}$ given by $\mathcal{K}^{i}(q)=\varphi(\|c_{i}-q\|)$ with $\varphi:\mathbb{R}_{+}\to\mathbb{R}_{+}$ are radial basis functions for a set of points $c_{i}$. This assumption is common in neural networks and justified as follows:

Theorem 1** ([25, 26]).**

For any continuous function $f:\mathbb{R}^{n}\to\mathbb{R}$ and any $\epsilon>0$, there is an RBF network with $p$ elements, a set of centers $\{c_{i}\}_{i=1}^{p}$, such that we can define

[TABLE]

with $\|f-\hat{f}\|_{L_{2}}^{2}\leq\epsilon=\mathcal{O}\left(p^{-\frac{1}{n}}\right)$.

The theorem tells us that we can approximate a continuous function to an arbitrary accuracy by using a network of RBF elements. An example of positive definite radial kernel is the Gaussian kernel,

[TABLE]

where $c_{i}$ are the centres of the Gaussian kernels. The main problem studied in this work is to accurately determine the parameters $a^{i}$ so that the scalar field $\phi(.)$ may be accurately reconstructed. We make the following assumption:

Assumption 1**.**

The centres $c_{i}$ of the radial functions are known to all the mobile agents.

The strengths $a_{i}$ of individual radial functions are unknown and need to be estimated. To proceed, we require the following theorem:

Theorem 2** (Micchelli’s Theorem [27]).**

Given $p$ distinct points $c_{1},c_{2},\dots,c_{p}$ in $\mathbb{R}^{q}$, the $p\times p$ matrix $\mathrm{K}$, whose elements are $\mathrm{K}_{ij}=\mathcal{K}^{i}(c_{j})=\varphi(\|c_{i}-c_{j}\|)$ is non-singular.

The theorem says that for positive definite radial kernels, the $p\times p$ matrix formed by evaluating the radial functions at each of the centres is non-singular. In what follows, we assume that $\phi(.)$ can be exactly parameterized by the RBF kernels. A consequence of theorem 2 is given below:

Lemma 1**.**

The matrix $S$ given by

[TABLE]

where $\mathcal{K}(q)=\left[\mathcal{K}^{1}(q)\,\,\,\mathcal{K}^{2}(q)\,\,\,\dots\,\,\,\mathcal{K}^{p}(q)\right]^{\scriptscriptstyle\top}$ and $\phi$ is parameterized as in (1a), is positive definite.

Proof.

From the definition of $S$, we know it is atleast positive semi-definite. Therefore for any $v\neq 0$, $v^{\scriptscriptstyle\top}Sv\geq 0$ or

[TABLE]

Now, since $\mathcal{K}(q)$ consists of positive definite radial kernels, we have from theorem 2 that

[TABLE]

is positive definite. This implies that the vectors $\mathcal{K}(c_{j});\,\,j=1,2,\dots,p$ are linearly independent. Thus, given any $v\neq 0,v\in\mathbb{R}^{p}$, there exists some $j\in\{1,2,\dots,p\}$ such that $\mathcal{K}(c_{j})^{\scriptscriptstyle\top}v$ is non-zero. This along with the fact that $\mathcal{K}(\cdot)$ is continuous allows us to conclude that

[TABLE]

Hence, $S$ is positive definite. ∎

III Single Mobile Robot Sensor

In this section, we consider the case of a single mobile sensor ($N=1$) with position $x(t)$ at time $t$ deployed in the region $\mathcal{Q}$ to estimate the scalar field parameter $a$ (as given by equation (1a)). The estimate of $a$ is denoted by $\hat{a}$. Then we can state the following corollary to lemma 1.

Corollary 1**.**

Suppose the mobile sensor moves continuously within the domain $Q$, such that in time $T$, it passes through each of the RBF centres $c_{i}\,;\,\,i=1,2,\dots,p$, then

[TABLE]

is positive definite.

Proof.

The proof is essentially the same and follows from lemma 1. ∎

Now consider the following integrators running on the mobile sensor:

[TABLE]

where $\mathcal{K}(t):=\mathcal{K}(x(t))$ denotes the value of function $\mathcal{K}(\cdot)$ at the point where the robot is at time $t$ and $\phi(t)$ is the measured value of the density function $\phi(\cdot)$ by the robot at time $t$.

Proposition 1**.**

Suppose the mobile sensor moves such that it passes through each of the centres $c_{i};\,\,i=1,2,\dots,p$ in some finite time $T>0$, and during this motion updates its estimate $\hat{a}$ of $a$ by

[TABLE]

where $\Gamma$ is a positive definite gain matrix, then the estimate $\hat{a}$ is bounded and converges asymtotically to the true value $a$.

Proof.

Under the assumptions of the proposition 1 and corollary 1,

[TABLE]

is positive definite. This implies that

[TABLE]

is positive definite for all $t\geq T$.

Now consider the positive definite candidate Lyapunov function,

[TABLE]

where $\tilde{a}=\hat{a}-a$ is the estimation error. Taking the derivative of $V$, we obtain

[TABLE]

Substituting the update law from (6) and simplifying, we get

[TABLE]

where $\alpha=\frac{\lambda_{\min}(S(T))}{\lambda_{\max}(\Gamma^{-1})}>0$, $\lambda_{min}(\cdot)$ and $\lambda_{\max}(\cdot)$ denoting the minimum and maximum eigenvalues of their argument. Since $V$ is always non-increasing and bounded from below, $\tilde{a}(t)$ is bounded for all $t>0$. Since $\dot{V}<0$ for all $t\geq T$, then we have $V(t)\to 0$ as $t\to\infty$. This implies that $\tilde{a}\to 0$ as $t\to\infty$. ∎

Remark 1**.**

The matrix $S(t)$ being positive definite for all $t\geq T$ is a sufficient excitation condition, similar to (but weaker than) the persistency of excitation condition, on the robot trajectories which ensures parameter convergence. See [21] for more information.

III-A Relaxing the condition in corollary 1

In corollary 1, it was required that the mobile sensor passes through the centres $c_{i}$ of the radial kernels. This can be relaxed so that the mobile sensor need only move through a sufficiently small neighbourhood of each of the centres $c_{i}$, as described in [28]. Consider the vector $\mathcal{X}(q):=\mathrm{K}^{-1}\mathcal{K}(q)$ where $\mathrm{K}$ is the matrix specified in theorem 2. Then $\mathcal{X}(q)$ has the property that $\mathcal{X}^{j}(c_{k})=\delta_{jk}$ where $\delta_{jk}$ is the Kronecker delta function and $\mathcal{X}^{j}(c_{k})$ is the $j$-th component of $\mathcal{X}(c_{k})$. Now consider the diagonal dominance sets defined by ($0<\varepsilon<1$)

[TABLE]

It can be easily seen that $\mathcal{A}_{j}^{\varepsilon}$ contains the centre $c_{j}$ and thus $\mathcal{A}_{j}^{\varepsilon}$ is an open subset containing $c_{j}$. The following lemma is an adaptation of theorem $1$ in [28]:

Lemma 2**.**

Suppose that the mobile sensor moves continuously throughout the domain $\mathcal{Q}$ such that in time $T$, the trajectory traverses through each of the neighbourhoods $\mathcal{A}_{j}^{\varepsilon},\,\,j=1,2,\dots,p$, then the matrix $\mathcal{S}_{T}$ given by equation (4) is positive definite.

Proof.

$\mathcal{S}_{T}$ can be written as $\mathcal{S}_{T}=\mathrm{K}\bar{\mathcal{S}}_{T}\mathrm{K}^{\scriptscriptstyle\top}$ where

[TABLE]

Since $\mathrm{K}$ is invertible, $\mathcal{S}_{T}$ is positive definite iff $\bar{\mathcal{S}}_{T}$ is positive definite. $\bar{\mathcal{S}}_{T}$ is positive definite iff there exists some $\delta>0$ such that \underaccent{\bar}{\sigma}(\bar{\mathcal{S}}_{T})\geq\delta where \underaccent{\bar}{\sigma}(A) denotes the minimum singular value of $A$. Suppose $\bar{\mathcal{S}}_{T}$ is not positive definite under the conditions of the theorem. Then there exists no $\delta>0$ such that \underaccent{\bar}{\sigma}(\bar{\mathcal{S}}_{T})\geq\delta. This implies that for any $\delta>0$, there exists $u\neq 0,\|u\|=1$ such that $u^{\scriptscriptstyle\top}\bar{\mathcal{S}}_{T}u<\delta$, i.e.,

[TABLE]

Let $i$ be the index of the components of $u$ which has the largest absolute value. i.e., $|u^{i}|\geq|u^{j}|\,\,\forall j$. Also let $[t_{i1},t_{i2}]\subset[0,T]$ be the subinterval during which the mobile sensor trajectory is contained in the set $\mathcal{A}_{i}^{\varepsilon}$. Clearly since the set $\mathcal{A}_{i}^{\varepsilon}$ is open and the trajectory is continuous, $[t_{i1},t_{i2}]$ has finite positive length. Then,

[TABLE]

Choosing $\delta<(t_{i2}-t_{i1})\varepsilon^{2}|u^{i}|^{2}$ leads to a contradiction. Therefore, $\bar{\mathcal{S}}_{T}$ is positive definite and hence $\mathcal{S}_{T}$ is positive definite. ∎

A sufficient condition for satisfaction of lemma

2’s assumptions:

Since checking the condition of the mobile sensor traversing througn the sets $\mathcal{A}_{j}^{\varepsilon}$ in lemma 2 involves transforming the vector $\mathcal{K}(q)$ at each instant which can be cumbersome if the number of parameters are large, we present a simpler sufficient condition which ensures that a given point $q$ is inside the set $\mathcal{A}_{j}^{\varepsilon}$. Note that the conditions derived are not equivalent to the conditions of the lemma, but only sufficient and thus can be conservative. However it is beneficial during implementations.

Lemma 3**.**

Given the mobile sensor position $x$, if

[TABLE]

then $x\in\mathcal{A}_{j}^{\varepsilon}$.

Proof.

We have the $i$-th component of $\mathcal{X}(x)$, $\mathcal{X}^{i}(x)=\left[\mathrm{K}^{-1}\mathcal{K}(x)\right]^{i}$. Then

[TABLE]

Now consider the mapping

[TABLE]

where

[TABLE]

The $1$ in the first row and the [math] in the second row occurs at the $j$-th column. If the infinity-norm of $y=[y_{1},y_{2}]^{\scriptscriptstyle\top}$, $\|y\|_{\infty}<(1-\varepsilon)/2$, then it is guaranteed that $x\in\mathcal{A}_{j}^{\varepsilon}$. We also have

[TABLE]

Requiring the above bound to be less than $\frac{(1-\epsilon)}{2}$ and noting that $\|B\|_{\infty}=(p-1)$ we have

[TABLE]

∎

Any point $p$ which satisfies the above condition will lie in the set $\mathcal{A}_{j}^{\varepsilon}$ although all points in $\mathcal{A}_{j}^{\varepsilon}$ are not characterized by the above condition.

IV Mobile Sensor Network

Suppose that we have $N$ mobile sensors deployed in the region $\mathcal{Q}$, with the position of the $i$-th mobile sensor denoted by $x_{i}$. We want to estimate the function $\phi:\mathcal{Q}\to\mathbb{R}_{+}$ collectively. We assume that equation (1a) holds so that we can linearly parameterize $\phi(\cdot)$ in terms of radial basis functions. We partition the region into $N$ components $\mathcal{Q}_{i}\,\,(i=1,2,\dots,N)$. Correspondingly we partition the basis function vector $\mathcal{K}(q)$ and the parameter vector $a$ as

[TABLE]

Each region $\mathcal{Q}_{i}$ contains the centres of the basis functions in the sub-vector $\mathcal{K}^{(i)}$. We assign each region $\mathcal{Q}_{i}$ to one of the mobile sensors where the sensor operates. This assignment is permanent and each mobile sensor starts within its region $\mathcal{Q}_{i}$ and moves in $\mathcal{Q}_{i}$. The algorithms presented below do not depend on any particular partition or assignment of mobile sensors, and this can be done arbitrarily. One particular method to divide the region and assign the sensors will be discussed in section VI. Assuming the region $\mathcal{Q}$ is partitioned and the mobile sensors are assigned to each partition, we consider the graph $\mathcal{G}$ with the vertices representing the mobile sensors and an edge existing between two sensors if they belong to adjacent partitions. By adjacent partitions, we mean two partitions which share a subset of their boundary with each other that is of non-zero length. See figure 1 for an illustration.

Now we consider two cases: (1) each mobile sensor estimates the entire parameter vector, and (2) each mobile sensor estimates only part of the parameter vector.

IV-A Each mobile sensor estimates the full parameter vector

In this subsection, we consider the case where each mobile sensor estimates the entire parameter vector, the estimate of sensor $i$ being denoted by $\hat{a}_{i}$. To proceed, we consider the following integrators running on mobile sensor $i$:

[TABLE]

where $\mathcal{K}_{i}(t)=\mathcal{K}(x_{i}(t))$ and $\phi_{i}(t)=\phi(x_{i}(t))$ is the measurement of $\phi(.)$ obtained by sensor $i$ at its location at time $t$.

We consider the following update law for the parameter estimate of mobile sensor $i$:

[TABLE]

with $\hat{a}_{i}(0)$ being arbitrary; where $\zeta$ is a positive constant, $l_{ij}$ is the weight of the edge between sensors $i$ and $j$. The weight $l_{ij}$ is zero if there is no edge between sensor $i$ and $j$ and positive otherwise. The first term corresponds to the measurement update of mobile sensor $i$ and the second term is a consensus term to ensure that the estimates of all the mobile sensors asymptotically agree or come close to each other. This is critical in establishing the convergence of the estimation error as will be shown below.

Lemma 4**.**

Suppose the mobile sensors translate continuously such that in some time $T>0$, each sensor $i$ passes through each of the centres in the region $\mathcal{Q}_{i}$ so that

[TABLE]

where $\mathcal{K}_{i}^{(i)}(t)$ denotes part of the vector $\mathcal{K}_{i}(t)$ corresponding to the partition (21). Then, we have

[TABLE]

Proof.

Since each mobile sensor $i$ passes through the centres in the region $\mathcal{Q}_{i}$, the union of the trajectories of all mobile sensors cover all the centres, which implies that the matrix

[TABLE]

is positive definite using the same arguments as in proof of corrollary 1 and lemma 1. ∎

Remark 2**.**

Lemma 4 states that each agent passing through the centres in its partition $\mathcal{Q}_{i}$ is sufficient to ensure that the total sum matrix (25) is positive definite.

Now we have the following result:

Theorem 3**.**

Suppose the $N$ mobile sensors adopt the parameter adaptation law (24). Further assume that each mobile sensor $i$ traverses a trajectory going through all the basis function centres in $\mathcal{Q}_{i}$. Then

[TABLE]

for each $i\in\{1,2,\dots,N\}$, i.e. the mobile sensors arrive at a common value for the parameters, the common value being the true parameter value.

Proof.

Consider the function

[TABLE]

Taking the derivative of $V$,

[TABLE]

Substituting for the variables $\Lambda_{i}$, $\lambda_{i}$ and rearranging the second term,

[TABLE]

where $\hat{a}^{\alpha}=[a_{1}^{\alpha}\,a_{2}^{\alpha}\,\dots\,a_{N}^{\alpha}]^{\scriptscriptstyle\top}$ is the vector of the estimate of parameter $\alpha$ of all the sensors stacked together. The function $V$ is lower bounded and non-increasing, and therefore tends to a limit. This implies that $\dot{V}$ is integrable and also that the estimates $\hat{a}_{i}$ are bounded. $\dot{V}$ is also uniformly continuous since the derivative of each term in $\dot{V}$ is bounded. Using Barbalat’s lemma, we conclude that $\dot{V}$ tends to zero as $t\to\infty$. From the second term in $\dot{V}$, noting that $L$ is the laplacian matrix of the connected graph $\mathcal{G}$ with nullspace $k\mathbf{1}$ where $\mathbf{1}$ is the vector of ones and $k\in\mathbb{R}$, we see that as $t\to\infty$, $\hat{a}^{\alpha}\to k_{\alpha}\mathbf{1}$ for some $k_{\alpha}$. Then,

[TABLE]

since $\hat{a}_{i}=[a_{i}^{1}\,a_{i}^{2}\,\dots\,a_{i}^{p}]^{\scriptscriptstyle\top}$. Now from the first term of $\dot{V}$ we have, as $t\to\infty$,

[TABLE]

where $\tilde{a}$ is the common value to which the mobile sensor parameter estimation errors $\tilde{a}_{i}$ converge. Then using lemma 4, it follows that $\lim_{t\to\infty}\tilde{a}=0$ and the parameter estimates converge to the true parameter values. ∎

Remark 3**.**

Although lemma 4 and theorem 3 requires that the mobile sensors move through the centres, the relaxation given in section III-A (requiring that the mobile sensors move only through the neighbourhoods $\mathcal{A}_{j}^{\varepsilon}$ of the centres) also applies here, as well as in all the following results which requires the sensors to move through the centres.

IV-B *Each mobile sensor estimates only part of the

parameter vector*

If the number of parameters $p$ is large as could be the case when the density function is completely unknown, each mobile sensor estimating the entire parameter vector could be computationally intensive, as it would require computing $\left(\frac{p(p+1)}{2}+p\right)$ filter variables in addition to the $p$ parameter estimates. In such cases it would be beneficial to have each mobile sensor estimate only part of the parameters. Suppose each mobile sensor $i$ is to estimate only part of the $a$-vector $a^{(i)}$ given by (21). Now we use $\hat{a}_{i}$ to denote the estimate of $a^{(i)}$ by sensor $i$. We write

[TABLE]

where $\mathcal{K}(q)$ and the parameter $a$ are partitioned appropriately. Since the mobile sensor $i$’s measurement is denoted by $\phi_{i}(t):=\phi(x_{i}(t))$, we have

[TABLE]

where $\mathcal{K}_{i}(t):=\mathcal{K}(x_{i}(t))$ and $\Delta\phi_{i}(t):=\bar{\mathcal{K}}_{i}^{(i)}(t)^{\scriptscriptstyle\top}\bar{a}^{(i)}$. The basis functions in $\bar{\mathcal{K}}_{i}^{(i)}(t)$ are centred outside the region $\mathcal{Q}_{i}$ and thus their values at the points $p_{i}(t)$ are assumed to be small. Under this condition, we consider the contribution to $\phi(.)$ from these terms as a disturbance $\Delta\phi_{i}(t)$.

Let $C=\{c_{1},c_{2},\dots,c_{p}\}$ be the set of centres of the basis functions, $C_{i}\subset C$ be its subset which belongs to $\mathcal{Q}_{i}$. We can then bound $\Delta\phi_{i}(t)$ as follows:

Lemma 5**.**

For each mobile sensor $i$, $i\in\{1,2,\dots,N\}$,

[TABLE]

*where $\delta_{i}:=\max\limits_{j\in\{1,\dots,p\}}\exp\left\{-\frac{d_{i}^{2}}{\sigma_{j}^{2}}\right\}$, $d_{i}:=\mbox{dist}(C_{i},C\setminus C_{i})$, $\mbox{dist}(A,B)=\min\limits_{a\in A,b\in B}\|a-b\|$, and $a_{\max}$ is an upper bound for the parameters, i.e., $|a^{i}|\leq a_{\max}\,\,\forall i\in\{1,2,\dots,p\}$.

Further the bound can be made independent of $i$ as follows,*

[TABLE]

where $\delta=\max\limits_{j\in\{1,\dots,N\}}\delta_{i}$.

Proof.

The lemma follows from the definition of $\Delta\phi_{i}(t)$ using Cauchy-Schwartz inequality. ∎

We again define the following integrators:

[TABLE]

where $s$ is a switching signal which takes values in the set $\{0,1\}$. Consider the following adaptation law:

[TABLE]

Then we have the following result:

Theorem 4**.**

Suppose the $N$ mobile sensors implement the parameter adaptation law (38) with each sensor $i$ only estimating part of the full parameter vector $a^{(i)}$. Further assume that each mobile sensor $i$ produces a trajectory going through all the basis function centres in $\mathcal{Q}_{i}$ in time $T>0$. Then

[TABLE]

where $r_{i}=\frac{Tp\delta_{i}a_{\max}}{\alpha\eta_{i}}$, $a_{\max}$ is the upper bound on the parameter values in $a^{(i)}$, $\alpha\in(0,1)$ and $\eta_{i}$ is the smallest eigen-vlaue of the matrix $\int_{0}^{T}\mathcal{K}_{i}^{(i)}\mathcal{K}_{i}^{{(i)}^{\scriptscriptstyle\top}}d\tau$.

Proof.

Consider

[TABLE]

Taking derivative,

[TABLE]

For $t\geq T$, the first term becomes negative definite (assuming $s>0$). Setting $s=1$ for $t\leq T$ and $s=0$ for $t>T$, we have

[TABLE]

for $t>T$. Then

[TABLE]

where $\kappa=\frac{\eta_{\min}}{\lambda_{\max}(\Gamma^{-1})}$ and $\alpha\in(0,1)$. Thus for $\|\tilde{a}_{i}\|>r_{i}$, we have $\dot{V}\leq-\kappa V$ and $V$ decays exponentially. Therefore the statement of the theorem holds. ∎

IV-C Improving the steady state error

In this section, we propose a strategy to improve the steady state error with the strategy in theorem 4. Note that the strategy in theorem 4 is completely decentralized in that there is no real-time communication required between the mobile sensors to implement the estimation strategy. On the other hand, we can get better parameter estimates at the cost of exchanging information about parameter estimates with other mobile sensors.

The term $\Delta\phi_{i}(t)$ depends on the true value of parameters corresponding to the other mobile sensors (denoted $\bar{a}^{(i)}$). Since we do not know the true values, we cannot cancel this term and treat it as a disturbance. However we know that the other mobile sensors have estimates for the true values of $\bar{a}^{(i)}$. We can use these parameter estimates to reduce the effect of the $\Delta\phi_{i}(t)$ term on the estimation algorithm. Note that the vector $\bar{a}^{(i)}$ consists of the sub-vectors $a^{(j)}$ for all $j\neq i$. Now, corresponding to each $a^{(i)}$, we construct a directed graph with a rooted outbranching (see [29]), denoted $\mathcal{G}_{i}$ which is a sub-graph of the undirected graph $\mathcal{G}$ with mobile sensor $i$ as the root node. An illustration is shown in figure 2.

For each mobile sensor $i$, we introduce additional states $b_{i}^{j}$ for each $j\in\{1,2,\dots,N\}$ and $j\neq i$, which evolves according to the equation

[TABLE]

where we define $b_{i}^{i}:=\hat{a}_{i}$ for ease of notation and $l^{d}_{ik}$ is zero if there is no directed path from node $i$ to $k$ in graph $\mathcal{G}_{j}$, and non-zero constant value otherwise. This implements a directed consensus protocol on the variables $b_{i}^{j}$ with $i=1,2,\dots,N$ (see [29]) converging to the root value $b_{j}^{j}=\hat{a}_{j}$ for each $j$. Thus $b_{i}^{j}$ is an estimator of $\hat{a}_{j}$ by mobile sensor $i$. We now use the modified integrators:

[TABLE]

where $b_{i}$ is the concatenated vector given by $b_{i}=\left[b_{i}^{1^{\scriptscriptstyle\top}}\dots b_{i}^{j^{\scriptscriptstyle\top}}\dots b_{i}^{N^{\scriptscriptstyle\top}}\right]^{\scriptscriptstyle\top}$ ($j=i$ not included). Using the adaptation law (38) we can see that the disturbance term now becomes

[TABLE]

which is expected to be smaller than $\Delta\phi_{i}(t)$, although we cannot put a theoretical bound better than $r_{i}$ in theorem 4. The stability and convergence in case of the above modification is not proved here as it is essentially a similar exercise to that in the previous section. We will investigate the effect of the above modification in section VI.

V Unknown Centres

In this section, we assume as before that the scalar field is a finite linear combination of radial basis functions. We further assume that the centres are not exactly known, but known to within an accuracy of $\epsilon_{c}$, i.e., $\|\hat{c}_{i}-c_{i}\|\leq\epsilon_{c}$. We will evaluate the quality of parameter estimates in this case. Define

[TABLE]

where $\hat{\mathcal{K}}(q)$ is the RBF evaluated at the known values of the centres and $\mathcal{K}(q)$ corresponds to the true values of the centres.

V-A *Each mobile sensor estimates only a part of

the parameter vector*

As in section IV-A, we assume that each mobile sensor estimates part of the parameter vector $a^{(i)}$ corresponding to the partition $\mathcal{Q}_{i}$. In this case we propose the following modified filters,

[TABLE]

with equation (38) as the adaptation law. Then we have the following result.

Proposition 2**.**

Assuming the centres are only known to within an accuracy of $\epsilon_{c}$ ($\|\hat{c}_{i}-c_{i}<\epsilon_{c}\|$), let each mobile sensor pass through the set of known (inaccurate) centres $\hat{c}_{i}$ in $\mathcal{Q}_{i}$. If each mobile sensor implements the adaptation law (38) along with (51)- (52), the estimation error $\tilde{a}_{i}$ converges to within a bound $r_{i}$ of the origin, where $r_{i}=\frac{Tpa_{\max}(\sqrt{p}k\epsilon_{c}+\delta_{i})}{\alpha\eta_{i}}$.

Proof.

Consider the same Lyapunov function as before,

[TABLE]

Taking the time derivative,

[TABLE]

Also note that $|\hat{\mathcal{K}}^{i}(q)|\leq 1\implies\|\hat{\mathcal{K}}(q)\|\leq\sqrt{p}$, and $|\tilde{\mathcal{K}}^{i}(q)|\leq k\epsilon_{c}\implies\|\tilde{\mathcal{K}}(q)\|\leq\sqrt{p}k\epsilon_{c}$ for some $k$ (lipschitz constant), Setting $s=1$ for $t\leq T$ and $s=0$ for $t>T$ as before and, assuming the first term becomes negative definite at time $T$, we now have

[TABLE]

for $t\geq T$. Therefore, the statement of the theorem follows. ∎

V-B *Each mobile sensor estimates the entire parameter

vector*

We define the following filter equations,

[TABLE]

The adaptation law is given by equation (24). In this case, we have the following proposition.

Proposition 3**.**

Suppose the $N$ mobile sensors adopt the parameter adaptation law (24) with the integrators (53)- (54). Also assume that each mobile sensor $i$ produces a trajectory going through all the approximate basis function centres $\hat{c}_{i}$ in $\mathcal{Q}_{i}$. Then the parameter estimation errors of the mobile sensors converge to within a bound $r_{i}$ of origin, where $r_{i}=\frac{Tp\sqrt{p}k\epsilon_{c}a_{\max}}{\alpha\eta_{\min}}$.

Proof.

Consider the lyapunov function

[TABLE]

Taking the derivative of $V$,

[TABLE]

Substituting for the variables $\Lambda_{i}$, $\lambda_{i}$ and rearranging the second term,

[TABLE]

Simplifying,

[TABLE]

for $t\geq T$. We can write the first and last terms in the above equation in terms of stacked vectors as

[TABLE]

where $\tilde{\underline{a}}=\left[\tilde{a}_{1}^{\scriptscriptstyle\top}\,\tilde{a}_{2}^{\scriptscriptstyle\top}\,\dots\,\tilde{a}_{N}^{\scriptscriptstyle\top}\right]^{\scriptscriptstyle\top}$,

[TABLE]

and $P$ is the permutation matrix

[TABLE]

of dimension $Np\times Np$. We show that the matrix $\left(\underline{Q}+P^{\scriptscriptstyle\top}\underline{L}P\right)$ is positive definite. Each of the terms are positive semi-definite. The nullspace of matrix $\underline{L}$ contains elements of the form

[TABLE]

Therefore $P^{\scriptscriptstyle\top}\underline{L}P$ has nullspace elements of the form

[TABLE]

i.e., elements of the form $\left[\,c_{1}\,c_{2}\,\dots\,c_{N}\,c_{1}\,c_{2}\,\dots\,c_{N}\right]^{\scriptscriptstyle\top}$. Correspondingly the $\underline{Q}$ term can be written as

[TABLE]

where $c=\left[\,c_{1}\,c_{2}\,\dots\,c_{N}\right]^{\scriptscriptstyle\top}$. Under the assumptions of the proposition, and lemma 4, the above term is strictly positive. Hence $\left(\underline{Q}+P^{\scriptscriptstyle\top}\underline{L}P\right)$ is positive definite. Let $\eta_{\min}$ be the smallest eigen-value of $\left(\underline{Q}+P^{\scriptscriptstyle\top}\underline{L}P\right)$. Then we have

[TABLE]

for some $\kappa>0$. Thus for $\|\tilde{a}_{i}\|>\frac{Tp\sqrt{p}k\epsilon_{c}a_{\max}}{\alpha\eta_{\min}}$, $V$ decreases exponentially and the result holds. ∎

VI Simulations

In this section, we verify the algorithms presented using simulations. First we consider the exact parameterization case where the true scalar field is a linear combination of RBFs with the centres of the RBFs being known. This case allows us to verify the correctness of the algorithms presented in the paper. Next we consider a scalar field which is completely unknown, and use the algorithms presented to reconstruct the scalar field. The mobile sensors in the simulations are assumed to be single integrators with dynamics given by $\dot{x}_{i}=u_{i}$ where $x_{i}$ is the position of sensor $i$ and $u_{i}$ is its control input. For ease of comparing various algorithms, we refer to the algorithm in section IV-A as Algorithm S1, the algorithm presented in section IV-B as Algorithm S2, and the modified version of algorithm $S2$ in section IV-C as Algorithm S3.

VI-A Exact parameterization

We consider the unit square region $\mathcal{Q}$ with four mobile sensors. The scalar field to be estimated is exactly parameterized in terms of Gaussian RBFs (given by equation (2)), the $x$ and $y$ coordinates of the RBF centres $c_{i}$ being given in table I. The standard deviation of each of the gaussians $\sigma_{i}$ is chosen to be $0.1$. The true parameter values $a^{i}$ are also given in table I.

The scalar field is shown in figure 3.

The initial positions of the mobile sensors were chosen randomly and shown in figure 4. The partition of the region was done by constructing the voronoi cells for each mobile sensor. The Voronoi cell of mobile sensor $i$ (denoted $\mathcal{Q}_{i}$) consists of those points which are closer to sensor $i$ as compared to all other sensors:

[TABLE]

For motion control of the sensors, we use a proportional control law $u_{i}=k(x_{i}-x_{gi})$ where $x_{gi}$ is made to switch between all the centres in the region $\mathcal{Q}_{i}$ making sure the condition in lemma 2 is satisfied. The control gain $k$ was chosen to be $5$. The simulation ran for $16.5$ seconds. The excitation condition was achieved in $T=1.5$ seconds. The reconstructed scalar field with algorithm S$1$ is shown in figure 4 on the right and the average (across all the mobile sensors) parameter estimation error is shown in figure 5. It can be seen that the parameters converge exactly to the true values and exact reconstruction is achieved.

The reconstructed field with algorithm S$2$ and algorithm S$3$ are shown in figure 6. The corresponding estimation errors are shown in figures 7 and 8 respectively.

The maximum parameter estimation error using algorithm S$2$ was found to be $0.030$ and using the algorithm S$3$ was found to be $0.017$. Thus the algorithm S$3$ is seen to give better parameter estimates in this case.

We also present simulation results where we do not know the exact value of the centres of the RBFs (as in section V). We assume we know the centres within an accuracy of $\epsilon_{c}=0.05$. For this, we add a random perturbation (bounded by $\epsilon_{c}$) to the true centre coordinates and use the perturbed centres in the estimation algorithm. The reconstructed fields with algorithms S$1$, S$2$ and S$3$ are shown in figures 9 and 10 respectively. Table 9(a) also compares the maximum steady state parameter errors in the three cases. As expected, algorithm S$1$ has much lower steady state error compared to algorithm S$2$ and algorithm S$3$ performs better than algorithm S$2$. It should be noted that all the algorithms identify the main features of the true field, as seen from the reconstructed field plots.

VI-B Fully unknown scalar field

Now we test the estimation algorithms on a more general scalar field which is not a linear combination of RBFs. For this we consider the continuous scalar field given by

[TABLE]

over the unit square region $\mathcal{Q}$. A plot of $\phi(\cdot)$ is shown in figure 11. We use $N=5$ mobile sensors with the partitions $\mathcal{Q}_{i}$ determined as follows: We first run a uniform coverage algorithm (coverage algorithm presented in [5] with a uniform density function $\phi(q)\equiv 1$). This makes the mobile sensors uniformly spread out in the region $\mathcal{Q}$. We then compute the voronoi partition (55) of the sensors and use it as the required partition $\mathcal{Q}_{i}$.

We first show the results for approximating the field $\phi(\cdot)$ with $p=100$ Gaussian RBFs. The centres of the Gaussian are arranged in a uniform grid over the region $\mathcal{Q}$. The reconstructed field plots for two values of $\sigma_{i}$ (standard deviation of the Gaussian RBFs) are shown in figures 12, 13 and 14 with the three algorithms.

To compare the various algorithms, we use the integral error (see theorem 1)

[TABLE]

where $\hat{a}$ is the final parameter estimate obtained from the given algorithm. The integral error for approximation of $\phi(\cdot)$ using $p=100$ parameters is shown in table II. The table also shows the time $T$ in seconds at which the excitation (positive definiteness) condition is achieved. The total runtime of the estimation algorithms was $T+20$ seconds.

The reconstructed field plots for $p=196$ parameters is shown in figures 15, 16 and 17 with the three algorithms. The comparison of various algorithms is given in table III.

We see that algorithm S$1$ gives better approximation compared to the others as expected. Also the algorithm S$3$ performs significantly better compared to algorithm S$2$. Increasing the number of parameters gives better approximation as expected for algorithm $1$, though for the other algorithms this is not guaranteed due to the extra error incurred (see theorem 4) which may increase with larger $p$ depending on other variables such as the location of centres. $\sigma_{i}$ also plays an important role in the reconstruction of the original field. For $p=100$, $\sigma_{i}=0.05$ seems to provide a better approximation compared to $\sigma_{i}=0.04$, and for $p=196$, $\sigma_{i}=0.04$ seems to provide a better approximation compared to $\sigma_{i}=0.03$. To summarize, algorithm S$1$ gives better approximation compared to the others though it is more computational and memory intensive. The algorithm S$3$ also gives a good approximation requiring much less memory. It may also be noted that in many applications, we may only be interested in identifying the main features of the original field which was successfully done in most of the cases discussed.

VII Conclusion

In this paper we consider the estimation of a scalar field motivated by tools from adaptive control theory and lyapunov analysis. We derived two estimation algorithms, one in which each mobile sensor estimates the entire parameter vector, and another in which each mobile sensor estimates only part of the parameter vector. We verified and tested the algorithms using simulations. Further work involves improving upon the proposed algorithms, and possibility of estimation of time-varying fields by persistent motion of the mobile sensors.