Various processes that exhibit spatial and temporal variations, such as the dispersion of pollutants in the atmosphere or water bodies, are often represented by two-dimensional (2D) PDEs over a domain R . A typical example is the 2D advection-diffusion equation, which models the transfer of substances via advection (the movement of substances through a fluid) and diffusion (the spreading of substances from areas of higher to lower concentration). This can be mathematically expressed as: ∂ z ∂ t r , t = θ ∇ 2 z r , t + v T ∇ z r , t , r ∈ R , (1) where z ( r , t ) denotes the concentration function of the field at position r at time instance t , ∇ and ∇ 2 are the gradient and Laplacian operators, respectively. The coefficient θ > 0 represents the rate of diffusion, and v is the two-dimensional advection coefficient, representing the speed at which a quantity such as heat, concentration, or pollutant is transported by the bulk movement of a fluid. Both θ and v are considered constant but possibly unknown over a fixed interval.

In this work, we consider the field as a linear superposition of multiple non-overlapping advection diffusion phenomena, each governing a spatial-temporal region R i , i = 1 , … , M where M is the number of the advection-diffusion processes. The concentration z_i(r,t) in each region satisfies the advection-diffusion equation: ∂ z i ∂ t ( r , t ) = θ ∇ 2 z i ( r , t ) + v T ∇ z i ( r , t ) ,   r ∈ R i . The global concentration field is then represented as z r , t = ∑ i χ i r z i r , t , r ∈ Ω , (2) where χ_i (r) is an indicator function defined as χ i ( r ) = 1 if r ∈ R i and 0 otherwise, and Ω = ⋃ R i , i = 1 , … , M . Moreover, in various real-world environmental monitoring scenarios, the domain Ω is significantly larger than the dimensions of the operational robots, enabling the approximation of the boundary as essentially flat. Under these conditions, we apply the initial and Dirichlet boundary conditions as shown in Equation 3 at the boundary ∂ Ω (Demetriou et al., 2013): z r , 0 = z 0 r , z r , t = 0 , r ∈ ∂ Ω . (3)

In this work, we consider a group of N mobile sensing robots moving in a coordinated formation in the field Ω represented in Equation 2. The algorithm proposed in this work commands the formation to travel on planned paths to efficiently reconstruct the unknown field. We make the following assumption regarding these mobile sensing robots.

Assumption 2.1

Each sensing robot is equipped with sensors to localize itself and to measure the field concentration value at its current location at each discrete time step k .

The measurement of the i th sensing robot at time step k is modeled as follows: p r i k , k = z r i k , k + n i , i = 1 , … N , (4) where r i k represents the location of the i th robot at the discrete time step k and n i is assumed to be i.i.d Gaussian noise. Additionally, using the locations of all the robots in the formation at time step k , we can determine the location of the formation center r c k at time step k as r c k = 1 N ∑ i N r i k .

When the robots move in a desired formation, it covers a time-varying view-scope Γ ( k ) , which is the area of the field domain Ω that lies inside the polygon formed by sensing robot locations. As shown in Figure 1, the shaded region illustrates the time-varying view-scope Γ ( k ) at discrete time step k , the blue circles represent the four mobile robots, and the red circle represents the formation center. At any given time, the mobile sensing robots can measure and exchange concentration values at their specific locations as shown in Equation 4 and the field values z ( r k , k ) , r k ∈ Γ ( k ) can be estimated by interpolating the measured values from the robots. Consequently, it is reasonable to assume that the estimated field values, z ( r k , k ) , r k ∈ Γ ( k ) , are available to us at all times.

In this work, to facilitate the implementation of the PPO algorithm for source mapping in Section 4.1, we discretize the global field Ω to a E × F grid, where each grid cell represents a single location r k in the map and associates with a concentration value z ( r k , k ) . The following assumption holds for the formation center.

Assumption 2.2

Robots travel in a coordinate formation and the formation center moves along the eight possible directions “up”, “down”, “left”, “right”, “up-left”, “up-right”, “down-left”, “down-right” in the discretized domain.

With the robots moving in a formation, a CKF developed in (You et al., 2016; Wu et al., 2020) is employed to output estimates of concentration z ( r c k , k ) and gradients ∇ z ( r c k , k ) at the formation center r c at time step k . These estimated values will play a major role in the developed algorithm in Section 4. Furthermore, we apply the parameter identification algorithms developed in (Wu et al., 2020) to estimate the unknown diffusion coefficient θ in real-time using the output from the CKF and the RLS algorithm. Thus, in the following discussions, we consider θ as a known value for field reconstruction.

