Distributed Multi-Robot Multi-Target Tracking Using Heterogeneous Limited-Range Sensors

In recent years, exploring robotic networks’ ability to adaptively adjust their sensor coverage, multi-robot multi-target tracking (MR-MTT) problems have been widely investigated. In this paper, we discuss distributed algorithms to find scalable and fault-tolerant approaches to the problems. Distributed MR-MTT mainly consists of cooperative multi-target state estimation and multi-robot active motion planning. The former can be implemented using either Bayesian methods [15, 16, 17, zhu2020fully, 19], which require each robot to recursively estimate the statistical distribution about target states over the entire task space, or consensus algorithms [20, 22, 23], which enable robots to continuously exchange information until they reach an agreement. The latter is frequently formulated as an optimal control problem aiming to maximize the tracking performance, which is defined as the degree of information acquisition [24, 25, 26], the distance between a robot and a target [27], or some general measure of the target observability [28, 19]. However, the above-mentioned schemes could be susceptible to local extrema, where multiple robots track the same target while others are left unobserved, and provide no guarantee to cover the entire task space with the robots’ limited-range sensors. As a result, none of them is able to jointly detect, localize, and track an unknown and time varying number of multiple targets.

Perhaps the most suitable MR-MTT solution to pair space exploration with active target state estimation is through coverage control framework. Coverage control is a popular way to optimally gather target information while guaranteeing complete space coverage. It aims to maximize a performance measure such as the robots’ distances to targets. The most widely used distributed coverage control strategies are Voronoi-based methods [29, 30, 31, 32] that leverage Lloyd’s algorithm [33]. The basic idea is to partition an assigned task space using the Voronoi tessellation and then drive each robot towards the weighted centroid of its assigned Voronoi partition. This process will enable the network to reach a local optimum in its coverage performance.

In coverage control, the mission space is typically convex, and robots are only required to communicate with their neighbors and are assumed to have perfect knowledge of their locations. A number of recent works have provided various solutions to accommodate uncertainty in sensor localization [34, 35, 36, 37] and to allow for non-convex spaces [38, 39, 40, 41]. In particular, leveraging the coverage control framework, the work of Dames [42] proposes to couple a distributed multi-target state estimation scheme and a distributed planning algorithm in order for each robot to estimate the density of targets and move to the centroid of the estimated target density. Consequently, networked robots with a limited-range isotropic sensor are able to actively move to maximize the total detection probability of all targets despite measurement errors and false alarms, while maintaining coverage of the entire task space.

Nevertheless, the above-mentioned work exploits MR-MTT with coverage control framework using fixed types of sensor. Yet, cooperative heterogeneous sensor networks composed of arbitrary types of sensors expand the use of robotic networks to complex scenarios where a variety of different types of robots and sensors are needed to complete the task. A variety of approaches have been proposed to incorporate measurements from multiple heterogeneous resources to improve the target tracking quality [43, 44, 45] or to enhance the robustness and resilience of a robotic network [46, 47]. However, none of these methods aims to attain the globally optimal tracking performance by taking advantage of heterogeneous sensing capability while maintaining coverage of the task space.

Weighted Voronoi diagrams, which include power diagrams, have recently been used to account for heterogeneity in the radii of circular sensor footprints [41, 48, 49, 50] and energy levels [51]. Voronoi-based methods have also been used for wedge-shaped sensor field of views (FoV) [52]. Stergiopoulos and Tzes [53, 54, 55] used novel partitioning and distribution algorithms for sensing networks where all sensor FoVs share the same arbitrary shape but may be of a different size. Hexsel et al. [56] introduced a gradient ascent-based coverage control algorithm to maximize the joint probability of detection in anisotropic sensing networks. Nevertheless, as these approaches use gradient ascent, robots can prematurely converge to local extrema of an objective function they aim to maximize.

Building on our previous conference version [57], this paper proposes a novel multi-robot multi-target joint state estimation and planning approach for heterogeneous limited-FoV robotic networks. The proposed method allows a team of robots to search for targets over a task space and actively maintain coverage of a majority of detected targets in a distributed manner. We summarize our main contributions as follows.

