Abstract
Mobile sensor networks are an important part of modern robotics systems and are widely used in robotics applications. Therefore, sensor deployment is a key issue in current robotics systems research. Since it is one of the most popular deployment methods, in recent years the virtual force algorithm has been studied in detail by many scientists. In this paper, we focus on the virtual force algorithm and present a corresponding parameter investigation for mobile sensor deployment. We introduce an optimized virtual force algorithm based on the exchange force, in which a new shielding rule grounded in Delaunay triangulation is adopted. The algorithm employs a new performance metric called ‘pair-correlation diversion', designed to evaluate the uniformity and topology of the sensor distribution. We also discuss the implementation of the algorithm's computation and analyse the influence of experimental parameters on the algorithm. Our results indicate that the area ratio, φs, and the exchange force constant, G, influence the final performance of the sensor deployment in terms of the coverage rate, the convergence time and topology uniformity. Using simulations, we were able to verify the effectiveness of our algorithm and we obtained an optimal region for the (φ s , G)-parameter space which, in the future, could be utilized as an aid for experiments in robotic sensor deployment.
Keywords
Introduction
With the development of information technology, nowadays robots play a key role in industrial manufacturing, intelligence services and military assistance. Many scientists are currently conducting research and development in robotics applications. Due to the fact that informational sensors have become more efficient and inexpensive over time, mobile wireless sensor networks (WSNs) have been assigned an important role in autonomous robots with sensing and communication abilities. Mobile sensors collect a variety of data based on a distributed network that have the advantages of low power consumption, low cost, high reliability and self-organization. Such networks are widely used in robotics technology and in applications such as range sensing, target tracking, robot control, robot vision and environment monitoring, etc. [1].
Sensor deployment algorithms can be categorized as either deterministic or movement-assisted [2]. In special environments (e.g., remote or unmanned cases), sensor nodes are randomly distributed and only movement-assisted deployment can be applied. Since a good deployment method can offer robotics applications higher flexibility, practicability and adaptivity, deploying mobile sensors effectively is a key issue in current robotics WSN research. In general, in deployment algorithms for mobile sensor networks, each sensor knows its position, can communicate with other nodes, and can move to a better position according to collected or computed information in order to improve the coverage rate. Good intelligent deployment algorithms mainly include the following aspects: (1) processes to save sensor energy to extend the lifetime of network systems; (2) processes to obtain a uniform network topology in order to enlarge the coverage area; and (3) processes for accelerating the entire deployment procedure. Many approaches have been proposed for robotics sensor network deployment that have their own advantages, including: (1) computational geometry and grid-based methods [3–5]; (2) virtual force- or potential-based strategies [6–15]; (3) swarm intelligence and bionics algorithms [16–23]; and/or (4) a combination of the above approaches [24–26].
Among the above methods, virtual force-based strategies have been deemed effective, are popular and have become a hot topic in WSN research for robotics systems over the last decade. Recently, Zou and Chakrabarty presented a virtual force algorithm (VFA) using attractive and repulsive forces in order to determine the virtual motion paths and movement rates of sensor nodes. They also proposed a probabilistic target localization algorithm that can maximize the sensor field coverage [7]. To increase the robustness of the robot sensor network, Kribi et al. reported a serialized VFA with good coverage, connectivity and fault-tolerance [9]. Garetto et al. investigated a distributed sensor relocation algorithm based on virtual force for environmental control that can provide regular tessellation for a given geographical area. The algorithm provides an effective distributed method as compared with an optimal centralized solution [10]. Yu et al. proposed a new sensor deployment algorithm based on Van der Waals [12] and the virtual spring force [14]. The spring force-based algorithm resulted in proper damping for the deployment of sensor nodes in mobile WSNs [14]. Yuαŕs methods are capable of supplying a better coverage rate, more uniformity for the network configuration, and a moderate convergence time. Moreover, a novel VFA based on physical laws for a dusty plasma system was presented by Wen et al. The method is capable of forming an hexagonal topology in wireless sensor networks. These authors also introduced a performance metric based on the pair-correlation function in a crystal structure in order to evaluate the similarity between sensor network topologies resulting from the VFA and a perfect hexagonal structure [15]. The algorithm provides maximum field coverage with a minimum system cost.
Here, based on the work of Garetto, we report modifications and improvements for the basic VFA [10]. We first describe changes to the shielding rule using Delaunay triangulation in order to replace the original shielding rule that determines the adjacent nodes of a given node. Second, to evaluate the uniformity of the node distribution, a new performance evaluation metric (improved from paper [14, 15]) is introduced. Our simulation results indicate that the approach exhibits good performance in terms of the convergence time and uniformity for the network configuration. We also provide a detailed discussion of the factors influencing our algorithm and obtain an optimal region for the (φ s ,G)-parameter space.
The paper is organized as follows. In Section 2, we introduce the basic concept of the VFA and present our optimized node deployment solution, where we propose a new shielding rule for our algorithm and improve the metrics in order to evaluate it. In Section 3, we simply describe the implementation routine for our algorithm. Simulation results and discussions are provided in Section 4. We also verify the effectiveness of our algorithm and provide a detailed analysis of the corresponding parameter effects. Concluding remarks are provided in Section 5.
The Optimized Node Deployment Algorithm
Review of the VFA
Since each sensor node is regarded as a particle in the robotic sensor network, the virtual force algorithm is influenced by interactions among the particles of the physical system. In general, two assumptions are made in the physical system [7], as follows: (1) Each sensor node has a sensing radius, R s , and a communication radius, R c , (R s ≪ R c ). Within the sensing range, the sensor can detect local environmental conditions and can communicate with other nodes that fall into its communication range. (2) The position of all the nodes can be determined using various methods (e.g., GPS). All of the nodes can effectively move according to the calculation results of the algorithm. Thus, all of the nodes should eventually reach equilibrium through interactions of the virtual forces and form a stable lattice structure that corresponds to a sensor network topology.
Here, we adopt a second-order differential equation for node motion in the VFA [10, 14], whose nodes are subjected to variable motion because the second-order differential equation is normally easier to determine in dynamic equilibrium in a physical system. The equation of a given node i is as follows,
where m is the virtual mass of a sensor node,
When the distance d ik (t) between two given nodes i and k at time t is less than the communication radius R c , there is an exchange force between the two nodes. If A(i,t) is the set of nodes acting with exchange forces on node i at time t, then k ∊ A(i,t). We define the formula for total exchange forces on node i at time t as
where
The friction forces acting on a node i at time t can be written as
where
At the beginning of our algorithm, when the initial nodes are randomly deployed, all of the nodes are stationary. We then compare
In general, the issue of how to define adjacent nodes of a given node will determine the computational complexity of the total virtual force in a VFA simulation. Thus, a good shielding rule used to restrict adjacent nodes plays a key role in improving and simplifying the entire algorithm.
In a perfect hexagonal network topology, each node should be surrounded by six adjacentnodes. When the equilibrium distance D
m
= √3R
s
, the sensor network ideally forms a triangular tessellation and can yield the optimal hexagonal coverage [10, 15]. Garetto et al. introduced a shielding rule as follows: If k and k’ are two nodes within the communication range of node i, node k’ shields k only if k’ is closer to i and the angle formed by the vectors
Therefore, we present a new method for optimizing the shielding rule. We generate Delaunay triangulation in order to accurately ascribe all of the neighbours of each node. A Delaunay triangulation for a set of points P in a plane is a triangulation DT(P) such that no point in P is inside the circumcircle of any triangle in DT(P). The method maximizes the minimum angle for all of the angles of the triangles in the triangulation, which tends to avoid skinny triangles [27]. Yu et al. introduced this method into WSN applications and presented its effectiveness in adjacent node analysis [11].
We let A nei (i,t) be the set of first-round adjacent nodes surrounding node i generated by Delaunay triangulation at time t. We then computed the total virtual force from these accurate nodes in the set A(i,t)∩A nei (i,t). Formula (2) can then be rewritten as:
Since Delaunay triangulation can accurately determine the adjacent relationship of all the nodes, it is a useful tool for investigating the total force acting from second- or third-round adjacent nodes and - later on - enlarging the shielding range.
To objectively evaluate the performance of various VFAs, assessing performance from different perspectives is essential. Based on research on crystallization in a dusty plasma [28, 29], Wen et al. introduced a novel performance metric, known as the pair-correlation diversion (PCD), in order to evaluate the proximity between any network coverage topology and a perfect hexagonal topology [15]. The method is based on the pair-correlation function, g(r) = (N(r,Δ)S)/(2πrΔN), in a crystalline structure and represents the probability of finding two nodes separated by a distance r. In the equation, S is the area occupied by a 2D system, N is the number of nodes contained in S, and N(r,Δ) indicates the number of nodes located between r − Δ/2 and r + Δ/2 (see references [14, 15, 29] for details). The PCD works well when the node density of the sensor network is close to the node density of a perfect hexagonal configuration, but not so well when the node density is too large or too small.
Here, we make improvements to the PCD in order to evaluate the uniformity of the network configuration. Based on the definition of g(r), we further improve the performance metric formula as follows,
where Ω is the network topology to be analysed, Ω H is the perfect hexagonal topology, and r T is the bound of the radial distance. For a sensor network, when the value of δ(Ω,Ω H ) is small, the network coverage topology will approximate a perfect hexagonal structure. A perfect network configuration has a δ(Ω,Ω H ) value close to zero. Moreover, δ(Ω,Ω H ) can be applied in order to characterize the convergence rate of a VFA. Beginning with an initial node distribution, δ(Ω,Ω H ) decreases while the deployment progresses over time. When the network reaches the equilibrium state, the value becomes stable. We can then use the slope of the PCD function over the convergence time in order to specify the convergence rate.
In this section, we describe our implementation routine for the optimized VFA. At the beginning of the simulation, we pre-configure values for the initial parameters required for the configuration. Considering that N nodes are deployed in a circular area, S, the main process of the algorithm is as follows:
Set the centre of S at (0,0). Initially, N sensor nodes are randomly distributed within a circular area, S
init
, and the location of node i is
Use the Bowyer-Watson algorithm [30] in order to generate a Delaunay triangulation and obtain a relational structure among the nodes; then, choose the first-round of adjacent nodes for each node i.
Cyclically execute the following three steps for all of the sensor nodes i (i = 1,2, ⋅⋅⋅, N), as follows:
Calculate the distance between node i and the adjacent node in A nei (i,t), such as d ik . If d ik > R c , abandon the value. Rest nodes are located in A(i,t) ∩ A nei (i,t).
Combine formulae (4) and (3) in order to calculate the virtual force and the acceleration, expressed as
Use leapfrog integration in order to calculate the new velocity and the displacement of node i after time step dt; then, update the original velocity and the displacement. If the velocity of node i exceeds V max ’ set it to V max .
Repeat Steps 2 and 3 in order to calculate the acceleration, velocity and displacement for all of the nodes after each time step dt until the system reaches the equilibrium state.
Simulation Results and Discussion
In this section, to verify the effectiveness of the algorithm and to give details regarding the influence of experimental parameters on the optimized VFA, we provide our simulation results.
Effectiveness Verification
Using the parameters employed by Garetto et al. [10], we modified R s = 1, R c = 3 and D m = √3. Theoretically, the optimal coverage range is approximate to S perf = N D m 2√3/2 (about N hexagonal area, ignoring the approximate border). We used G = 0.1, k0 = 1.0e − 5, k v = 0.1, V max = 0.01 and dt = 1.0 as the initial parameters. For N = 500, all of the sensor nodes were randomly distributed around the centre of the deployment area. We choose an initial distribution area of S init = 0.7S perf = 525√3.
In Figure 1, we provide the initial distribution and the final distribution when the network reached an equilibrium state. The black point indicates a single sensor node and the blue circle indicates its sensing range. Obviously, the final network achieved a uniform configuration and a good coverage.

