Accurate Distributed Time Synchronization in Mobile Wireless Sensor Networks from Noisy Difference Measurements

# Accurate Distributed Time Synchronization in Mobile Wireless Sensor Networks from Noisy Difference Measurements

Chenda Liao, and Prabir Barooah,  This work has been supported by the National Science Foundation by Grants CNS-0931885 and ECCS-0955023.C. Liao is with Nuance Communications,Inc., Burlington, MA 01803 USA (e-mail: cdliao84@gmail.com, ph: 352 870 5251); he was formerly with the Dept. of Mechanical and Aerospace Engineering, University of Florida.P. Barooah is with the Department of Mechanical and Aerospace Engineering, University of Florida, Gainesville, FL 32611 USA (e-mail: pbarooah@ufl.edu, ph: 352 392 0614)
###### Abstract

We propose a distributed algorithm for time synchronization in mobile wireless sensor networks. Each node can employ the algorithm to estimate the global time based on its local clock time. The problem of time synchronization is formulated as nodes estimating their skews and offsets from noisy difference measurements of offsets and logarithm of skews; the measurements acquired by time-stamped message exchanges between neighbors. A distributed stochastic approximation based algorithm is proposed to ensure that the estimation error is mean square convergent (variance converging to 0) under certain conditions. A sequence of scheduled update instants is used to meet the requirement of decreasing time-varying gains that need to be synchronized across nodes with unsynchronized clocks. Moreover, a modification on the algorithm is also presented to improve the initial convergence speed. Simulations indicate that highly accurate global time estimates can be achieved with the proposed algorithm for long time durations, while the errors in competing algorithms increase over time.

Time synchronization, Stochastic approximation, Wireless sensor networks, Difference measurement, mobile ad-hoc networks, distributed estimation.

## I Introduction

Time synchronization is critical for the functionality and performance of wireless sensor networks. For example, in TDMA-based communication schemes, accurate time synchronization ensures that each node communicates with others in the correct time slots without interfering with others. Furthermore, operation on a pre-scheduled sleep-wake cycle for energy conservation also requires a common notion of time among nodes. However, clocks in sensor nodes run at different speeds due to the imperfectness of quartz crystal oscillators, and a tiny difference on the oscillators of two clocks will cause time readings drift apart over time.

In the common clock model, the local time of node , , as a function of the global time is , is defined as:

 τu(t)=αut+βu, (1)

where the scalar parameters are called its skew (the relative speed of the clock with respect to ) and offset (the time reading when ), respectively [1]. In practice, skews are time-varying due to temperature change, aging etc. However, it is common to model the skew of a clock as a constant since its variation is negligible during time intervals of interest [2]. A node that knows the global time is called a reference node. The global time can be Coordinated Universal Time (UTC) if the reference node(s) can access it through GPS. If none of the nodes can access the UTC, one node in the network is elected as a reference node so that the local time in the node becomes the global time. A node can determine the global time from its local time if it knows its skew and offset. Specifically, if has estimates of its true skew and offset , it can estimate the absolute time as

 ^tu=τu(t)−^βu^αu (2)

Hence the problem of the time synchronization can be alternatively posed as the problem of skews and offsets estimation among nodes. A node can use a pairwise synchronization method, such as those in [3, 4, 5], to estimate and if the paired neighboring node is a reference node. However, most nodes in sensor networks are not connected to the reference nodes directly due to the limited communication range. It is therefore not possible for all nodes to obtain their skews and offsets directly.

Network-wide time synchronization in sensor networks has been intensely studied in recent years. Work in this area can be grouped into three categories: cluster-based protocols, tree-based protocols and distributed protocols. In cluster-based [6] and tree-based protocols [7, 8], synchronization relies on establishing a pre-specified network structure. In mobile networks, however, network topology continually changes, which results in frequent re-computation of a cluster and spanning tree, or re-election of a root node. This introduces considerable communication overhead to the networks, therefore the above cluster-based and tree-based protocols are primarily targeted to networks of static or quasi-static nodes.