1. 1.

   We propose a new measure called the normalized unused sensing capacity that quantifies the difference between the current information that a sensor gathers and the theoretical maximum. This can be computed using entirely local information and does not require any assumptions about the type of sensor or the shape of the sensor FoV.

2. 2.

   Leveraging this measure, we first replace the standard Voronoi diagram in Lloyd’s algorithm with the power diagram, with the goal of balancing the task load across the team. In particular, we assign robots that have more accurate sensors, larger sensor FoVs, and/or are not currently tracking any objects to cover larger areas.

3. 3.

   While power diagram implementation provides fast and near optimal space allocation for heterogeneous robotic teams, it does not yield the best tracking accuracy. Therefore, we propose CCVD, a closed-form optimal space partitioning algorithm to further improve the tracking performance and to evaluate the power diagram method by comparing it with theoretical optimum.

4. 4.

   We demonstrate the efficacy of our approach through a series of experiments in simulated environments to show that the approach yields higher quality and more reliable tracking than the standard Voronoi diagram-based approach and the zigzag coverage path, especially when the robots are highly heterogeneous. To the best of our knowledge, our work presents the first distributed algorithm that allows heterogeneous sensors to jointly detect, localize and track unknown and time-varying number of targets. Therefore, our results are only comparable to baseline algorithms [42][torres2016coverage].

Distinct from our conference paper [57], this paper includes the following novel contributions.

1. 1.

   Whereas only the power diagram was used in [57], in this paper, we introduce a new space partitioning algorithm based on CCVD. The area of the task space allocated to each robot by CCVD is strictly proportional to the sensing capacity of that robot, which will further optimize the space partitioning by leveraging the heterogeneity of the robots, but at the cost of greater computing resource requirements.

2. 2.

   We analyze the computational complexity of the proposed algorithms, and discuss the communication cost incurred by the algorithms.

3. 3.

   We present ROS simulations with visualization in Gazebo and RViz. We conduct additional Matlab simulations for quantitative analysis of the new methods.

4. 4.

   We compare our methods with a coverage path planning approach quantitatively. Meanwhile, we conduct an ablation study to demonstrate the usage of both the centroid of detection and the partitioning algorithms.

The rest of this paper is organized as follows. In Section II, we formulate the MR-MTT problem and introduce the prerequisites of our proposed methods. Section III presents our proposed space partitioning algorithms using both the power diagram and CCVD along with a distributed control algorithm based on these partitions. In Section IV we analyze the performance of our algorithms in terms of computational complexity and communication loac and introduce the quantitative metrics used to validate our algorithms. Finally, we show simulation results in Section V and draw conclusions in Section VI.

Given a density function $\phi(x)$ defined over the 2-dimensional Euclidean space, the objective of Lloyd’s algorithm is to compute $Q={q_{1},\cdots,q_{n}}$ and $\mathcal{W}=(\mathcal{W}_{1},\cdots,\mathcal{W}_{n})$ minimizing the following functional:

where $\mathcal{W}_{i}$ is the dominance region of robot $s_{i}$ (i.e., the region that is assigned to robot $s_{i}$ for coverage), $\|\cdot\|$ is the Euclidean distance, $x\in E$, and $f$ is a monotonically increasing function. The function $f$ defines how the cost of robots’ sensing depends on its distance to each sensing location $x$. The dominance regions $\mathcal{W}$ partition $E$, meaning the regions have disjoint interiors (i.e., $\mathrm{int}(\mathcal{W}_{i})\cap\mathrm{int}(\mathcal{W}_{j})=\emptyset\;\forall i\neq j$) and the union of all regions is $E$ (i.e., $\cup_{i=1}^{n}\mathcal{W}_{i}=E$) [29].