The initial node distribution locations (left) and the final distribution locations (right)
In Figure 2, we plot the Delaunay triangulation and the dual graph Voronoi diagram for the final network configuration. The green lines clearly indicate the triangular tessellation generated by Delaunay triangulation, and the blue-line structure generated by the Voronoi diagram provides a very good hexagonal coverage topology.

The final network coverage topology represented by the Delaunay triangulation method
In Figure 3, we present the PCD curve along the calculation time step. We determined that VFA-DT reached an equilibrium state after 1,150 time steps and that the corresponding δ(Ω,Ω H ) = 0.158. However, VFA-LJ required 1,900 time steps and the corresponding δ(Ω,Ω H ) = 0.206. The results indicate that VFA-DT converged faster than VFA-LJ and that the final distribution of VFA-DT was more approximate to a perfect hexagonal topology, which also means that our algorithm can obtain a better convergence time and uniformity compared with the general VFA of Garetto et al. [10].

A performance evaluation curve for δ(Ω,Ω H ) as a function of the time step, as determined from our optimized algorithm (VFA-DT) and the algorithm (VFA-LJ) of Garetto et al.[10]
Previously, we arbitrarily chose S init = 0.7S perf to verify our algorithm, although we know that S init also indicates the node deployment density, according to ρ = N / S, for a fixed node number, N. In an actual experiment, the node deployment area is easy to identify, so below we adjust the initial distribution area, S init , in order to discuss the influence of different node densities on sensor deployment. In general, ρ init = N / S init reflects the degree of sparseness for the sensor nodes. If ρ init is too large, the nodes should move a longer distance in order to be steady, and the convergence time should decrease with high energy consumption. If ρ init is too small, blank monitoring areas may appear after the sensors are deployed. The area is only a local optimum but not a global optimum.
By setting the ratio φ s = S init / S perf and keeping the other parameters the same, as in Section 4.1, in Figure 4 we present typical profiles for the PCD value δ(Ω,Ω H ) as a function of the time step for cases φ s = 0.1,0.2, ⋅⋅⋅, 1.2.

