CVT-based swarming techniques for convex and non-convex cases

Abstract

In the multi-robot swarm control algorithm, the traditional rigid formation algorithm has low formation redundancy, and it is difficult to flexibly adapt to the external environment. In this study, centroidal Voronoi tessellation (CVT) was used to control a group of robots for forming the required flexible formations and increasing the structural redundancy and internal carrying capacity of the multi-robot system. When the CVT area is divided in a non-convex environment, the Voronoi and non-convex areas will overlap, which will cause the robot to collide with the obstacles and hinder the optimization of the cost function. In this study, constrained CVT was used to solve the problem of region overlap in non-convex environments. The proposed path planning algorithm ensured the safe movement of the robot and explored the applicability of the algorithm in different scenarios. Finally, the performance of the algorithm in the non-convex area was verified and discussed by analyzing the simulation and experimental results for the convex and non-convex environments. The results showed that the constrained CVT could control the robot to generate the formation as quickly as possible in the non-convex area, and it was easy to dynamically switch the formation according to the external environment.

Keywords

Multi-robot swarm centroidal Voronoi tessellation (CVT)convex non-convex

Introduction

Swarm robotics has attracted wide attention in the field of intelligent control and has replaced humans in natural resource prospecting, forest fire monitoring,¹ spacecraft exploration,² multi-robot mapping, disaster relief, and other areas. Multi-robot formation control is an essential branch of multi-intelligent systems wherein multiple robots are required to maintain a particular formation and avoid collision during a mission. Over the years, extensive research results have been reported on multi-robot formation control and various formation control methods have been developed based on the leader–follower method,³ behavior-based control method,⁴ virtual structure method,⁵ artificial potential field method,⁶ and self-organizing coverage method based on gradient and grouping.⁷ These formation methods have different advantages and disadvantages and are often used in combination to achieve the desired formation effect in various application scenarios. The advantage of the leader-following method is that the formation structure is simple and easy to implement, but the leader's mistake will be propagated through the entire system. The behavior- based method decomposes the robot motion into basic behaviors, but it cannot clearly define the overall behavior of the formation system, which is not conducive to the stability analysis of the system. The advantage of an artificial potential field is that it can enable the robot to avoid obstacles, but it can easily fall into local optimization. Nevertheless, it can reduce the system cost, improve the system's robustness and efficiency, and provide redundancy, reconfigurability, and structural decomposition for the system simultaneously.^8,9

Formation control mostly aims for a fixed formation. All robots would implement their mission well even if they do not have to maintain the desired distance between them, provided they do not collide with each other. The centroidal Voronoi tessellation (CVT) algorithm used in this study ensures collision avoidance and facilitates the proper formation of multiple robots.¹⁰ The CVT algorithm implements the distributed calculation to optimize the cost function by using a Voronoi partition.¹¹ Further, according to the properties of the Voronoi diagram, this method can ensure that the agents move only within their own Voronoi areas to avoid collisions. A careful review of the literature shows that most researchers have used this algorithm for dealing with robot swarm control in convex areas.^12,13,14 Corte et al.¹⁵ proposed a Voronoi partition and Lloyd's descent algorithms. It was assumed that the sensor moves according to the first-order dynamics. They proved that the sensor position converges to the Voronoi configuration of the center of mass, which is also the local optimal coverage for the cost function related to the sensing performance.

However, their algorithm cannot directly adapt to a non- convex environment with obstacles owing to the following two reasons. First, the sensing ability may be attenuated because of obstacles, thus invalidating the sensing model based on the Euclidean distance. Second, when the center of mass is outside the area, it cannot be guaranteed that the center of mass of the Voronoi area is set as the desired sensor target position.

In the case of robots with relatively simple functions in the swarm robot system, the limitations of computing power and communication load prevent them from adapting to a larger-scale complex task environment. Therefore, complex environments mainly have the following two challenges: how can swarm robots create a formation more flexibly and improve the internal robustness and flexibility to adapt to large-scale expansions? How can swarm robots adapt to a variety of non-convex environments of greater complexity to ensure the safety of the robots during their movements while maintaining their performance.

In response to the above challenges, this paper proposes a swarm control algorithm for a convex environment, and swarm control and path planning algorithms for a non-convex environment. It is ensured that the swarm robot system can adapt to a more complex non-convex environment while flexibly forming a swarm robot system. Simultaneously, it can optimize the cost function without reducing the performance of the swarm, and guarantee the safe movement and expansion ability of the swarm robots. The CVT-based swarm robot control can be applied to more practical scenes. These scenes are non-convex, it contains randomly distributed obstacles and holes.

The organization of this article is as follows. Section “State of the art” presents the state of the art. Section “CVT-based swarm control algorithm” presents the swarm control algorithm based on CVT. It then explains the difficulty in CVT tessellation in a non-convex environment, and describes the constrained CVT and path planning algorithm developed in this study. Subsequently, it presents the control law and related proofs. Section “Simulation of convex and non-convex areas” describes and analyzes the simulation experiments and results for different scenarios. Section “Swarm robotics experiment” presents the experimental results of this study. Finally, Section “Conclusion and future work” presents the conclusion and discusses the scope of future work.

State of the art

Most groundbreaking studies have considered only barrier-free environments. As is well known, the existing results for such an ideal environments are difficult to generalize to more complex obstacle areas. To overcome these limitations for a non-convex area, Boardman et al.,¹⁶ Prabowo1 et al.,¹⁷ and Collins et al.¹⁸ solved the global coverage control problem of multirobots in a non-convex environment by combining the Voronoi tessellation and corresponding path planning algorithms. Among them, Prabowo et al. used the greedy best-first search algorithm to calculate the best center of mass outside the obstacle according to whether the center of mass is inside the obstacle to increase the coverage time performance. Although these methods are advantageous in terms of scalability and calculation time, they cannot guarantee the optimal trajectory of the robot, and the cost function converges slowly.

Stergiopoulos et al.¹⁹ and Nair et al.²⁰ used the geodesic distance instead of the traditional Euclidean distance, allowing multiple robots to be deployed in a non-convex environment. However, it is difficult to avoid collisions when the center of mass and current position of the robot pass over obstacles. Finally, owing to the high cost of geodetic distance calculation, it is easy to cause the robot to fall into the extreme point owing to obstacles.