When $f(x)=x^{2}$, each $\mathcal{W}_{i}$ of $\mathcal{W}$ is given by $\mathcal{W}_{i}=\{x\mid\|x-q_{i}\|=\min_{k=1,\ldots,n}\|x-q_{k}\|\}$ [29]. In other words, $\mathcal{W}_{i}$ is the collection of all points that are the nearest neighbors of $s_{i}$. This is widely known as the Voronoi tessellation, which is illustrated in Figure 1. We refer to $\mathcal{W}_{i}$ as a Voronoi partition, which is convex by construction, and $q_{i}$ as the generation point of $\mathcal{W}_{i}$. When each $q_{i}$ is the weighted centroid of the $i$-th Voronoi partition given by

the robot positions $Q$ minimize $\mathcal{H}$ [29].

By applying Lloyd’s algorithm, coverage control sets the control input for robot $s_{i}$ to

where $k_{\rm prop}>0$ is a positive gain. In this paper, we use the control law in a discrete time manner, following this direction at the maximum speed of the robot. The angular velocity uses a bang-bang strategy, maximizing the angular velocity until the robot is facing directly towards the goal $q_{i}^{\ast}$. By following this control law, the robots asymptotically converge to the weighted centroids of their associated Voronoi partitions. This still holds even when $\phi$ varies with time.

In an MTT setting, a natural choice for $\phi(x)$ is to capture the target density at each location $x$. This time-varying density function can be estimated using any standard MTT algorithms. In this paper, we use the PHD filter, as it does not require any explicit data association.²²2We note that the overall framework presented in this paper can be easily adapted to any other choice of Bayesian MTT tracker.

The Probability Hypothesis Density (PHD), denoted by $v(x)$, is the first order moment of a random finite set (RFS), i.e., a set with a random number of random elements such as a time-varying set of moving targets, and takes the form of a density function over the state space. The integral of the PHD over a region is equal to the expected number of targets in the region at a given time. By assuming that the RFSs are Poisson, i.e., the number of targets follows a Poisson distribution and the spatial distribution of targets is i.i.d., the PHD filter recursively propagates the PHD in order to track the distribution over target sets [11]. Similar to any other Bayesian state estimators, the PHD filter recursively predicts the target state using state transition probability (i.e., target motion), and updates the prediction using the probabilistic models of the sensor measurements.

The PHD filter uses three models to describe the motion of targets: 1) The motion model, $f(x\mid\xi)$, describes the likelihood of an individual target transitioning from an initial state $\xi$ to a new state $x$. 2) The survival probability model, $p_{s}(x)$, describes the likelihood that a target with state $x$ will continue to exist from one time step to the next. 3) The birth PHD, $b(x)$, encodes both the number and locations of the new targets that may appear in the environment.

The PHD filter uses another three models to describe the ability of robots to detect targets: 1) The detection model, $p_{d}(x\mid q)$, gives the probability of a robot with state $q$ successfully detecting a target with state $x$. For instance, we can set the detection probability to zero for all $x$ outside the sensor FoV. 2) The measurement model, $g(z\mid x,\,q)$, gives the likelihood of a robot with state $q$ receiving a measurement $z$ from a target with state $x$. 3) The false positive (or clutter) PHD, $c(z\mid q)$, describes both the number and locations of the clutter measurements.

Using these target and sensor models, the PHD filter prediction (5a) and update (5b) equations are as follows:

where $\psi_{z,q}(x)$ is the probability of a sensor at $q$ receiving measurement $z$ from a target with state $x$ and (5c) is used within the prediction step (5b). The subscript $t$ refers to time. Note that $v(x)$ is not a probability density function (PDF), but can converted into one by normalizing it with the expected number of targets in a sub-space of $E$,

Dames [42] developed a distributed PHD filter in which each robot maintains the PHD within a unique subset, $\mathcal{W}_{i}$, of the environment while ensuring the distributed filtering scheme yields the same target estimation performance as its centralized counterpart. Three algorithms then account for motion of the robots (to update the subsets $\mathcal{W}_{i}$), motion of the targets (in (5a)), and measurement updates (in (5b)). In this paper, we adopt the same strategy for each robot to locally maintain its portion of the PHD.