Recently, a number of fully distributed algorithms that do not require the establishment of clusters or trees have been proposed. These typically perform synchronization by estimating skews and/or offsets and then computing the global time from them. The algorithms proposed in [9, 10, 11, 12] belong to this category. In these algorithms, estimates of a log-skews (and offsets) are obtained from noisy measurements of the difference of log-skews (and offsets) between pairs of neighbors. These distributed algorithms for skew and offset estimation are more readily applicable to mobile networks than the previous two categories of algorithms, though their convergence analyses were provided only for static networks. Convergence of these algorithms in mobile networks is analyzed in [13]. The algorithm analyzed in [13] is applicable to mobile networks and it subsumes the algorithms in [9, 10, 11, 12]. The algorithm in [13] is called JaT algorithm (Jacobi-type) due to its similarity to the Jacobi algorithm first proposed in [9]. The time-varying topology of a network of mobile nodes is modeled as the state of a Markov chain. Under certain conditions, it was shown that the variances of estimation errors of log-skews and offsets converge to positive values. However, even a small error in the skew estimate leads to poor absolute time estimate over long time periods, cf. (2). Thus, even a small steady variance of the skew estimates may lead to poor time estimates over time, requiring frequent restarting of the synchronization process.

In this paper, we revisit the problem of distributed estimation of clock skews and offsets from noisy difference measurements. The main contribution is an algorithm (called DiSync) that achieves steady-state variance of the skews and offsets under mild assumptions on the pairwise measurement noise. Mean square convergence of the algorithm is proved for both random (Markovian) as well as deterministic switching of graphs. Time varying gains in the proposed algorithm that make the variances converge to are adopted from stochastic approximation, which is also used in [14, 15] to attenuate noise. The gains need to vary in a specific manner with time, which poses a challenge in implementation in a network of unsynchronized clocks. We address this issue by using a novel approach: an iteration schedule is pre-specified to the nodes so that they can effectively perform a synchronous update without having synchronized clocks. This makes DiSync fully distributed and asynchronous. Furthermore, we propose a DiSync-I algorithm in which the effect of slow convergence rate of the DiSync is ameliorated while retaining theoretical convergence guarantees. We evaluate the accuracy of the DiSync and DiSync-I algorithms when applying to global time estimation through Monte Carlo simulations. Time estimation accuracy of the proposed algorithms are compared with that of the JaT algorithm. Simulations indicate the global time estimation error in the proposed algorithms stay close to for long time intervals, while that in JaT increases over time.

The simulation studies in papers [9, 10, 11, 12, 16, 13] use noisy difference measurements generated by adding noise on true values, while in practice these difference measurements are supposed to be obtained from processing multiple time-stamps by using an existing pairwise synchronization protocol, such as the ones proposed in [3, 4, 5]. In contrast, here we generate noisy difference measurement by simulating the pairwise synchronization protocol of [3] with random delays in packet reception. The noise in the difference measurements obtained are thus likely to have more realistic characteristics.

A new type of virtual time-synchronization protocols has been proposed recently. They let nodes estimate a common virtual global time that is not the local time of any clock [17, 18]. These algorithms are potentially applicable to mobile networks, though not useful when knowing an absolute global time is critical. Still, we compare the proposed algorithm to the ATS algorithm proposed in [17]. Although ATS does not provide estimate of an absolute global time, we compare the algorithms in terms of the maximum synchronization error - the maximum deviation in the estimates of global time (virtual or absolute) over two arbitrary nodes. It turns out that the proposed DiSync and DiSync-I algorithm outperforms ATS under this metric.

