Scalable Gas Sensing, Mapping, and Path Planning via Decentralized Hilbert Maps

A Hilbert map is a continuous probability map developed by formulating the mapping problem as a binary classification task. Let $\mathbf{x}\in \mathcal{W}$ be any point in $\mathcal{W}$ and $Y\in \{0,1\}$ be defined as a categorical random variable such that where the realization is denoted by y. The Hilbert map describes the probabilities $P(Y=1|\mathbf{x})=p(\mathbf{x})$ and $P(Y=0|\mathbf{x})=1-p(\mathbf{x})$ at the position $\mathbf{x}$.

Consider the training measurement data set, $\mathcal{M}={\left\{({\mathbf{x}}_m,y_m)\right\}}_{m=1}^M$, where ${\mathbf{x}}_m\in \mathcal{W}$ indicates the measurement position, the Boolean variable $y_m\in \{1,0\}$ indicates whether the concentration measurement, $\widehat{c}\left({\mathbf{x}}_m\right)$, belongs to the range, $\mathcal{R}$ (where $1=Yes$ and $0=No$), and M is the number of measurements. The probability $P(Y|\mathbf{x})$ is defined as, where, and $f\in \mathcal{H}$ is a Hilbert function defined on $\mathcal{W}$. $\mathcal{H}$ is a reproducing kernel Hilbert space (RKHS) associated with a kernel $k(·,·)$. The kernel mapping is denoted by $k(\mathbf{x},·)=\phi \left(\mathbf{x}\right),\phi \left(\mathbf{x}\right)\in \mathcal{H}$, and is an injective map from $\mathcal{W}$ to $\mathcal{H}$ [28,29,30,31,32]. According to the kernel trick [28,29,30,31,32], the evaluation of the Hilbert function can be expressed in the form of inner product, where ${〈·,·〉}_{\mathcal{H}}$ indicates the inner product in the RKHS. To learn the Hilbert map, the loss function $J_H$ is defined as, where ${\ell }_m=-ln\left[P(Y=y_m|{\mathbf{x}}_m)\right]$ is the negative log-likelihood (NLL) of the data $({\mathbf{x}}_m,y_m)$. Here, the term $\lambda$ is the regularization term, which is a small user-defined positive scalar, $\lambda \ll 1$. Then the gradient of the loss function with respective to f is expressed as, where $\mathsf{\Phi }=\left[\phi \left({\mathbf{x}}_1\right),\dots ,\phi \left({\mathbf{x}}_M\right)\right]$, $\mathbf{p}={\left[p\left({\mathbf{x}}_1\right),\dots ,p\left({\mathbf{x}}_M\right)\right]}^T={\left[\frac{e^{f\left({\mathbf{x}}_1\right)}}{1+e^{f\left({\mathbf{x}}_1\right)}},\dots ,\frac{e^{f\left({\mathbf{x}}_M\right)}}{1+e^{f\left({\mathbf{x}}_M\right)}}\right]}^T$, and $\mathit{y}={\left[y_1,\dots ,y_M\right]}^T$. In addition, the Hessian operator is expressed as, where $\mathbf{I}$ is an identity operator (or matrix) defined in the domain of $\mathcal{H}\times \mathcal{H}$, so that for any function $h\in \mathcal{H}$, and $\mathbf{I}h=h$. Here, $\Lambda$ is an $M\times M$ diagonal matrix defined as, and ${\mathbf{1}}_M={[1,\dots 1]}^T$ is an $M\times 1$ vector with all elements equal to one.

Using the Newton-Raphson method, the Hilbert function, is updated iteratively, where ${\mathbf{R}}_i={\left[{\Lambda }_i{\mathsf{\Phi }}^T\mathsf{\Phi }+\lambda {\mathbf{I}}_M\right]}^{-1}$, ${\mathbf{I}}_M$ is an $M\times M$ identity matrix, and the subscript “i” indicates the ith iteration for learning function f. The Hilbert function, $f\left(\mathbf{x}\right)$, is evaluated iteratively as follows, where $\mathbf{k}={\mathsf{\Phi }}^T\phi \left(\mathbf{x}\right)={\left[k({\mathbf{x}}_1,\mathbf{x}),\dots ,k({\mathbf{x}}_M,\mathbf{x})\right]}^T$.