We assume that each robot $s_{i}$ has a finite field of view (FoV) $F_{i}$, which could be different across the sensors $i$. Examples of $F_{i}$ include a wedge shape for a camera (e.g., Figure 2) or a circle for a lidar. Let $p_{d}(x|q_{i},\theta_{i})$ denote the probability of robot $s_{i}$ at position $q_{i}$ and with orientation $\theta_{i}$ detecting a target with state $x\in E$, which is the same model as in the PHD filter. The robot cannot detect targets outside its FoV $F_{i}$, i.e., $p_{d}(x|q_{i},\theta_{i})=0\,\forall x\notin F_{i}$. We assume that $p_{d}$ is time-invariant in the robot’s local coordinate frame, but we make no other assumptions about the shape of $F_{i}$ or the functional form of $p_{d}$. In practice, $p_{d}$ can be estimated using data-driven approaches [59] and/or prior knowledge of the sensor model.

We define the total detecting capability of robot $s_{i}$ as

which depends on the size of $F_{i}$ and the sensor’s ability to detect information within $F_{i}$.

In this paper, we consider sensor heterogeneity in terms of not only the sensing capability, such as FoV and detection accuracy, but also its current usage to track targets. When a sensor detects a target and starts to track it, the sensor’s sensing capability will be reduced. We assume that we are tracking targets of finite size that cannot overlap, and that the target size is much smaller than the sensor’s FoV. Therefore, for a robot $i$, there exists some area $\mathcal{B}$ such that only 1 object can be in that area, i.e., $\max\int_{\mathcal{B}}v(x)\,dx=1$. Then

where the last approximate equality holds for small $\mathcal{B}$, such that $p_{d}$ is approximately constant over $\mathcal{B}$. We can then write robot $s_{i}$’s FoV as the union of disjoint regions, $\mathcal{B}_{k}$, with the above property, so that $F_{i}=\cup_{k}\mathcal{B}_{k}$ and $\mathcal{B}_{i}\cap\mathcal{B}_{j}=\emptyset,\,\forall i\neq j$. Letting $p_{d,k}\approxeq\max_{x\in\mathcal{B}_{k}}p_{d}(x|q_{i},\theta_{i})$, we then define robot $s_{i}$’s maximum sensing capacity as

where $D_{i}$ comes from (7) and $\mu$ is a tuning parameter associated with the maximum target density in the task space. For instance, as the expected maximum distance between the targets becomes larger, a smaller $\mu$ is selected to discount the maximum sensing capacity of a sensor. Remark 1 discusses the selection of $\mu$ in more detail.

The expected number of target detections at $x$ is given by

where $p_{t}(x)$ is defined in (6). Thus, the total expected target detection probability is given by

Note that (12) can be easily modified by replacing $v(x)$ with a density function propagated via a different Bayesian filter. The normalized unused sensing capacity quantifies current capacity for a sensor to track targets. The larger the sensing capability of a robot, the higher the number of targets it can track. However, as the robot tracks an increasing number of targets, its remaining sensing capability begins to decay.

To maximize the total detection probability of targets, we control the robots, considering their heterogeneity in spatial deployment. To this end, we optimize both space partitioning and each sensor’s location within its assigned partition. Unlike the Voronoi diagram, which is suitable for the sensors with a homogeneous and isotropic sensing model, power diagrams [60] are often used to compute the optimal dominance regions when the team has heterogeneous sensing models.

The power diagram is a variant of the standard Voronoi diagram that uses the power distance,

where $\rho(s_{i})$ is the weight or power radius of $s_{i}$, and $p_{i}$ is the generation point. Figure 1 illustrates an example of the power diagram. Existing power diagram-based approaches utilized sensor positions as the generation points, $p_{i}$, and the radii of the sensor FoVs as the weights, $\rho(s_{i})$, to account for sensor heterogeneity [50, 41]. However, these approaches are limited to isotropic sensors. On the other hand, our approach extends to heterogeneous anisotropic sensors.

We utilize the normalized unused sensing capacity $U_{i}$ to set the power radius in (13). This is a novel strategy to account for the heterogeneity in computing the dominance regions. To achieve this, we proceed by expressing the optimization functional (2) as

