Fully Distributed Informative Planning for Environmental Learning with Multi-Robot Systems

The first part of our work is multi-robot environmental learning in a distributed manner. For environmental learning, some useful techniques exist such as Gaussian mixture model (GMM) [8, 16, 17], finite element method (FEM) [9], and Gaussian process (GP) regression [4, 5, 7]. In particular, GP is a popular approach that derives a spatial relationship between sampled data using a kernel and performs Bayesian inference for prediction at an unknown region.

However, most GP-related researches focus on centralized systems, making it difficult to expand to large-scale multi-robot systems due to network resource limitations such as channel bandwidth and transmit power. Distributed multi-agent Gaussian regression is introduced in [12], which designs a finite-dimensional GP estimator by using Karhunen–Loève (KL) expansion [13]. In contrast to the decentralized GP presented in [10], the distributed GP provides a common copy of the global estimate to all agents by exchanging the estimated information with their neighbors. This paper extends [14], which shows that distributed GPs can construct environmental models using mobile robots in order to take the distributed path planning into account.

The second part of our work is informative path planning in a distributed multi-robot system. As an initial study of informative path planning, the problem of optimal sensor placement has been investigated to create an environmental map in a given space by properly placing a finite number of sensors [15, 16, 17]. Since then, by applying GPs and information theory, the research of optimal sensor placement has grown into the informative path planning research as presented in [5, 8, 18, 19, 20, 21, 22]. Some studies have combined GP with conventional planning algorithms such as rapidly-exploring random tree (RRT) [23], dynamic programming (DP) [20], or Monte Carlo tree search (MCTS) [24].

Besides the above approaches that mainly focus on informative path planning for single agents, many studies have applied informative planning for multi-robot systems. In [10, 7], although both studies deal with decentralized multi-robot exploration using GP, these algorithms are not scalable as they consider only two robots. [33] introduced the combination of the Kalman filter (KF) and the reduced value iteration (RVI) method for the parallelized active information gathering. While this technique is scalable to a large number of robots, it is noted that the environmental model has to be known, and only discrete environments can be represented since the model is expressed in KF. Considering the scalability for multi-robot systems, we extend the MCTS path planning in a distributed manner to be compatible with the distributed GP.

To achieve our goal of fully distributed multi-robot informative planning, we divide the whole process into two phases: environmental learning and path planning. During these phases, we focus on three main contributions as follows.

• We develop an online distributed GP algorithm for environmental learning through Karhunen–Loève expansion and an infinite impulse response filter. This algorithm is capable of learning a dynamic environment.

• We propose a distributed informative path planning algorithm using a distributed MCTS combined with GP. In addition, we introduce the trajectory merging method to consider predicted trajectories of other agents.

• We build a fully distributed exploration and learning architecture using only local peer-to-peer communication for system scalability, as shown in Table I.

We perform a multi-robot exploration simulation with a virtual environment setting and real-world dataset [34] provided by the National Climate Data Center (NCDC) in South Korea as shown in Fig. 1.

The outline of this paper is as follows. Section \@slowromancapii@ briefly describes a multi-robot system setup and preliminaries. Section \@slowromancapiii@ presents a method for online distributed environmental learning. Section \@slowromancapiv@ combines environmental learning and MCTS in the distributed system. Simulations for the synthetic environment and real-world dataset are presented in Section \@slowromancapv@. Section \@slowromancapvi@ concludes the paper.

As depicted in Fig. 1, we consider $n$ robot agents exploring the environment. Each robot $i$ takes the measurement $y_{k}^{i}$ of an unknown environmental process $f(\cdot)$ in its position $\mathbf{x}_{k}^{i}\in\mathbb{X}_{w}$ ($i=1,\cdots,n$) at time $k$ which has the following relationship: {ceqn}

where the measurement of $f(\mathbf{x}_{k}^{i})$ is corrupted by the additive white Gaussian noise $v_{k}^{i}\sim\mathcal{N}(0,\sigma_{v}^{2})$.

Each robot has its process modules, Distributed GP and Distributed MCTS, for the distributed monitoring task. During these processes, they share GP variables and predicted trajectories through a peer-to-peer communication network. This operation process is summarized in Fig. 2. The controller design process is not covered in this work.