The algorithm we are proposing bears a close resemblance to average-consensus (leaderless) algorithms in [17, 18]. The estimation error dynamics in our problem turns out to be a leader-follower consensus algorithm, where the leader states - corresponding to the estimation error of the reference nodes - are always 0. The convergence analysis in this paper is inspired from [17, 18]. There are some technical differences since our scenario is leader-follower consensus while those in [17, 18] are leaderless consensus.

A preliminary version of this paper has been published in [19]. Compared to [19] this paper contains several major extensions. First of all, the switching topology was assumed to be deterministic for the convergence analysis in [19], while in this paper we extend the analysis of convergence to Markovian switching. It has been shown in [20] that switching topology due to random node motion can be modeled as Markovian. The modified algorithm DiSync-I, which ameliorates the slow convergence rate of theDiSync algorithm, is another novel aspect of this paper compared to [19]. Moreover, practical implementation details of the algorithm, including extensive simulation comparisons with competing algorithms, is provided here which was lacking in [19].

## Ii Problem formulation

The time synchronization problem is formulated as nodes estimating their skews and offsets. It is possible for a pair of nodes , who can communicate with each other to estimate their relative skew and relative offset . The reason for this terminology is the following relationship , that can be derived from (1). The estimation of relative skews and offsets is called “pairwise synchronization”. Several protocols for pairwise synchronization from time-stamped messages are available; see [3, 4, 5] and references therein. We assume nodes can estimate relative skews and offsets by using one of these existing protocols. Only those nodes that can communicate directly with the reference nodes can estimate their (absolute) skews and offsets, since they can employ pairwise synchronization with reference nodes. Most of the nodes cannot estimate their skews and offsets due to limited communication range.

Suppose between a pair and , node obtains noisy estimates of the parameters by using a pairwise synchronization protocol. We model the noisy estimate as , where is the estimation error. Therefore, by ,

 log^αu,v =logαu−logαv+ξsu,v, (3)

where . The quantity obtained from pairwise synchronization is therefore a noisy difference measurement of log-skews. If , which is usually the case, and is small, then the measurement noise is small. Similarly, the noisy estimate of relative offset is modeled as , where is the error. Again, by ,

 ^βu,v =βu−βv+ξou,v, (4)

which is a noisy difference measurement of the offsets between the two nodes, with measurement noise . Due to the nonzero , the measurement error is biased even if is zero mean. Since is close to for most clocks, the bias is usually small.

We see from (3) and (4) that and are the noisy measurements of log-skew difference and offset difference , respectively. We now seek to estimate the log-skews and offsets of all the nodes in a distributed manner from these noisy pairwise difference measurements. Note that once a node estimates its log-skew, it can recover the skew, and then compute the global time from its local time using estimated skew and offset.

To facilitate further discussion, in this section we only consider the estimation of scalar valued node variables from noisy difference measurements. If an algorithm of solving this problem is available, two copies of the algorithm can be executed in parallel to obtain both log-skews and offsets. Let -th node in a -node network have an associated constant scalar node variable , . Nodes in do not know their node variables, while the reference nodes are the remaining nodes in , who know the values of their own node variables. Here represents for skew estimation and for offset estimation. Without loss of generality, we assume node variables of reference nodes are all 0, i.e. skews are 1 and offsets are 0. Time is measured by a discrete time-index . The mobile nodes define a time-varying undirected measurement graph , where if and only if and can obtain a difference measurement of the form

 ζu,v(k)=xu−xv+ξu,v(k), (5)

during the time interval between and , where is measurement error. We assume that between and , whoever obtains the measurement first shares it with the other so that it is available to both and . We also follow the convention that the difference measurement between and that is obtained by the node is always of while that used by is always of . Since the same measurement is shared by a pair of neighboring nodes, if receives the measurement from , then it converts the measurement to by assigning . For similar reasons, between a pair and , the node who computes in node pair and is fixed for all time . This can be achieved by comparing the magnitude of the index of nodes. For example, if , then computes first and then sends it to . The neighbors of at , denoted by , is the set of nodes that has an edge with in the measurement graph . We assume that if , then and can also exchange information through wireless communication at time .