Aiming at the robot collision problem that may occur under the geodesic distance, Breitenmoser et al.²¹ introduced the geodesic distance and combined the classic Lloyd algorithm with a local path planning algorithm, TangentBug, to calculate the robot's movement near obstacles and corners. The unreachable centroid is mapped to the reachable area, which proves the convergence and optimality of the algorithm. Palacios-Gasós et al.²² realized the Voronoi tessellation in a non-convex area by utilizing the geodesic distance combined with the weighted Voronoi diagram. In addition, the coverage area scan of the non- convex area was realized through regional gridding and an online path planning algorithm.

Caicedo-Nunez et al.²³ and Xu et al.²⁴ mapped the connected regions in a class of non-convex environments to the convex regions through diffeomorphism and used the Lloyd algorithm to effectively solve the coverage problem in simple non-convex environments. However, this method is not effective in coverage problems under multiple obstacles or more complex non-convex environments.

Howard et al. et al.²⁵ proposed a combination of Voronoi tessellation and the potential field method. The robot is used as a virtual particle subject to make the robot and the obstacle repel each other so that multiple robots can cover the area to the greatest extent. On this basis, Renzaglia et al.²⁶ made the robot repel the obstacles while being attracted to the Voronoi centroid, thus completing the coverage of the robot in a non-convex environment.

Lu et al.²⁷ and Mahboubi et al.²⁸ defined a visibility- based Voronoi diagram based on the attenuation of sensor capabilities caused by boundary walls or obstacles. Among them, Lu et al. used the centroid projection method to solve the problem of the centroid lying outside the Voronoi region, and introduced smoothness constraints during the deployment process to avoid the phenomenon of sharp turns of the robot.

In addition, unmanned aerial vehicles (UAVs) are important in the research field of optimal coverage control. Some core features and control methods of UAVs and the related future work have been presented by Tahir et al.²⁹ Du et al.³⁰ reviewed the scheduling problem and expounded on the intelligent algorithm of the UAV swarm. They also discussed the constraints in various tasks. Chen et al.³¹ applied swarm drones to agricultural applications such as spraying of pesticides. CVT was used to identify the best position of the drone during the spraying process, and its combination with the proposed intelligent health balance system extended the overall service life of the drone. In another research, Chen et al.³² eliminated the toxic substance diffusion process through a set of statically connected grid sensor networks, which was modeled as a distributed parameter system controlled by parabolic partial differential equations. Govindaraju et al.³³ used drone swarming to monitor the forest and proposed a probabilistic sensor model to alleviate visibility degradation caused by vegetation. CVT was used to calculate the way points, and the linear combination of the terrain contour and vegetation cover was used as the density distribution function. By compromising on the average visibility, the strategy of pushing the observation points to dense vegetation areas can provide near-uniform visibility across the area. The strategy of pushing the observation points toward the local peak can provide higher average visibility across the entire area. CVT solved the problem of optimal planning of the movement trajectory of multiple UAVs in the air. However, the research was limited to the convex environment. As the areas in most practical application are non-convex, this largely limits the application of CVT-based swarm robots.

CVT-based swarm control algorithm

CVT swarm control

In the deffnition of the Voronoi tessellation algorithm, a region containing n generated points is partitioned into n cells. Each cell contains a generated point, and the other points in the cell are close to the generated point of that cell, thus ensuring that the cost function of the system reaches the minimum value.³⁴ CVT is a particular type of Voronoi tessellation, which sets the density function of the target region, and its generated point is the centroid of each Voronoi cell.

Suppose Q represents a convex polytope. A single robot in a multi-robot system has the position $p_{i} \in Q, i = 1, 2, \dots n$ , where n robots move within the bounded interval $Q \in S^{2}$ . $S^{2}$ is a two-dimensional grid map of the target area. The position of the wireless sensor is represented as $q_{j}, q_{j} \in S^{2}$ . Then, the Voronoi cell $V = \sum_{i = 1}^{n} v_{i}$ in the region is Q given as follows:

v_{i} = {q \in Q | ‖ q - p_{i} ‖ \leq ‖ q - p_{j} ‖, i, j = 1, 2, 3, \dots n, j \neq i}

(1)

Set

σ (q)

as the density function deffned over the target region Q. A higher value of

σ (q)

implies that some properties that define an attribute are denser or more speciffc to a location. The properties are usually represented by the measurements within the sensor network, and the cost function within the bounded region Q is defined as follows:

E (p) = \sum_{i = 1}^{n} E^{i} = \sum_{i = 1}^{n} \int_{v_{i}} | | q_{j} - p_{i} | |^{2} σ (q_{j}) d q .

(2)

The number of wireless sensors is usually limited and non-continuous, assuming that the sensor network is uniformly distributed over the target area

S^{2}

. Then, the cost function within the bounded region Q can be expressed as Eq. (3). Du et al.²⁹ presented a detailed derivation of the CVT algorithm to minimize the cost function.

E (p) = \sum_{i = 1}^{n} E^{i} = \sum_{i = 1}^{n} \sum_{j \in ν_{i}} | | q_{j} - p_{i} | |^{2} σ (q_{j}) .

(3)

We explain how the system ensures that multiple robots maintain the proper formation during the mission and how the formation reduces the overall energy loss. We transform the formation problem into the tracking problem of centroids in the Voronoi interval. The density function is adjusted so that the robot maintains the Voronoi cells of similar masses. The mass of the cell can be expressed as a continuous version and discrete version, as shown in Eqs. (4) and (5), respectively. The calculation of the target point (centroid) is shown in Eq. (6), and Eq. (7) presents the discrete version of the centroid.

M_{v_{i}} = \int_{v_{i}} σ (q) d q .

(4)

M_{v_{i}} = \sum_{j \in v_{i}} σ (q_{j}) .

(5)

p_{i} * = \frac{\int_{v_{i}} σ (q) q d q}{\int_{v_{i}} σ (q) d q}, i = 1, 2, 3 \dots, n .

(6)

p_{i} * = \frac{\sum_{v_{i}} \sum_{j = 1}^{n} σ (q_{j}) q_{j}}{\sum_{v_{i}} \sum_{j = 1}^{n} σ (q_{j})}, i = 1, 2, 3 \dots, n .

(7)

Assume the average mass of .. of all cells in the region Q. The difference between the Voronoi cells is expressed as the difference in the average mass. When the difference reaches a minimum, the default mass of each cell is approximately equal, that is the workload of a single robot is approximated. The Voronoi mass difference is expressed as follows:

E_{2} (V) = \sum_{i = 1}^{n} E_{2}^{i} = \sum_{i = 1}^{n} | | M_{v_{i}} - \bar{M} | |^{2} .

(8)

In summary, keeping

E (p)

and

E_{2} (V)

to their minimum values during the task execution of a multi-robot system minimizes the information interaction energy loss and maintains the target formation, that is, a reasonable formation control is transformed into the multi-objective optimization problem.

CVT algorithm under constraints

According to the geometric shape, the task environment can be divided into convex and non-convex areas. If any two points in a certain area can be connected by a line segment that has no common point with the graph, then this area is called a convex area, as shown in Figure 1(a). The convex area is defined algebraically as shown below.

Figure 1.

Task area classification. (a) Convex area. (b) Non-convex area.

In an area $Q \subset S^{2}$ , for any two points $x, y \in Q$ and a real number $a (0 \leq a \leq 1)$ , if

a x + (1 - a) y \in Q,

(9)

the area Q is a convex area.

In the Euclidean space, a convex set is a set of points in which the points on a straight line between every two points are included in the set of points. For convex areas, the controller based on the Lloyd algorithm can ensure that the generated points asymptotically converge to the centroid of their respective Voronoi areas, as shown in Figure 1(a). However, this result has some problems for non-convex areas. Figure 1(b) depicts the centroid position generated by the CVT algorithm at the initial time, but there is an obstacle (big black disk) between the robot (brown circle) and the target point (blue cross). This will lead to errors in the Euclidean distance calculation. Therefore, the CVT partition in a non-convex environment needs improvement.

To solve the problem of robot Voronoi tessellation in different areas under different density functions, the CVT algorithm calculates the position of the center of mass in each Voronoi cell. First, it is necessary to ensure that the Voronoi cells cannot contain obstacles inside the region.

Assuming that the Voronoi area corresponding to each robot in the swarm robotics is $v_{i}, i = 1, 2, \dots n$ , and the obstacle area in the environment is $w_{i}, i = 1, 2, \dots n$ , to ensure that the Voronoi area does not contain obstacles, the corresponding Voronoi area of the robot should be

V_{i} = v_{i} - w_{i}, i, j = 1, 2 \dots n .

(10)

If there are no obstacles or convex areas in the environment, then

w = 0

and the Voronoi area will not change. When there is an obstacle in the Voronoi area, the original Voronoi area will regard the boundary of the obstacle area as the inner boundary to redefine the Voronoi area.

After ensuring that each Voronoi cell is within the region, the resulting centroid position may appear outside the region,considering that the region is non-convex.Hence, the centroid position needs to be restricted.This is accomplished by strictly constraining the centroid of the Voronoi region within the Voronoi region, as shown in Eq. (11).

p_{i} * \subset V_{i}, i = 1, 2 \dots n .

(11)

If the constraint solution is not obtained, it indicates that the position of the centroid is not in the Voronoi cell. Then, the center of mass of the Voronoi region will not be updated until the other Voronoi cells change, which again affects the updating of the center of mass.

The target position and robot position relative to the obstacles are divided into the four cases shown in Figure 2.

Figure 2.

Relationship between obstacles and paths.

The brown circle represents the current position of the robots, the blue cross represents the target position of the robots, the blue shaded area represents the obstacle, and the dotted line represents the safe path of the robots relative to the obstacle. The trajectory of the robot is marked with a red line. r is the radius of the obstacle, and $r^{'}$ is the radius of the robot. R is the distance between the center of the robot and the center of the obstacle, where $R = r + r^{'}$ , and $θ$ is the central angle of the arc path.

In Figure 2(a), the path between the robot and the target point does not pass through the obstacle, and the robot can directly track its target. At this time, $θ = 0$ .

In Figure 2(b), the path between the robot and the target point is tangential to the obstacle; hence, the robot can still track the target directly, and $θ = 0$ .

In Figure 2 (c), the straight path between the robot and the target point is blocked by the obstacle. The robot moves along the arc with radius R to avoid collision with the obstacle, and $0 < θ < π$ .

In Figure 2 (d), the line between the robot and the target point passes through the center of the obstacle, that is, $θ = π$ . The algorithm steps for this case are shown in Table 1.

Table 1.

Constrained CVT algorithm.

	Constrained CVT algorithm
1	Randomly generate n points
2	for $(i \leq n)$
3	Calculate the Voronoi cell of the i point;
4	if there are no obstacles in the Voronoi cell；
5	integrate the Voronoi cell to find the centroid；
6	else replace the inner boundary of the Voronoi cell with the boundary of the obstacle;
7	calculate the centroid
8	End
9	if the straight line connecting the robot and the target point is blocked by obstacles;
10	plan a new path; // shown in Figure 2
11	drive the robot to the centroid;
12	else drive the robot to the centroid;
13	end
14	if the algorithm convergence;
15	Break; // Stop the iteration
16	else Jump to step 2, calculate the Voronoi cell again;
17	End
18	End

$n$ points are randomly generated in the region. First, the Voronoi cells of each generated point are calculated, and the presence of obstacles in the Voronoi cells is examined. If there is no obstacle, the centroid of the cell is calculated. Otherwise, the boundary of the obstacle is replaced by the inner boundary of the Voronoi cell. Next, the presence of obstacles on the straight line connecting the robot and the target point is examined. If the path is blocked by an obstacle, plan a new path and drive the robot to the centroid. If not, drive the robot to the centroid directly. The last step is to judge whether the algorithm convergence meets the termination condition, and then interrupt or jump to the second step to calculate the Voronoi cells again .

Proof of convergence

The robot control rate in CVT is expressed as

{\ddot{p}}_{i} = u_{i} = - k_{p} (p_{i} - p_{i} *) - k_{d} {\dot{p}}_{i},

(12)

where

p_{i} *

is the center of mass of

V_{i}

, and both

k_{p}

and

k_{d}

are positive constants.

k_{p}

is a coefficient of proportionality,

k_{d}

is the damping coefficient, and

{\dot{p}}_{i}

is the velocity of the robot i. These two terms eliminate the oscillations as the robot approaches its destination.Consider the following Lyapunov function for the control law:

V (p_{i}) = \frac{1}{2} k_{p} (| | p_{i} - p_{i} * | | - 2 R + π R)^{2} + \frac{1}{2} {\dot{p}}_{i}^{2},

(13)

where r is the radius of the obstacle and

