Testing a Firefly-Inspired Synchronization Algorithm in a Complex Wireless Sensor Network

Data acquisition is the foundation of soft sensor and data fusion. Distributed data acquisition and its synchronization are the important technologies to ensure the accuracy of soft sensors. As a research topic in bionic science, the firefly-inspired algorithm has attracted widespread attention as a new synchronization method. Aiming at reducing the design difficulty of firefly-inspired synchronization algorithms for Wireless Sensor Networks (WSNs) with complex topologies, this paper presents a firefly-inspired synchronization algorithm based on a multiscale discrete phase model that can optimize the performance tradeoff between the network scalability and synchronization capability in a complex wireless sensor network. The synchronization process can be regarded as a Markov state transition, which ensures the stability of this algorithm. Compared with the Miroll and Steven model and Reachback Firefly Algorithm, the proposed algorithm obtains better stability and performance. Finally, its practicality has been experimentally confirmed using 30 nodes in a real multi-hop topology with low quality links.


Introduction
Data acquisition is the foundation of soft sensor and data fusion. In industrial processes, data often comes from multiple sources, which may be separated physically or virtually, and need to be synchronized. Therefore, distributed data acquisition and its synchronization are the important technologies to ensure the accuracy of soft sensors [1,2]. At present, Wireless Sensor Networks (WSNs), as an important distributed data acquisition technology, have been widely used in industry [3][4][5]. Synchronization is a key feature of WSNs, and it has been used in many wireless protocols. In general, the synchronization algorithm can be classified according to the centrality. On one hand, the centralized algorithm (e.g., Reference Broadcasts Synchronization [6], Timing-sync Protocol [7], and Flooding Time Synchronization Protocol [8]) sets a base time for the coordinator and the nodes modify their time correspondingly, thus the accuracy of this algorithm is high for structures with simple topology [9]. However, if the topological structure is complex, the centrality algorithm is inapplicable. On the other hand, compared with the centralized algorithm for synchronization, the distributed algorithm, which avoids single points of failure and is more robust, is more suitable for large-scale complex topological structures [10].
The distributed algorithm can be implemented with the classical packet-coupling synchronization or the firefly-inspired synchronization. The packet-coupling synchronization has accuracies ranging from microseconds to milliseconds [11]. However, a large number of packets need to be exchanged, which leads to increased computational complexity, energy expenditure, and poor scalability. Development of synchronization technology has led to the firefly-inspired synchronization [12][13][14].

1.
By using a discrete state vector instead of a continuous phase, the algorithm can be easily implemented on a Microcontroller Unit (MCU) in a wireless module. 2.
The step size can be adjusted, which ensures the stability according to the global convergence property of Markov chains without adjusting the coupling strength.

3.
A better performance emerges from the stochastic coupling at multiscale quantitative levels and the channel congestion is alleviated significantly.
The rest of this paper is organized as follows: in Section 2, the discrete phase dynamics at multiscale quantitative levels is proposed. Section 3 presents the design of a stochastic coupling algorithm, including the compensation for frequency drift and communication delay. The stability of this algorithm is verified in Section 4. Simulation and experiment results are discussed in Sections 5 and 6. Finally, conclusions and future research directions are presented in Section 7.

Discrete Phase Dynamics
In this section, we generalize the integrating and coupling dynamics. According to the M&S phase model [26], the phase is discretized at a single-scale quantitative level, making it is easier to implement on a WSN node. The single-scale quantitative level is expanded to multiscale quantitative levels for better performance.

Discretization of a Single-Scale Quantitative Level
In the discrete phase, each node is characterized by a cycle-count timer k (the period is T). Assume that ∆t is a short duration of time which denotes the resolution of timer and T can be divided exactly by ∆t. The timer k increases by one for each ∆t. Thus, in each node, the local timer k can be expressed as a discrete phase, which consists of a series of moments t and: where Z denotes the set of integers, rather than the continuous phase contained in the M&S model. The dynamics of each node are composed of the integrating and coupling dynamics.

Integrating Dynamics
The integrating dynamics are determined by the properties of the node, which does not consider the network coupling. The node's phase would not be influenced by other nodes in integrating dynamics. The integrating dynamics represent the process of increasing the discrete phase k independently in each node as a counter, shown in Equation (2): As time increases from t n to t n+1 with ∆t, the phase k becomes k + 1. When the phase reaches its phase threshold T/∆t, k(t n+1 ) will return to zero immediately. From Equation (2), the period of the phase is T and the update interval is ∆t (which denotes the quantization resolution). Here, the quantitative level is defined as the threshold of the timer, as shown in Equation (3): The timer k may be any natural number less than the quantitative level, and T is the threshold time corresponding to the phase threshold.