To implement the communication network of $n$ robots, we define a set of neighbors for robot $i$ as $\mathcal{N}_{k}^{i}=\{j|\ {||\mathbf{x}_{k}^{i}-\mathbf{x}_{k}^{j}||}<d_{comm},j\in\mathcal{N}/i\}$, where $\mathcal{N}=\{1,2,\cdots,n\}$ is the index set of agents and $d_{comm}$ is the communication range. $\mathcal{N}^{i+}$ means $\mathcal{N}^{i}\cup\{i\}$. For arbitrary variable $A$, $A^{i+}$ means $\{A^{j}\}_{j\in\mathcal{N}^{i+}}$, and $A^{i-}$ means $\{A^{j}\}_{j\in\mathcal{N}^{i}}$ for brevity.

GP regression, which is data-driven non-parametric learning, can provide Bayesian inference over the set $\mathbb{X}_{w}$, taking into account joint Gaussian probability distribution between the sampled dataset [27]. In (1), the unknown process model $f(\cdot)$ is assumed to follow a zero-mean Gaussian process as

$\kappa(\mathbf{x^{\prime}},\mathbf{x^{\prime\prime}})$ is a kernel or covariance function for positions $\mathbf{x^{\prime}},\mathbf{x^{\prime\prime}}\in\mathbb{X}_{w}$, and the original squared exponential (SE) kernel is defined as

where $\sigma_{s}^{2}$ is the signal variance of $f(\cdot)$, and $\Sigma_{l}$ is the length scale. The hyper parameters $\sigma_{s}^{2}$ and $\Sigma_{l}$ can be determined by maximizaing the marginal likelihood [27].

Formally, let $D_{k}^{i}=\{(\mathbf{x}_{t}^{i},y_{t}^{i})\}_{t\in M_{k}}$ be the training dataset sampled by the robot $i$, where $t$ is the sampling time. $M_{k}$ is the set of sampling time indices up to time $k$. With the dataset $D_{k}^{i}$ of size $m_{k}^{i}$, we can simply define the input data matrix as $\mathbf{X}_{k}^{i}=[\mathbf{\bar{x}}_{1}^{i},\cdots,\mathbf{\bar{x}}_{m_{k}^{i}}^{i}]^{\top}\in\mathbb{R}^{m_{k}^{i}\times 3}$ and the output data vector as $\mathbf{y}_{k}^{i}=[\bar{y}_{1}^{i},\cdots,\bar{y}_{m_{k}^{i}}^{i}]^{\top}\in\mathbb{R}^{m_{k}^{i}\times 1}$. According to the test point $\mathbf{x}\in\mathbb{X}_{w}$, the posterior distribution over $f(\mathbf{x})$ by robot $i$ is derived as follows: {ceqn}

where

$\mathbf{K}(\mathbf{X}_{k}^{i},\mathbf{X}_{k}^{i})$ is the $m_{k}^{i}\times m_{k}^{i}$ kernel matrix whose $(u,v)$-th element is $\kappa(\mathbf{\bar{x}}_{u}^{i},\mathbf{\bar{x}}_{v}^{i})$ for $\mathbf{\bar{x}}_{u}^{i},\mathbf{\bar{x}}_{v}^{i}\in\mathbf{X}_{k}^{i}$. $\mathbf{k}(\mathbf{X}_{k}^{i},\mathbf{x})$ is the $m_{k}^{i}\times 1$ column vector that is also obtained in the same way.

To obtain the better description of a spatial process model, robots perform informative path planning. It maximizes the information gain $\mathbb{I}(;)$, which is the mutual information between the process $f$ and measurements $\mathbf{y}$:

where $\textrm{H}(\cdot)$ is the entropy of a random variable. Let $X_{1:k}^{i}$ be the possible trajectory of robot $i$ and $X_{1:k}=\{X_{1:k}^{1},\cdots,X_{1:k}^{n}\}$ be the possible trajectories of all robots. Then, the multi-robot team’s global objective function $J(X_{1:k})$ is defined as follows:

$\mathbf{y}_{1:k}$ is the measurements corresponding to $X_{1:k}$. As a result, the optimal trajectories for all agents are defined as follows:

In this study, the entropy $\textrm{H}(\cdot)$ is obtained using GP. With the result of GP estimation (5), the optimal trajectory generation for each agent will be addressed in Section IV.

Let the usual GP consider $n$ robots. We can simply define the input data matrix for $n$ robots as $\mathbf{X}_{k}=[\mathbf{X}_{k}^{1\top},\cdots,\mathbf{X}_{k}^{n\top}]^{\top}\in\mathbb{R}^{mn\times 3}$. For simplicity, it is assumed that $m_{k}^{i}$’s are same for all robots, and we omit the subscript $k$, so $m^{i}_{k}=m$ hereafter. With the matrix $\mathbf{X}_{k}$, the usual GP requires all the sampled data $\mathbf{X}_{k}$ and inversion of $K(\mathbf{X}_{k},\mathbf{X}_{k})$ with $O((mn)^{3})$ operations. These requirements are impractical when peer-to-peer communication is only used, and the computational burden also increases depending on the data size. For this reason, a new kernel method is needed. The kernel (3) can be expanded in terms of eigenfunctions $\phi_{e}$ and corresponding eigenvalues $\lambda_{e}$ as follows [13]:

where $\lambda_{e}\phi_{e}(\mathbf{x}^{\prime})=\int_{\mathbb{X}_{w}}\kappa(\mathbf{x}^{\prime},\mathbf{x}^{\prime\prime})\phi_{e}(\mathbf{x}^{\prime\prime})d\mu(\mathbf{x}^{\prime\prime})$. It is difficult to derive the kernel eigenfunctions in a closed-form, but the SE kernel expansion has already been obtained via Hermite polynomials, as mentioned in [28]. Then, the process model $f$ for the position $\mathbf{x}\in\mathbb{X}_{w}$ is expanded as

where $a_{e}\sim\mathcal{N}(0,\lambda_{e})$ for $e=1,2,\cdots,\infty$. $f_{E}(\mathbf{x})$ is the $E$-dimensional model of $f(\mathbf{x})$ where $E$ is a constant design parameter. This parameter can be tuned by the SURE strategies [12]. As shown in [28], the optimal $E$-dimensional models can be obtained by a convex combination of the first $E$-kernel eigenfunctions as the size of sampled dataset increases to infinity.

We apply $E$-dimensional approximation to the GP estimator in (5) to derive the estimation of $f_{E}(\mathbf{x})$. According to $E$-dimensional approximation, the kernel function (9) can be described as $\kappa(\mathbf{x}^{\prime},\mathbf{x}^{\prime\prime})\approx\sum\limits_{e=1}^{E}\lambda_{e}\phi_{e}(\mathbf{x}^{\prime})\phi_{e}(\mathbf{x}^{\prime\prime})$. For the input data matrix $\mathbf{X}_{k}$, kernel matrices included in (5) are defined by

where $\Phi(\mathbf{x}):=\begin{bmatrix}\phi_{1}(\mathbf{x}),\cdots,\phi_{E}(\mathbf{x})\end{bmatrix}^{\top}$ and $G:=\begin{bmatrix}\Phi(\mathbf{\bar{x}}_{1}^{1})\cdots\Phi(\mathbf{\bar{x}}_{m}^{1}),\cdots,\Phi(\mathbf{\bar{x}}_{1}^{n})\cdots\Phi(\mathbf{\bar{x}}_{m}^{n})\end{bmatrix}^{\top}$. $\Lambda_{E}$ is the diagonal matrix of kernel eigenvalues.

With (5a) and (13a), the $E$-dimensional estimator for GP is expressed as follows [12]:

where

Because each agent cannot obtain $G$ and $\mathbf{y}$ in (15) without a fully connected network, we decompose the associated terms included in (16) as follows:

where $\alpha_{m}^{i}:=\sum_{t=1}^{m}\Phi(\mathbf{\bar{x}}_{t}^{i})\Phi^{\top}(\mathbf{\bar{x}}_{t}^{i})/m$ and $\beta_{m}^{i}:=\sum_{t=1}^{m}\Phi(\mathbf{x}_{t}^{i})\bar{y}_{t}^{i}/m$ are GP states after the $m$-th sensor measurements. Now (15) is reformulated in the following distributed form:

As the results of average consensus protocol [29], (18) converges to (15) after iterative communication. Similarly, the distributed form of $\Sigma(\mathbf{x})$ in (5b) is expressed as

When we compare (5) with (18) and (20), the computational complexity of the distributed algorithm is $O(E^{3})$, whereas that of the centralized algorithm is $O((mn)^{3})$ due to the matrix inversion [12]. Therefore, the distributed algorithm is more scalable since $E\ll mn$ in general.

If the $(m+1)$-th new training dataset $\{(\mathbf{\bar{x}}_{m+1}^{i},\bar{y}_{m+1}^{i})\}_{i=1}^{n}$ are obtained, $\{\alpha_{m}^{i}\}_{i=1}^{n}$ and $\{\beta_{m}^{i}\}_{i=1}^{n}$ have to be discarded to include new data so that the consensus process must be restarted from scratch. To avoid repeated restarts and keep the continuity of environmental estimate, we introduce the online information fusion algorithm.

Let us assume that the sensor measurement frequencies of all agents are same for convenience. The update rule of $\alpha_{m}^{i}$ and $\beta_{m}^{i}$ is defined as follows:

where $\alpha_{0}^{i}=0$ and $\beta_{0}^{i}=0$. This rule is an infinite impulse response (IIR) filter. If $r=(m+1)^{-1}$, The update rule reflects all dataset equally, so it is suitable for static environmental learning. If $r>(m+1)^{-1}$, this rule reflects more of the recent data, so it is suitable for dynamic environmental learning. Simulation results for each environmental learning are shown in Chapter V.

For each agent $i$, $X_{k+1:k+T}^{i}=(\mathbf{\hat{x}}_{k+1}^{i},\cdots,\mathbf{\hat{x}}_{k+T}^{i})$ denotes the predicted trajectory with the prediction length $T$, or it can be represented by $X_{T}^{i}$ for brevity. Assuming that the agent $i$ receives the predicted trajectories of neighboring agents $X_{T}^{i-}=\{{X_{T}^{j}}\}_{j\in\mathcal{N}_{i}}$, we modify the GP sate $\alpha_{m}^{i}$ as follows: {ceqn}

$\mathrm{n}(X_{T}^{i-})$ is the number of sensing points included in $X_{T}^{i-}$. With (22) and the $E$-dimensional estimator in (20), we define the trajectory-merged GP estimator as follows:

In this way, predicted trajectories of neighboring agents are temporarily included in the acquired data set of the GP estimator. Because (22) and (23) are temporary values for tree search in distributed MCTS, they do not affect $\alpha_{m}^{i}$ and disappear after getting new predicted trajectories. This process is summarized in the Distributed MCTS block of Fig. 2. With this result, the path planning process will be explained in the next section.

As we mentioned in (7), the information gain $\mathbb{I}$ is the objective function we have to maximize. With the definition in (8), the optimal trajectory considering neighboring paths is defined as follows.

$\mathbf{y}_{T}^{i}$ and $\mathbf{y}_{T}^{i-}$ are the measurements corresponding to $X_{T}^{i}$ and $X_{T}^{i-}$, respectively. $\mathbb{X}_{k}^{i}$ is the domain of possible trajectories for agent $i$. As shown in (6), information gain is represented with entropies as follows:

Using the mesuarement model (1) and the entropy calculation for the normal distribution [32], conditional entropy becomes

$\textrm{H}(\mathbf{y}_{T}^{i}\cup\mathbf{y}_{T}^{i-}\cup\mathbf{y}_{1:k})$ is decomposed using conditional entropy, and it is obtained by calculating the entropy of GP as follows:

As a result, with the trajectory-merged GP estimator in (23), the optimal trajectory for agent $i$ is defined as follows:

We call $\mathcal{R}^{i}(\cdot)$ the informational reward function, which is utilized in the tree search algorithm.