r^{'}

is the radius of the robot. To avoid collision between the robot and obstacle, the definition of the obstacle radius R in the algorithm includes the robot radius r, that is,

R = r + r^{'}

. Find the derivative of the Lyapunov function as

\dot{V} (p_{i}) = k_{p} (p_{i} - p_{i} *) {\dot{p}}_{i} + {\dot{p}}_{i} {\ddot{p}}_{i} .

(14)

Factor out

{\dot{p}}_{i}

and substitute the control law for

{\ddot{p}}_{i}

\dot{V} (p_{i}) = {\dot{p}}_{i} (k_{p} (p_{i} - p_{i} *) - k_{p} (p_{i} - p_{i} *) - k_{d} {\dot{p}}_{i}) .

(15)

Cancel the terms and simplify the expression as

\dot{V} (p_{i}) = - k_{d} {\dot{p}}_{i}^{2} \leq 0.

(16)

Therefore, the Lyapunov derivative is negative semi-definite. However,

V (p_{i}) = 0

if and only if

p_{i} = 0

. Therefore, by LaSalle's principle, the robots asymptotically converge to a specific CVT

(p_{i} = p_{i} *)

if the tessellation set is finite.

u_{i} = - k_{p} (p_{i} - p_{i} *) \Rightarrow p_{i} = p_{i} * .

(17)

Simulation of convex and non-convex areas

Formation control in convex (obstacle-free) areas

The simulation experiments in this paper were performed on a MATLAB simulation platform with Windows OS. In MATLAB, the robot was assumed as a first-order dynamics model, and the specified area size was $10 \times 10 m^{2}$ . It was proposed to control four robots in this area to complete the square and line formations. The area was equipped with 1000 × 1000 uniformly placed virtual sensor. As shown in Figure 3(a) and Figure 4(a), the x-axis and y-axis are divided into 1000 grids, each with a unit of 1 cm. The sensor readings were fitted to the density function of the CVT algorithm by modifying the density function to guide the centroid toward the desired position, thereby controlling the robot to track the centroid to form the desired formation. As shown in Figure 5, the initial positions of all four robots are near the origin.

Square formation

Figure 3.

Multi-robot square formation simulation. (a) Density function. (b) Square formation. (c) Position error. (d) Cost function value.

Figure 4.

Multi-robot line formation simulation. (a1) Density function. (b) Square formation. (c) Position error. (d) Cost function value.

Figure 5.

Initial positions of the robots in the multi-robot system.

The desired formation determined the value of the density function. For a square formation, the density function was initialized to a constant value of 1. As shown in Figure 3(a), when the initial density function in the region Q was constant, the region was partitioned using the Voronoi algorithm, and the robot tracked the location of the generated point to complete the square formation. The unit of each grid is 1 cm, and the lengths of the x-axis and y-axis are both 10 m. In Figure 3(b), the solid lines of different colors represent the robot's tracking trajectory.

In the initial phase of the experiment, the positions of the robots were randomly distributed near the origin, where the value of the cost function was the maximum. As the CVT algorithm was iterated, the generated point tended to move toward the centroid of each Voronoi cell, and the value of the cost function gradually decreased. When the generated points in the Voronoi cell reached their respective centroid positions, the cost function value was zero and the CVT algorithm converged, as shown in Figure 3(d). Figure 3(c) shows the position error between the robot and its corresponding centroid. At the beginning, the sum of the errors between the robots’ positions and the Voronoi generated points was the maximum. As the iterations of the CVT algorithm progressed, the error fluctuations gradually decreased owing to the PID controller, and the distance between the generated point and the centroid tended to the gradient descent. When the CVT algorithm converged, the error between the robots’ positions and the generated points eventually tended to zero.

Line formation

Based on the density function as a constant, the line formation of the robots can be completed by adjusting the density function to a “V”-shaped function. The expression for the “V”-shaped function is shown in Eq. (10). The robots will remain in the high-density areas so that they can create the line formation. This density function distribution is shown in Figure 4(a), with a region size of 10 × 10 m². The line formation is shown in Figure 4(b).

In the line formation process, the robots need to first complete the Voronoi partition under the constant density function distribution; hence, the cost function value at the initial moment is the maximum. During the iteration of the CVT algorithm, the controller overshoots due to integral saturation. The value of the cost function decreases significantly when the number of iterations is approximately 80. The optimal solution is reached again when the number of iterations is approximately 270, as shown in Figure 4(d). Figure 4(c) shows the position error between the robot and its corresponding centroid. As with the square formation, the robot's position error and the generated point are most significant during the initial stage of control. During the iteration of the CVT algorithm, the error fluctuation decreases gradually due to the gradient descent of the distance between the generated point and the centroid.

σ (x, y) = e x p^{- δ (a | x - x_{c} | + b (y - y_{c}) + c)}

(18)

Non-convex (obstacle) swarm control

The task environment of swarm robotics usually contains irregular non-convex areas. For strictly convex areas, the Lloyd algorithm's controller can ensure that the generated points asymptotically converge to the centroid of their individual Voronoi cells. Unfortunately, this result is inapplicable to non-convex areas and convex areas with obstacles. As shown in Figure 6(a), ten points were randomly generated in a convex environment. Based on the Lloyd algorithm, the generated points finally converge to the centroid of their respective regions to obtain the distribution. It is assumed that there are circular obstacles in the environment and the area becomes non-convex. Under this condition, if the same method as in the convex environment is applied here, the center of mass will occur at the obstacle, as shown in Figure 6(b). This distribution method is unreasonable, as the robot can move only outside the obstacle, as shown in Figure 6(c).

Figure 6.

Effect of obstacles on CVT. (a) Convex area. (b) CVT area over. (c) Non-convex area laps with obstacle.

To conclude, we consider that the following two fundamental issues must be satisfied to allow the use of swarm robotics in a non-convex environment:

(i) The divided Voronoi area cannot be within an obstacle;

(ii) The center of mass of each Voronoi cell cannot be inside the obstacle.

Only by strictly limiting the position of the Voronoi cell and the center of mass can it be ensured that when the robot follows the corresponding center of mass to cover the environment, it remains inside the environment area.

A set of simulations was designed for the distribution of robots in a non-convex environment with different parameters. By adjusting the size, distribution, and number of obstacles, the influence of the obstacles on the value of the cost function was examined, and the effects of different parameters on the distribution efficiency of swarm robotics were observed to identify the most energy-efficient swarm control scheme. Furthermore, in order to determine the optimal number of robots, simulations were conducted with fixed obstacle sizes. By adjusting the number of robots, the impact of robot quantity on coverage rate and coverage efficiency was observed. The criterion for judging the stability of the swarm robotics is that the average value of the sum of errors between the position of each robot and its center of mass is less than 0.001, that is, the swarm reaches the target position.

Radius changes of obstacle

The time consumed, that is, the total iteration time, refers to the time required for the system to reach stability. In this study, the experiments were conducted with Microsoft Windows 10 operating system, Intel(R) Core (TM) i7-10750H, and 512 GB memory.

The sizes of the obstacles were changed to verify the performance of the algorithm proposed in this paper under obstacles of different sizes. The 10 × 10 environment was again used in this experiment. Figure 7(a), (b), (c), and (d) shows the final distributions of multiple robots for obstacles of different sizes when the obstacle center is at (5,5) and the radius r is 0, 1, 2, 3. The corresponding cost function is shown in Figure 8. The minimum values were obtained when the number of iterations was 124, 113, 82, and 102, respectively. It is considered that the swarm reached the stable state. The experimental data are shown in Table 2.

Figure 7.

Different sizes of obstacles. (a) r = 0. (b) r = 1. (c) r = 2. (d) r = 3.

Figure 8.

Cost function values for different sizes of obstacles. (a) r = 0. (b) r = 1. (c) r = 2. (d) r = 3.

Table 2.

Different size of obstacles.

Radius	No. of robots	Center	Iterations	Time consumed (s)
r = 0	10	(5,5)	124	43.48
r = 1	10	(5,5)	113	37.27
r = 2	10	(5,5)	82	28.97
r = 3	10	(5,5)	102	36.80

From Figure 7 and Table 2, it can be seen that the size of the obstacle affects the degree of non-convexity of the environmental area. The algorithm can ensure a reasonable partition of the environment for obstacles of different sizes. However, a greater number of obstacles will increase the computation effort required for the algorithm. When the obstacle radius gradually increases, the area that the robot can effectively divide becomes increasingly smaller, which decreases the number of iterations and time in the steady state. When $r = 2$ , the time required is the smallest. Because the area is fixed in the simulation, the obstacle area in the case of $r = 3$ is nearly 1/3 of the total area, and the amount of calculation increases dramatically, resulting in increased overall iteration time.

Different locations of obstacle

In the virtual environment, in addition to obstacles being of different sizes, their positions are usually random. The experiment reported here verified the effectiveness of the proposed algorithm when obstacles occur at different positions. The 10 × 10 environment was again used in this experiment. For

r = 2

, the positions of the center of the obstacle were selected as (3,3), (5,5), (7,7). The experimental results are shown in Figure 9 and the cost function is shown in Figure 10. The minimum values were obtained when the number of iterations was 107, 82, and 90, respectively. The experimental data are presented in Table 3.

Figure 9.

Different location of obstacles. (a) Position (3,3). (b) Position (5,5). (c) Position (7,7).

Figure 10.

Cost function values for different obstacle locations. (a) Position (3,3). (b) Position (5,5). (c) Position (7,7).

Table 3.

Different location of obstacles.

Center	Radius	No. of robots	Iterations	Time consumed (s)
(3,3)	r = 2	10	107	39.24
(5,5)	r = 2	10	82	28.97
(7,7)	r = 2	10	90	33.02

Figure 9 and Table 3 show that the algorithm can ensure a reasonable partition of the area for different locations of an obstacle. Figure 10 shows the graphs of the cost function for different locations of the obstacle. The iterations were the smallest when the obstacle was located at the center of the area. When the obstacle was at (3,3), which was close to the robots, the partition conditions of the Voronoi cell were not maintained; hence, the algorithm recalculated the centroids of the Voronoi cell, which increased the computation cost. The robots moved in the opposite direction when the obstacle was at (7,7). The obstacle was encountered toward the end of the area. Therefore, the fluctuation was not apparent for a few iterations at the beginning. In the range of 20 to 30 iterations, the cost function fluctuated wildly but eventually stabilized. However, in these two cases, the robots overcame several situations where there were obstacles in the environment. When the obstacle was at (5,5), the iterations were smaller than in the other cases because the initial positions of the robots were significantly different from the desired locations. It can be seen that as the obstacle position moved closer to the center, the iteration time required was longer.

Different quantities of robots

This section verifies the robustness of the algorithm to variation in the number of robots. Figure 11 shows a graph of the experimental results and Figure 12 shows the graphs of the cost function for different numbers of robots in the swarm. The minimum values were obtained when the number of iterations was 82, 392, and 619 for 10, 50, and 100 robots in the swarm, respectively.

Figure 11.

CVT for different numbers of robots in the swarm. (a) n = 10. (b) n = 50. (c) n = 100.

Figure 12.

Cost function value for different numbers of robots in the swarm. (a) n = 10. (b) n = 50. (c) n = 100.

Figure 11 and Table 4 show that the algorithm can ensure a profitable partition of the environment for different numbers of robots in the swarm. If the other parameters remain unchanged, as the number of robots increases, the amount of calculation inevitably increases, and the number of iterations required for system area stabilization will increase accordingly. It can be concluded that the greater the number of robots in the swarm, the greater will be the amount of calculation and iteration time.

Different quantities of obstacles

Table 4.

Different numbers of robots in the swarm.

No. of robots	Center	Radius	Iterations	Time consumed (s)
10	(5,5)	r = 2	82	28.97
50	(5,5)	r = 2	392	213.81
100	(5,5)	r = 2	619	407.97

It can be seen from Figure 13, Figure 14, and Table 5 that as the number of obstacles increased, the number of iterations for achieving the steady state increased. Figure 14(c) shows that for iterations from 1 to 90, owing to the obstacle position and the initial position of the swarm, the swarm robotics system repeatedly calculated the centroid positions in the lower right corner area, resulting in the robots staying in this area for a long time. Therefore, the cost function between 30 and 90 iterations was zero. When the 90th iteration recalculated the robot's Voronoi cell and the center of mass, the position of the center of mass was iterated again, resulting in a large error. The system stabilized at 188 iterations. The iteration time was the shortest in the presence of one obstacle, and increased as the number of obstacles increased.

Figure 13.

CVT for different numbers of obstacles. (a) n = 10. (b) n = 50. (c) n = 100.

Figure 14.

Cost function values for different numbers of obstacles. (a) N = 1. (b) N = 2. (c) N = 4.

Table 5.

Different numbers of obstacles.

Starting point	No. of robots	No. of obstacles	Radius	Iterations	Time consumed (s)
（1,1）	10	1	r = 2	82	28.97
（1,1）	10	2	r = 2	147	47.38
（1,1）	10	4	r = 1	180	52.68

The size, number, and distribution of the obstacles and the number of robots were varied in the above four simulations. In the non-convex area, the CVT algorithm has a good control effect on the swarm robotics, and its efficiency is not affected by the size of the obstacles. The position of the objects has little effect on the formation, and its efficiency is affected only by the number of robots. Thus, the constrained CVT algorithm can realize the formation control of swarm robotics.

Optimal Number of Robots

This section verifies how to determine the optimal number of robots given the size of the task area. If the number of robots is too small, it is difficult to achieve full coverage of the task area. Conversely, if the number of robots is too large, it will lead to resource wastage. Therefore, it is necessary to determine the optimal number of robots to achieve coverage of the entire target area.

Determining the initial number of robots

Firstly, consider the basic perception model of the robots, as illustrated in Figure 15. Here, $R_{s}$ represents the radius of the robot communication range, meaning the entire disk within this range is the robot's perception area. If a point within the task area lies within the robot's perception range, the robot's perception probability for that area is recorded as 1. Conversely, if it lies outside the perception range, the robot will not perceive it, and the perception probability is recorded as 0.

Figure 15.

Perception mode.

Therefore, for any robot node $p_{i}$ and task area point q, the perception probability of the robot to the environment is given by:

P (q, p_{i}) = {\begin{matrix} 1, & ‖ q - p_{i} ‖ \leq R_{s} \\ 0, & otherwise \end{matrix}

(19)

After defining the coverage range of robots, the problem of determining the optimal number of robots needs to be addressed. Based on existing literature,³⁵ the number of robot nodes can be estimated by dividing the task area into multiple coverage sets. This method requires that any point within the task area can be perceived by at least one robot node, as illustrated in Figure 16.

Figure 16.

Area coverage.

The communication diameter of a robot is defined as

d = 2 R_{s} .

(20)

The length of each side of the square is

L = d \sin 45 \circ .

(21)

If the coverage area is a rectangular region with side lengths x and y, then the number of robots can be determined by Equation (22).

n = x y / L^{2} .

(22)

Determining the optimal number of robots

Typically, the methods mentioned in the previous section are only used to determine the initial number of robots, while the specific optimal number of robots needs to be determined based on different task requirements and performance metrics. Two common coverage performance evaluation metrics are defined as follows:

Coverage rate

The node coverage rate is used to evaluate the extent to which robots cover the task area. It is determined by the ratio of the sum of the coverage areas of robot nodes to the total area of the task area, calculated as:

R = \frac{\sum_{j = 1}^{x y} P (q, p_{i})}{x y} .

(23)

Coverage efficiency

Robot node coverage efficiency is mainly used to measure the utilization rate of robot node coverage. A higher value indicates a higher utilization rate of robot node coverage. It is determined by the ratio of the sum of effective coverage areas of robot nodes to the sum of all robot perception range areas, calculated as:

T = \frac{R x y}{n (π R_{s}^{2})} .

(24)

In addressing the question of determining the optimal coverage range and quantity of robots, varying the number of robot nodes will yield different coverage effects. The initial number of robots calculated from Equation (22) is 45. To validate the issue of the optimal number of robots, this section conducts separate comparative analyses on coverage rate and coverage efficiency using different numbers of robots. Table 6 serves as the experimental parameter settings table.

Table 6.

Experimental parameter settings.

Parameters	Parameter values
Area size	88 m²
Perception radius	1 m
Number of robots	30, 45, 60

In this experiment, the quantities of robots selected are 30, 45, and 60 respectively, with the obstacle's central position set at (5,5). The experimental results are illustrated in Figure 17. Figure 18 presents the variations in the cost function for different numbers of robots. The cost function reaches its minimum values at iterations 128, 228, and 231 respectively. Experimental data are detailed in Table 7.

Figure 17.

CVT for different numbers of robots. (a) n = 128. (b) n = 228. (c) n = 231.

Figure 18.

Cost function value for different numbers of robots. (a) n = 128. (b) n = 228. (c) n = 231.

Table 7.

Performance metrics under different numbers of robots.

Number of robots	30	45	60
Coverage rate	93.06%	99.26%	99.73%
Coverage efficiency	86.34%	69.06%	46.27%

According to Table 7 and Figure 17, it can be observed that as the number of robots continues to increase, the overall coverage efficiency decreases continuously.

However, based on the principle of robot coverage efficiency formula, it is understood that with a constant task area, a smaller number of robots ensures a smaller overlapped coverage area, resulting in higher node coverage efficiency. As the number of robots increases, the overlapped area grows, leading to a decrease in node coverage efficiency while the coverage rate of robots increases. Therefore, the selection of the optimal number of robots needs to be dynamically adjusted based on the requirements of the task. Starting from the initial number of robots obtained, a higher coverage rate requires an increase in the number of robots, while a higher coverage efficiency necessitates a reduction in the number of robots.

Swarm robotics experiment

ROS robot formation control

The experiments conducted in this study used the Turtlebot3 robot, which is a low cost, development friendly mobile robot platform based on the robot operating system (ROS). Each Turtlebot3 can independently build a local map of its surroundings through the mapping function package of ROS, which provides SLAM (simultaneous localization and mapping) for LiDAR and a grid-based 2D map based on the LiDAR input and attitude data. In the actual experiment, the robot's motion model is a second-order dynamics model. The lab map was selected as the designated area for the robot formation Q. The area was free of obstacles, as shown in Figure 19(a). Four Turtlebot3 robots were randomly placed on the map, and the robots used the SLAM technology to complete their positioning. The results of the experiment are shown in Figure 19(b). The experiments were set such that when the Voronoi cost function reached the optimal solution, the robots reached the position required for the desired formation under that density function, that is, they achieved the target formation. The difference between the simulation and actual experiment was that the density function was not initialized as a constant at the initial stage of the experiment, but rather, the target expected formation was determined before the commencement of the experiment.

Figure 19.

Experimental environment for multi-robot formation. (a) Experimental. (b) SLAM map (RVIZ).

When the density function was set to a constant, the map was evenly divided into four generated points, as shown in Figure 20(a). Figure 20(c) shows the initial position of the robot. The CVT algorithm converged when the number of iterations was approximately 200, as shown in Figure 20(b), at which point the cost function reached its minimum value. The robots completed a square formation on the map, as shown in Figure 20(d).

Figure 20.

Square formation experiment. (a) CVT area. (b) Cost function. (c) Initial position. (d) Final position.

Figure 21(a) depicts the Voronoi partition under the “V”-shaped density function. The initial positions of therobots are shown in Figure 21(c). Figure 21(b) depicts the cost function for the “V”-shaped distribution of the density function. The initial value of the cost function was approximately 6.3 × 104. When the CVT algorithm converged after approximately 120 iterations, the robots completed a straight formation, as shown in Figure 21(d).

Figure 21.

Line formation experiment. (a) CVT area. (b) Cost function. (c) Initial position. (d) Final position.

The algorithm was also validated on the dynamic formations achieved with four Turtlebot3 robots. The experiments involved dynamic formations in the obstacle-free and random obstacle scenarios. The robots were initially positioned in the upper left corner of the map, marked by stars in Figure 22(a) and (b). The target point was located in the lower right corner of the map, marked by red flags in Figure 22(a) and (b), and the black line denotes the motion trajectory of the multi-robot system. We modified the map created by SLAM. The obstacles of irregular shapes in the original map were transformed into circular obstacles, as shown in Figure 22 (c). The modified map was used as the environmental information for the robots.

Figure 22.

SLAM map. (a) Obstacle-free. (b) Random obstacles. (c) Circular obstacles.

Obstacle-free scene (convex)

As shown in Figure 23, the experimental scene can be regarded as a rectangle and there were no obstacles in the area, which represented a convex environment. The robots were in the upper left of the area at the instant t₁. The Voronoi tessellation algorithm converged when the robots reached the centroid of each Voronoi cell, as shown is Figure 23(a). Figure 23(b) and (a) shows that the robots move to the target point at the instants t3 and t4, respectively. Finally, they reach the target at t₄ as shown in Figure 23(d). In comparison with the traditional CVT, our improved CVT can divide the area in a non-convex environment as well, which allows the robots to move to the target point in a suitable formation.

Figure 23.

Convex (obstacle-free) scene. (a) t₁ = 0 s. (b) t₂ = 64 s. (c) t₃ = 75 s. (d) t₄ = 120 s.

Randomly distributed obstacles scene (non-convex)

To verify the feasibility of the constrained CVT, we added obstacles to the scenario shown in Figure 23 to transform it to a non-convex environment, as shown in Figure 24 The robots were positioned in the upper left of the map at the instance t₁. In Figure 24(a), all robots reached the centroid of each Voronoi cell. Then, the constrained area Q moved along the desired direction while obtaining the two- dimensional map information as the region density. For a non-zero threshold, we set the density value to zero so that the robots could avoid the obstacles, as shown in Figure 24(b) and (c). In Figure 24(d), the robots reached the desired position at t₄. The experiment showed that the constrained CVT could separate the non-convex environment from the obstacles to ensure that the robots could bypass the obstacles and reach the target position. This proved the feasibility of the algorithm for a non-convex environment.

Figure 24.

Random obstacles. (a) t₁ = 0 s. (b) t₂ = 32 s. (c) t₃ = 80 s. (d) t₄ = 102 s.

In summary, the experimental results showed that the constrained CVT algorithm could be applied to both convex and non-convex regions. The time consumed by the algorithm was 120 s for the experiment on the convex area and 102 s for the experiment on the non-convex area. The results are consistent with the simulation results in Section 4.2-A. It can be inferred from the experiments that for regions of the same size, obstacles in the non-convex region reduce the area that can be effectively divided as well as the number of iterations; hence, the process is less time consuming when compared with the convex region.

Conclusion and future work

This study successfully applied the CVT algorithm to the ROS, saving a lot of research work. Through the CVT algorithm, the central processing unit processes position information and sends particle information from the CVT partition to each robot as its target point, thereby achieving control over multiple robots. Compared to other multi-robot system designs, the CVT algorithm converges to the centroid of each Voronoi cell, and each robot moves within its Voronoi cell, which can minimize collisions. The performance of the algorithm was validated through 5 different simulations, and the results showed that the constrained CVT algorithm is not sensitive to the size and position of obstacles, but an increase in the number of robots will increase the total number of iterations. The optimal number of robots needs to be dynamically selected based on different performance indicators. The experimental results show that using the CVT algorithm can ensure that multiple robots form a formation according to the optimal distribution in the designated area, while increasing the redundancy and robustness of the system. Easy to dynamically switch formations based on external environment, improving the flexibility and scalability of formations. However, the current method assumes that all robots are in good condition and have the same power consumption, and future research will consider more practical scenarios, such as irregular and dynamic obstacle environments.

Footnotes

Acknowledgements

This study was co-supported by the funds of National Natural Science Foundation of China Youth Fund (62103315, 62303368); Key R & D Projects in Shaanxi Province (2022GY-238); Special Project for Guiding Technological Innovation in Shaanxi Province (2022QFY01-16, 2023-ZDLNY-61); Ministry of Education Key Laboratory of Information Fusion Technology Open Fund Project (Project No. 202312-IFTKFKT-007).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Key R & D Projects in Shaanxi Province, Ministry of Education Key Laboratory of Information Fusion Technology Open Fund Project, National Natural Science Foundation of China Youth Fund, Special Project for Guiding Technological Innovation in Shaanxi Province (grant numbers 2022GY-238, 202312-IFTKFKT-007, 62103315, 62303368, 2022QFY01-16, 2023-ZDLNY-61).

ORCID iDs

Kai Cao

Yangquan Chen

Song Gao

Manying Cheng

References

Zhang

. A solution for searching and monitoring forest fires based on multiple UAVs. In: Proceedings of the inter-national conference on unmanned aircraft systems (ICUAS), Atlanta, USA, 11–14 June 2019, pp.661–666. IEEE

Luo

, et al. Target-enclosing affine formation control of two layer networked spacecraft with collision avoidance. Chin J Aeronaut 2019; 32: 2679–2693.

Qiu

Duan

. Multiple UAV distributed close formation control based on in-flight leadership hierarchies of pigeon flocks. Aerosp Sci Technol 2017; 70: 471–486.

Lee

Chwa

. Decentralized behavior-based formation control of multiple robots considering obstacle avoidance. Intell Serv Rob 2018; 11: 127–138.

Zhou

Wang

Schwager

. Agile coordination and assistive collision avoidance for quadrotor swarms using virtual structures. IEEE Transactions on Robotics 2018; 34: 916–923.

Chen

Meng

Xing

. Shape control of spacecraft formation using a virtual spring-damper mesh. Chin J Aeronaut 2016; 29: 1730–1739.

Wang

Zhang

. A self-organizing area coverage method for swarm robots based on gradient and grouping. Symmetry 2021; 13: 680.

Chen

Wang

. Formation control: a review and a new consideration. In: IEEE/RSJ international conference on intelligent robots and systems, Edmonton, AB, Canada, 2–6 August 2005, pp.3181–3186.

Ahn

H-S

. Formation control-approaches for distributed agents. Switzerland: Springer, 2019; ISBN :978-3-030-15186- 7.

10.

Moon

Frew

. Distributed cooperative control for joint optimization of sensor coverage and target tracking. In: Proceedings of the 2017 international conference on unmanned aircraft systems (ICUAS), Miami, USA, 13–16 June 2017, pp.759–766.

11.

Sun

Cassandras

Meng

. Exploiting submod- ularity to quantify near-optimality in multi-agent coverage problems. Automatica 2019; 100: 349–359.

12.

Teruel

Aragues

López-Nicolás

. A distributed robot swarm control for dynamic region coverage. Rob Auton Syst 2019; 119: 51–63.

13.

Xiong

Chen

, et al. Path planning of multiple autonomous marine vehicles for adaptive sampling using Voronoi-based ant colony optimization. Rob Auton Syst 2019; 115: 90–103.

14.

Pierson

Schwarting

Karaman

. Weighted buffered voronoi cells for distributed semi-cooperative behavior. In: Proceedings of the 2020 IEEE international conference on robotics and automation (ICRA), Paris, France, 2020 May 31–June 4, pp.5611–5617: IEEE.

15.

Cortes

Martinez

Karatas

, et al. Coverage control for mobile sensing networks. IEEE Trans Rob Autom 2004; 20: 243–255.

16.

Boardman

Harden

Martínez

. Limited range spatial load balancing in non-convex environments using sampling-based motion planners. Auton Rob 2018; 42: 1731–1748.

17.

Prabowo

Trilaksono

. Collision-free coverage con- trol of swarm robotics based on Gaussian process regression to estimate sensory function in non-convex environ- ment. Int J Electr Eng Inf 2019; 11: 125–143.

18.

Collins

Ghassemi

Esfahani

, et al. Scalable coverage path planning of multi-robot teams for monitoring non- convex areas. In: 2021 IEEE international conference on robotics and automation (ICRA), Xian China, 30 May–5 June 2021, pp.7393–7399.

19.

Stergiopoulos

Thanou

Tzes

. Distributed collaborative coverage-control schemes for non-convex domains. IEEE Trans Autom Control 2015; 60: 2422–2427.

20.

Nair

Guruprasad

. GM-VPC: an algorithm for multi- robot coverage of known spaces using generalized voronoi partition. Robotica 2019; 38: 1–16.

21.

Breitenmoser

Schwager

Metzger

, et al. Voronoi coverage of non-convex environments with a group of networked robots. In: 2010 IEEE international conference on robotics and automation, Anchorage, Alaska, 3-8 May 2010, pp.4982–4989.

22.

Palacios-Gasos

Tardioli

Montijano

, et al. Equitable persistent coverage of non-convex environments with graph- based planning. Int J Rob Res 2019; 38: 1674–1694.

23.

Caicedo-Nunez

Zefran

. A coverage algorithm for a class of non-convex regions. In: 2008 47th IEEE conference on decision and control(CDC), Cancun, Mexico, USA, 9–11 December 2008, pp.4244–4249.

24.

Diaz-Mercado

. Multi-Robot control using coverage over time-varying non-convex domains. In: 2020 IEEE international conference on robotics and automation (ICRA), Paris, France, 31 May–31 August 2020, pp.4536–4542.

25.

Howard

Mataric

Sukhatme

. Mobile sensor network deployment using potential fields: a distributed, scalable solution to the area coverage problem. Distrib Auton Rob Syst 2002; 5: 299–308.

26.

Renzaglia

Martinelli

. Distributed coverage control for a multi-robot team in a non-convex environment. In: IEEE IROS09 3rd workshop on planning, perception and navigation for intelligent vehicles, 2009, pp.1–6.

27.

Choi

Wang

. Visibility-based coverage of mobile sensors in non-convex domains. In: 2011 eighth international symposium on voronoi diagrams in science and engineering, 2011, pp.105–111.

28.

Mahboubi

Sharifi

Aghdam

, et al. Distributed control of multi-agent systems with limited communication range in the fixed obstacle environments. IEEE Access 2019; 7: 118259–118268.

29.

Tahir

Bling

Haghbayan

, et al. Swarms of unmanned aerial vehicles – a survey. J Ind Inf Integr 2019; 16: 100–106.

30.

Xing

Cai

. Survey of intelligent dispatching technology of swarm unmanned aerial vehicle. Acta Autom Sin 2020; 46.

31.

Chen

Wang

Liang

. Automatic dynamic flocking in mobile actuator sensor networks by central voronoi tessellations. In: Proceedings of the 2005 IEEE international conference mechatronics and automation, Niagara Falls, Ont, Canada, July 29–August 1 2005, pp.1630–1635.

32.

Chen

Wang

Moore K

. Optimal spraying control of a diffusion process using mobile actuator networks with fractional potential field based dynamic obstacle avoidance. In: Proceedings of the 2006 international conference on networking, sensing, and control, Piscataway, NJ, USA, 23–25 April 2006, pp.107–112.

33.

Govindaraju

Leng

Qian

. Multi-UAV surveillance over forested regions. Photogramm Eng Remote Sens 2014; 80: 1129–1137.

34.

Qiang

Vance

Max

. Centroidal Voronoi tes- sellations: applications and algorithms. SIAM Rev 1999; 41: 637–676.

35.

El Khamlichi

Tahiri

Abtoy

, et al. A hybrid algorithm for optimal wireless sensor network deployment with the minimum number of sensor nodes. Algorithms 2017; 10: 80.