Remark 2.1

Multi-robot formation control is a well-studied topic and researchers have developed numerous formation control algorithms (Zhang and Leonard, 2010; Ren and Beard, 2008; Wu and Zhang, 2012). In this work, we employ the formation control strategy developed in (Zhang and Leonard, 2010) and applied in (Wu et al., 2020). The strategy uses the Jacobi transform to decouple the formation control from the motion control of the multi-robot formation, which enables us to only plan the path and design the controller for the formation center r c . The individual robot controllers are then derived using the formation controller.

In real-world scenarios, the task of mapping complex dynamic fields for cases like gas-leaking and wildfires is important and is often time-critical. It is essential for the robot formation to explore and detect diffusion sources in unknown areas and generate a map as quickly as possible. With the field defined in Section 2.1 and the multi-robot formation defined in Section 2.2, the goal of this study is to design a path for the multi-robot formation so that the formation can identify the multiple non-overlapping diffusion sources in the dynamic field and reconstruct the field in real-time with the limited concentration measurements collected by the multi-robot formation along its trajectory. To achieve the goal, we will introduce a two-mode strategy in Section 4, which consists of a Source Mapping mode and a Field Exploration mode. In the Source Mapping mode, we employ the RL-based algorithm and train a PPO model to guide the multi-robot formation toward a diffusion source in the field and reach a stationary state, where the formation arrives at the source and moves with the field at the same speed as the advection flow. In the Field Exploration mode, we develop a K-means clustering-based exploration strategy to enable efficient exploration of unknown areas.

PPO (Schulman et al., 2017) is a significant development in reinforcement learning, introduced as a more efficient and simpler alternative to Trust Region Policy Optimization (TRPO) (Schulman et al., 2015). PPO is based on the policy gradient approach, a class of algorithms that optimize policies by directly computing gradients of expected rewards for policy parameters. This approach allows the learning agent to improve its policy iteratively by following the gradient of expected rewards. PPO enhances this process by addressing the complexities of earlier methods while retaining their benefits, particularly in maintaining stable and reliable policy updates.

PPO operates as an on-policy method. Unlike traditional policy gradient methods that apply a single update after each interaction with the environment, PPO refines the policy by using multiple updates on the same batch of data. The core of PPO is its surrogate objective function, designed to prevent large, potentially destructive policy updates. This is achieved through a probability ratio r k between the new and old policies, which is clipped to keep updates within a safe range. The surrogate objective function L CLIP ( θ ) is expressed as: L C L I P θ = E ̂ k min r k θ A ^ k , clip r k θ , 1 − ε , 1 + ε A ^ k . (5)

In Equation 5, E ̂ k denotes the expectation over timestep k , θ represents the policy parameters, A ̂ k is an estimate of the advantage function at time step k , and ε is a small hyperparameter that controls the clipping range. It is important to note that while r k and θ follow the conventional notations used in literature, they differ from the notations in other sections of this paper, where r refers to the locations in the field and θ refers to the diffusion coefficient. The clipping mechanism ensures that if the probability ratio deviates outside the predefined range [ 1 − ε , 1 + ε ] , the function applies the clipped values to prevent excessively large updates, thereby maintaining the stability of the learning process. By constraining the probability ratio, PPO effectively controls the size of policy updates, balancing stability and performance. PPO is particularly well-suited for discrete action spaces, which makes it an ideal choice for our environment setup. The PPO algorithm is summarized in Algorithm 1.

Algorithm 1

PPO-Clip Algorithm.

Input: Initial policy parameters  θ 0 , initial value function parameters ϕ 0