It can be seen from (18) that the evaluation of $f_{i+1}\left(\mathbf{x}\right)$ only depends on the evaluations of $f_i\left(\mathbf{x}\right)$ and $k({\mathbf{x}}_m,\mathbf{x})$, $m=1,\dots ,M$. Therefore, the evaluation $f_i\left(\mathbf{x}\right)$ for each iteration is not needed if $\mathbf{x}\notin \{{\mathbf{x}}_1,\dots ,{\mathbf{x}}_M\}$. Instead, $f_i\left(\mathbf{x}\right)$ can be calculated at the last iteration directly from the evaluations of $f_i\left({\mathbf{x}}_m\right)$, such that where $f_0\left(\mathbf{x}\right)$ is the initial function evaluation at $\mathbf{x}$, prior to learning the function from $\mathcal{M}$. The following matrix only depends on the measurement data, $\mathcal{M}$, and can be updated iteratively.

Furthermore, consider Q collocation points, ${\mathcal{X}}^c={\{{\mathbf{x}}_q^c\in \mathcal{W}\}}_{q=1}^Q$, in the ROI, characterized by the same spatial interval, labeled by the superscript “c”. Let ${\mathsf{\Phi }}^c=\left[\phi \left({\mathbf{x}}_1^c\right),\dots ,\phi \left({\mathbf{x}}_Q^c\right)\right]$ denote the feature matrix of the collocation points. The evaluations of function $f\left({\mathbf{x}}_q^c\right)$ at all collocation points can be updated by, where ${\mathbf{f}}_i={{\mathsf{\Phi }}^c}^T{f}_i={\left[f_i\left({\mathbf{x}}_1^c\right),\dots ,f_i({\mathbf{x}}_Q^c)\right]}^T$. According to (19), the evaluations, ${\mathbf{f}}_i={{\mathsf{\Phi }}^c}^T{f}_i$, can be updated directly by, where ${\mathbf{f}}_0$ comprises the initial function evaluations at the collocation points prior to learning the function from the measurement data $\mathcal{M}$.

In summary, instead of learning the function f or its coefficients as in traditional KLR or GPR [26], we update the evaluations of $f\left(\mathbf{x}\right)$ at the collocation points, ${\mathcal{X}}^c$. Given the Hilbert function f, the evaluations at the collocation points provide the Hilbert map, $\mathbf{f}$, defined as

The previous subsection develops a method for learning the Hilbert map from the measurement data set $\mathcal{M}$, obtained by the sensors during one time interval. This subsection considers a Hilbert map updated according to the next data set, ${\mathcal{M}}_k={\left\{({\mathbf{x}}_{k,m},y_{k,m})\right\}}_{m=1}^{M_k}$, obtained during the kth time interval. Here, $M_k$ is the number of measurements in the kth time interval, and ${\mathcal{X}}_k={\left\{{\mathbf{x}}_{k,m}\right\}}_{m=1}^{M_k}$ and ${\mathcal{Y}}_k={\left\{y_{k,m}\right\}}_{m=1}^{M_k}$ are the measurement positions and the corresponding classification estimates, respectively. To learn the next Hilbert map, the loss function over a period T can be expressed as, where $\lambda$ is the regularization term as in (13), $\gamma (T-k)$ is the “forgetting factor”, and ${\ell }_{k,m}=-ln\left[P(Y=y_{k,m}|{\mathbf{x}}_{k,m})\right]$ is the NLL of the data $({\mathbf{x}}_{k,m},y_{k,m})$. Similarly to (14), the gradient is expressed as, where ${\mathsf{\Phi }}_k=\left[\phi \left({\mathbf{x}}_{k,1}\right),\dots ,\phi \left({\mathbf{x}}_{k,M_k}\right)\right]$, ${\mathbf{p}}_k={\left[p\left({\mathbf{x}}_{k,1}\right),\dots ,p\left({\mathbf{x}}_{k,M_k}\right)\right]}^T={\left[\frac{e^{f\left({\mathbf{x}}_{k,1}\right)}}{1+e^{f\left({\mathbf{x}}_{k,1}\right)}},\dots ,\frac{e^{f\left({\mathbf{x}}_{k,M_k}\right)}}{1+e^{f\left({\mathbf{x}}_{k,M_k}\right)}}\right]}^T$, ${\mathbf{y}}_k={\left[y_{k,1},\dots ,y_{k,M_k}\right]}^T$ for $k=1,\dots ,T$, and ${\tilde{\mathsf{\Phi }}}_T=\left[{\mathsf{\Phi }}_1,\dots ,{\mathsf{\Phi }}_T\right]$. Furthermore, ${\Gamma }_T$ is a diagonal matrix obtained by placing the vector, ${\left[\begin{array}{ccccc}\gamma (T-1){\mathbf{1}}_{M_1}^T& \dots & \gamma (T-k){\mathbf{1}}_{M_k}^T& \dots & \gamma \left(0\right){\mathbf{1}}_{M_T}^T\end{array}\right]}^T$ on the diagonal of a zero matrix. In addition, as with (15), the Hessian operator can be expressed as, where, and ${\Lambda }_k={\mathbf{p}}_k{({\mathbf{1}}_{M_k}-{\mathbf{p}}_k)}^T$.