Now the reformulated problem is to estimate the node variables , , by using the difference measurements that become available over time. We assume (i.e., there exists a least one reference node), otherwise the problem is indeterminate up to a constant.

## Iii The DiSync algorithm

We present an iterative algorithm that nodes can use to solve the problem of node variable estimation from noisy difference measurements in a distributed manner. Since nodes do not have synchronized clocks, iterative updates have to be performed asynchronously. Each node keeps its local iteration index and maintains an estimate of its node variable in its local memory. The estimates can be initialized to arbitrary values. In executing the algorithm, node starts its -th iteration at a pre-specified local time , for , which will be described in Section III-A. Then, node obtains current estimates along with the measurements from its current neighbors . After a fixed time length (measured in local time), node increments its local iteration index by 1, and updates its new estimate based on current measurements and neighbors’ estimates by using the following update law:

 ^xu(ku+1)= ^xu(ku)+m(ku)∑v∈Nu(ku)auv(ku)(^xv(ku) +ζu,v(ku)−^xu(ku)), (6)

where the time varying gain has to be specified to all nodes a-priori. Note that when , . The choice of will play a crucial role in the convergence of the algorithm and will be described in Section IV. In this paper, we let weight if . Recall that the reference nodes take part by helping their neighbors obtain difference measurements, but they do not update their own node variables. The algorithm is summarized in Algorithm 1. Note that, since obtaining difference measurements requires exchanging time-stamped messages, current estimates can be easily exchanged during the process of obtaining new measurements.

### Iii-a Iteration schedule and synchronous view

We will later describe that the gains is chosen to be a decreasing function of time, which helps reduce the effect of measurement noise. This is a well-known idea in stochastic approximation. However, using this idea in a network of unsynchronized clocks presents an unique challenge since no node has a notion of a common time index, at least in the initial phase when they do not have good estimates. If nodes waits for a constant length of time (measured in their local clocks) before starting a new iteration, a node with faster skew might finish the -th iteration while a node with slower skew hasn’t even finished the -th iteration. Therefore, specifying a function to all the nodes does not ensure that nodes use the same gain at the same (global) interval, which is required by stochastic approximation.

We address the problem by providing the nodes a priori the sequence of local time instants , mentioned at the beginning of the Section III. This sequence is called an iteration schedule, and the formula for computing it is described below. Let the skews and offsets of all clocks be lower and upper bounded by those in two fictitious clocks and , such that , . Recalling (2), therefore for all . The formula for calculating is

 τ(i+1)=αcHαcL(τ(i)+δt−βcL)+βcH, (7)

where has to satisfy . This schedule ensures that nodes operating on their unsynchronized local clocks still perform updates in an effectively synchronous manner over time. To see this, define as a global interval and as the global time interval with respect to -th local iteration of node . Eq. (7) guarantees that, at each , for all . In other words, there exists a sequence of global time intervals such that each -th global interval contains, and only contains, the -th local iteration (in global time) of all . Figure 1 shows the relationship between intervals of local iterations and the corresponding global intervals. In Figure 1, we pick the 3rd global interval from Figure 1, and show the global time intervals when local iteration updates occur. We emphasize that is the same for all nodes and every node starts and ends its -th iteration at the same local time instants and . Each node is provided the values of the parameters , ,, and ahead of time, which are design variables.

We next address the issue of how to pick values for , and without knowing a real bounds on skews and offsets of all clocks in a network. In wireless sensor nodes, a pair of clocks in sensor nodes usually drift apart up to  [2]. Therefore, we can pick . To pick reasonable values of the offset bounds, the following procedure should be used to initialize the synchronization. The reference node first broadcasts a message (to indicate the beginning of synchronization) and sets its local clock time to zero simultaneously. A node that receives this message sets its own clock to zero and broadcast such message again. The nodes that hear this message also set their local clocks to , and so forth. Since nodes (except the reference node) start their local clocks after – but close to – the instant of , their offsets are negative and small. Therefore, can be chosen as zero and can be picked as an estimate of the time it takes for all active nodes to receive the “synchronization start” signal. For a node who was out of communication range at the beginning but joins the networks later, it can set the local time to the current local time of a neighbor that has already started the time synchronization, and record neighbor’s iteration index as well. In this way, the newly joined node can take part in the synchronization process as if it started at the very beginning.