Coupling Dynamics
The coupling dynamics adjust the phase k of a node coupled with another node in the WSN. For the purposes of the following description, suppose that node j receives the synchronization message of node i. i−j Diff ( i k(t), j k(t)) represents the phase difference between i and j at moment t, where the phase in i and j is i k(t) and j k(t), respectively. The phase difference can be calculated from Equation (4), and k can be obtained by processing the synchronization packet. A further explanation is provided in Sections 3.2 and 3.4.
From Equation (4), i−j Diff is an integer and meets −T/2∆t ≤ i−j Diff ≤ T/2∆t . Considering the discretization, the phase is non-derivable and the adjustment of the phase must be divided exactly by ∆t. The adjustment of the discrete phase is determined by a step-increasing function for the discrete phase difference in order to retain the concave property of the phase mapping function [26]. In order to avoid cumbersome floating-point operations, the adjustment is determined according to a step-increasing function f ( i−j Diff ) as shown in Equation (5): In Equation (5), r denotes the refractory period which usually meets r < 0.4T/∆t [28]. Here, the refractory period is shown in terms of the phase difference. From Equation (5), the coupling dynamics transform the continuous phase adjustment into a phase state transition, similar to the state transition in a Markov chain. The details are shown in Section 4.1. The conversion conditions are shown in Equation (5), which are key to achieving synchronization. The stability is ensured by the property of Markov chains, which will be explained in Section 4. Finally, the coupling dynamics are presented in Equation (6): Figure 1 shows an example. In this figure, i−j Diff ( i k(0), j k(0)) = 2 while the refractory period is 1, and node j is coupling with node i. Thus, f ( i−j Diff ) = 1. The coupling occurs when the node j reaches its phase coupling (t = T). These details will be explained in Section 3.
The coupling dynamics adjust the phase k of a node coupled with another node in the WSN. For the purposes of the following description, suppose that node j receives the synchronization message of node i. i−jDiff(ik(t), jk(t)) represents the phase difference between i and j at moment t, where the phase in i and j is ik(t) and jk(t), respectively. The phase difference can be calculated from Equation (4), and k can be obtained by processing the synchronization packet. A further explanation is provided in Sections 3.2 and 3.4.
From Equation (4), i−jDiff is an integer and meets −T/2Δt ≤ i−jDiff ≤ T/2Δt. Considering the discretization, the phase is non-derivable and the adjustment of the phase must be divided exactly by Δt. The adjustment of the discrete phase is determined by a step-increasing function for the discrete phase difference in order to retain the concave property of the phase mapping function [26]. In order to avoid cumbersome floating-point operations, the adjustment is determined according to a step-increasing function f(i−jDiff ) as shown in Equation (5): In Equation (5), r denotes the refractory period which usually meets r < 0.4T/Δt [28]. Here, the refractory period is shown in terms of the phase difference. From Equation (5), the coupling dynamics transform the continuous phase adjustment into a phase state transition, similar to the state transition in a Markov chain. The details are shown in Section 4.1. The conversion conditions are shown in Equation (5), which are key to achieving synchronization. The stability is ensured by the property of Markov chains, which will be explained in Section 4. Finally, the coupling dynamics are presented in Equation (6): Figure 1 shows an example. In this figure, i−jDiff(ik(0),jk(0)) = 2 while the refractory period is 1, and node j is coupling with node i. Thus, f(i−jDiff ) = 1. The coupling occurs when the node j reaches its phase coupling (t = T). These details will be explained in Section 3.
Send the synchronization packet at a random moment from node i.
Couple with node i.