where $g:\mathbb{R}\rightarrow\mathbb{R}$ is a mapping from the unused sensing capacity to the power radius. Since the normalized unused sensing capacity has units of area, we choose $g$ such that the resulting power radius, $g(U_{i})$, is equal to the radius of a perfect (i.e., $p_{d}=1$) isotropic sensor with the same total sensing capability, $D_{i}$, i.e., $\pi g(U_{i})^{2}=U_{i}$. Therefore we have

Existing methods based on the power diagram, such as [50, 41, 51], assume that the sensor’s detection probability at $x$ in its assigned power partition is a non-increasing function of the distance from $x$ to the sensor. In other words, the location of a sensor is the location that maximizes its detection probability of targets in its power partition. Thus, it makes sense that they use the sensor location as the generation point, i.e., let $p_{i}=q_{i}$. However, this no longer holds true for anisotropic sensors.

Instead, we find the weighted centroid of the detection probability as

which we call the centroid of detection (COD). We use COD as the generation points for our power diagram, i.e., $p_{i}=q_{\textrm{cod},i}$.³³3Note that for an isotropic sensor the COD will be the same as the sensor position. Thus, the power partition of each robot becomes

The capacity-constraint Voronoi diagram (CCVD) [61], visualized in Figure 1, computes the optimal task space assignment with the weight of each generation point as a hard constraint and yields closed-form optimal space partition in a discrete space. This is done by two steps: initial cell assignment and cell swapping. 1) Firstly, the task space is segmented into a finite set of regular grid cells $X=\{x_{1},\cdots,x_{|X|}\}$, where each cell is identified by its center $x_{i}$. The $i$th cell is indexed by $X[i]$. Then, we conduct an initial cell assignment satisfying a capacity constraint which specifies the maximum number of cells that can be assigned to each generation point. 2) Secondly, we iteratively revise the cell assignment to minimize the total cost defined by

where $A(x)$ is the generation point assigned to $x\in E$, $f(\cdot)$ is a monotonically increasing function as introduced in (2), and $\rho(\cdot)$ denotes the weight of a generation point as introduced in (13). The right-most term in (18) is a constant for all assignments, and can therefore be omitted. The cell assignment minimizing (18) results in a discrete power diagram where the number of cells in each power partition is equivalent to the capacity of its generation point.

Once the optimized partition is retrieved, each robot must move to its optimized pose to improve the detection probability. By adopting the same analytical approach presented in [29], we can compute the partial derivative of $\mathcal{H}_{p}(Q,\mathcal{W})$ with respect to $q_{\textrm{cod},i}$ as

where $M_{\mathcal{W}_{i}}$ and $C_{\mathcal{W}_{i}}$ are the mass and center of mass associated with the region $\mathcal{W}_{i}$, respectively, defined as

Thus, as $q_{\textrm{cod},i}$ is recursively driven to $C_{\mathcal{W}_{i}}$, the partial derivative approaches to $0$ and thus sensing capability of $s_{i}$ in $\mathcal{W}_{i}$ is optimized. The control inputs in (1) for $s_{i}$ are then given by

where $\textrm{ang}(\cdot)$ denotes the angle of a position vector in global frame.

Algorithm 3 outlines the distributed control algorithm for each robot. At each time step, the robot first finds the optimized partition using one of the approaches explained in Sections III-B and III-C, and maintains the PHD within its partition using the same strategy as in [42], which requires the exchange of local PHD with its neighbors. Then, the robot computes its control using (24).

Our distributed algorithms for the space partitioning using a power diagram, as described in Sections III-B and III-C, are similar to distributed computation methods for constructing a Voronoi diagram. Resorting to the computational complexity analysis reported in [64], the complexity of such methods is $\mathcal{O}{\left(\log{|\mathcal{N}_{i}|}\right)}$ and, hence, so is ours, where $|\mathcal{N}_{i}|$ is the number of neighbors of each robot $s_{i}$.