Another practical issue is the unbounded growth of the inter-synchronization intervals . For example, if is chosen as second, the choice of ensures that the time interval between two successive iterations, , will increase from second to seconds after iterations, or hours. If the updates are to be done more slowly, a larger will be used, which will slow down the growth of the iteration interval . If it is desired that the inter-synchronization interval does not increase beyond a certain pre-specified value, a reference node to restart the synchronization after a certain time to maintain within the desired bound. When the restart should occur can be computed from the recursion (7).

###### Remark 1.

Nodes perform skew and offset estimation simultaneously in a distributed and iterative fashion by using two copies of the DiSync algorithm, one for skew estimation and one for offset estimation. With current estimated skew and offset , a node can compute the global time using (2), i.e. , which is the final step of the time synchronization. The entire time synchronization procedure is illustrated in Figure 2.

## Iv Convergence analysis of DiSync algorithm

Since there exists a common iteration counter that can be used to describe the local updates in the nodes by using the iteration schedule (even though none of the nodes has access to it), we consider only the synchronous version of the algorithm using global index from now on. We rewrite (III) as

 ^xu(k+1)= ^xu(k)+m(k)∑v∈Nu(k)auv(k)(^xv(k) +ζu,v(k)−^xu(k)). (8)

Defining the estimation error as , Eq. (IV) reduces to the following using (5):

 eu(k+1)= eu(k)+m(k)∑v∈Nu(k)auv(k)(ev(k) −eu(k)+ϵu,v(k)). (9)

We introduce the following stipulations and notations to pursue subsequent analysis. First,let for . Secondly, the Laplacian matrix of the graph is defined as if , and if . By removing the rows and columns of with respect to reference nodes, we get the principle submatrix (so called grounded or Dirichlet Laplacian matrix [21]). Let , the corresponding state space form of the estimation error is

 e(k+1)=(I−m(k)Lb(k))e(k)+m(k)D(k)ϵ(k), (10)

where

 ϵ(k):=[¯ϵ1(k)T,…,¯ϵnb(k)T]T, ¯ϵu(k):=[ϵu,1(k),…ϵu,n(k)]T, D(k):=diag(¯a1(k),…,¯anb(k)), ¯au(k):=[au,1(k),…,au,n(k)],

where and . When , is a pseudo random variable with the same mean and variance as the measurement noise on any existing edge. Moreover, as a node computes measurement and sends to , . Now, we introduce the following assumptions:

###### Assumption 1.

Measurement noise vector is with mean and bounded second moment, i.e. , where denotes 2-norm. Furthermore, and are independent for . In addition, is independent of , where .

In practice, may be time-varying even if all the underlying processes are wide sense stationary. For instance, the bias in offset difference measurement computed from node is different from that computed from node , as . Therefore, depends on which node initializes pairwise synchronization at time . To meet this requirement in Assumption 1, we can stipulate that the node who computes between a pair and is fixed for all time . This can be achieved by comparing the magnitude of the index of nodes. For example, if , then computes first and then sends it to . Indeed, the purpose of this requirement is to provide formula to compute the steady state value of estimation error, and the system may still achieve convergence without it.

###### Assumption 2.

The non-increasing positive sequence (step size of the stochastic approximation) is chosen as , where are constant real numbers.

Note that Assumption 2 is a special case of the standard requirement in stochastic approximation: and . The assumption is made to simplify the subsequent analysis, but we believe it is not necessary.