Discretization of Multiscale Quantitative Levels
Discretization divides the period time into pieces. Quantitative level denotes the amount of pieces. According to [32], a larger quantitative level will lead to a higher synchronization accuracy. However, a large quantitative level will lead to a smaller ∆t and a larger number of state elements, which may lead to partial synchronization and slow the convergence. Therefore, multiscale quantization is adopted in this section. Utilizing multiscale quantitative levels, the synchronization process is accelerated by using large quantitative levels, and the accuracy is improved by using small quantitative levels.
Consider the quantization process with different levels. Each level has a different quantization resolution. After quantization, the discrete phase value is converted to a multiscale phase space. In order to decouple each level, we sort them hierarchically by quantitative level and define the threshold time T (T can be modified as needed) in one layer as the quantization resolution in next lower layer, as shown in Equations (7) and (8): where l denotes the layer index of quantitative level; m is the amount of layers and k lmax represents the quantitative level of layer l. Thus, when layer l is adjusted, the phase in other layers will not be affected unless the threshold is reached. The value of period and quantization resolution in each layer can be selected as needed. The smallest ∆t determines the time resolution of the node, which is usually valued according to the accuracy requirement and hardware limitations. A smaller ∆t may get a higher time resolution and accuracy but require a more tedious calculation. Here we use the smallest ∆t to form a resolution vector RES, as shown in Equation (9): Because the time resolution at the highest layer of quantitative level (the mth layer) is the smallest ∆t. Thus RES m = 1. In multiscale dynamics, the real phase value can be expressed by the discrete phase value in each layer. Combining the different phase values of m layers forms the state vector i K(t) of the node i, as shown in Equation (10): where i k l (t) represents the phase in i at the lth layer at the moment t.
The real normalized phase in node i can be rewritten as i K in Equation (11): In each layer, the integrating and coupling dynamics are the same as the dynamics of the single-scale quantitative level. According to Equation (2), the integrating dynamics at the lth layer are shown in Equation (12).
As with the single-scale dynamics, it is necessary to define the multiscale phase differences in a matrix before determining the multiscale coupling dynamics, shown in Equation (14): where i−j Diff l denotes the phase difference at the lth layer of quantitative levels between node i and node j.
The step-increasing function can be determined by mapping the phase difference matrix into the adjustment step vector as shown in Equation (15): . . .
where f can be found by Equation (5). Thus, the coupling dynamics at multiscale quantitative levels are shown in Equation (16): Node j coupling with Node i ⇒ ∀j = i, j K(t + ) = K(t) + F( i−j Diff) (16) Sensors 2017, 17, 544 6 of 19 As with the single-scale dynamics, it is necessary to define the multiscale phase differences in a matrix before determining the multiscale coupling dynamics, shown in Equation (14): where i−jDiffl denotes the phase difference at the lth layer of quantitative levels between node i and node j.
The step-increasing function can be determined by mapping the phase difference matrix into the adjustment step vector as shown in Equation (15): where f can be found by Equation (5). Thus, the coupling dynamics at multiscale quantitative levels are shown in Equation (16):

Stochastic Coupling Synchronization Algorithm
Based on the dynamics presented in Section 2, a stochastic coupling synchronization algorithm is proposed in this section, as shown in Figure 2. It consists of five processing tasks: 1. Self-increasing the state vector. 2. Sending the synchronization packet at a random moment. 3. Delay and frequency drift compensation.

Stochastic Coupling Synchronization Algorithm
Based on the dynamics presented in Section 2, a stochastic coupling synchronization algorithm is proposed in this section, as shown in Figure 2. It consists of five processing tasks: 1.
Self-increasing the state vector.

2.
Sending the synchronization packet at a random moment.

3.
Delay and frequency drift compensation.

Self-Increasing the State Vector
According to the integrating dynamics, the state vector should be updated to meet Equation (12). Unlike the general continuous phase model, the phase growth in each layer is step-like. As described in Section 2, the discrete phase value should increase by one after time ∆t l in the lth layer. Thus, the state vector should be updated at each time point ∆t m , which is easy to achieve with a timer in the MCU.

Sending the Synchronization Packet at a Random Moment
The algorithm adopts a stochastic coupling method. The triggering synchronization packet includes the regular data format of the wireless network and the state vector information is added at the end of the MAC Payload. It is compatible with existing mainstream network formats (e.g., IEEE802.11 and IEEE802.15.4). Figure 3 illustrates the synchronous message structure of IEEE802.15.4 as an example. The phase information in this figure presents the state vector of the trigger node, and the coupling can be realized in each layer. Unlike the traditional triggering method, the synchronous message sending is carried out stochastically during the synchronization. In addition, a synchronous message can be combined with other network transmission packets to reduce the network load. In addition, we employ the Carrier Sense Multiple Access/Collision Avoidance (CSMA/CA) for the idle listening which is compatible with IEEE802.11 or IEEE802.15.4. The Maximum Transmission Unit (MTU) depends on the protocol and the application communication requirements. 4. Synchronization packet processing. 5. Reachback state vector adjustment.

Self-Increasing the State Vector
According to the integrating dynamics, the state vector should be updated to meet Equation (12). Unlike the general continuous phase model, the phase growth in each layer is step-like. As described in Section 2, the discrete phase value should increase by one after time Δtl in the lth layer. Thus, the state vector should be updated at each time point Δtm, which is easy to achieve with a timer in the MCU.

Sending the Synchronization Packet at a Random Moment
The algorithm adopts a stochastic coupling method. The triggering synchronization packet includes the regular data format of the wireless network and the state vector information is added at the end of the MAC Payload. It is compatible with existing mainstream network formats (e.g., IEEE802.11 and IEEE802.15.4). Figure 3 illustrates the synchronous message structure of IEEE802.15.4 as an example. The phase information in this figure presents the state vector of the trigger node, and the coupling can be realized in each layer. Unlike the traditional triggering method, the synchronous message sending is carried out stochastically during the synchronization. In addition, a synchronous message can be combined with other network transmission packets to reduce the network load. In addition, we employ the Carrier Sense Multiple Access/Collision Avoidance (CSMA/CA) for the idle listening which is compatible with IEEE802.11 or IEEE802.15.4. The Maximum Transmission Unit (MTU) depends on the protocol and the application communication requirements.