Cooperative Kalman Filter (CKF) is a collaborative state estimation scheme first developed in (Zhang and Leonard, 2010), then used in later studies (Wu and Zhang, 2012; You et al., 2016; Wu et al., 2020; You et al., 2022), by combining live sensor data collected by the network of multiple mobile robots to collaboratively improve the accuracy of the state estimation process. In particular, when applied to the state estimation in dynamic fields, the authors incorporated the dynamics of the mobile robot formation and the diffusion equation into the formulation of the state equation of the CKF. This integration facilitates reliable and accurate state estimation, taking into account how changes in diffusion fields and the formation trajectory over time affect sensor data measurements. More specifically, the state vector X ( k ) at each time step k is defined as: X ( k ) = z r c k , k , ∇ z r c k , k , z r c k , k − 1 , ∇ z r c k , k − 1 T , where z r c k , k and z r c k , k − 1 denote the field concentration values at location r c k at two consecutive time steps k and k − 1 , respectively, and ∇ z r c k , k and ∇ z r c k , k − 1 denote the field gradient at location r c k at two consecutive time steps k and k − 1 , respectively. Given that the mobile robots maintain a symmetrical formation while traversing the environment, CKF estimates the state vector X ( k ) along the trajectory of the formation center. Note that since the field is spatial-temporal varying, the field concentration values and gradients are different at time steps k and k − 1 even at the same location r c k . This fact is critical in the construction of the CKF to provide reliable estimates of the state vector. The estimated state vector is subsequently employed to iteratively identify the diffusion coefficients of the field over time, based on the RLS algorithm. The estimated diffusion coefficients are vital for the task of identifying and reconstructing spatial diffusion fields in our paper. To save time and avoid excessive length in this paper, we will not provide the complete derivation of the CKF. For additional details, interested readers may refer to the original papers.

As discussed in the high-level overview, the goal of the formation in Source Mapping mode is to move toward the source of a diffusion field and facilitate diffusion field reconstruction by estimating advection and diffusion parameters. For this purpose, we train a PPO model that takes the field information vector state to predict the optimal action. In this section, we describe the setup of our training environment and the architecture of our PPO model.

We define the observation input state S ( k ) for our PPO model as follows: S k = z r c k , k , ∇ z x r c k , k , ∇ z y r c k , k , (6) where ∇ z x r c k , k and ∇ z y r c k , k represent the estimated concentration gradients at the formation center at time step k , in the x and y directions, respectively. As discussed in the previous section, we rely on the CKF to provide estimates of the field concentration and gradients at the formation center, which form the complete input state of our model. By incorporating concentration gradients in the observation state, we provide the model with the direction of the largest concentration value change at the current location, which can be useful for heading toward the source center. Additionally, our definition of the state vector S ( k ) in Equation 6 limits the model to learning state-action values solely based on the characteristics of the field. With the formation controller, the mobile sensor robots can maintain a constant desired formation while traversing the environment; thus; we only plan the path for the formation center instead of individual robots.

For every time step, our robot formation can move to any adjacent cells (including diagonal) or stay at the current location. Thus, we can define the action space A ( k ) consisting of nine actions as follows: A ⁡ k = “ up ” , “ down ” , “ left ” , “ right ” , “ up - left ” , “ up - right ” , “ down - left ” , “ down - right ” , “ stay ” . (7)

Since our goal is to train a PPO model that can guide the formation toward the center of a diffusion field and maintain stationary state as long as possible, it is crucial to develop a reward function that incentivizes this behavior. For this reason, we model the reward function based on the concentration values inside the formation view-scope as follows: R k = α ∑ r ∈ Γ k z r k , k , (8) where α is a rescale constant. As regions with high concentration values play a significant role in field reconstruction, a reward proportional to the total concentration values within the view-scope as shown in Equation 8 motivates the model to learn to navigate towards areas with higher concentration values. This, in turn, greatly reduces the error in reconstruction and prioritizes information-rich trajectory. We chose PPO as our model due to its long-standing role as a crucial component in various state-of-the-art solutions in reinforcement learning. Furthermore, PPO can be used for environments with discrete action space, as in our case. Figure 3 shows the architecture of our PPO model. Since PPO is a type of actor-critic algorithm, it has two neural networks - the actor network and critic network. In our case, we employed the same architecture for both. This network consists of three dense layers of size 128 each, an input layer of size 3 and an output layer of size 9. We added random dropout layers and regularization between the hidden layers, with ReLu as the non-linear activation function (Agarap, 2018). The network is trained using Adams optimizer (Kingma and Ba, 2015).

Using the trained PPO model, the robot formation takes actions chosen from the action space in Equation 7, and is guided toward the source of the diffusion field until it reaches the stationary state. This stationary state is achieved when the estimated field concentration reaches a local maximum and the estimated field gradient approaches zero, i.e., z r c k , k > z r c k − 1 , k − 1 , z r c k , k > z r c k + 1 , k + 1 , and ∇ z r c k , k ≈ 0 . At this point, the formation moves at the same speed as the advection flow for a designated period of T steps before switching to the Field Exploration mode to search for other diffusion fields in unexplored areas.