Using the informational reward function defined in IV-B, the tree search algorithm iteratively explores and evaluates predictive path candidates according to the discounted upper confidence bound (D-UCB) rule to find the optimal path. D-UCB rule assigns the probabilistic search priority to the action candidates.

The tree structure consists of nodes $s$ and edges $(s,a)$ for all legal actions $a\in\mathcal{A}(s)$. Each edge contains a set of variables $\{N_{s}^{a},W_{s}^{a},\tau_{s}^{a},C_{s}^{a}\}$ where $N_{s}^{a}$ is the visit count, $W_{s}^{a}$ is the total action value, $\tau_{s}^{a}$ is the number of iterations for tree search (shown in line 5 of Algorithm 2), and $C_{s}^{a}$ is a closing variable which will be discussed. We follow the tree search process in Algorithm 1 and the distributed MCTS with GP in Algorithm 2. The MCTS process can be divided into four main steps as follows.

We perform the fully distributed informative planning simulation for multiple agents. They conduct exploration to obtain an estimate of the environmental map, considering collision avoidance and coordination. They can communicate only with neighbors within a range of 10 m (the map size is $20$ m $\times$ $20$ m) and move at $1$ m/s constantly. We set $\sigma_{s}^{2}=1$ and $\Sigma_{l}=\text{diag}([0.02,0.02])$ for the Gaussian kernel (3), and we set $E=80$ for $E$-dimensional estimator (15) and (19).

The progress over time from $0$ to $50$ seconds is shown in Fig. 5. As shown in Figs. 5(a)-(b), twelve agents search the map together and generate the GP estimate presenting the ground truth model in 5(c). The agents scatter naturally and find the next locations to be updated based on the variance map. Also, as they avoid the places where the estimate is already reliable, they can minimize the redundant actions that can decrease the exploration efficiency. Through Fig. 5(b), it can be confirmed that the information of all agents is diffused through a communication link.

Although Figs. 5(a)-(b) show results from agent $\#1$ only, all the distributed GP estimates of each agent converge to the same by the average consensus as shown in Fig. 5(d). In other words, all agents do not simply use local information only in the exploration process but construct a global GP estimation map in a distributed manner.

We conduct more simulations in various environments to investigate the factors that affect search performance. In the tree search algorithm, the search depth $T$ and the number of iterations $N_{iterations}$ are the factors that directly affect the search result. The simulation results in Figs. 6(a)-(b) show that the search performance is proportional to both $T$ and $N_{iterations}$. The deeper the search, the more distant paths are considered. As the number of searching iterations increases, the probability of finding an optimal route increases. Fig. 6(c) shows that the search performance can be improved as the number of agents increases through a distributed algorithm. From these results, we can see that our algorithm is scalable for a large number of robots as well.

Fig. 7 compares the search performance of distributed systems, centralized systems, and independent systems. In the independent system, agents explore the area without communication. The centralized system has a central server that gathers all the information regardless of the communication range, and the server calculates paths for all agents. Because the action space of centralized system $(\mathrm{n}(\mathcal{A})^{n})$ is much bigger than that of distributed system $(\mathrm{n}(\mathcal{A}))$, we set about 22 times more $N_{iterations}$ for the centralized system than the distributed system. Even with limited communication and much fewer iterations, the distributed system performs similarly to the centralized system.

This section presents simulation result for the exploration in a dynamic environment (Fig. 1), using a real-world meteorological dataset. The simulation uses temperature data collected from weather stations in Seoul, South Korea [34]. The reason we choose the weather data for Seoul is that weather stations are densely distributed (the area of Seoul is $605.25\ \text{km}^{2}$). We create a heat map for a ground truth based on the data measured from 0 am to 5 am on July 16th, 2020.

In this scenario, a team of UAVs flies over the search area and gathers temperature data from the current UAV’s location. They fly at a constant velocity of 20 km/h. Because their communication distance is limited to 20 km, sometimes they can be disconnected from one another.

Some snapshots taken during the simulation are shown in Fig. 8. The ground truth over time is shown in Fig. 8(a), and the GP estimation is shown in Fig. 8(b). After 1.7 hours, the estimation result is similar to the ground truth. After that, the results track the true value continuously even when the actual environment changes.