In the CCVD implementation described in Section  III-C, Algorithm 2 requires each of the $n$ robots to find the optimal assignment of the total $|X|$ cells. For the worst-case complexity analysis, we assume that each robot is a neighbor of all the other robots, resulting in $n-1$ execution of the loop in Lines 3-21. The most time-consuming part of the algorithm is in finding the cell $\tau$, i.e., the median point in the unordered array $P_{ij}$. In our implementation, we employ Hoare’s quick-select algorithm [65] to find the median value $\tau$ by choosing a pivot randomly and dividing the array into two parts, i.e., one with values less than or equal to the pivot and the other with the values greater than the pivot. The median value is either included in the former part if this array has $(n+1)/2$ elements or more, where $n$ is the number of elements in the array, or the latter part otherwise. This divide-and-conquer algorithm allows the algorithm to run in linear time on average, with the time complexity of $\mathcal{O}{\left(\frac{|X|}{n}\right)}$ [63]. Hence, the total time complexity for each single neighbor iteration is $\mathcal{O}{\left(\frac{|X|}{n}\right)}$. Considering that the robot has $n-1$ neighbors by assumption, the complexity for a single robot iteration is $\mathcal{O}{\left(|X|\right)}$.

Therefore, the distributed construction of power diagrams has lower computational complexity compared to the construction of CCVDs. Additionally, the complexity of constructing power diagrams depends only on the number of neighbors, making it a better choice when the computational power of the robots is limited, despite the fact that it provides only a near-optimal solution for space partitioning. One may choose to use the most suitable algorithm depending on the available computational resources and trade-off between computational time and tracking accuracy.

We assume that the communication between each robot and its neighbors is lossless and delay-free. The bandwidth requirement for the robot network is mainly determined by the number of neighbors each robot has and both the size and the frequency of data exchange. On the one hand, given a fixed data size and communication frequency, the worst-case bandwidth requirement occur when all robots are neighbors of one another , in which case $n\left(n-1\right)$ non-overlapping communication links are required, where $n$ is the number of robots. In fact, this worst case assumption is unlikely to happen as the robots are expected to spread out over a large task space.

On the other hand, given a fixed number of communication channels, each robot only needs to exchange information about its space assignment (either power partition vertices or CCVD cell sets), location, and local PHD to its neighbors. The Voronoi methods, at least with constant weighting function, typically result in each robot having 6 neighbors in a hexagonal packing structure [66]. Hence, the bandwidth requirement is low, and will not increase as the environmental dimension (number of robots and targets, size of the task space, etc) increases.

To analyze the communication load of the space partitioning algorithms, we consider the data size of the information exchanged between robots. For the power diagram, the data size is determined by the number of vertices in the power partition, which is proportional to the number of cells in the task space. Each vertex is represented by a 2D position vector, which requires 8 bytes per vertex for a double-precision floating point number. For the CCVD, each robot shares its current cell assignments, capacity, and local information about detected targets. The data size is proportional to the number of cells in the task space and the number of detected targets. Each cell is represented by a 2D position vector, which requires 8 bytes per cell. Each detected target is represented by a 2D position vector and a detection probability, which requires 12 bytes per target. Hence, for an average of 6 neighbors, 10 targets per robot, and a frequency of 1 Hz, the bandwidth requirement for the power diagram is $8\times 6\times 1=48$ bytes per second, and for the CCVD is $\left(8+12\times 10\right)\times 6\times 1=768$ bytes per second.

Currently, robots share all information with their neighbors in every iteration. However, the communication frequency and bandwidth requirements can be further reduced through selective sharing approaches, where re-partitioning only happens when robot positions change significantly.

First, we apply CCVD implementation and show the result from a single run using 5 TurtleBot3 Burger differential drive robots, namely $s_{1},\ldots,s_{5}$, tracking 40 holonomic moving targets in a $10\,\mathrm{m}\times 10\,\mathrm{m}$ obstacle-free square task space as visualized in Figure 3.

The simulation was implemented using Ubuntu 18.04 with ROS Melodic. The robots move at a maximum linear velocity of $0.2\,\mathrm{m/s}$ and a maximum angular velocity of $1\,\mathrm{rad/s}$, and they are able to localize themselves using position data retrieved from ROS topic /odom , i.e., the ground truth state from Gazebo. Our MR-MTT algorithm recursively generates next waypoints for the robots, and they navigate through the waypoints using the move_base ROS package which uses the dynamic window approach (DWA) as a local planner and Dijkstra’s algorithm as a global planner [69].