Whenever the robot formation reaches a stationary state, indicating the detection of the source of a diffusion field, the robot formation switches back to Field Exploration mode and moves toward unexplored area in the map to look for the remaining diffusion fields. During this phase, it is essential for the robot formation to come up with a new destination that is away from the already explored locations to avoid revisiting the same source, but also not too far to cause the formation to go back and forth when scanning the whole map. In short, our main objective is to generate a path that allows our robot formation to scan the whole field as quickly as possible without leaving any diffusion field undetected. Lawn mowing is a good example of generating a deterministic trajectory for scanning an unknown map. However, since our mobile sensing robots are often deployed in time-critical missions. It is necessary to opt for a more aggressive exploration strategy that allows the formation to discover all sources as quickly as possible. In our algorithm, we partition the unvisited cells in the entire map into multiple clusters using the K-Means clustering algorithm (Na et al., 2010). The selection of K directly affects how aggressive or cautious the scanning behavior of our robot formation will be. The initial value of K is selected based on the initial estimate of the minimum size of a diffusion field. In this work, we use the term “size” to refer to the bounded area of a diffusion field where the concentration value exceeds a certain threshold. For different scenarios (such as wildfires and gas leaks), we are often able to come up with a rough estimation of the average size of a diffusion field. Let us denote this as S ̃ field and the global map size as S map = E × F . Then, K is estimated as: K = max 2 , S map S ̃ field . (9)

After partitioning the map, the robot formation selects the centroid of the nearest cluster as the new destination and moves toward that destination to explore the field. It continues to visit the centroids of other clusters, prioritizing nearby clusters, as long as it is in the Field Exploration mode. The formation switches to Source Mapping mode when the estimated field concentration value exceeds a chosen threshold, i.e., z ( r c k , k ) > δ . At all times, the formation maintains a record of visited clusters and diffusion fields, effectively creating a mask that distinguishes between explored and unexplored regions.

Whenever a switch occurs from Field Exploration mode to Source Mapping mode and the formation reaches the stationary state in the Source Mapping mode (indicating a new diffusion field is detected), the robot formation evaluates the field and computes a new estimated K value to repartition the remaining unexplored areas in the map. The “size” of the newly discovered diffusion field can be roughly estimated based on the circular area with radius R dist extending from the source center to the location where the concentration first surpasses the threshold. This is also where the formation switches from Field Exploration to Source Mapping mode. Denote the size as S new represented as: S new ≈ π ∗ R dist 2 . (10)

With the size of the latest detected diffusion field calculated based on Equation 10, we can update the estimated average size of the diffusion field S ̃ field as follows: S ̃ field ′ = 1 − β ⋅ S ̃ field + β ⋅ S new , (11) where β ∈ [ 0,1 ] is the weighting factor that determines the influence of the newly discovered field on the updated estimate S ̃ field ′ . This approach allows for increasing accurate field size estimation with new discovery. Note that K will only decrease or maintain unchanged with each new source found. When the size of an unexplored area in the map falls below a certain threshold, we consider the whole map has been adequately covered. At this point, The robot formation can either be set to idle mode or be directed to transition to a new map (potentially a neighboring global field) to initiate its operations from a different starting location. The exploration strategy can be summarized in Algorithm 2.

Algorithm 2

K-Means Clustering Based Exploration Mode.

1: Input: Unexplored region M 0 as an N × N grid map.

The field reconstruction begins when the formation detects the source of a diffusion field within the global domain. Given that our environment is modeled as a 2D grid, we discretize the diffusion Equation 1 to enable field reconstruction. Assuming the domain of interest Ω is divided into square cells of size Δ r x = Δ r y , as illustrated in Figure 4, where a 3 × 3 grid is demonstrated. Here, r 1 , r 2 , r 3 , r 4 are the neighboring cells of grid cell r 0 . Let z ( r 0 , k ) , z ( r 1 , k ) , z ( r 2 , k ) , z ( r 3 , k ) and z ( r 4 , k ) denote the concentration values at grid cells r 0 to r 4 at discrete time step k . Using the finite difference method, the discretized advection-diffusion equation can be expressed as: z r 0 , k + 1 − z r 0 , k t s = θ z r 2 , k + z r 4 , k − 2 z r 0 , k Δ r x 2 + z r 1 , k + z r 3 , k − 2 z r 0 , k Δ r y 2 + v T ∇ z r 0 , k + e r 0 , k , (12) where t s denotes the sampling interval and e ( r 0 , k ) represents the approximation error. With the symmetric property, Equation 12 can be further simplified to: z r 0 , k + 1 − z r 0 , k t s = θ ∑ i = 1 4 z r i , k − 4 z r 0 , k Δ r x 2 + v T ∇ z r 0 , k + e r 0 , k . (13)

