Understanding natural phenomena is crucial in many fields of science and technology. However, collecting data with stationary sensors is often costly and time-consuming. Mobile robotic sensor networks (MRSN) offer a new way to generate such spatio-temporal (ST) data since MRSN can be quickly deployed and target data collection in areas of high information value. With the development of unmanned vehicles in all fields (ground, surface water, underwater, and air), this approach can solve many monitoring and observation tasks [1]- [3]. Compared to sensing by stationary sensor nodes, the main challenge of using the data collected by MRSN to model ST phenomena is that the number of observation positions, hence the number of measurements taken at a sampling instant, is restricted by the limited number of robots and their mobility constraints. Thus, successful ST mapping solutions with mobile robots must take into account these limitations.
Recently, Gaussian process regression (GPR) has received significant interest as a technique for discovering ST data correlations. GPR provides a fundamental framework for nonlinear non-parametric Bayesian inference widely used in soil organic matter mapping [4], temperature mapping [5] and leakage detection [6]. The use of non-parametric models opens possibilities for mapping solutions to remain generic and flexible, since hyperparameters are able to be adjusted to create more accurate practical models for some specific applications. GPR also provide an estimation of forecast uncertainties and provides opportunities for future planning algorithms focusing on uncertainties. This spatial and temporal mapping technique GPR can be used in any situation in which a mobile robot is faced with a phenomenon that differs in time and space. For example, as an typical exploration task, GPR for spatial temporal maps can be established for understanding temperature, chemical concentration and water flow. In the field of robotics, it may be extremely valuable for precise control to have the ability to model and predict environmental disturbances and then to make appropriate strategies against the impacts of disturbances.
With the help of MRSN, ST mapping has been intensively investigated to observe and model temporal changes of unknown environments. The authors in [7] take advantage of mobile sensors to build a map of spatio-temporal phenomena via GPR; however, the restrictions on the movement of mobile sensors were not considered. The study in [8] attempted to capture a slow-changing phenomenon in real-time operation, but it is assumed that the phenomenon is static during robot measurements. Recently, the authors in [9] presented spatialtemporal mapping with observations from a single robot traversing on a fixed-path design.
Among robot planning methods in exploration tasks, informative path planning (IPP) [10] has excellent performance, as future paths are generated by estimated environment models. The core idea of the IPP is based on minimizing prediction uncertainties, which leads to designing optimal routes to collect measurements. In the literature, the centralized IPP can be found in [11], [12], where the observed data is collected in a central unit to update a surrogate model. Then, optimal paths are computed and sent to each robot. These works meet inherent restrictions since a tremendous amount of data collected by many robots possibly results in congestion in both communication and computation. In recent years, several studies have been devoted to distributed IPP [13], [14] with regard to GPR. However, none of them takes into account the mapping of spatial-temporal phenomena and the dynamics of actual robots. Furthermore, the cost functions of the IPP used in these studies are separated, i.e., each robot has its own cost function related only to its future path. This setup ignores the cross-relation between the future paths of neighbor robots.
Motivated by the above discussion, this article presents a new distributed IPP approach for mapping spatial-temporal fields by using multiple robots. In other words, we propose a spatially and temporally distributed prediction scheme based on GPR while the connectivity of the robot team is preserved during their movements. Another key contribution of this paper is the novel local cost functions for the IPP optimization problem with respect to the future paths of neighbor robots.
The proposed spatially temporally distributed IPP approach was validated by mapping spatio-temporal temperature using a real-world dataset. The organization of this paper is as follows. Section II briefly presents models of mobile robots for monitoring a spatio-temporal fields. The IPP optimality with connectivity preservation is then presented in Section III. Next, Section IV describes the distributed implementation of the proposed IPP algorithm. Finally, simulation results obtained by implementing the proposed approach using the real-world dataset in synthetic environments are discussed in Section V. The conclusions are described in Section VI.
Notations: Let us denote N and R as the sets of natural and real numbers, respectively, → as the Kronecker product, and 1 n ↑ R n as a vector in which each element is 1. With a set of integers (index set)
↑ as a block of matrices M i with appropriate dimension.
SPATIO-TEMPORAL FIELDS Consider a convex set Q ↑ R ω standing for an operation space of all robots. Let us define n as a number of robots working in Q, and we assume that the communication area of each robot i at anytime is a ball (or a circle in 2D) centered at p i,k with radius R. Here, p i,k is the location of robot i at time t k . In this paper, the spatio-temporal field of interest is considered as a latent relationship z : (Q, R + ) ↓ R mapping a location of measurement in Q and its current time t k to a spatio-temporal phenomenon. The robot i observes a noisy measurement y i,k ↑ R of the spatio-temporal field z at its current position for every time step. In addition, robot movements are described as
where u i,k stands for the bounded control input (↔u i,k ↔ ↓ ↗ ū) of robot i between two consecutive time steps with ū being the maximum magnitude of each element in control input vector u i,k . Additionally, A k and B k represent the matrices obtained by linearizing the dynamics of the robot at time step k.
Based on the above setups, let us define that robots i and j are connected at step k if ↔p i,k ↘ p j,k ↔ 2 ↗ R. Accordingly, let E k be the set of pairs of connected robots (i, j) at time step k, that is,
, . . . , M} as a set of indexed vertices in which each robot represents a vertex. Then, let G k be an undirected graph established by set of vertices V and edges E k . Note that G k varies over time. The undirected graph G k is connected if there exists at least a path between any pair of robots. Robot i is considered a neighbor of robot j if they are connected.
Denote N i,k as a set of neighbors of robot i and let y i,k be a measured value of the robot i at the time t k at location p i,k . Denote D i,k as a dataset of robot i collected up to time t k . The local data set D i,k can be decomposed from the sets D y i,k , D p i,k , D ε i,k of all measurements, locations, and timestamps. Consequently, the data exchange in robot
In this setup, each robot has a measurement model as follow
) is an independent and identically distributed zero-mean Gaussian noise with standard deviation
is the random/latent variable with covariance funcion K and mean µ which can be set as a deterministic function (constant, polynomial or periodic) or determined by observed data such as neural network [15], [16]. After moving to new locations, p 1,k+1 , . . . , p n,k+1 , at the next step t k+1 , the robots collect additional measurements from the spatio-temporal phenomenon. The Gaussian process model, z ⇒ GP (•, •), is then updated with the newly acquired data. The central challenge in this context is determining the optimal new locations for measurement collection.
The movements of robots possibly disrupt the connectivity of sensor network. Thus, this section proposes a distributed algorithm to ensure that the connectivity of robot network is preserved in the next step. To be specific, if the network is currently connected, then by maintaining some edges, the network will be connected in the next step. We then formulate the IPP optimization problem given the robot dynamics and connectivity constraints.
At the beginning, we recall previous results in [14]. Let us define
). In this paper, we assume that the graph G 0 is connected at the initial time t 0 .
With the help of Lemma 3.1, at current time k, define a set S i,k of neighbors of robot i that robot i attempts to preserve its connectivity to robots in S i,k at time k + 1. Then, we nail down Algorithm 1 for finding S i,k such that the robot network at the next time step k + 1 is connected. The following Theorem presents a distributed mechanism to preserve the connectivity of robot network.
end if 6: end for determined by Algorithm 1 for all i ↑ V. Then G k+1 is connected.
Proof: Based on Lemma 3.1, the proof of Theorem 3.2 follows the same steps with similar arguments as in the proof of [14,Lemma 3.6].
Apart from this, let y D i,k be a vector of all measurements that the robot i took up to time step k, then the vector of local measurements y D i,k follows a multivariate Gaussian distribution in the following form
where µ D i,k is a mean vector with regard to the local dataset
k I denotes a covariance matrix with noise term. In what follows, for unobserved locations of interest pi,H = [p ↑ i,k+1 , . . . , p↑ i,k+H ] ↑ ↑ R ω↘H (4) corresponds to specific time instants t H = [t k+1 , . . . , t k+H ] ↑ and an unobserved vector of measurements ŷi,H = [ŷ(p i,k+1 , t k+1 ), . . . , ŷ(p i,k+H , t k+H )] ↑ . Here, H = {0, 1, . . . , H ↘ 1} with 0 < H ↑ N represents the predictive horizon. Additionally, let us use (•) for an unobserved vector of latent variables, locations, and their corresponding covariance matrices. Then, following a multivariate Gaussian distribution, it has
in which matrices !D i,k H , !i,H are obtained from a spatiotemporal covariance function K(p, p ↗ , t, t ↗ ) with respect to the data set D i,k and unobserved locations pi,H . In addition, μi,H denotes mean vectors with respect to unobserved locations pi,H . According to [17], the conditional distribution of unobserved positions is ŷi,
where matrices μH|D i,k , !i,H|D i,k are given by μi,H|D
) Let P k = {s 1 , s 2 , . . . , s q } ↑ Q q be a set of distinguished locations in which all robots can measure in H consecutive sampling instants t k+1 , t k+2 , . . . , t k+H . With a very large q, P k can be chosen such that for arbitrary s ↑ Q and small ϖ > 0 (depends on q), we have dist(s, P k ) = min e→P k ↔e ↘ s↔ 2 < ϖ. Let T k,H = {t k+1 , t k+2 , . . . , t k+H } be a set of the corresponding time stamps. Based on such the setups, we present a IPP for a robot.
where S(•) is the conditional entropy [18] with respect to the vector of measurements U k,H = {ŷ(s, ϱ)|s ↑ P i,k , ϱ ↑ T k,H } corresponding to P i,k and T k,H . By using the chain rule for conditional entropy [18], we have
It should be noted that ŷi,H is a vector of latent variables at locations pi,H and time t k+1 , . . . , t k+H . We assume that pi,H ↑ P i,k then ŷi,H is contained in U k,H (t k+1 , . . . , t k+H ↑ T k,H ). Then,
) is constant. Therefore, it can be clearly seen that ( 8) is concerted to pi,H = arg max S(ŷ i,H |D y i,k ).
The conditional entropy of a multivariate Gaussian distribution of random variables ŷi,H at unobserved locations pi,H at time t H is given by a closed form [18]:
When robot i determines its future path pi,H independently (by solving (9)), the informative correlation between the optimal paths will be ignored, i.e. , for the neighbor j of robot i, pj,H is not considered for obtaining pi,H . Based on Problem 3.3, we modify (8) to set up an informative planning problem for multiple robots in a distributed manner. To begin with, define p+ i,H = pj,H j→N + i,k and ŷ+ i,H = ŷj,H j→N + i,k . Problem 3.4 (Multiple robot informative planning): Find optimal paths p+ i,H of robots i = 1, 2, . . . , n in the mobile robot network at time k, leading to the lowest uncertainties at all unobserved locations of interest p1,H , . . . , pn,H = arg min n i=1
Unlike ( 8), we use future measurements of neighbors ŷ+ i,H in (10) to evaluate uncertainties. Consequently, the conditional uncertainties at unobserved locations will have better local evaluations. In the same manner as derivation of (9), the optimization problem (10) is equivalent to p1,H , . . . , pn,H = arg max
where !+ i,H|D i,k is the covariance matrix associated to conditional distribution of ẑ+ i,H determined via (7b) by using p+ i,H instead of pi,H .
In the light of ( 11), Problem 3.4 for optimally mapping a spatio-temporal field is formulated in the following optimization problem: for all h ↑ H max
(12d) In the above optimization, constraints (12b) and (12c) denotes system dynamics and their physical limitations, and (12d) stands for connectivity conditions. Remark 3.5: As seen from (12a), log det !+ i,H|D i,k stands for the local cost function of robot i with respect to its optimal future locations pi,H and its neighbors pj,H (see Fig. 1). Compared to previous work [4], [14], the correlation between the data collected by the robot i and its neighbor's future paths is exploited. Thus, each robot has an awareness of the future movements of its neighbors for planning.
It should be noted that the objective function (12a) and the connectivity constraint (12d) are involved in at least two robots. To solve the optimization in a distributed way, each robot should create copies of its neighbors. Let
↑ (for all j ↑ N i,k ) represent the virtual positions of the robot j computed by the robot i via solving the optimization problem (12a). Denote a vector
where ω ii = pi,H . Our objective here is to achieve ω ij,H = pj,H (for all j ↑ N i,k ). With the above definitions, the optimization problem (12a) is equivalent to:
) for all h ↑ H. For the sake of simplicity, we define
The optimization problem (13a) can be solved in a distributed fashion by using proximal alternating direction method of multiplier (proximal ADMM) presented in [12], [19]. However, the method requires gradient updates ∝f i for every iteration of the ADMM in a single time step. It should be noted that the computational complexity of ∝f i is O(d|D i,k | 2 ) where d the input dimension and |D i,k | is the number of local data, therefore updating ∝f i at every iteration accounts for many computing resources.
The cost function (13a) is highly nonconvex, resulting in a great computational burden. Thus, let us approximate (convexify) the cost function around the previous value
2 where r > 0 denotes a regulation parameter and ∝f i (ω k↔1 i ) represent the gradient of f i at ω k↔1 i . In the initial step, ω k↔1 i is selected from the initial position Algorithm 2 Distributed solving of optimization Input: Number of robots n, set of neighbors N i,k , and a small tolerate error ϖ.
neighbors and receives ε (n) ji from them 4: Compute ω (n+1) i by (17) 5: i sends ω (n) ij to neighbors and receives ω (n) ji from them 6: Compute ε (n+1) i by (16) 7:
end if 10: end loop of the robots, that is, ω 0 ii,h = p i,0 and ω 0 ij,h = p j,0 for all h ↑ H, j ↑ N i,0 . To simplify local constraints, let us define
as a set of local constraints including robot dynamics and network connectivity. At the initial step, positions
↑ B i,0 for all i because we already assumed that G 0 is connected. Consequently, B i,0 is a non-empty set, and (12a) is feasible at the initial time step, and then
Sequentially with the next steps k + 1, B i,k+1 is also nonempty. It can be observed that B i,k is a convex set. The optimization (13a) can be rewritten as
for all i ↑ V and j ↑ N i,k . Next, the Lagrangian function of ( 14) is defined by
is the dual variable. Then, using the dual decomposition method [20], the optimization ( 14) is handled by
The optimization problem (15) can be distributively solved. Indeed, the Lagrangian is written as
= arg min
where "
↑ . Note that ( 17) is formulated as a convex quadratic programming quadratic constraints (QCQP) that is solved effectively in polynomial time by solvers such as OSQP [21] or SOCP [22]. The distributed algorithm 2 is tailored to describe the steps to solve the optimization problem (14).
To demonstrate the effectiveness of the proposed approach, we implemented it in a synthetic environment by using the real-world temperature dataset [23]. It is noted that the temperature dataset was spatially and temporally collected by 12 fixed-location sensors during 24 hours in a crop area of 20 m ≃ 100 m, which resulted in 756 measurements in total. To verify our IPP algorithm, we first simulated a spatio-temporal field of the real-world temperature dataset by building a model from all 756 measurements. This model is called ground truth (GT). Then whenever a robot moves to a particular location in an unknown area and takes a virtual measurement at a particular time, the GT model would estimate that virtual measurement for the robot. In other words, the robots in our simulation virtually took the "real" measurements when they were exploring the field.
In the simulations, 6 robots with a communication range of R = 20 [m] were chosen to conduct a task of mapping the spatio-temporal temperature field in an unknown area with the same dimensions of 20 m ≃ 100 m. At the beginning, none of the robots knew anything about the field. They could only gather temperature information over their navigation. In addition, let us take
The covariance function was selected as a serial combination between square exponential and Matérn ( 12 ) functions as
, where ϑ s > 0 and ϑ t > 0 are spatial and temporal length scales, respectively; and ε > 0 is a scalar scale. We also considered wheeled mobile robots with the dynamics ṗi = v i cos φ i sin φ i ↑ (18) where v i is the longitude velocity and φ i is the heading angle of robot i. The dynamics of the robots are linearized by using first-order approximation and discretizing (18) by the Euler method as pi,k+h = pi,k+h↔1 + ϱ cos φ i,k↔1 , sin φ i,k↔1 v i,k+h↔1
+ ϱ v i,k↔1 ↘ sin φ i,k↔1 , cos φ i,k↔1 #φ i,k+h↔1 , where #φ i,k = φ i,k+h↔1 ↘ φ i,k↔1 . With respect to (1), we have
In the IPP context, a group of robots aims to map a spatio-temporal field in an unknown environment. The robots navigate the environment while taking measurements at their moving steps. And our proposed distributed IPP algorithm provides the robots with optimal navigation in terms of gaining maximal information of the field in both space and time. In other words, the measurements taken by the robots run by our IPP approach carry most informative content of the spatio-temporal temperature field. To validate this fact, we exploited the measurements collected by the robots along their navigation paths and learned a GP model. It is noticed that this GP model was updated after every moving step of the robots as the temperature field was varying over both space and time. Since the robots could traverse to only a limited number of locations in the environment, we utilized the learned GP model to predict the temperature in the whole space at any expected time. Mappings of the spatio-temporal temperature field in the whole environment over moving steps of the robots are illustrated in the right column of Fig. 2. For the comparison purposes, we also generated the mappings of the field by using the GT model, which are depicted in the left column of Fig. 2. As can be seen from Fig. 2, our method provides the comparative results that the robots could build the spatio-temporal maps intensively comparable to the ground truth. It is also demonstrated in the right column of Fig. 2 that connectivity of the robots was well maintained overtime. Code for the simulation is written in Julia and can be found in github.com/AACLab/SpaTemIPP.git.
In practice, apart from mapping a spatio-temporal field in a whole environment, one may be interested in values of the field at some specific locations. Of course, these specific locations are not accessible by robots; hence no measurement can be made. In that case, we can use the learned GP model to temporally predict the field at those locations. To verify efficacy of our algorithm in location level, we chose 21 testing points on a grid with X = 20, 30, 40, 50, 60, 70, 80 and Y = 0, ↘5 ↘10 . We exploited our learned GP model to predict the temperature at these 21 locations over 80 time steps. The prediction uncertainties at all 21 locations were summarized in a box plot. All the box plots over 80 time steps are demonstrated in Fig. 3. Apparently, in the first few steps when the robots did not have much information about the field, the prediction uncertainties are high. However, after about 16 time steps when the robots learned well about the field in both space and time, the uncertainties significantly reduce. Though the prediction uncertainties are considerably small from 20 time steps onwards, there are still some minor variations among them since the temperature kept changing.
This paper has addressed the problem of mapping spatiotemporal environmental field using multiple robots based on Gaussian process regression. The IPP problem has been formulated in terms of multiple prediction steps with crosscorrelation cost functions that guarantee the connectivity of the robot network during the exploration time. By using the dual composition method, we have solved the IPP problem in a distributed manner. The efficacy of the proposed approach was verified in a synthetic experiment utilizing a real-life dataset. In future work, we will consider the synchronous update of measurements in the robot network.