Then, the update rule for learning the Hilbert function, f, at the ith iteration is given by, where ${\tilde{\mathbf{R}}}_i={\left[{\Gamma }_T{\tilde{\Lambda }}_{T,i}{\tilde{\mathsf{\Phi }}}_T\tilde{\mathsf{\Phi }}+\lambda {\mathbf{I}}_{M_{tol}}\right]}^{-1}$ and $M_{tol}={\sum }_{k=1}^T{M}_k$. According to the above equation, the function, f, is expressed by using the set of all the measurement positions, $\left\{{\bigcup }_{k=1}^T{\mathcal{X}}_k\right\}$, and the corresponding coefficients.

To reduce the computational load, the forgetting factor, $\gamma (T-k)$, is modeled by, such that each robot stores only the measurements obtained during the past $\tau$ time steps. By setting $\tau =1$, the update rule (28) for learning the function f at the ith iteration is expressed as, where ${\mathbf{R}}_{T,i}={\left[{\Lambda }_i{\mathsf{\Phi }}_T^T{\mathsf{\Phi }}_T+\lambda {\mathbf{I}}_{M_T}\right]}^{-1}$.

The previous subsections present a method for learning the Hilbert function, f, from the available measurement data, comprising the measurement position data set ${\mathcal{X}}_k={\left\{{\mathbf{x}}_{k,m}\right\}}_{m=1}^{M_k}$ and corresponding classification estimates ${\mathcal{Y}}_k={\left\{y_{k,m}\right\}}_{m=1}^{M_k}$ for $k=1,\dots ,T$ that indicate whether the concentration measurement at the measurement position belongs to the range, $\mathcal{R}$, defined in (9). If L classification estimates are obtained from the same concentration measurements, indicating that the concentration measurement belongs to the L concentration ranges, such as ${\mathcal{R}}_1,\dots ,{\mathcal{R}}_L$, respectively, using the learning method described in Section 3.2, one can approximate L Hilbert functions, $f_1,\dots ,f_L$. Furthermore, the probabilities in (6) can be evaluated at all the collocation points, ${\mathcal{X}}^c={\left\{{\mathbf{x}}_q\right\}}_{q=1}^Q$, such that

Now, consider Hilbert functions, ${\left\{f_l\right\}}_{l=1}^L$, approximated locally by each robot. Although all the concentration intervals are mutually exclusive and complete, i.e., it is not guaranteed that the sum of all the learned probabilities, ${\left\{p_l\left({\mathbf{x}}_q^c\right)\right\}}_{l=1}^L$, is equal to one, or ${\sum }_{l=1}^L{p}_l\left({\mathbf{x}}_q^c\right)=1$. Therefore, the learned probabilities must be normalized to obtain the discrete distribution, $\mathit{\pi }\left({\mathbf{x}}_q^c\right)=[{\pi }_1\left({\mathbf{x}}_q^c\right),\dots ,{\pi }_L\left({\mathbf{x}}_q^c\right)]$. Each component of the discrete distribution is calculated by, where $p_0\left({\mathbf{x}}_q^c\right)$ is a normalization term. Let $\overline{\mathbf{c}}=[{\overline{c}}_1,\dots ,{\overline{c}}_L]$ denote the medians of these intervals, where ${\overline{c}}_l=(H_l-L_l)/2$. Then, the expectation of the random variable $C\left({\mathbf{x}}_q^c\right)$ can be expressed as, where $\mathbb{E}(·)$ is the expectation operator.