Delay and Frequency Drift Compensation
Due to the communication delay and in order to reduce the instability of the algorithm, it is necessary to compensate by optimizing and calibrating the delay. The total communication delay tdly satisfies the relationship in Equation (17): where tuc denotes the uncertain time delay, which includes the sending time, access time, and receiving time. The send time is that of the sender constructing the time message to transmit on the network. The access time is that of the MAC layer delay in accessing the network. Finally, the receive time is the receiving node processing the message and transferring it to the host. tuc can be nearly eliminated with the Start Frame Delimiter (SFD) timestamp technique [28]. In IEEE802.11, chips may not provide the SFD interrupt nor the SFD output pin. Instead, Transmit (TX) start interrupt and Receive (RX) end interrupt can be utilized. Then calculate the timestamp by the TX or RX timing. tc denotes the deterministic delay, which is the difference between the running and response times, which is caused by the interrupted service that is a function of sending and receiving SFD, and can be obtained by calibration. After obtaining the communication delay, it is necessary to convert the delay into the form of a state vector. The conversion method is shown in Equations (18) and (19):

Delay and Frequency Drift Compensation
Due to the communication delay and in order to reduce the instability of the algorithm, it is necessary to compensate by optimizing and calibrating the delay. The total communication delay t dly satisfies the relationship in Equation (17): where t uc denotes the uncertain time delay, which includes the sending time, access time, and receiving time. The send time is that of the sender constructing the time message to transmit on the network. The access time is that of the MAC layer delay in accessing the network. Finally, the receive time is the receiving node processing the message and transferring it to the host. t uc can be nearly eliminated with the Start Frame Delimiter (SFD) timestamp technique [28]. In IEEE802.11, chips may not provide the SFD interrupt nor the SFD output pin. Instead, Transmit (TX) start interrupt and Receive (RX) end interrupt can be utilized. Then calculate the timestamp by the TX or RX timing. t c denotes the deterministic delay, which is the difference between the running and response times, which is caused by the interrupted service that is a function of sending and receiving SFD, and can be obtained by  (18) and (19): Since the deterministic delay can be calculated in advance, the state vector of the delay can be calculated in advance, and there is no need to calculate it again while the process is running. The modified state vector meets the relationship in Equation (20). Thus, a more accurate state vector K can be obtained: A low-performance RC oscillator is often used as the clock driver in WSNs for cost savings, which introduces a large frequency drift error. To reduce this error and improve the synchronization accuracy, this algorithm employs the correction algorithm presented in [25]. First, the real frequency f RC of the RC oscillator is calibrated and compared with the theoretical frequency f ideal . The adjusted minimum quantization resolution is obtained on the basis of the theoretical minimum quantization resolution ∆t m , shown in Equation (21):

Synchronization Packet Processing
In the synchronization process, the node receives the synchronization packet sent by the neighboring node and couples with each layer. During the processing, there is a buffer that only stores the difference of the state vector, which is the smallest in a period of the refractory interval. If we assume that node j has received a synchronization packet sent by node i, then the process is as follows: 1.
Node j should record the local state vector j K exactly when receiving the synchronization packet of node i and get the corrected state vector by compensation as described in Section 3.3.

2.
Analyze the synchronization packet to obtain node i's state vector i K.

3.
According to Equation (14) calculate the difference i−j Diff (iK,jK). Due to the discretization of phase, when the phase difference equal to 1 in a certain quantitative layer, the real continuous difference value should be less than the time resolution in this layer. So the difference should be reflected in the former layer which has a more precise time resolution. Thus, in order to prevent an oscillatory response between two adjacent layers when the synchronization is almost achieved, it should be equivalently converted to the former layer as shown in Equation (22): 4. If i−j Diff (iK,jK) in Step 3 is out of the refractory interval, update the buffer i−j Diff buf by the smaller difference between the one just obtained and the one stored in the buffer. Thus, after a period, the buffer only stores the difference of state vector that is the smallest in the period. That is to satisfy Equation (23): where RES is the resolution vector in Equation (8) and r is the length of the refractory interval in [28]. The phase difference in the buffer is used to provide input for the reachback state vector adjustment.

Reachback State Vector Adjustment
As the node's hardware constraints and receiving and adjusting the phase right after receiving the packet may cause hysteresis and instability, the adjustment uses the reachback mechanism of RFA [20]. Under that mechanism, after the longest time threshold, the node obtains the state vector difference in the buffer and calculates F( i−j Diff buf ) as the starting state vector of the next cycle, shown in Equation (24). This can reduce the amount of computing and avoid the impact of disordering reception:

Stability
As the adjustment in one layer does not affect the other layers, the stability of the target algorithm can be determined by verifying the stability in a single layer. Therefore, this section focuses on proving the stability of the model with a single-scale quantitative level. The stability of the two-node network is demonstrated, and the conclusions are extended to a multi-node network.