### Iv-a Deterministic Topology Switching

In this section, we analyze the convergence of (IV) when the topology switching is deterministic.

###### Assumption 3.

There exists s.t. for any , is connected, where is set of edges in .

###### Assumption 4.

The limits exist: , , .

Assumption 3 implies that information can go from any node to the rest nodes within a uniformly bounded length of time. Furthermore, as is bidirectional, another equivalent assumption is that contains a spanning tree. The proposed algorithm is also robust to permanently adding or deleting nodes in case the new resulting graph satisfies the assumption on connectivity. To understand the Assumption 4, we define the finite state space as the set of graphs that can occur over time. If the sequence of can be divided into a sequence of finite intervals , , such that the percentage of times that each state occurs is fixed in all except finitely many such intervals , then , and exist. Another example is that the state occurs according to a sample path of a stationary ergodic process. In addition, we denote sets of matrices and , where and correspond to . If the percentage of all states occurring is , then

 ¯Lb:=N∑i=1πiLbi, ¯D:=N∑i=1πiDi (11)
###### Theorem 1.

Under Assumption 1-4, the Algorithm 1 ensures that in (10) converges to in mean square, i.e., .

The theorem states that under the assumptions, the variance of the estimation error decays to . If additionally all the difference measurements are unbiased (), then the bias of the estimates converge to as well. Proof of the theorem uses the next two lemmas. The proof of the first lemma, which is inspired by [14, 15], can be found in [22], while the second is from[23].

###### Lemma 1.

If difference measurement is unbiased, i.e., , under assumption 1-3, the Algorithm 1 ensures that in (10) converges to in mean square, i.e., .

When , (10) can be regarded as a leader-following consensus problem with time-varying topology and zero-mean noisy input. The leaders are reference nodes , which hold their variables as zero. Then, for is driven to zero by the reference nodes as goes to in mean square sense.

###### Lemma 2.

[23]Denote by A an unknown symmetric and positive semi-definite matrix in , and we have to solve the equation Ax=y for an unknown . Assume that exists. We are given a sequence of matrices and a sequence , where . In addition, suppose that , , exists. Consider the sequence : is arbitrary,

 xk+1=xk+c1k+c2(yk−Akxk), (12)

where and are constant real numbers. Then, .

###### Proof of Theorem 1.

Taking expectations on both sides of (10) with respect to measurement noise , we obtain

 η(k+1)=(I−m(k)Lb(k))η(k)+m(k)D(k)γ, (13)

where . By substituting (13) in (10), we get

 ~e(k+1)=(I−m(k)Lb(k))~e(k)+m(k)D(k)ξ(k), (14)

where and . Note that is zero mean and satisfies Assumption 1. By Lemma 1, converges to 0 in mean square. Therefore, is mean square convergent to . Now, we examine the convergence of . From the definition of and the symmetry of , is a symmetric grounded Laplacian of . Since is connected, by Lemma 1 in [21], is positive definite. Consequently, . Now, it follows from Lemma 2 that . Consequently, converges to in mean square, i.e., .

### Iv-B Markovian Topology Switching

In this section, we analyze convergence when network topology switches randomly. We model the switching of the network topologies as a Markov chain; the reasonableness of this model for mobile networks has been established in [13].

###### Assumption 5.

The temporal evolution of the measurement graph is governed by an N-state homogeneous ergodic Markov chain with state space , which is the set of graphs that can occur over time. Furthermore is connected, where is set of edges in . In addition, the processes and are independent for all and .

In Assumption 5, the Markovian switch on the graphs means that where . The requirement for ergodicity of the Markov chain ensures that there is an unique steady state distribution with non-zero entries. This means every graph in the state space of the chain occurs infinitely often. Since is connected, ergodicity implies that information from the reference node(s) will flow to each of the nodes over time. Again, note that none of the graphs that ever occur is required to be a connected graph. Since the Markov chain is ergodic, the steady state distribution of the chain is defined as . Recalling that and correspond to , we can use the same formula as that in (11) to define and .