Consider two robots that have each learned their respective Hilbert functions, $f_1$ and $f_2$, from two different sets of measurement data, ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$, respectively. Then, the fused Hilbert function, $f_F$, can be obtained based on the following theorem.

The proof of Theorem 1 is provided in the Appendix A.1. According to Theorem 1, the following corollary can be obtained.

Because the prior concentration distribution is unknown, the Hilbert functions are fused according to (41) in Corollary 1, unless otherwise stated. Furthermore, assume that all the robots have the same collocation points, ${\mathcal{X}}^c$. The information fusion can be implemented by fusing the Hilbert maps among neighboring robots, where ${\mathbf{f}}_{H,1}$ and ${\mathbf{f}}_{H,2}$, and ${\mathbf{f}}_{H,F}$ are the Hilbert maps associated with the Hilbert functions $f_1\left(\mathbf{x}\right)$, $f_2\left(\mathbf{x}\right)$, and $f_F\left(\mathbf{x}\right)$, respectively.

Any pair of robots in the network can share and update their Hilbert maps efficiently according to (42). To implement data fusion for a very large number of robots, however, a communication strategy is also required to determine how to share data between robots characterized by active communication links. In this paper, gas measurement data are communicated at every time step $T_c$. The communication protocol requires four steps at every time k, as follows. In the first, the N robots are partitioned into smaller communication networks according to their positions and communication range $r_c$, and ${\mathcal{G}}_{ı}$ denotes the index set of the robot in the ıth communication network. Then, for any robot, $a_n$, with $n\in {\mathcal{G}}_{ı}$, there exists another robot, $a_{n^{\prime }}$, such that where $d_{n,n^{\prime }}$ is the distance between robots $a_n$ and $a_{n^{\prime }}$, and $\parallel ·\parallel$ indicates the Euclidean norm.

As a second step, one robot in every communication network is selected as a temporary data center (TDC) denoted by $a_{n_{ı}^{*}}$, $n_{ı}^{*}\in {\mathcal{G}}_{ı}$. The other robots in ${\mathcal{G}}_{ı}$ send the Hilbert map change $\Delta {\mathbf{f}}_{n,k}$ to robot $a_{n_{ı}^{*}}$. The Hilbert map change, $\Delta {\mathbf{f}}_{n,k}$, is defined as the change between the nth robot’s Hilbert map at the current time step k and its Hilbert map at the previous communication time step, $k-T_c$, such that where ${\mathbf{f}}_{n,k}$ and ${\mathbf{f}}_{n,k-T_c}$ denote the Hilbert maps of the nth robot at times kth and $(k-T_c)$th, respectively.

In the third step, the sum of the Hilbert map changes obtained from all robots in ${\mathcal{G}}_{ı}$, defined as, is communicated back to the other robots, except the TDC robot, $a_{n_{ı}^{*}}$. Finally, in the fourth step, all the robots update their own Hilbert maps by adding the received total Hilbert map changes, $\Delta {\mathbf{f}}_{G_{ı},k}$, to the current Hilbert maps, such that where ${\mathbf{f}}_{n,k^{+}}$ represents the Hilbert map after the data fusion process.

An information-driven approach is developed by planning the path of the robots such that they move towards collocation points with higher entropy. The collocation points are treated as goal positions characterized by attractive artificial potential fields defined as, where the superscript “att” indicates the attractive field. The corresponding attractive gradient is expressed as

Similarly to classic artificial potential field methods, the repulsive potential functions $U^{rep}$ generated from the other robots are also considered, such that where $r_{rep}$ is the distance threshold to create a repulsion effect on the robot. The repulsive gradient is expressed as

Using (52) and (54), the total potential gradient for the nth robot is expressed as, where $\xi$ and $\eta$ are user-defined coefficients.

Based on the total gradient, ${\mathbf{g}}^{tol}\left({\mathbf{x}}_n\right)$, the nth robot can be controlled to visit the collocation points with higher uncertainty and avoid collisions with other robots. The algorithm is developed based on the entropy attraction, which is named as EAPF algorithm.

The particle swarm optimization (PSO) algorithm proposed by Clerc and Kennedy in [33] and its variants use a “constriction coefficient” to prevent the “explosion behavior” of the particles, and have been successfully applied to GDM and gas source localization problems [12,13]. In the original PSO algorithm and its variants, the concentration measurements are used to update the robot controls. Considering the different objective function in (5), an entropy-based PSO (EPSO) is proposed.