State Space of the System
As shown in Section 3, all possible values of the discrete phase in node i can be gathered in a set as shown in Equation (25): To explain this process clearly, suppose T is an even multiple of ∆t. Therefore, i X is a finite state space of the node itself and each node can be regarded as a finite state machine (FSM), which constitutes a Markov chain.
Similarly, for a two-node network, the absolute value of i−j Diff can be selected as the state of the network. All possible values can be gathered in a set as shown in Equation (26): Hence, the system state can be modeled by a discrete-time Markov chain. Synchronization is achieved when the discrete phase difference between the two nodes decreases to zero.

One-Step Transition Probability Matrix
For nodes i and j, suppose they satisfy i k > j k and: Consider the next adjustment of the phase in a two-node network. From Equation (5), we can see that the adjustment step is −f ( i−j Diff ), which always has an opposite sign of the input i−j Diff. Hence, the adjustment is a negative feedback and always reduces the absolute value of i−j Diff.
Therefore, if the synchronization has not yet been achieved, the probability of state transition from the adjustment l to (l + 1) can be found from Equations (28) and (29): For the two-node system that has a state space shown in Equation (30) and the Markov chain is shown as Figure 4:

One-Step Transition Probability Matrix
For nodes i and j, suppose they satisfy ik > jk and: Consider the next adjustment of the phase in a two-node network. From Equation (5), we can see that the adjustment step is −f(i−jDiff), which always has an opposite sign of the input i−jDiff. Hence, the adjustment is a negative feedback and always reduces the absolute value of i−jDiff.
Therefore, if the synchronization has not yet been achieved, the probability of state transition from the adjustment l to (l + 1) can be found from Equations (28) and (29): For the two-node system that has a state space shown in Equation (30) and the Markov chain is shown as Figure 4:

Proof of Stability
Suppose the state is n (n ≥ 1). After one step transition, the state will be n − 1 according to the Equation (31). And if the state is 0, after one step transition, the state will always be 0. Since the max of state is T/2∆t in Equation (30), after T/2∆t steps, the state will be 0 constantly and the probability of it will be 1.
From Equation (30), we know that this Markov chain is a finite irreducible chain. The synchronization state is an absorptive state, and has global reachability, producing Equation (32): The two-node system reaches a steady state after T/2∆t steps starting from any initial state. Because the steady state is an absorptive state, the algorithm is converged and stable.

Stability of the Multi-Node Network
The stability of the two-node network can be extended to verify the stability of a multi-node network. In a multi-node network, a specific node can be selected as a reference node, and the corresponding discrete phase differences of the other nodes can be calculated. These phase differences constitute a vector, and all the possible vectors constitute a state space of the multi-node system. The process of the algorithm is transformed into the state transition process of the Markov chain.
According to Equation (31), for any two adjacent nodes, the steady state is an absorbing state and has global reachability while the initial state is a non-absorbing state. As a multi-node network is composed of multiple pairs of two-node networks, when all nodes stochastically match each other, the steady state of the state space also has global reachability [31]. Through deeply connecting and limiting the number of nodes in the network, the composition of the Markov chain is a finite irreducible chain. Therefore, the convergence of the continuous-time Markov process in [31] is also applicable to the discrete-phase Markov chain model of the network. First, the state space in a Markov chain is represented by a set of state vectors S = {s 1 ,s 2 , . . . ,}. Each element in S represents a possible state vector of the system. Suppose the steady state in the Markov chain constitutes the set C. The remaining states are non-persistent states, constituting a set N. P represents the Markov chain transition probability matrix, and Equation (33) can be deduced: Then: For ∀x,y ∈ N, because N is a non-persistent state set: According to the theory of Markov chains, we know that there always exists a value of transfer times n that can realize P n > 0 and P n=1 p T , p ∈ C. According to Equation (36), since the steady state is the absorptive state, as the algorithm iterates, other nodes tend to be synchronized and finally reach stability.

Simulation Verification
In this section, we verify the stability of the synchronization algorithm proposed in Section 3 using Matlab. We divide the simulations into two groups. One group employs the typical RFA mechanism based on the M&S model and the other one employs the target algorithm to perform synchronization experiments in a non-fully connected mesh network. By recording the phase information of the nodes during the process, the influence of different algorithms on the synchronization can be analyzed.

Simulation Parameters
In order to enhance the contrast between the algorithms, there are several common simulation parameters in both groups. In the simulation experiment, 30 or 50 nodes are randomly arranged in a 1000 m × 1000 m area as shown in Figure 5. The communication distance of the node is 350 m, and the simulation time is 300 s with a 4 µs simulation step. The initial phase of the nodes is randomly distributed between 0 and T. The period of the cycle is 1.048576 s. The refractory period of the node is set to 16 µs. The communication delay and frequency drift of the node both influence the synchronization performance [11]. To simulate an extremely poor communication environment in spite of the compensation, we set a 50 ppm frequency drift and 100 µs communication delay. Additionally, a packet loss rate up to 20% is set which is uniformly distributed.