To reconstruct a diffusion field using the measurements taken by the robot formation with Equation 13, we need the estimated field concentration values z ̂ ( r , k ) within the view-scope of the robot formation at each time step k , as well as the estimated diffusion coefficient θ ̂ and the advection vector v ̂ . The former can be obtained through interpolation at each time step using the measurements taken by the robots. As mentioned previously, we employ the strategy developed in (Wu et al., 2020) to identify the diffusion coefficient θ . In many scenarios, the advection coefficient v is assumed to be a known constant. When the advection coefficient is unknown, we estimate it when the robot formation reaches a stationary state, where the robot’s velocity matches that of the advection flow. Let v = ( v x , v y ) T and k S > 0 represent the time step when a stationary state is detected. Assuming the formation stays at the stationary state for T steps, ( v x , v y ) can be estimated as v ̂ x = r c , x k S + T − r c , x k S T , v ̂ y = r c , y k S + T − r c , y k S T . (14)

With these values determined, the field values across the diffusion field can be propagated through Equation 13. This approach enables field reconstruction using only the sparse measurements gathered along the robots’ paths.

To assess the overall solution, we developed a low-fidelity simulation environment within a discrete space. This environment is structured as a 100 × 100 square grid, incorporating between one to four non-overlapping spatial diffusion fields of varying sizes and configurations, as depicted in Figure 5. We generated a total of 15 different environments, with increasing level of complexity, to investigate and assess the model’s efficiency in mapping unknown environments. The color of each cell in the grid denotes the concentration value of the diffusion field. Each diffusion field possesses distinct and independent advection and diffusion coefficients. Figure 6 shows the evolution of a sample diffusion field over time. The environment was simulated for up to 400 time steps with the discrete interval Δ t = 0.1 , Δ x = 0.1 , and Δ y = 0.1 . Table 1 lists the configurations of the 15 diffusion field environments with advection terms ( v x , v y ) and diffusion term θ . The center of the field is denoted as p o s ( x , y ) . Note that s i z e is only used internally by the generator as a scale factor to control how large a diffusion field appears on the map.

Since the map is discretized in our approach, selecting an appropriate grid cell size also plays a role in both the accuracy of data collection and computational efficiency. Each grid cell should be small enough to capture meaningful concentration gradients but large enough to reduce computational demands. Ideally, the grid cell size should reflect both the overall map dimensions and the characteristics of the environment being monitored. For example, when studying gas leaks, where subtle concentration changes are significant, a finer grid may be required. On the other hand, wildfire propagation fields, which tend to cover larger areas, can accommodate slightly larger cells. Adjusting the grid cell size based on the specific characteristics of the environment allows us to find the right balance between resolution and computational cost, facilitating effective and efficient exploration and reconstruction.

For the training of the PPO model, we follow a curriculum learning approach (Wang et al., 2023; Wang et al., 2022) that involves gradually increasing the complexity of the training environment. We created a 100 × 100 training environment featuring a single diffusion field at the center of the map p = ( 50,50 ) with a constant diffusion coefficient θ = 1 . The spread of this diffusion field at the beginning of an episode is drawn randomly based on the parameter s i z e , which controls how large the area with non-zero concentration is due to the presence of the diffusion field. This environment has two variants: a “static” environment with the advection term set to zero, and a “dynamic” environment with a fixed non-zero advection term.

We simulate a group of four mobile robots in a symmetric formation as shown in Figure 1 to move in the environments, with the formation controller running to maintain the desired formation. With the CKF providing the state S ( k ) which includes the estimated field value z ( r c k , k ) and gradient along the trajectory of the formation center ∇ z x ( r c k , k ) and ∇ z y ( r c k , k ) , the PPO model underwent training on static maps first before advancing to training on dynamic maps. After training, our PPO model outputs the policy to direct the robot formation to move toward the source of a diffusion field and maintain the stationary state, which is required for the estimation of advection coefficients.