###### Theorem 2.

Under Assumption 1,2 and 5, in (10) is mean square convergent, i.e., , where is expectation of w.r.t. measurement noise , and converges to almost surely.

The proof of the theorem uses the next two lemmas.

###### Lemma 3.

If relative measurements are unbiased, i.e., , under 1,2 and 5, in (10) converges to in mean square, i.e., .

###### Lemma 4.

[Proposition 1 in [24]] Assume , , is stochastic process on , where is an symmetric positive semidefinite matrix and is . Consider the following iteration with arbitrary :

 xk+1=xk+c1k+c2(yk−Akxk), (15)

where and are real constants. If , , and is positive definite, then almost surely.

Proof of Lemma 3 is omitted due to its length; it can be found in [22]. Lemma 4 follows in a straightforward manner from the results in [24], so its formal proof is omitted.

###### Proof of Theorem 2.

The proof is similar to that for Theorem 1. Define , where the expectation is taken with respect to measurement noise . Take the expectation on both sides of (10), we get

 η(k+1)=(I−m(k)Lb(k))η(k)+m(k)D(k)γ. (16)

By substituting (16) in (10), we get

 ~e(k+1)=(I−m(k)Lb(k))~e(k)+m(k)D(k)ξ(k), (17)

where and . Note that is zero mean and satisfies Assumption 1. By Lemma 3, converges to 0 in mean square. Now, we examine the convergence of . We rewrite (16) as

 η(k+1)=η(k)+m(k)(D(k)γ−Lb(k)η(k)). (18)

It follows from Assumption 5,

 limk→∞1kk∑i=1E[Lb(k)]=N∑i=1πiLbi=¯Lb, limk→∞1kk∑i=1E[D(k)]=N∑i=1πiDi=¯D (19)

From the definition of and the symmetry of , is a symmetric grounded Laplacian of . Since is connected, by Lemma 1 in [21], is positive definite. Consequently, . Furthermore, is positive-semi definite. Then, by Lemma 4 that converges to almost surely.

## V Ameliorating slow convergence: DiSync-I

The proposed DiSync algorithm ensures mean square convergence by attenuating the measurement noise using ideas from stochastic approximation. In particular, the gain that decays slowly with time is instrumental in driving the variances of the estimation error to zero. This slowly decaying gain, however, also makes the convergence rate slow. We’ll see numerical evidence of this in Section VI. In practice, performance of the DiSync algorithm can be further improved by modifying the update law, which we describe next. The modified algorithm is called the DiSync-I(DiSync-Improved) algorithm.

First, we define a scalar state for each node , which is a surrogate for the distance (number of hops) of node from the reference nodes. The state is called average distance of from the reference nodes at time . Furthermore, we define the set , the subset of neighbors of node whose average distances are smaller than that of node at . The average distance is updated according to the iterative law described in Algorithm 2. In brief, is updated by averaging all from its neighbors who have smaller average distance. If none of its neighbors have smaller average distance than itself, increases. Both and are maintained in at index .

The node variable update law for the DiSync-I algorithm is given below; where the subscript on the local iteration counter is suppressed:

 ^xu(k+1)= ^xu(k)+h(k)∑v∈Hu(k)auv(k)(^xv(k) +ζu,v(k)−^xu(k)), (20)