Simulation Verification
In this section, we verify the stability of the synchronization algorithm proposed in Section 3 using Matlab. We divide the simulations into two groups. One group employs the typical RFA mechanism based on the M&S model and the other one employs the target algorithm to perform synchronization experiments in a non-fully connected mesh network. By recording the phase information of the nodes during the process, the influence of different algorithms on the synchronization can be analyzed.

Simulation Parameters
In order to enhance the contrast between the algorithms, there are several common simulation parameters in both groups. In the simulation experiment, 30 or 50 nodes are randomly arranged in a 1000 m × 1000 m area as shown in Figure 5. The communication distance of the node is 350 m, and the simulation time is 300 s with a 4 μs simulation step. The initial phase of the nodes is randomly distributed between 0 and T. The period of the cycle is 1.048576 s. The refractory period of the node is set to 16 μs. The communication delay and frequency drift of the node both influence the synchronization performance [11]. To simulate an extremely poor communication environment in spite of the compensation, we set a 50 ppm frequency drift and 100 μs communication delay. Additionally, a packet loss rate up to 20% is set which is uniformly distributed. The state vector of each node is recorded at intervals of 1 s and the phase difference of any two points in the network can be calculated according to Equation (3). The maximum value of the difference is taken as the synchronization error. In the typical RFA synchronization mechanism, the coupling coefficient is set as 0.5, 0.01, and 0.005. In the multiscale phase dynamics, the number of multiscale layers is set to three, and the quantization levels are 64, 32, and 32. The minimum quantization resolution is 16 μs and the coupling meets Equation (16) while the coupling coefficient is calculated by Equation (15) dynamically. Figure 6 shows the initial phases where 20 nodes are in the non-synchronized state. The synchronization error in the 20-node network is shown in Figure 7. It can be seen that the initial phase of the node set is dispersed in both algorithms. Only the multiscale phase dynamics stimulation achieves phase convergence, and the convergence time is 44 s with an 80 μs error while the standard The state vector of each node is recorded at intervals of 1 s and the phase difference of any two points in the network can be calculated according to Equation (3). The maximum value of the difference is taken as the synchronization error. In the typical RFA synchronization mechanism, the coupling coefficient is set as 0.5, 0.01, and 0.005. In the multiscale phase dynamics, the number of multiscale layers is set to three, and the quantization levels are 64, 32, and 32. The minimum quantization resolution is 16 µs and the coupling meets Equation (16) while the coupling coefficient is calculated by Equation (15) dynamically. Figure 6 shows the initial phases where 20 nodes are in the non-synchronized state. The synchronization error in the 20-node network is shown in Figure 7. It can be seen that the initial phase of the node set is dispersed in both algorithms. Only the multiscale phase dynamics stimulation achieves phase convergence, and the convergence time is 44 s with an 80 µs error while the standard  Figure 8 shows the initial phase where 50 nodes are in the non-synchronized state. The synchronization error in the 50-node network is shown in Figure 9. We can see that the initial phase of the node set is dispersed in both algorithms. The RFA is only able to achieve phase convergence with a coupling strength of 0.01, and the target algorithm achieves phase convergence, while the rest are unstable. The RFA converges in a spread interval of 0.11 s while the target algorithm converges in a spread interval of 0.032 ms with the standard deviation well under 0.00876 ms. These results indicate that the synchronization algorithm based on the proposed mechanism is superior to the traditional RFA in both convergence speed and precision.  Figure 8 shows the initial phase where 50 nodes are in the non-synchronized state. The synchronization error in the 50-node network is shown in Figure 9. We can see that the initial phase of the node set is dispersed in both algorithms. The RFA is only able to achieve phase convergence with a coupling strength of 0.01, and the target algorithm achieves phase convergence, while the rest are unstable. The RFA converges in a spread interval of 0.11 s while the target algorithm converges in a spread interval of 0.032 ms with the standard deviation well under 0.00876 ms. These results indicate that the synchronization algorithm based on the proposed mechanism is superior to the traditional RFA in both convergence speed and precision.  Figure 8 shows the initial phase where 50 nodes are in the non-synchronized state. The synchronization error in the 50-node network is shown in Figure 9. We can see that the initial phase of the node set is dispersed in both algorithms. The RFA is only able to achieve phase convergence with a coupling strength of 0.01, and the target algorithm achieves phase convergence, while the rest are unstable. The RFA converges in a spread interval of 0.11 s while the target algorithm converges in a spread interval of 0.032 ms with the standard deviation well under 0.00876 ms. These results indicate that the synchronization algorithm based on the proposed mechanism is superior to the traditional RFA in both convergence speed and precision.