We provide the training results in Figure 7. The PPO model was initially trained for 400,000 time steps on the static environment where advection terms are set to zero, as shown in Figure 7A. After that, we proceeded to train the PPO model on the dynamic environment for an additional 2,000,000 time steps, as shown in Figure 7B The model was trained multiple times with orthogonal random weight initialization using the Stable-Baselines3 framework (Raffin et al., 2021), and the best-performing model was selected for use in our experiments. In both environments, the formation is initially placed in a low-concentration region at the start of each episode to avoid starting too close to the source. Additionally, in our dynamic environment, the source is assigned random, non-zero advection and diffusion coefficients, causing it to move in a different direction in each episode. This setup encourages the model to adapt to various scenarios but also introduces some fluctuations in performance early in training, as the formation may take suboptimal actions initially and struggle to catch up to the moving source. In both training phases, however, the average reward per episode increases steadily, indicating that the model successfully improves its given task over time. Additionally, Figure 7C provide some samples of source-heading operation performed by our PPO model post-training. In these samples, the 4-robot formation, shown as the red square in the map, are tasked with locating the source while maintaining its formation. The yellow region denotes the area with high concentration - where the source is located, the red square represents the robot formation, and the red dot is the starting position of the formation center. The robot formation is directed towards the center of the source, which has maximum concentration, per the PPO’s objective of reward maximization. PPO’s role in this task is to update the policy per iteration to make an informed decision on where to go next. The simulation results show the PPO algorithm’s efficacy in source mapping.

The PPO algorithm’s training process effectively learns a policy that guides the robot formation to explore the field, detect diffusion sources, and assist in reconstruction while adapting to dynamic changes in the environment. The algorithm’s capability in handling complex and dynamic environments indicates its potential in real-life scenarios where rapid environmental change happens, and accurate detection of diffusion sources is crucial. This capability is precious in pollution tracking or gas leak detection scenarios, where time-sensitive and precise localization is essential.

When the PPO model leads the robot formation to move toward a diffusion source and the formation center reaches the stationary state, the advection coefficient can be estimated based on Equation 14, and the field reconstruction process can start using Equation 13. Along with the field reconstruction process, the robot formation switches to the field exploration mode, where the K-means clustering-based exploration algorithm 2 plays a role.

Figure 8 provides an example of how the K-means clustering-based exploration works in a 100 × 100 grid map with two diffusion fields. The formation started in Field Exploration mode with the starting position shown as the red dot. Based on our initial assumption about the average size of diffusion fields, the map is partitioned into six clusters with K = 6 as illustrated in Figure 8B. The formation moves toward the centroid of the nearest cluster and continues to other neighboring centroids until an area with a high concentration is detected. When that happens, the formation transitions to Source Mapping mode and attempts to move toward the source center to reconstruct the field. When the new diffusion field is fully reconstructed, the formation switches back to Field Exploration mode, re-calculates a new K as described in Equations 9-11 and re-partitions the unexplored area, as shown in Figure 8C. Note that the unexplored area has excluded visited clusters and any detected diffusion fields. This process continues to repeat until the map is fully covered.

With the trained PPO model and the K-means clustering-based exploration algorithm, we now ready to implement the overall multi-robot Source Mapping and Field Exploration strategy to reconstruct a dynamic field. Figure 9 shows different trajectories of the robot formation obtained from the simulation across multiple spatial-diffusion environments. In this setup, the robot formation started at the same initial position r c 0 = ( 90,10 ) (top-right corner, shown as a red marker). The red square displays the final position of the robot formation at the end of the episode ( k = 300 ) . Despite variations in spatial-diffusion field distributions and characteristics, we can see that the robot formation managed to detect all sources in the map while efficiently covering the entire map before the episode ends.

In Figure 10, we compare the trajectory and mapping errors of our solution against other alternative path-planning algorithms, including Random-Walking (shown in black) and Lawn-Mowing (shown in green). With the Random-Walking strategy, the formation takes a random action in every time-step, resulting in an arbitrary trajectory. On the other hand, with Lawn-Mowing strategy, the formation attempts to scan the field row-by-row until the entire field is fully covered. Compared to these approaches, the trajectory generated by our solution is faster at detecting and mapping diffusion fields, resulting in a much lower mapping error. In this map, our K-Mean-based approach is the only one that detects all diffusion fields before the episode ends ( k = 300 ) .