where

 h(k) ={11+|Hu(k)|fork

where and are constants - that satisfy - that are pre-specified to all nodes.

If , it means that node has been closer to the reference nodes than node on average (even if node may be father from the reference node than at current index ). This indicates that node is likely to contain better estimates that . If and are neighbors, should use the estimate from to perform update, while should not use estimates from .

Note that when and , (V) becomes the JaT algorithm of [16, 13]. As shown in [16, 13], JaT algorithm ensures that the mean of the estimation error converges to a constant (zero if measurement is unbiased) and variance to a constant. In addition, when and , the resulting update law (V) is a modified version of JaT algorithm: it now uses Algorithm 2 during the initial phase up to . We will call it the JaT-I algorithm.

Table I shows how the update law (V) can produce different algorithms depending on the values of the parameters . For simulation studies in this paper on the DiSync-I algorithm, we pick somewhat arbitrarily. In this case, DiSync-I first adopts JaT-I when and then becomes DiSync when . The reason behind the modification (V) over (IV) is the following. First, JaT has better convergence speed during initial phase than that of DiSync. Second, JaT-I has even better convergence rate than JaT due to the use of only those neighbors that have been closer to the reference node(s).

Note that the DiSync-I differs from DiSync only during the initial phase (up to ), otherwise it is the same. As a result, the mean square convergence results of DiSync holds for DiSync-I as well.

## Vi Simulation evaluation

We now examine through simulations the time synchronization performance of the DiSync and DiSync-I algorithms, as well as that of JaTand JaT-I algorithms. Finally, they are compared with the virtual time synchronization algorithm ATS [17] in terms of pairwise synchronization error. Simulations are performed in a 10-node mobile network within a square field. Nodes’ motions are generated according to the widely used random waypoint (RWP) mobility model [25]. It has been shown in [13] that when nodes move according to the RWP mobility model, the switching of the graphs can be modeled as a Markov chain. A pair of nodes can communicate when distance between them is less than . The true skews and offsets of 9 nodes are picked uniformly from and respectively according to [8]. The single reference node (10th) has skew and offset . Denote as update intervals (also called synchronization periods), and as the global instant of the beginning of -th interval. In this simulation, the interval is chosen as 1 sec, therefore . For the sake of convenience, simulations are carried out in a synchronous fashion.

### Vi-a Implementation of pairwise synchronization

The simulation of the DiSync, DiSync-I, JaT and JaT-I algorithms requires pairs of nodes to obtain difference measurements by exchanging time stamped messages when they are within communication range. In order to evaluate the entire time synchronization procedure, unlike [10, 11, 12, 16], in which difference measurements are generated by adding random noise to true difference of log-skew/offset, we select the pairwise synchronization algorithm proposed in [3] to compute the relative skew and relative offset , and the difference measurements and are then obtained from these as described in Section II. According to [3], at the beginning of the th interval, node sends a message to that contains the value of the local time at when the message is sent: . When node receives this message, it records the local time of reception: . After a waiting period, node sends a message back to that contains both and . When it arrives at , node again records the local time of reception: . Two nodes and in communication range performs this procedure, called two-way time-stamped message exchange, twice - at the beginning and in the middle of each synchronization period. At the end of the synchronization period, node uses the obtained eight time stamps for to estimate and via the formula provided in [3]. Finally, node sends back to these estimates.

There is a random delay between the time a node sends a message and the other node receives the message. This delay directly induces errors in the estimated and , and thus determines the level of noise in the resulting difference measurement. Therefore, the random delay ultimately affects the time synchronization accuracy. In employing the pairwise synchronization protocol of [3], we subject the message exchanges to a random delay that is Gaussian distributed with mean and standard deviation , as these values are considered realistic for wireless sensors networks with current hardware and communication protocols with uncertain-delay elimination (e.g., MAC-layer time-stamping) [8]. Although the relation between the statistics of the random delays and the noise of difference measurements of log-skew and offsets defined in Section II are complex, the noise levels in the difference measurements used are likely to be realistic due to realistic choice of delays.

### Vi-B Performance in estimating global time

We conduct 1000 Monte Carlo simulations of running algorithms for 800 second (iterations). Figure 3 shows two snapshots of the network during a simulation. As we can see, only a limited number of nodes can communicate with each other. Recall that if for all . The step size function is chosen as . We pick somewhat arbi