Accurate Distributed Time Synchronization in Mobile Wireless Sensor Networks from Noisy Difference Measurements
aa r X i v : . [ c s . S Y ] D ec Accurate Distributed Time Synchronization inMobile Wireless Sensor Networks from NoisyDifference Measurements
Chenda Liao, and Prabir Barooah,
Member, IEEE
Abstract
We propose a distributed algorithm for time synchronization in mobile wireless sensor networks. Each node canemploy the algorithm to estimate the global time based on its local clock time. The problem of time synchronizationis formulated as nodes estimating their skews and offsets from noisy difference measurements of offsets and logarithmof skews; the measurements acquired by time-stamped message exchanges between neighbors. A distributed stochasticapproximation based algorithm is proposed to ensure that the estimation error is mean square convergent (varianceconverging to 0) under certain conditions. A sequence of scheduled update instants is used to meet the requirementof 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 thathighly accurate global time estimates can be achieved with the proposed algorithm for long time durations, while theerrors in competing algorithms increase over time.
Index Terms
Time synchronization, Stochastic approximation, Wireless sensor networks, Difference measurement, mobile ad-hoc networks, distributed estimation.
I. I
NTRODUCTION
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 withothers 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 sensornodes run at different speeds due to the imperfectness of quartz crystal oscillators, and a tiny difference on theoscillators of two clocks will cause time readings drift apart over time.
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: [email protected], ph: 352 870 5251); he was formerlywith 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)
In the common clock model, the local time of node u , τ u ( t ) , as a function of the global time is t , is defined as: τ u ( t ) = α u t + β u , (1)where the scalar parameters α u , β u are called its skew (the relative speed of the clock with respect to t ) and offset (the time reading when t = 0 ), 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 duringtime intervals of interest [2]. A node that knows the global time is called a reference node . The global time canbe Coordinated Universal Time (UTC) if the reference node(s) can access it through GPS. If none of the nodescan access the UTC, one node in the network is elected as a reference node so that the local time in the nodebecomes the global time. A node u can determine the global time t from its local time if it knows its skew andoffset. Specifically, if u has estimates ˆ α u , ˆ β u of its true skew and offset α u , β u , it can estimate the absolute time as ˆ t u = τ u ( t ) − ˆ β u ˆ α u (2)Hence the problem of the time synchronization can be alternatively posed as the problem of skews and offsetsestimation among nodes. A node u can use a pairwise synchronization method, such as those in [3–5], to estimate α u and β u if the paired neighboring node is a reference node. However, most nodes in sensor networks are notconnected to the reference nodes directly due to the limited communication range. It is therefore not possible forall nodes to obtain their skews and offsets directly.Network-wide time synchronization in sensor networks has been intensely studied in recent years. Work in thisarea can be grouped into three categories: cluster-based protocols, tree-based protocols and distributed protocols. Incluster-based [6] and tree-based protocols [7, 8], synchronization relies on establishing a pre-specified network struc-ture. In mobile networks, however, network topology continually changes, which results in frequent re-computationof a cluster and spanning tree, or re-election of a root node. This introduces considerable communication overheadto the networks, therefore the above cluster-based and tree-based protocols are primarily targeted to networks ofstatic or quasi-static nodes.Recently, a number of fully distributed algorithms that do not require the establishment of clusters or trees havebeen proposed. These typically perform synchronization by estimating skews and/or offsets and then computing theglobal time from them. The algorithms proposed in [9–12] belong to this category. In these algorithms, estimates ofa log-skews (and offsets) are obtained from noisy measurements of the difference of log-skews (and offsets) betweenpairs of neighbors. These distributed algorithms for skew and offset estimation are more readily applicable to mobilenetworks than the previous two categories of algorithms, though their convergence analyses were provided only forstatic networks. Convergence of these algorithms in mobile networks is analyzed in [13]. The algorithm analyzedin [13] is applicable to mobile networks and it subsumes the algorithms in [9–12]. The algorithm in [13] is calledJaT algorithm (Jacobi-type) due to its similarity to the Jacobi algorithm first proposed in [9]. The time-varyingtopology of a network of mobile nodes is modeled as the state of a Markov chain. Under certain conditions, itwas 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 frequentrestarting of the synchronization process.In this paper, we revisit the problem of distributed estimation of clock skews and offsets from noisy differencemeasurements. The main contribution is an algorithm (called DiSync) that achieves steady-state variance of theskews and offsets under mild assumptions on the pairwise measurement noise. Mean square convergence of thealgorithm is proved for both random (Markovian) as well as deterministic switching of graphs. Time varying gainsin the proposed algorithm that make the variances converge to are adopted from stochastic approximation, whichis also used in [14, 15] to attenuate noise. The gains need to vary in a specific manner with time, which posesa challenge in implementation in a network of unsynchronized clocks. We address this issue by using a novelapproach: an iteration schedule is pre-specified to the nodes so that they can effectively perform a synchronousupdate 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 whileretaining theoretical convergence guarantees. We evaluate the accuracy of the DiSync and DiSync-I algorithms whenapplying to global time estimation through Monte Carlo simulations. Time estimation accuracy of the proposedalgorithms are compared with that of the JaT algorithm. Simulations indicate the global time estimation error inthe proposed algorithms stay close to for long time intervals, while that in JaT increases over time.The simulation studies in papers [9–13, 16] use noisy difference measurements generated by adding noise ontrue values, while in practice these difference measurements are supposed to be obtained from processing multipletime-stamps by using an existing pairwise synchronization protocol, such as the ones proposed in [3–5]. In contrast,here we generate noisy difference measurement by simulating the pairwise synchronization protocol of [3] withrandom delays in packet reception. The noise in the difference measurements obtained are thus likely to have morerealistic characteristics.A new type of virtual time-synchronization protocols has been proposed recently. They let nodes estimate acommon virtual global time that is not the local time of any clock [17, 18]. These algorithms are potentiallyapplicable to mobile networks, though not useful when knowing an absolute global time is critical. Still, wecompare the proposed algorithm to the ATS algorithm proposed in [17]. Although ATS does not provide estimateof an absolute global time, we compare the algorithms in terms of the maximum synchronization error - themaximum deviation in the estimates of global time (virtual or absolute) over two arbitrary nodes. It turns out thatthe 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 theleader states - corresponding to the estimation error of the reference nodes - are always 0. The convergence analysisin this paper is inspired from [17, 18]. There are some technical differences since our scenario is leader-followerconsensus 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] thatswitching 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 comparedto [19]. Moreover, practical implementation details of the algorithm, including extensive simulation comparisonswith competing algorithms, is provided here which was lacking in [19].II. P ROBLEM FORMULATION
The time synchronization problem is formulated as nodes estimating their skews and offsets. It is possible for a pairof nodes u, v , who can communicate with each other to estimate their relative skew α u,v := α u α v and relative offset β u,v := β u − β v α u α v . The reason for this terminology is the following relationship τ u ( t ) = α u α v τ v ( t )+ β u − β v α u α v , thatcan be derived from (1). The estimation of relative skews and offsets is called “pairwise synchronization”. Severalprotocols for pairwise synchronization from time-stamped messages are available; see [3–5] and references therein.We assume nodes can estimate relative skews and offsets by using one of these existing protocols. Only thosenodes that can communicate directly with the reference nodes can estimate their (absolute) skews and offsets, sincethey can employ pairwise synchronization with reference nodes. Most of the nodes cannot estimate their skews andoffsets due to limited communication range.Suppose between a pair u and v , node u obtains noisy estimates ˆ α u,v , ˆ β u,v of the parameters α u,v , β u,v byusing a pairwise synchronization protocol. We model the noisy estimate as ˆ α u,v = α u,v + e su,v , where e su,v is theestimation error. Therefore, by α u,v = α u α v , log ˆ α u,v = log α u − log α v + ξ su,v , (3)where ξ su,v = log(1 + e su,v α v α u ) . The quantity obtained from pairwise synchronization is therefore a noisy differencemeasurement of log-skews. If α v /α u ≈ , which is usually the case, and e su,v is small, then the measurement noise ξ su,v is small. Similarly, the noisy estimate of relative offset is modeled as ˆ β u,v = β u,v + e ou,v , where e ou,v is theerror. Again, by β u,v = β u − β v α u α v , ˆ β u,v = β u − β v + ξ ou,v , (4)which is a noisy difference measurement of the offsets between the two nodes, with measurement noise ξ ou,v = β v (1 − α u α v ) + e ou,v . Due to the nonzero β v (1 − α u α v ) , the measurement error is biased even if e ou,v is zero mean.Since α u α v is close to for most clocks, the bias is usually small.We see from (3) and (4) that log ˆ α u,v and ˆ β u,v are the noisy measurements of log-skew difference log α u − log α v and offset difference β u − β v , respectively. We now seek to estimate the log-skews and offsets of all the nodesin a distributed manner from these noisy pairwise difference measurements. Note that once a node estimates itslog-skew, it can recover the skew, and then compute the global time from its local time using estimated skew andoffset. To facilitate further discussion, in this section we only consider the estimation of scalar valued node variables fromnoisy difference measurements. If an algorithm of solving this problem is available, two copies of the algorithm canbe executed in parallel to obtain both log-skews and offsets. Let u -th node in a n -node network have an associatedconstant scalar node variable x u ∈ R , u ∈ V = V b ∪ V r = { , . . . , n } . Nodes in V b = { , . . . , n b } do not knowtheir node variables, while the reference nodes are the remaining n r nodes in V r = { n b + 1 , . . . , n } , who know thevalues of their own node variables. Here x u represents log( α u ) for skew estimation and β u for offset estimation.Without loss of generality, we assume node variables of reference nodes are all 0, i.e. skews are 1 and offsets are0. Time is measured by a discrete time-index k = 0 , , . . . . The mobile nodes define a time-varying undirected measurement graph G ( k ) = ( V , E ( k )) , where ( u, v ) ∈ E ( k ) if and only if u and v can obtain a differencemeasurement of the form ζ u,v ( k ) = x u − x v + ξ u,v ( k ) , (5)during the time interval between k and k + 1 , where ξ u,v ( k ) is measurement error. We assume that between u and v , whoever obtains the measurement first shares it with the other so that it is available to both u and v . We alsofollow the convention that the difference measurement between u and v that is obtained by the node u is always of x u − x v while that used by v is always of x v − x u . Since the same measurement is shared by a pair of neighboringnodes, if v receives the measurement ζ u,v ( k ) from u , then it converts the measurement to ζ v,u ( k ) by assigning ζ v,u ( k ) := − ζ u,v ( k ) . For similar reasons, between a pair u and v , the node who computes ζ u,v ( k ) in node pair u and v is fixed for all time k . This can be achieved by comparing the magnitude of the index of nodes. For example,if u > v , then u computes ζ u,v ( k ) first and then sends it to v . The neighbors of u at k , denoted by N u ( k ) , is theset of nodes that u has an edge with in the measurement graph G ( k ) . We assume that if v ∈ N u ( k ) , then u and v can also exchange information through wireless communication at time k .Now the reformulated problem is to estimate the node variables x u , u ∈ V b , by using the difference measurements ζ u,v ( k ) , ( u, v ) ∈ E ( k ) that become available over time. We assume n r ≥ (i.e., there exists a least one referencenode), otherwise the problem is indeterminate up to a constant.III. T HE D I S YNC ALGORITHM
We present an iterative algorithm that nodes can use to solve the problem of node variable estimation from noisydifference measurements in a distributed manner. Since nodes do not have synchronized clocks, iterative updateshave to be performed asynchronously. Each node u ∈ V b keeps its local iteration index k u and maintains an estimate ˆ x u ( k u ) ∈ R of its node variable x u in its local memory. The estimates can be initialized to arbitrary values. Inexecuting the algorithm, node u starts its i -th iteration at a pre-specified local time τ ( i ) , for i = 0 , , . . . , which willbe described in Section III-A. Then, node u obtains current estimates ˆ x v ( k u ) along with the measurements ζ u,v ( k u ) from its current neighbors v ∈ N u ( k u ) . After a fixed time length δt (measured in local time), node u increments itslocal iteration index k u by 1, and updates its new estimate based on current measurements and neighbors’ estimatesby using the following update law: ˆ x u ( k u + 1) =ˆ x u ( k u ) + m ( k u ) X v ∈N u ( k u ) a uv ( k u )(ˆ x v ( k u )+ ζ u,v ( k u ) − ˆ x u ( k u )) , (6)where the time varying gain m ( · ) : Z + → R + has to be specified to all nodes a-priori. Note that when N u ( k u ) = ∅ , ˆ x u ( k u + 1) = ˆ x u ( k u ) . The choice of m ( · ) will play a crucial role in the convergence of the algorithm and willbe described in Section IV. In this paper, we let weight a uv ( k u ) = 1 if ( u, v ) ∈ E ( k u ) . Recall that the referencenodes take part by helping their neighbors obtain difference measurements, but they do not update their own nodevariables. The algorithm is summarized in Algorithm 1. Note that, since obtaining difference measurements requiresexchanging time-stamped messages, current estimates can be easily exchanged during the process of obtaining newmeasurements. Algorithm 1
DiSync algorithm at node u while u is performing time synchronization do if Local time τ u = τ ( i ) , i = 0 , , . . . then u collects current local indices k v from neighbors v ∈ N u ( k u ) . for all v ∈ N u ( k u ) do if k u = k v and u does not have ζ u,v ( k u ) then u and v perform pairwise communication; u saves ζ u,v ( k u ) and ˆ x v ( k u ) ; v saves ζ v,u ( k v ) and ˆ x u ( k v ) ; else u and v stop the communication; end if end for end if if τ u = τ ( i ) + δt , i = 0 , , . . . then if N u ( k u ) = ∅ then u updates ˆ x u ( k u + 1) using (6); else ˆ x u ( k u + 1) = ˆ x u ( k u ) ; end if u updates, k u = k u +1; end if end while A. Iteration schedule and synchronous view
We will later describe that the gains m ( · ) is chosen to be a decreasing function of time, which helps reduce theeffect of measurement noise. This is a well-known idea in stochastic approximation. However, using this idea ina network of unsynchronized clocks presents an unique challenge since no node has a notion of a common timeindex, 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 ( i + 1) -thiteration while a node with slower skew hasn’t even finished the i -th iteration. Therefore, specifying a function m ( · ) to all the nodes does not ensure that nodes use the same gain at the same (global) interval, which is requiredby stochastic approximation.We address the problem by providing the nodes a priori the sequence of local time instants τ ( i ) , i = 0 , . . . mentioned at the beginning of the Section III. This sequence is called an iteration schedule , and the formula forcomputing it is described below. Let the skews and offsets of all clocks be lower and upper bounded by thosein two fictitious clocks c L and c H , such that α c L ≤ α u ≤ α c H , β c L ≤ β u ≤ β c H . Recalling (2), therefore τ c L ( t ) ≤ τ u ( t ) ≤ τ c H ( t ) for all u ∈ V . The formula for calculating τ ( i ) is τ ( i +1) = α c H α c L ( τ ( i ) + δt − β c L ) + β c H , (7)where τ (0) has to satisfy τ (0) > β c H . This schedule ensures that nodes operating on their unsynchronized lo-cal clocks still perform updates in an effectively synchronous manner over time. To see this, define I ( i ) :=( τ ( i ) − β cH α cH , τ ( i +1) − β cH α cH ) as a global interval and I ( i ) u := ( τ ( i ) − β u α u , τ ( i ) + δt − β u α u ) as the global time interval withrespect to i -th local iteration of node u . Eq. (7) guarantees that, at each i , I ( i ) u ⊂ I ( i ) for all u ∈ V . In other words,there exists a sequence of global time intervals such that each i -th global interval contains, and only contains, the i -th local iteration (in global time) of all u ∈ V . Figure 1(a) shows the relationship between intervals of localiterations and the corresponding global intervals. In Figure 1(b), we pick the 3rd global interval from Figure 1(a),and show the global time intervals when local iteration updates occur. We emphasize that τ ( i ) is the same for allnodes and every node u starts and ends its i -th iteration at the same local time instants τ ( i ) and τ ( i ) + δt . Eachnode is provided the values of the parameters α cH α cL , β c L , β c H , and δt ahead of time, which are design variables.We next address the issue of how to pick values for α cH α cL , β c L and β c H without knowing a real bounds on skewsand offsets of all clocks in a network. In wireless sensor nodes, a pair of clocks in sensor nodes usually drift apartup to µsec/sec [2]. Therefore, we can pick α c H /α c L ≈ × − . To pick reasonable values of the offsetbounds, the following procedure should be used to initialize the synchronization. The reference node first broadcastsa message (to indicate the beginning of synchronization) and sets its local clock time t to zero simultaneously. Anode that receives this message sets its own clock to zero and broadcast such message again. The nodes that hearthis message also set their local clocks to , and so forth. Since nodes (except the reference node) start their localclocks after – but close to – the instant of t = 0 , their offsets are negative and small. Therefore, β c H can be chosenas zero and β c L can be picked as an estimate of the time it takes for all active nodes to receive the “synchronizationstart” signal. For a node who was out of communication range at the beginning but joins the networks later, it canset the local time to the current local time of a neighbor that has already started the time synchronization, andrecord neighbor’s iteration index as well. In this way, the newly joined node can take part in the synchronizationprocess as if it started at the very beginning.Another practical issue is the unbounded growth of the inter-synchronization intervals τ ( i ) . For example, if δt ischosen as second, the choice of α c H /α c L = 1 + 4 × − ensures that the time interval between two successive iterations, τ ( i +1) − τ ( i ) , will increase from second to seconds after . × iterations, or . hours. Ifthe updates are to be done more slowly, a larger δt will be used, which will slow down the growth of the iterationinterval τ ( i +1) − τ ( i ) . If it is desired that the inter-synchronization interval does not increase beyond a certainpre-specified value, a reference node to restart the synchronization after a certain time to maintain τ ( i +1) − τ ( i ) within the desired bound. When the restart should occur can be computed from the recursion (7). Lo c a l t i m e Localinterval Globalinterval τ max (t)=1.1t+0.3 τ min (t)=0.9t−0.3 PSfrag replacements I (3) cH I (3) u I (3) v I (3) cL I (3) (a) PSfrag replacements I (3) cH I (3) u I (3) v I (3) cL I (3) (b)Fig. 1. In (a), the X-axis is labeled by the time instants of the beginning of global intervals and the Y-axis is labeled by the sequence of τ ( i ) . The red solid slanted lines represent two fictitious clocks c L and c H as the bounds for local clocks in nodes. The black solid verticallines divide global time into a sequence of I ( i ) . Each i -th interval from black solid vertical to black dotted-dash vertical line is the interval ( τ ( i ) − β cH α cH , τ ( i ) − β cL α cL ) , which contains the global time instants of τ ( i ) for all u ∈ V . In (b), the 3rd global interval is picked as also circledin (a). The segments in the second and fifth rows correspond to I (3) c H and I (3) c L of two fictitious clocks respectively. The segments in the thirdand fourth rows present I (3) u and I (3) v of any two nodes u, v ∈ V accordingly. Remark 1.
Nodes perform skew and offset estimation simultaneously in a distributed and iterative fashion by usingtwo copies of the DiSync algorithm, one for skew estimation and one for offset estimation. With current estimatedskew ˆ α u ( k ) and offset ˆ β u ( k ) , a node u can compute the global time using (2) , i.e. ˆ t u ( k ) = ( τ u − ˆ β u ( k )) / ˆ α u ( k ) ,which is the final step of the time synchronization. The entire time synchronization procedure is illustrated inFigure 2. IV. C
ONVERGENCE ANALYSIS OF D I S YNC ALGORITHM
Since there exists a common iteration counter k that can be used to describe the local updates in the nodes byusing the iteration schedule (even though none of the nodes has access to it), we consider only the synchronousversion of the algorithm using global index k from now on. We rewrite (6) as ˆ x u ( k + 1) =ˆ x u ( k ) + m ( k ) X v ∈N u ( k ) a uv ( k )(ˆ x v ( k )+ ζ u,v ( k ) − ˆ x u ( k )) . (8) PairwiseSynchronization DiSyncAlgorithm
PSfrag replacements log ˆ α u,v ( k )ˆ β u,v ( k ) ˆ α u ( k )ˆ β u ( k ) ˆ t u ( k ) k ← k + 1 Fig. 2. Time synchronization by using DiSync.
Defining the estimation error as e u ( k ) := ˆ x u ( k ) − x u , Eq. (8) reduces to the following using (5): e u ( k + 1) = e u ( k ) + m ( k ) X v ∈N u ( k ) a uv ( k )( e v ( k ) − e u ( k ) + ǫ u,v ( k )) . (9)We introduce the following stipulations and notations to pursue subsequent analysis. First,let a uv ( k ) = 0 for v / ∈ N u ( k ) . Secondly, the n × n Laplacian matrix L ( k ) of the graph G ( k ) is defined as L uv ( k ) = P nv =1 a uv ( k ) if u = v , and L uv ( k ) = − a uv ( k ) if u = v . By removing the rows and columns of L ( k ) with respect to referencenodes, we get the n b × n b principle submatrix L b ( k ) (so called grounded or Dirichlet Laplacian matrix [21]). Let e ( k ) := [ e ( k ) , . . . , e n b ( k )] T , the corresponding state space form of the estimation error is e ( k + 1) = ( I − m ( k ) L b ( k )) e ( k ) + m ( k ) D ( k ) ǫ ( k ) , (10)where ǫ ( k ) := [¯ ǫ ( k ) T , . . . , ¯ ǫ n b ( k ) T ] T , ¯ ǫ u ( k ) := [ ǫ u, ( k ) , . . . ǫ u,n ( k )] T ,D ( k ) := diag (¯ a ( k ) , . . . , ¯ a n b ( k )) , ¯ a u ( k ) := [ a u, ( k ) , . . . , a u,n ( k )] , where ǫ u,u ( k ) / ∈ ¯ ǫ u ( k ) and a u,u ( k ) / ∈ ¯ a u ( k ) . When a uv ( k ) = 0 , ǫ u,v ( k ) is a pseudo random variable withthe same mean and variance as the measurement noise on any existing edge. Moreover, as a node u computesmeasurement ζ u,v ( k ) and sends − ζ u,v ( k ) to v , ǫ u,v ( k ) = − ǫ v,u ( k ) . Now, we introduce the following assumptions: Assumption 1.
Measurement noise vector ǫ ( k ) is with mean E [ ǫ ( k )] = γ and bounded second moment, i.e. E [ k ǫ ( k ) k ] < ∞ , where k · k denotes 2-norm. Furthermore, ǫ ( k ) and ǫ ( j ) are independent for k = j . In addition, { ǫ ( k ) } is independent of e (0) , where E k e (0) k < ∞ . In practice, E [ ǫ ( k )] may be time-varying even if all the underlying processes are wide sense stationary. Forinstance, the bias in offset difference measurement computed from node u is different from that computed fromnode v , as β v (1 − α u α v ) = β u (1 − α v α u ) . Therefore, E [ ǫ ( k )] depends on which node initializes pairwise synchronizationat time k . To meet this requirement E [ ǫ ( k )] = γ in Assumption 1, we can stipulate that the node who computes ζ u,v ( k ) between a pair u and v is fixed for all time k . This can be achieved by comparing the magnitude of the index of nodes. For example, if u > v , then u computes ζ u,v ( k ) first and then sends it to v . Indeed, the purpose ofthis requirement is to provide formula to compute the steady state value of estimation error, and the system maystill achieve convergence without it. Assumption 2.
The non-increasing positive sequence { m ( k ) } (step size of the stochastic approximation) is chosenas m ( k ) = c k + c , where c , c are constant real numbers. Note that Assumption 2 is a special case of the standard requirement in stochastic approximation: P ∞ k =0 m ( k ) = ∞ and P ∞ k =0 m ( k ) < ∞ . The assumption is made to simplify the subsequent analysis, but we believe it is notnecessary. A. Deterministic Topology Switching
In this section, we analyze the convergence of (8) when the topology switching is deterministic.
Assumption 3.
There exists d ∈ N s.t. for any t ≥ , ˆ G dt := S t + d − k = t G ( k ) = ( V , S t + d − k = t E ( k )) is connected, where E ( k ) is set of edges in G ( k ) . Assumption 4.
The limits ¯ L, ¯ L b , ¯ D exist: ¯ L := lim t →∞ t P tk =0 L ( k ) , ¯ L b := lim t →∞ t P tk =0 L b ( k ) , ¯ D :=lim t →∞ t P tk =0 D ( k ) . Assumption 3 implies that information can go from any node to the rest nodes within a uniformly boundedlength of time. Furthermore, as G ( k ) is bidirectional, another equivalent assumption is that ˆ G dt contains a spanningtree. The proposed algorithm is also robust to permanently adding or deleting nodes in case the new resultinggraph satisfies the assumption on connectivity. To understand the Assumption 4, we define the finite state space G = {G , . . . , G N } as the set of graphs that can occur over time. If the sequence of G ( k ) can be divided into asequence of finite intervals I j , j = 1 , , . . . , such that the percentage of times that each state G k occurs is fixed in allexcept finitely many such intervals I j , then ¯ L , ¯ L b and ¯ D exist. Another example is that the state G i occurs accordingto a sample path of a stationary ergodic process. In addition, we denote sets of matrices L b = { L b , . . . , L bN } and D = { D , . . . , D N } , where L bi and D i correspond to G i ∈ G . If the percentage of all states occurring is π = { π , π , . . . , π N } , then ¯ L b := N X i =1 π i L bi , ¯ D := N X i =1 π i D i (11) Theorem 1.
Under Assumption 1-4, the Algorithm 1 ensures that e ( k ) in (10) converges to ¯ L − b ¯ Dγ in mean square,i.e., lim k →∞ E ( k e ( k ) − ¯ L − b ¯ Dγ k ) = 0 . (cid:3) The theorem states that under the assumptions, the variance of the estimation error decays to . If additionallyall the difference measurements are unbiased ( γ = 0 ), 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 foundin [22], while the second is from[23]. Lemma 1.
If difference measurement is unbiased, i.e., γ = 0 , under assumption 1-3, the Algorithm 1 ensures that e ( k ) in (10) converges to in mean square, i.e., lim k →∞ E ( k e ( k ) k ) = 0 . (cid:3) When γ = 0 , (10) can be regarded as a leader-following consensus problem with time-varying topology andzero-mean noisy input. The leaders are reference nodes u ∈ V r , which hold their variables as zero. Then, e u ( k ) for u ∈ V b is driven to zero by the reference nodes as k goes to ∞ in mean square sense. Lemma 2. [23]Denote by A an unknown symmetric and positive semi-definite matrix in R n × n , and we have to solvethe equation Ax=y for an unknown y ∈ R n . Assume that A − exists. We are given a sequence of matrices A k and asequence y k , where k = 0 , , . . . . In addition, suppose that lim k →∞ k k P ki =1 y i − y k = 0 , lim k →∞ k k P ki =1 A i − A k = 0 , lim k →∞ k P ki =1 k A i k exists. Consider the sequence x k : x is arbitrary, x k +1 = x k + c k + c ( y k − A k x k ) , (12) where c and c are constant real numbers. Then, lim k →∞ x k = A − y . (cid:3) Proof of Theorem 1.
Taking expectations on both sides of (10) with respect to measurement noise ǫ ( k ) , we obtain η ( k + 1) = ( I − m ( k ) L b ( k )) η ( k ) + m ( k ) D ( k ) γ, (13)where η ( k ) = E [ e ( k )] . By substituting (13) in (10), we get ˜ e ( k + 1) = ( I − m ( k ) L b ( k ))˜ e ( k ) + m ( k ) D ( k ) ξ ( k ) , (14)where ˜ e ( k ) = e ( k ) − η ( k ) and ξ ( k ) = ǫ ( k ) − γ . Note that ξ ( k ) is zero mean and satisfies Assumption 1. ByLemma 1, ˜ e ( k ) converges to 0 in mean square. Therefore, e ( k ) is mean square convergent to η ( k ) . Now, weexamine the convergence of η ( k ) . From the definition of ¯ L b and the symmetry of L bi , ¯ L b is a symmetric groundedLaplacian of ˆ G . Since ˆ G is connected, by Lemma 1 in [21], ¯ L b is positive definite. Consequently, λ m ( ¯ L b ) > .Now, it follows from Lemma 2 that lim k →∞ η ( k ) = ¯ L − b ¯ Dγ . Consequently, e ( k ) converges to ¯ L − b ¯ Dγ in meansquare, i.e., lim k →∞ E ( k e ( k ) − ¯ L − b ¯ Dγ k ) = 0 . (cid:4) B. Markovian Topology Switching
In this section, we analyze convergence when network topology switches randomly. We model the switchingof the network topologies as a Markov chain; the reasonableness of this model for mobile networks has beenestablished in [13].
Assumption 5.
The temporal evolution of the measurement graph G ( k ) is governed by an N-state homogeneousergodic Markov chain with state space G = {G , . . . , G N } , which is the set of graphs that can occur over time. Furthermore ˆ G := S Nk =1 G k = ( V , S Nk =1 E k ) is connected, where E k is set of edges in G k . In addition, theprocesses G ( k ) and ǫ ( j ) are independent for all k and j . In Assumption 5, the Markovian switch on the graphs means that P ( G ( k + 1) = G i |G ( k ) = G j ) = P ( G ( k + 1) = G i |G ( k ) = G j , G ( k −
1) = G ℓ , . . . , G (0) = G p ) where G i , G j , G ℓ , . . . , G p ∈ G . The requirement for ergodicity of theMarkov chain ensures that there is an unique steady state distribution with non-zero entries. This means every graphin the state space of the chain occurs infinitely often. Since ˆ G is connected, ergodicity implies that informationfrom the reference node(s) will flow to each of the nodes over time. Again, note that none of the graphs that everoccur is required to be a connected graph. Since the Markov chain is ergodic, the steady state distribution of thechain is defined as π = { π , π , . . . , π N } . Recalling that L bi and D i correspond to G i ∈ G , we can use the sameformula as that in (11) to define ¯ L b and ¯ D . Theorem 2.
Under Assumption 1,2 and 5, e ( k ) in (10) is mean square convergent, i.e., lim k →∞ E ( k e ( k ) − E ǫ ( e ( k )) k ) = 0 , where E ǫ ( e ( k )) is expectation of e ( k ) w.r.t. measurement noise ǫ ( k ) , and E ǫ ( e ( k )) converges to ¯ L − b ¯ Dγ almost surely. The proof of the theorem uses the next two lemmas.
Lemma 3.
If relative measurements are unbiased, i.e., γ = 0 , under 1,2 and 5, e ( k ) in (10) converges to inmean square, i.e., lim k →∞ E ( k e ( k ) k ) = 0 . Lemma 4. [Proposition 1 in [24]] Assume { A k , y k } , k = 0 , , . . . , is stochastic process on (Σ , F , P ) , where A k is an n × n symmetric positive semidefinite matrix and y k is n × . Consider the following iteration with arbitrary x : x k +1 = x k + c k + c ( y k − A k x k ) , (15) where c and c are real constants. If A := lim k →∞ k P ki =1 E [ A i ] , y := lim k →∞ k P ki =1 E [ b i ] , and A is positivedefinite, then lim k →∞ x k = A − y almost surely. Proof of Lemma 3 is omitted due to its length; it can be found in [22]. Lemma 4 follows in a straightforwardmanner 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 η ( k ) = E ǫ [ e ( k )] , where the expectation istaken with respect to measurement noise ǫ ( k ) . Take the expectation on both sides of (10), we get η ( k + 1) = ( I − m ( k ) L b ( k )) η ( k ) + m ( k ) D ( k ) γ. (16)By substituting (16) in (10), we get ˜ e ( k + 1) = ( I − m ( k ) L b ( k ))˜ e ( k ) + m ( k ) D ( k ) ξ ( k ) , (17) where ˜ e ( k ) = e ( k ) − η ( k ) and ξ ( k ) = ǫ ( k ) − γ . Note that ξ ( k ) is zero mean and satisfies Assumption 1. ByLemma 3, ˜ e ( k ) converges to 0 in mean square. Now, we examine the convergence of η ( k ) . We rewrite (16) as η ( k + 1) = η ( k ) + m ( k )( D ( k ) γ − L b ( k ) η ( k )) . (18)It follows from Assumption 5, lim k →∞ k k X i =1 E [ L b ( k )] = N X i =1 π i L bi = ¯ L b , lim k →∞ k k X i =1 E [ D ( k )] = N X i =1 π i D i = ¯ D (19)From the definition of ¯ L b and the symmetry of L bi , ¯ L b is a symmetric grounded Laplacian of ˆ G . Since ˆ G isconnected, by Lemma 1 in [21], ¯ L b is positive definite. Consequently, λ m ( ¯ L b ) > . Furthermore, L bi is positive-semi definite. Then, by Lemma 4 that η ( k ) converges to ¯ L − b ¯ Dγ almost surely.V. A MELIORATING SLOW CONVERGENCE : D I S YNC -IThe proposed DiSync algorithm ensures mean square convergence by attenuating the measurement noise usingideas from stochastic approximation. In particular, the gain m ( k ) that decays slowly with time is instrumental indriving the variances of the estimation error to zero. This slowly decaying gain, however, also makes the convergencerate slow. We’ll see numerical evidence of this in Section VI. In practice, performance of the DiSync algorithmcan be further improved by modifying the update law, which we describe next. The modified algorithm is calledthe DiSync-I(DiSync-Improved) algorithm.First, we define a scalar state y u ( k ) for each node u , which is a surrogate for the distance (number of hops)of node u from the reference nodes. The state y u ( k ) is called average distance of u from the reference nodes attime k . Furthermore, we define the set S u ( k ) := { v | y u ( k ) ≥ y v ( k ) , v ∈ N u ( k ) } , the subset of neighbors of node u whose average distances are smaller than that of node u at k . The average distance is updated according to theiterative law described in Algorithm 2. In brief, y u ( k ) is updated by averaging all y v ( k ) from its neighbors whohave smaller average distance. If none of its neighbors have smaller average distance than itself, y u ( k ) increases.Both S u ( k ) and y u ( k ) are maintained in u at index k .The node variable update law for the DiSync-I algorithm is given below; where the subscript on the local iterationcounter is suppressed: ˆ x u ( k + 1) =ˆ x u ( k ) + h ( k ) X v ∈H u ( k ) a uv ( k )(ˆ x v ( k )+ ζ u,v ( k ) − ˆ x u ( k )) , (20) where h ( k ) = |H u ( k ) | for k < k h m ( k − k h ) for k ≥ k h , (21) H u ( k ) = S u ( k ) for k < k H N u ( k ) for k ≥ k H , (22)where k h and k H are constants - that satisfy k H ≤ k h - that are pre-specified to all nodes.If y u ( k ) ≥ y v ( k ) , it means that node v has been closer to the reference nodes than node u on average (evenif node v may be father from the reference node than u at current index k ). This indicates that node v is likelyto contain better estimates that u . If u and v are neighbors, u should use the estimate from v to perform update,while v should not use estimates from u . Algorithm 2
Average distance algorithm at node u ∈ V b Initialize y u (0) = ∞ while u is performing iteration do for v ∈ N u ( k ) do if y u ( k ) ≥ y v ( k ) and y v ( k ) = ∞ then S u ( k ) ← v ; end if end for if S u ( k ) = ∅ then y u ( k + 1) = P v ∈S u ( k ) y v ( k ) |S u ( k ) | ; else y u ( k + 1) = y u ( k ) + 0 . ; end if k = k + 1 ; end while Note that when k h = ∞ and k H = 0 , (20) becomes the JaT algorithm of [13, 16]. As shown in [13, 16],JaT algorithm ensures that the mean of the estimation error converges to a constant (zero if measurement isunbiased) and variance to a constant. In addition, when k h = ∞ and k H < ∞ , the resulting update law (20) is amodified version of JaT algorithm: it now uses Algorithm 2 during the initial phase up to k ≤ k H . We will call itthe JaT-I algorithm.Table I shows how the update law (20) can produce different algorithms depending on the values of the parameters k h , k H . For simulation studies in this paper on the DiSync-I algorithm, we pick k h = k H somewhat arbitrarily. Inthis case, DiSync-I first adopts JaT-I when k < k h and then becomes DiSync when k > k h . The reason behind themodification (20) over (8) is the following. First, JaT has better convergence speed during initial phase than that ofDiSync. Second, JaT-I has even better convergence rate than JaT due to the use of only those neighbors that havebeen closer to the reference node(s). Fig. 3. Two graphs that occur during one simulation with nodes moving according to the random direction mobility model. Note that the DiSync-I differs from DiSync only during the initial phase (up to k h ), otherwise it is the same. Asa result, the mean square convergence results of DiSync holds for DiSync-I as well. TABLE IC
OMPARISON OF DIFFERENT ALGORITHMS k H = 0 k H < ∞ k H = ∞ k h = 0 DiSync — — k h < ∞ — DiSync-I — k h = ∞ JaT JaT-I —
VI. S
IMULATION EVALUATION
We now examine through simulations the time synchronization performance of the DiSync and DiSync-I algo-rithms, as well as that of JaTand JaT-I algorithms. Finally, they are compared with the virtual time synchronizationalgorithm ATS [17] in terms of pairwise synchronization error. Simulations are performed in a 10-node mobilenetwork within a m × m square field. Nodes’ motions are generated according to the widely used randomwaypoint (RWP) mobility model [25]. It has been shown in [13] that when nodes move according to the RWPmobility model, the switching of the graphs can be modeled as a Markov chain. A pair of nodes can communicatewhen distance between them is less than m . The true skews and offsets of 9 nodes are picked uniformly from [1 − × − , × − ] and [ − − , − ] sec respectively according to [8]. The single reference node (10th)has skew and offset . Denote k = 1 , , . . . as update intervals (also called synchronization periods), and t k asthe global instant of the beginning of k -th interval. In this simulation, the interval is chosen as 1 sec, therefore t k +1 − t k = 1 . For the sake of convenience, simulations are carried out in a synchronous fashion. A. Implementation of pairwise synchronization
The simulation of the DiSync, DiSync-I, JaT and JaT-I algorithms requires pairs of nodes to obtain differencemeasurements by exchanging time stamped messages when they are within communication range. In order toevaluate the entire time synchronization procedure, unlike [10–12, 16], in which difference measurements aregenerated 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 α u,v and relative offset β u,v , and the difference measurements β u − β v and log( α u ) − log( α v ) are then obtained from these as described in Section II. According to [3], at thebeginning of the k th interval, node u sends a message to v that contains the value of the local time at u when themessage is sent: τ (1) u . When node v receives this message, it records the local time of reception: τ (1) v . After a waitingperiod, node v sends a message back to u that contains both τ (2) v and τ (1) v . When it arrives at u , node u againrecords the local time of reception: τ (2) u . Two nodes u and v in communication range performs this procedure, calledtwo-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 u uses the obtained eight time stamps { τ ( i ) u , τ ( i ) v } for i = 1 , . . . to estimate α u,v ( k ) and β u,v ( k ) via the formula provided in [3]. Finally, node u sends back to v 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 α u,v ( k ) and β u,v ( k ) , and thus determines the level of noisein the resulting difference measurement. Therefore, the random delay ultimately affects the time synchronizationaccuracy. In employing the pairwise synchronization protocol of [3], we subject the message exchanges to a randomdelay that is Gaussian distributed with mean µsec and standard deviation µsec , as these values are consideredrealistic for wireless sensors networks with current hardware and communication protocols with uncertain-delayelimination (e.g., MAC-layer time-stamping) [8]. Although the relation between the statistics of the random delaysand the noise of difference measurements of log-skew and offsets defined in Section II are complex, the noise levelsin the difference measurements used are likely to be realistic due to realistic choice of delays. B. Performance in estimating global time
We conduct 1000 Monte Carlo simulations of running algorithms for 800 second (iterations). Figure 3 shows twosnapshots of the network during a simulation. As we can see, only a limited number of nodes can communicatewith each other. Recall that a uv ( k ) = 1 if ( u, v ) ∈ E ( k ) for all k . The step size function is chosen as m ( x ) = x +3 . We pick k h = k H = 40 somewhat arbitrarily. Moreover, to evaluate the algorithms under sleep-wake cycleimplementation for energy conservation, we force nodes to pause the updates when k ∈ [400 , . When k > ,the step size function is changed to m ( k − k h − , i.e. it resumes the value before the pause.Figure 4(a) and 4(b) show the mean and variance of estimation error of the skew of node . As expected, both thevariances of DiSync and DiSync-I are seen to converge to zero, while that of JaT and JaT-I converge to constantvalue. In addition, DiSync-I improves the accuracy of the mean of estimation error at the expense of slightlyincreasing the variance of estimation error.Figure 5 shows the global time estimation error in one experiment, i.e. ˆ t u ( t k ) − t k as a function of t k for node . Both DiSync and DiSync-I show much higher accuracy of global time estimation than that of JaT and JaT-I.The accurate skew estimates is crucial in getting good global time estimates, since even a tiny error in the skewestimate leads to a large error in the prediction of global time t over time. In addition, the DiSync-I further reducesthe initial transient period that DiSync suffer from. −6 Time (sec) M ean JaTJaT−I DiSyncDisync−I True (a) Mean −13 −12 −11 −10 Time (sec) V a r i an c e JaTJaT−I DiSyncDisync−I (b) VarianceFig. 4. Empirically estimated mean and variance of the estimation error of skews in node 3. Note that in (b), y-axis is in logarithm scale. T i m e e s t. e rr o r ( s e c ) JaTJaT−I DiSyncDisync−I True
Fig. 5. The estimation error of global time in node 3
C. Comparison with ATS
In ATS [17], each node u estimates the virtual global time using ˆ t ru ( t ) = ˆ ̺ u ( k ) τ u ( t ) + ˆ o u ( k ) for t k ≤ t ≤ t k +1 ,where variables ˆ ̺ u ( k ) and ˆ o u ( k ) can be thought of as the skew and offset of a virtual global time respect tothe local time of u during interval k . Here, we present a synchronous version of ATS algorithm that we use insimulation in order to be consistent with the discussion. Each node u updates its ˆ ̺ u ( k ) and ˆ o u ( k ) using ˆ ̺ v ( k ) , ˆ t rv ( k ) , ˆ α uv ( k ) from its neighbors, where ˆ α uv ( k ) is the estimated relative skew. ˆ α uv ( k ) is obtained by pairwisecommunication between u and v during k -th update interval, as part of the ATS algorithm. This is done as follows.Two time-stamped messages are sent from node v to node u : one at the beginning of k -th interval and the otherone in the middle of the interval. Note that no return messages from u to v is required. The computation of ˆ α uv ( k ) is performed by a low-pass filter as provided in ATS: ˆ α u,v ( k ) = ρ ˆ α u,v ( k −
1) + (1 − ρ ) τ (1) v − τ (2) v τ (1) u − τ (2) u , where ρ is a −4 −3 −2 −1 Time (sec) S y n c e rr o r ( s e c ) ATSJaT JaT−IDiSync Disync−I
Fig. 6. Maximum synchronization error along time in one experiment. tuning parameter and chosen as . (same value used in ATS). It has been shown that lim k →∞ α u ˆ ̺ u ( k ) = ¯ α , lim k →∞ ˆ o u ( k ) + β u ˆ ̺ ( k ) = ¯ β , where ¯ α and ¯ β is the skew and offset of the virtual clock with respect to t . The ATSalgorithm ensures that the estimated virtual global times in all nodes are eventually equal, i.e., lim t →∞ ˆ t ru ( t ) = ˆ t rv ( t ) for all u and v .The performance of ATS is guaranteed under the assumption that the time stamps are exchanged without randomdelay. To compare with other algorithms under identical conditions, we add random delay to τ ( i ) u for i = 1 , . Thedelay parameters are the same as those used during the simulation of the other algorithms. In addition, since ATSdoes not estimate the clock time at any of the nodes, we use the metric “maximum synchronization error” to compareATS with the other four algorithms. This is defined as max u,v | ˆ t ru ( t k ) − ˆ t rv ( t k ) | for ATS, and max u,v | ˆ t u ( t k ) − ˆ t v ( t k ) | for the other four algorithms.Figure 6 compares maximum synchronization error in one experiment for all five algorithms. Although themaximum synchronization error decreases faster in ATS, JaT and JaT-I at the beginning, the superior robustness tothe measurement noise of the DiSync algorithm helps it outperform them after about 100 sec. It can be seen fromthe figure that the DiSync-I achieves the lowest maximum synchronization error among all algorithms.VII. C ONCLUSION
We proposed a novel distributed asynchronous algorithm to estimate clock skews and offsets from pairwisedifference measurements of log-skews and offsets. The nodes measure log-skew difference and offset differencewith nearby neighbors by exchanging time stamped messages. A node fuses these measurements with currentestimates to iteratively update their estimates of skew and offset. The algorithm is inspired by the recent workon using stochastic approximation in consensus algorithms. The time varying gain from stochastic approximationensures that the variance of skew/offset estimation error asymptotically converges to zero under certain assumptions.Using estimated skews and offsets, nodes can estimate the time of global clock accurately which was demonstratedin numerical evaluation. Simulations also show that the proposed algorithms outperform competing algorithms.The fact that clocks are not synchronized poses an unique challenge in applying stochastic approximation ideasto log-skew and offset estimation problem, since all nodes need to reduce a gain in a synchronous manner. This was addressed by providing an iteration schedule to the nodes ahead of time so that all nodes can effectivelyperform updates synchronously. The scheduled iteration interval grows along time, though the growth rate is quiteslow. In many cases the growing interval between iterations schedules may be useful since it automatically slowsdown the rate of synchronization over time. This has the desired effect of more iterations and faster lowering ofsynchronization error initially at the cost of high communication overhead, while lowering the overhead as moreaccurate synchronization is achieved. In cases when the growth of synchronization interval is undesirable, thereference nodes could restart the synchronization process. Another possibility is to use the estimated skews aftersome time has passed instead of using pre-specified bounds at all times. Since the proposed algorithms providehighly accurate skew estimates, doing so will ensure synchronous updates while keeping the time interval betweenupdates from growing. The effectiveness of this strategy will be studied in future work. Another avenue of futurework is a study of the effect of design parameters k H and k h in the DiSync-I algorithm. Their values determinewhen the nodes switch out of the initial Jacobi-type iterations and start using time-varying gains to reduce thevariance of the estimates. Waiting long to make the switch may mean that steady state is reached before the switchis made, but then the steady state variance is determined by the JaT algorithm, which is non-zero, leading to poorglobal time estimates. Switching soon will make sure the variance keeps decreasing from early on, but the decay rateis now slow since it is governed by the underlying stochastic approximation algorithm. Determining the “optimal”values of these parameters is a topic of future research.R EFERENCES [1] B. M. Sadler and A. Swami, “synchronization in sensor networks: an overview,” in
IEEE MILCOM , October 2006, pp.1–6.[2] J. R. Vig, “Introduction to quartz frequency standards,” Army Research Laboratory, Tech. Rep., 1992.[3] K.-L. Noh, Q. M. Chaudhari, E. Serpedin, and B. W. Suter, “Novel clock phase offset and skew estimation using two-way timing message exchanges for wireless sensor networks,”
IEEE Transactions on Communications , vol. 55, no. 4, pp.766–777, Apr 2007.[4] S. Yoon, C. Veerarittiphan, and M. L. Sichitiu, “Tiny-sync: Tight time synchronization for wireless sensor networks,”
ACMTransactions on Sensor Networks , vol. 3, no. 2, pp. 1–34, Jun 2007.[5] M. Leng and Y.-C. Wu, “On clock synchronization algorithms for wireless sensor networks under unknown delay,”
IEEETransactions on Vehicular Technology , vol. 59, no. 1, pp. 182–190, Jan 2010.[6] J. Elson, L. Girod, and D. Estrin, “Fine-grained network time synchronization using reference broadcasts,” in the FifthSymposium on Operating Systems Design and Implementation (OSDI) , 2002.[7] S. Ganeriwal, R. Kumar, and M. B. Srivastava, “Timing-sync protocol for sensor networks,” in
ACM Conference onEmbedded Networked Sensor Systems (SenSys) , 2003.[8] M. Mar´oti, B. Kusy, G. Simon, and ´A. L´edeczi, “The flooding time synchronization protocol,” in
ACM Conference onEmbedded Networked Sensor Systems (SenSys) , 2004.[9] P. Barooah and J. P. Hespanha, “Estimation from relative measurements: Error bounds from electrical analogy,” in
Proc. ofthe 2nd International Conference on Intelligent Sensing and Information Processing(ICISIP) , January 2005, pp. 88–93.[10] P. Barooah, N. M. da Silva, and J. P. Hespanha, “Distributed optimal estimation from relative measurements for localization and time synchronization,” in International Conference on Distributed Computing in Sensor Systems DCOSS’06 , SanFrancisco, June 2006, pp. 266 – 281.[11] A. Giridhar and P. R. Kumar, “Distributed clock synchronization in wireless networks: Algorithms and analysis (I),” in , December 2006, pp. 4915 – 4920.[12] R. Solis, V. S. Borkar, and P. R. Kumar, “A new distributed time synchronization protocol for multihop wireless networks,”in
Proc. of the 45th IEEE Conference on Decison and Control , December 2006, pp. 2734–2739.[13] C. Liao and P. Barooah, “Distributed clock skew and offset estimation from relative measurements in mobile networkswith markovian switching topology,”
Automatica
Automatica , vol. 46, no. 10, pp. 1571–1583, Oct. 2010.[15] T. Li and J. Zhang, “Consensus conditions of multi-agent systems with time-varying topologies and stochastic communi-cation noises,”
Automatic Control, IEEE Transactions on , vol. 55, no. 9, pp. 2043–2057, 2010.[16] C. Liao and P. Barooah, “Time synchronization in mobile sensor networks from difference measurements,” in
In proceedingsof the 49th IEEE Conference on Decision and Control , December 2010, pp. 2118 – 2123.[17] L. Schenato and F. Fiorentin, “Average timesynch: A consensus-based protocol for clock synchronization in wireless sensornetwork,”
Automatica , vol. 47, no. 9, pp. 1878 – 1886, 2011.[18] R. Carli, E. D’Elia, and S. Zampieri, “A pi controller based on asymmetric gossip communications for clocks synchro-nization in wireless sensors networks,” in
Decision and Control and European Control Conference (CDC-ECC), 2011 50thIEEE Conference on , dec. 2011, pp. 7512 –7517.[19] C. Liao and P. Barooah, “Di-Sync: High-accuracy distributed clock synchronization in mobile ad-hoc networks from noisyrelative measurements,” in
American Control Conference , June 2013, pp. 3338–3343.[20] ——, “Estimation from relative measurements in mobile networks with markovian switching topology: Clock skew andoffset estimation for time synchronization,”
ArXiV preprint , 2013. [Online]. Available: http://arxiv.org/abs/1301.2218[21] P. Barooah and J. P. Hespanha, “Graph effective resistances and distributed control: Spectral properties and applications,”in
Proc. of the 45th IEEE Conference on Decision and Control , December 2006, pp. 3479–3485.[22] C. Liao, “Distributed time synchronization from relative measurement in mobile wireless sensor networks,” Ph.D.dissertation, University of Florida, April 2013.[23] L. Gyrfi, “Stochastic approximation from ergodic sample for linear regression,”
Probability Theory and Related Fields ,vol. 54, pp. 47–55, 1980.[24] M. Kouritzin, “On the convergence of linear stochastic approximation procedures,”
Information Theory, IEEE Transactionson , vol. 42, no. 4, pp. 1305 –1309, jul 1996.[25] T. Camp, J. Boleng, and V. Davies, “A survey of mobility models for ad hoc network research,”