At the kth time step, the update of the nth robot position, ${\mathbf{x}}_{n,k}$, can be described as where ${\mathit{\nu }}_{n,k}$ represents the velocity of the nth robot at the kth time step, ${\mathbf{b}}_{n,k}^{neig}$ and ${\mathbf{b}}_{n,k}^{glob}$ are the best neighboring and global collocation points, respectively. The best neighboring collocation point, ${\mathbf{b}}_{n,k}^{neig}$, is determined by, where $r_{neig}$ is a coefficient which specifies the neighbor area from ${\mathbf{x}}_{n,k}$. The best global collocation point, ${\mathbf{b}}_{n,k}^{glob}$, is determined by

The learning coefficients, ${\psi }_1\in [0,{\overline{\psi }}_1]$ and ${\psi }_2\in [0,{\overline{\psi }}_2]$, are two uniform random variables. The constant parameter $\chi >0$ prevents the explosion behavior. For efficient performance and prevention of the explosion behavior in (56), the parameter settings of the learning coefficients proposed in [12,13,33] is applied. The constriction parameter $\chi >0$ is calculated by (refer to [33]), where $\psi ={\overline{\psi }}_1+{\overline{\psi }}_2$ and $\kappa \in [0,1]$.

The decentralized sensing system consists of a network of $N=100$ robots characterized by single integrator dynamics, where ${\mathbf{x}}_i={[x,y]}^T$ is the robot state vector, $x,y$ are the inertial coordinates, and ${\mathbf{u}}_i={[u_1,u_2]}^T$ is the robot control vector comprised of the x- and y-velocity components. The robot state/position is assumed to be observable by a built-in GPS with zero-mean measure noise w, where w is Gaussian white noise $\mathcal{N}(\mathbf{0},\mathsf{\Sigma })$ with $\mathsf{\Sigma }=0.05\times {\mathbf{I}}_2$. The above distributed sensing system is tasked with mapping a gas distribution in an ROI $\mathcal{W}=[0,L_x]\times [0,L_y]$ where $L_x=200$ m and $L_y=160$ m, with an unknown gas distribution shown in Figure 1, and over a fixed time interval $[0,T_f]$, where $T_f=500$ min. The normalized gas concentration range $\mathcal{R}=[0,100]$ (Figure 1) is divided into $L=3$ intervals: ${\mathcal{R}}_1=[0,30)$, ${\mathcal{R}}_2=[30,70)$, and ${\mathcal{R}}_3=[70,100]$, representing low-hazard, medium-hazard, and high-hazard concentrations, respectively.

The robots are initially deployed at four corners in the ROI by sampling from a given Gaussian Mixed Model with 4 components where ${\mu }_1={[20,20]}^T$, ${\mu }_2={[20,140]}^T$, ${\mu }_3={[180,20]}^T$ and ${\mu }_4={[180,140]}^T$, and identical covariance matrices, $\mathsf{\Sigma }=\left[\begin{array}{cc}10& 0\\ 0& 10\end{array}\right]$. The initial robot positions are shown as red points in Figure 1. Each robot is equipped with a small metal oxide sensor array comprised by $M=5\times 5=25$ gas sensors, where the spatial intervals between two sensors are all 10 cm. The FOV of each sensor array covers an area of $40\times 40$ cm${}^2$ in the ROI, which is very small relative to the whole ROI. Assume that the measurement noise of the gas concentration $V\left(\mathbf{x}\right)$ is also characterized by Gaussian white noise, $\mathcal{N}(0,{\mathsf{\Sigma }}_c\left(\mathbf{x}\right))$ with the covariance matrix ${\mathsf{\Sigma }}_c\left(\mathbf{x}\right)=0.5\times {\mathbf{I}}_2$ everywhere in $\mathcal{W}$. All of robots communicate with each other at the same time with a communication period $T_c=10$ min. The communication range is $r_c=30$ m. In addition, $Q=100\times 100$ virtual collocation points are evenly deployed in the ROI to generate the Hilbert maps, where the Gaussian kernel is used for Hilbert function learning with a kernel size of $\sigma =2$ m, and the parameter $\tau$ is set to 3.