Typical profiles for PCD δ(Ω,Ω H ) as a function of the time step for specified values of φ s
In Figure 5, based on the results provided in Figure 4, we plot the convergence time and the PCD value δ(Ω,Ω H ) as a function of φ s using a cubic spline interpolation when the system first reaches the equilibrium state.

The convergence time and the PCD value δ(Ω,Ω H ) as a function of φ s when the system first reaches the equilibrium state
As can be seen in Figures 4 and 5, when φ s increases, S init increases (the initial node density ρ init deceases) for a fixed N and the convergence time decreases. When S init is in the interval of [0.1S perf , 0.9S perf ], the value of δ(Ω,Ω H ) is small and relatively stable, which indicates the final network topology is approximate to a hexagonal structure and that deployment is successful. When S init is in the interval of [0.9S perf ,1.2S perf ], the value of δ(Ω,Ω H ) significantly increases. The system converges quickly but the deployment result is not good, which may be caused by a local optimum with a blank monitoring area.

Typical profiles of PCD δ(Ω,Ω H ) as a function of the time step for the specified values of G
Based on Figure 5, we concluded that when φ s is in the interval of [0.7,0.9] or when S init is in the interval of [0.7S perf ,0.9S perf ], both the convergence rate and the deployment structure are relatively good. The corresponding optimal range of the initial node density, ρ init , is [(1.11ρ perf , 1.43ρ perf ].
From Equations (2) and (4), one can clearly see that the constant G significantly affects the total exchange force that determines the velocity of the node. Here, we discuss in detail the influence of different exchange force G values on the convergence time and the deployment results.
We fixed S init =0.7S perf and used the values for Section 4.2 for the other parameters. In Figure 6, we present profiles of the PCD δ(Ω,Ω H ) as a function of the time step for the specified cases of G=0.01,0.02, ⋅⋅⋅, 0.20. When G is small - say G < 0.03 - the system exhibits poor convergence all of the time. However, the convergence time decreases as the value of G increases. When G > 0.18, the sensor distribution ultimately obtains a bad network topology.
In Figure 7, we present the convergence time and the PCD value, δ(Ω,Ω H ), as a function of G using cubic spline interpolation according to the results of Figure 6. The system clearly exhibits a good convergence time when G is in the interval of [0.09,0.19] and when the value of δ(Ω,Ω H ) is small, while it is relatively stable when G is in the interval of [0.01,0.17]. Thus, the optimal range of G is [0.09,0.17] and both the convergence time and the network topology are relatively good.

The convergence time and the PCD value δ(Ω,Ω H ) as a function of G when the system first reaches the equilibrium state

The convergence time and the PCD value for δ(Ω,Ω H ) as a function of G for the specified values of φ s when the system first reaches the equilibrium state
Using the convergence time and the PCD δ(Ω,Ω H ) as performance indicators for sensor deployment, we obtained a detailed parameter investigation based on additional simulations and the optimal parameter space of φ s and G.
The optimal range of G for different φ
s
and the average convergence time and
under an optimal range of G
The optimal range of G for different φ
s
and the average convergence time and
Similar to the previous figures, in Figure 8 we plot the convergence time and the PCD value δ(Ω,Ω H ) as a function of G for the specified cases φ s =0.1,0.2, ⋅⋅⋅, 1.0 when the system first reaches the equilibrium state. Table 1 provides the optimal range of G for which both the convergence rate and the deployment results are relatively good for different φ s . We also provide the average δ(Ω,Ω H ) and the convergence time under the optimal range of G.
According to the data in Table 1, as shown in Figure 9, we draw a relationship diagram between φ s and the optimal ranges of G. The green curve displays the average convergence time under the optimal range of G for different φ s , and decreases as φ s increases. The red curve displays the PCD value δ(Ω,Ω H ) under this optimal range of G, which is relatively small. In general, sensor deployment requires a short convergence time in order to obtain a quick response for emergent situations and energy saving. In our simulations, when the convergence time is in the interval [500,1500], the system has the best performance. The purple shadow in Figure 9 shows areas of good (φs,G)-parameter space.

A relationship diagram between φ s and the optimal ranges of G
In this study, we introduced an optimized virtual force algorithm based on the exchange force, in which a new shielding rule based on a Delaunay triangulation which accurately defines adjacent nodes was adopted. We also used a new performance metric called ‘pair-correlation diversion’ to evaluate the uniformity of the distribution. Simulation results verified the effectiveness of the algorithm. We also provided a detailed discussion of the parameter effect and obtained an optimal region for the (φ s ,G)-parameter space. Our results have implications for choosing input parameters in actual experiments and also for the analysis of sensor data in robotics systems. In our future work, we will attempt to enlarge the shielding range and compare the shielding effects in different physical systems (e.g., a dusty plasma system).
Footnotes
6.
Our work was supported by the National Science Foundation of China (Grant No. 41464007), the Science Foundation of Jiangxi Province (Grant No. 2014BAB202008), the Specialized Research Fund for the Doctoral Program of Higher Education (Grant No. 20133601120005); and the Educational Foundation of Jiangxi Province (GJJ13050).