All robots are equipped with heterogeneous range-bearing sensors with specifications described in Table I. The standard deviation of range and bearing measurements for all sensors is 0.04 m and $0.1^{\circ}$, respectively. Since target density is low relative to the size of $E$ and the FoV of robots, we select $\mu=0.1$ in (9) for all robots. The task space is discretized into $50\times 50$ cell for both PHD representation and task space assignment. We assume that only one target can occupy each $0.2\,\mathrm{m}\times 0.2\,\mathrm{m}$ cell.

As explained in Section I, an important application of the MR-MTT is in tracking targets to study collective behaviours of animal groups and crowd dynamics of pedestrians. Motivated by related works in modeling their collective motion [70], we apply Boids algorithm, which has been used to study the emergent flocking behaviors arising in social animal groups from combinations of separation, alignment, and cohesion behaviors, to simulate the moving targets.

In our simulations, a target may exit the task space in which case another target will immediately enter into the space to maintain a constant number of targets. Targets move at a maximum speed of $0.2\,\mathrm{m/s}$.

The robots begin with a uniform PHD where we set the expected number of targets equal to $1$ over the entire task space, meaning that they have no prior knowledge of the target density distribution. At the beginning of the simulations, the robots are located alongside the lower edge of the task space where they cannot detect the majority of the targets. Four animation clips excerpted from a 12 minute 30 second long target tracking simulation are presented in Figure 4. Each column corresponds to one of the clips illustrating the following three behaviors exhibited by the robot network.

To quantitatively compare the efficacy of our power diagram (Section III-B) and CCVD (Section III-C) based method with baseline algorithms, we run batches of simulation trials in MATLAB using a point robot model with dynamic model given as (1). In order to demonstrate the effectiveness of two novel algorithms and to perform ablation experiments, we apply five methods in the following tests:

1. 1.

   Zigzag (Z) Method: This method utilizes a complete coverage path planning framework, a standard approach for target search. Firstly, a centroidal Voronoi diagram of $n$ Voronoi sites is computed using Lloyd’s algorithm. Each robot is assigned a Voronoi cell of the diagram. Then, a zigzag coverage path is planned [torres2016coverage] with the spacing of two adjacent parallel paths equaling to 1 m for each robot to exploit its cell. An example of paths planned for 60 robots is shown in Figure 9. Each robot tracks targets using PHD filter while traversing through the pre-planned coverage path.

2. 2.

   Voronoi (V) Method: Partitioning the space using Voronoi diagram with robot locations as generator points, similar to Dames’ method [42].

3. 3.

   Voronoi-COD (VC) Method: Similar to Voronoi method besides using robot CODs as generator points.

4. 4.

   Power-COD (PC) Method: Our first proposed method, i.e., partitioning the space using power diagram with robot CODs as generator points.

5. 5.

   CCVD-COD (CC) Method: Our second proposed method, i.e., partitioning the space using CCVD with robot CODs as generator points.

In these simulations, the size of the environment $E$ is $100\,\mathrm{m}\times 100\,\mathrm{m}$. Robots move at a maximum linear velocity of $2\,\mathrm{m/s}$ and a maximum angular velocity of $2\,\mathrm{rad/s}$. Each robot carries one of the five sensor types, as outlined in Table II. Assuming that the target density within the FoVs of sensors are completely unknown, we select $\mu=1$ in (9) for all robots in all trials. Based on these type definitions, we define six different network compositions, shown in Table III in which columns A-E list the number of robots with each sensor type. Distinguished from Section V-A, in this series of simulations, target motion has a higher degree of randomness. Targets move at a maximum velocity of $1\,\mathrm{m/s}$ with their heading directions randomly changing over time, and may enter or leave the environment by crossing over the boundary. This means that the number of targets varies over time, and the true number of targets is not known to the sensors.