Comparison with RFA
The simulation results show that the target algorithm is stable in both networks, while the RFA is stable only in the 50-node network. Additionally, the target algorithm is 50% faster and much more accurate. In terms of synchronization performance, the RFA is difficult to stabilize in the 20-node network because of low network connectivity and needs an appropriate coupling coefficient for convergence in the 50-node network. In contrast, the target algorithm is stable and has a better performance under all conditions. The proposed algorithm is proven to have high stability, high convergence speed, and high synchronization precision compared with RFA. Additionally, the target algorithm has similar advantages as the firefly-inspired algorithm.

Comparison with Traditional Synchronization Protocols
The target algorithm represents a radically different method. Compared with the traditional synchronization protocols, all of the nodes behave equally in our target algorithm. Unlike the traditional synchronization protocols, there are no root nor reference nodes. Our target algorithm is not affected by the network layer. As a result, it is more robust to the single points of failure and the flexible links, which is one of the main attractions. In addition, it is more conducive to network

Comparison with RFA
The simulation results show that the target algorithm is stable in both networks, while the RFA is stable only in the 50-node network. Additionally, the target algorithm is 50% faster and much more accurate. In terms of synchronization performance, the RFA is difficult to stabilize in the 20-node network because of low network connectivity and needs an appropriate coupling coefficient for convergence in the 50-node network. In contrast, the target algorithm is stable and has a better performance under all conditions. The proposed algorithm is proven to have high stability, high convergence speed, and high synchronization precision compared with RFA. Additionally, the target algorithm has similar advantages as the firefly-inspired algorithm.

Comparison with Traditional Synchronization Protocols
The target algorithm represents a radically different method. Compared with the traditional synchronization protocols, all of the nodes behave equally in our target algorithm. Unlike the traditional synchronization protocols, there are no root nor reference nodes. Our target algorithm is not affected by the network layer. As a result, it is more robust to the single points of failure and the flexible links, which is one of the main attractions. In addition, it is more conducive to network

Comparison with RFA
The simulation results show that the target algorithm is stable in both networks, while the RFA is stable only in the 50-node network. Additionally, the target algorithm is 50% faster and much more accurate. In terms of synchronization performance, the RFA is difficult to stabilize in the 20-node network because of low network connectivity and needs an appropriate coupling coefficient for convergence in the 50-node network. In contrast, the target algorithm is stable and has a better performance under all conditions. The proposed algorithm is proven to have high stability, high convergence speed, and high synchronization precision compared with RFA. Additionally, the target algorithm has similar advantages as the firefly-inspired algorithm.

Comparison with Traditional Synchronization Protocols
The target algorithm represents a radically different method. Compared with the traditional synchronization protocols, all of the nodes behave equally in our target algorithm. Unlike the traditional synchronization protocols, there are no root nor reference nodes. Our target algorithm is not affected by the network layer. As a result, it is more robust to the single points of failure and the flexible links, which is one of the main attractions. In addition, it is more conducive to network expansion. Nevertheless, comparing with traditional synchronization protocols, there is still a certain gap in terms of accuracy. For instance, as shown in Section 5.2, the target algorithm achieves 32 µs which is better than RFA, but still slightly less than the reported 14 µs accuracy of Flooding Time Synchronization Protocol (FTSP) [9]. We believe that this accuracy can be increased by adjusting the quantization layer and time resolution in each layer and much work remains to be done. Currently, the accuracy of our target algorithm is sufficient for sleep management and time-division multiplexing.

Discussion on Power Consumption
Although the algorithm does not achieve significant improvement in energy consumption, it is still more energy efficient than RFA because of its higher convergence speed. However, as described in Section 2.1.2, each synchronization packet cannot correct the local clock completely but by a little step. Therefore, comparing with the traditional synchronization protocols, more synchronization packets and more time will be needed for reaching the convergence, which means more power consumption. Once the time synchronization is achieved, the synchronization mainly depends on the internal accuracy of the components or clock drifts of the clock system given by the manufacturing process. The traditional synchronization protocols are adaptive. Thus, in this case, they can schedule the synchronization packets on demand. Similarly, if the clock drift is not too large, our target algorithm can reduce the energy consumption by operating at a much lower transmission frequency and combining the synchronization packet with other network transmission packets. The transmission frequency depends on the amount of out-of-phase synchronization packets received during a period of time. So it will be competitive to the traditional synchronization protocols during the maintaining period in terms of the energy consumption. In the future work, we consider the algorithm optimization as follows: When join in the network, each node will initialize its phase by traditional synchronization algorithm, and use our algorithm to tighten the accuracy and maintain synchronicity. It will greatly speed up the synchronization and reduce the energy consumption before convergence.

Hardware and Experiment Design
In order to verify the practicality and validity of the synchronization algorithm, a high-precision synchronization acquisition platform is designed with 30 wireless nodes, an acquisition board, and an industrial computer. The algorithm is implemented on the MCUs (JN5148) in wireless nodes [33]. The phase synchronization emerges from the interaction of the nodes by wireless communication, while the phase value is acquired by a wired National Instruments (NI) synchronization acquisition platform. The unified time reference is provided by the synchronous acquisition to verify the effectiveness of the algorithm.
In this experiment, the integrating dynamics are calculated by the timer on the MAC layer in the RF module, and the synchronization period is 1.0486 s. The minimum quantization resolution is 16 µs. In order to maintain a certain degree of consistency with the simulation, as in the simulation, the number of multiscale layers is three, and the quantization levels are 64, 32, and 32. The refractory period is set to 50 µs. To form a non-fully connected network in the laboratory, it is necessary to restrict each node to have at most eight neighbor nodes with which they can communicate directly, and this limits the network connectivity up to a factor of four. Meanwhile, a packet loss rate up to 20% is set by the software. When a node finishes a synchronization period, its MCU will pull up the level of a GPIO port for 300 ms. The GPIO is sampled by NI's PXIe6358 30-channel digital signal acquisition [34], and the sampling frequency is 100 ksps. For visual observation, the IO port is connected to an LED. The GPIO ports are connected via wires to the digital input channel of the acquisition card while the synchronization is implemented wirelessly. Labview is employed to acquire the real phase value of the diagram of the acquisition program (the Virtual Instrument block diagram in Labview) is shown as Figure 10.

Experimental Results and Discussion
It can be seen in Figure 11 that each node achieves phase synchronization. The maximum error spread of synchronization at each sampling point in time is plotted in Figure 11a. Due to the multiscale phase mechanism, synchronization of the node is divided into three phases. The maximum quantization resolution leads to early rapid decrease of the synchronization error, and the minimum quantization resolution reduces the jitter at later times. Figure 11b shows that the final error spread is in a 50 μs interval with the standard deviation well under 15.42 μs. And the average oscillation period is 1.0485904 s. However, interferences in network may lead to a synchronization error jitter after synchronization, but due to the synchronization algorithm, can still be limited to a certain range. Figure 12 shows the temporal dynamics of 30 nodes in an unsynchronized and synchronized state. In this figure, a vertical bar represents a fire occurrence. Figure 12a shows the chaotic fires in the unsynchronized network. By comparison, Figure 12b has the neat and regular fires in a synchronized state. Figure 13 shows the node synchronization effect. In Figure 13a, the network is in an unsynchronized state. The state at 50 s is shown in Figure 13b, and it can be seen that the network has reached synchronization which requires less than 1600 packets. Therefore, the algorithm achieves synchronization with stability and high performance. The effectiveness of the discrete phase dynamics and algorithm were proven to have a high practical value.

Experimental Results and Discussion
It can be seen in Figure 11 that each node achieves phase synchronization. The maximum error spread of synchronization at each sampling point in time is plotted in Figure 11a. Due to the multiscale phase mechanism, synchronization of the node is divided into three phases. The maximum quantization resolution leads to early rapid decrease of the synchronization error, and the minimum quantization resolution reduces the jitter at later times. Figure 11b shows that the final error spread is in a 50 µs interval with the standard deviation well under 15.42 µs. And the average oscillation period is 1.0485904 s. However, interferences in network may lead to a synchronization error jitter after synchronization, but due to the synchronization algorithm, can still be limited to a certain range. Figure 12 shows the temporal dynamics of 30 nodes in an unsynchronized and synchronized state. In this figure, a vertical bar represents a fire occurrence. Figure 12a shows the chaotic fires in the unsynchronized network. By comparison, Figure 12b has the neat and regular fires in a synchronized state. Figure 13 shows the node synchronization effect. In Figure 13a, the network is in an unsynchronized state. The state at 50 s is shown in Figure 13b, and it can be seen that the network has reached synchronization which requires less than 1600 packets. Therefore, the algorithm achieves synchronization with stability and high performance. The effectiveness of the discrete phase dynamics and algorithm were proven to have a high practical value.

Conclusions and Extensions
In this paper, discrete phase dynamics at multiscale quantitative levels is proposed. The stochastic coupling algorithm is employed to achieve synchronization in a WSN with a complex topology. The stability of the algorithm is verified in this paper. The simulation results show that the proposed algorithm is more stable in both a 20-node and 50-node network, while RFA is stable only in a 50-node network with a coupling strength of 0.01. Finally, the target algorithm is implemented in a realistic WSN with a complex topology for verification. The performance of the proposed algorithm is only demonstrated preliminarily. It is necessary to conduct further studies to determine the effect of the parameters in the proposed algorithm on the performance and the energy consumption will be taken into account. Because this algorithm is a software-based algorithm, it can and will be transplanted and tested on the ATmega256RFR2 chip for supporting sleep management functions. Nevertheless, the proposed algorithm improves the stability of the firefly synchronization algorithm and reduces the design difficulty of firefly-inspired synchronization algorithms for WSNs with complex topologies.