Clock Estimation for Long-Term Synchronization in Wireless Sensor Networks with Exponential Delays

Although the existing time synchronization protocols in wireless sensor networks (WSNs) are e ﬃ cient for short periods, many applications require long-term synchronization among the nodes, for example, coordinated sleep and wakeup modes, and synchronized sampling. In such applications, experiments have shown that even clock skew correction cannot maintain long-term clock synchronization and a quadratic model of clock variations can better capture the dynamics of the actual clock model involved, hence increasing the resynchronization period and conserving signiﬁcant energy. This paper derives the maximum likelihood (ML) estimator for all the clock parameters in a two-way timing exchange model with exponential delays. The same estimation procedure can be applied to one-way timing exchange models with little modiﬁcation.


INTRODUCTION AND RELATED WORK
A wireless sensor network (WSN) consists of several small scale miniature devices capable of onboard sensing, computing, and communications.WSNs are used in industrial and commercial applications to monitor data that would be difficult to monitor using wired sensors.Due to harsh operating conditions, they are mostly left unattended for their lifetimes without any maintenance or battery replacement.Therefore, limited energy and limited hardware are the most important constraints in the design of WSNs.
WSNs have numerous applications, for example, habitat monitoring, military surveillance, security, traffic monitoring, fire detection, object tracking, and industrial process monitoring.Time synchronization among the nodes in a WSN is important for various applications such as coordinated sleep and wakeup modes, object tracking, data fusion, security, and MAC protocols.Since energy is the scarcest resource in WSNs, a smart technique to conserve energy is to deploy coordinated turning on and off of radios in sensor nodes.If the nodes are time synchronized with each other, the efficient duty cycling operation of coordinated sleep and wakeup modes can be enabled which hugely boosts the lifetime of the network due to the minimal power consumption during the sleep mode.
Traditional clock synchronization techniques cannot be applied to WSNs because they were initially designed for general computer networks where communication comes for free and the computational resources are powerful.Such an environment comes in sharp contradiction with WSN requirements.As an example, network time protocol (NTP), the most widely used protocol for synchronizing computer networks [1], relies on a lot of communication messages between the server and the client, and hence is very energy expensive.Also, deploying global positioning system (GPS) on miniature nodes is very cost expensive, is not available indoors and underwater applications, can be easily jammed, and defeats the purpose of having a network of small-scale cheap nodes.
Recently, several efficient protocols have been proposed for short-term synchronization of WSNs such as timing synchronization protocol for sensor networks (TPSNs) [2], receiver broadcast synchronization (RBS) [3], and flooding time synchronization protocol (FTSP) [4], which can synchronize a pair of nodes within a few microseconds.TPSN [2] adjusts the clock offset between the two nodes only, while RBS [3] and FTSP [4] estimate both the clock offset and skew through linear regression.However, these schemes are only useful for short-term applications such as surveillance and object tracking and are unsuitable for efficient duty cycling and other applications that require continuous time synchronization such as synchronized sampling because they spend a lot of energy for resynchronization during a long time interval.To emphasize this fact, note that the most efficient time synchronization protocol reported thus far and implemented on real testbed, FTSP, has to resynchronize the nodes in the network every minute to achieve 90 microseconds synchronization error [4].This is the reason that even though RBS and FTSP estimate the clock skew alongside clock offset using linear regression, they are insufficient in practice for longterm synchronization, for example, the shooter localization system [5] uses FTSP to synchronize once every 45 seconds and the Center for embedded networked sensing (CENS) deployment at James Reserves [6] uses RBS to synchronize the nodes after every 5 minutes.Therefore, there is a need for a better model to estimate the clock parameters for achieving long-term synchronization in WSNs.And this paper targets the estimation of clock parameters by relating the clocks of two nodes through a quadratic model rather than a linear model used in previous research.
A detailed analysis of clock offset estimation assuming a symmetric exponential delay model was presented in [7].For a known fixed delay, the MLE of clock offset does not exist because the likelihood function does not possess a unique maximum with respect to the clock offset.However, [8] derived the MLE of the clock offset for an unknown fixed delay.This paper derives the MLE for the parameters which relate the clocks of two nodes through a model involving the clock offset, skew, and drift.Estimating the clock drift is important in light of the reasons mentioned above and finding the MLE is desirable due to its optimal properties for a large number of observations [9] (i.e., asymptotic unbiasedness, efficiency, and consistency).Although the estimation of clock parameters using a quadratic model is computationally more demanding than using the linear model, it helps in maintaining long-term synchronization among the nodes and subsequently less frequent communications and energy savings.Since it has been reported in [15] that the energy required to transmit 1 bit over 100 meters (3 Joules) is equivalent to the energy required to execute 3 millions of instructions, a tradeoff between spending reduced communication energy on the cost of more computational energy through estimating the long-term drift as well as the offset and the skew between clocks of two nodes is highly desirable.

THE MODEL
Figure 1 shows a model of a two-way timing message exchange mechanism between two nodes.When the two nodes start the synchronization process, Node 1 sends the timing message to Node 2 with its current time stamp T 1,r whose reception time is recorded as T 2,r at Node 2.Then, it sends at time T 3,r another synchronization message to Node 1 containing T 2,r and T 3,r , which time stamps the reception time of this message as T 4,r (see Figure 1).Hence, at the end of N such transmissions and receptions, Node 1 has access to a set of time stamps {T 1,r , T 2,r , T 3,r , T 4,r }, r = 1, . . ., N. Here, T 1,1 is the reference time, that is, every reading {T 1,r , T 2,r , T 3,r , T 4,r } is actually the difference between the actual recorded time and T 1,1 .Therefore, this model can be represented as where θ O , θ S , and θ D are the clock offset, skew, and drift of Node 2 with respect to Node 1 (the master node), respectively, d stands for the fixed portion of delay in the transmission of message from one node to another, for example, the sum of transmission time, propagation delay, reception time, X r and Y r denote the variable portion of delay and are assumed to be independent and exponentially distributed random variables with the same mean α.The justification of using the exponential model for the random delays is as follows.Several probability distribution function (PDF) models for random delays have been discussed in the literature where exponential, Gamma, lognormal, and Weibull distributions [10][11][12] have always been the most popular ones.As explained in [13], the exponential distribution fits quite well several applications.Also, a singleserver M/M/1 queue can fittingly represent the cumulative link delay for point-to-point hypothetical reference connection, where the random delays are independently modeled as exponential random variables [7].In addition, [7] not only justified the use of exponential distribution for modeling the delays but also presented several algorithms for estimating the clock offset amongst which the minimum link delay (MnLD) algorithm was experimentally demonstrated by [14] to perform the best.Using the minimum link delays to estimate the clock offset was mathematically proved by [8] as the ML estimator under the exponential delay model.This confirms that the exponential delay assumption for the random delays matches really well the experimental observations.

MAXIMUM LIKELIHOOD ESTIMATION
From (1), the general form of the likelihood function is given by where the indicator function I[•] is defined as From here onwards, without losing any generalization, we will assume that α is known.This is because even if α is unknown, due to the form of the reduced likelihood function L(d, θ O , θ S , θ D ), the MLE ( d, θ O , θ S , θ D ) remains the same (see [8]).Moreover, we assume that the clocks can neither stop nor run backwards, which implies that the clock skew θ S 0 and hence always positive.The actual values of practical clock skew is usually close to 1. Finally, for the simplification of derivation, θ D has been assumed to be positive.Following a similar procedure mentioned herein, a negative value of θ D will result in the same closed form expression of the MLE.
The constraints present in the likelihood function (2) can be expressed equivalently as These constraints can be viewed as 2N 4D curves due to the four unknowns.The 3D region where the two sets of N curves in ( 5) and ( 6) intersect each other yields θ O in terms of θ S and θ D as Plugging it back in (5), the sets of constraints can now be written as or equivalently, The above inequalities in (9) represent a 3D region due to three unknowns consisting of N 2 surfaces forming the boundary of the support region.To find this boundary of the support region as a function of θ D only, the intersection of these surfaces in (9) with each other is where and a is a simplified index notation as a function of the indices (i, j, k, l).Now plugging (10) into ( 9) yields the EURASIP Journal on Advances in Signal Processing support region in terms of d as a function of θ D only as where and b is again a simplified index notation as a function of the indices (m, n, p, q).Now the final form and uniqueness of the MLE can be proved by the following theorem.(1) From Figure 2, it is clear that the MLE lies on the boundary of the support region.This is because for any d lying somewhere within the support region, the likelihood function ( 2) can be further increased by increasing d until it reaches the boundary of the support region.
(2) Maximizing the likelihood function in ( 2) is equivalent to minimizing the exponential function argument ] in the likelihood function expression.Therefore, plugging (10) and ( 12) into the expression for Φ, it can be written in the form of a set φ a,b depending on indices a and b as (3) Starting from z b corresponding to min b {w b } (i.e., starting from the left with the first side of the semipolygon in Figure 2) and evaluating φ a,b on each subsequent z b on the boundary of the support region, observe that for each particular segment, φ a,b can be minimized by taking the largest possible (5) (e, f , g, h) = arg min m,n,p,q {θ m,n,p,q D }; where the indices {i, j, k, l, m, n, p, q} correspond to the two set of curves in (12) for which the sign of N r=1 {(T 2 4,r −T 2 1,r )+ v c (T 4,r −T 1,r )}−Nz b changes from negative to positive.Consequently, plugging θ D in ( 12), (10), and (7), we can write The complete procedure for finding this MLE ( θ D , d, θ S , θ O ) is explained in Algorithm 1. Algorithm 1 starts from the curve in (12) for which w has the least value.It selects the intersection of this curve with the neighboring curve intersecting it, and it checks the sign change condition of N r=1 {(T 2 4,r − T 2  1,r )+v c (T 4,r − T 1,r )}−Nz b .If the condition is not satisfied, the first curve is discarded and the same procedure is repeated for the second curve and so on until the same condition is satisfied.Figure 3 shows the MSE of θ D as a function of the number of timing messages N. It is evident that the MLE performs well even for small N and hence suitable for the power limited regime of WSNs.

CONCLUSIONS AND FUTURE WORK
Using a quadratic model for the relationship between the clocks of two nodes with a two-way timing message exchange mechanism, we have derived the MLE for the clock offset, skew, drift, and the fixed delay between the two nodes.In addition, complete steps for the algorithm required to find this MLE are also presented.Using the better model results in long-term synchronization between nodes and consequently saves a lot of energy in terms of considerably less frequent resynchronization through timing message communications.For future work, deriving the Cramer-Rao lower bound (CRLB) for the clock parameters derived in this paper is a very challenging open research problem.

Figure 1 :
Figure 1: Two-way timing message exchange mechanism involving N such pairs.

d 1 2 θDFigure 2 :
Figure 2: d as a function of θ D only.

Theorem 1 .
The MLE ( θ D , d, θ S , θ O ) is unique and is given by that intersection of two curves on the boundary of the support region in(12) where the termN r=1 {(T 2 4,r − T 2 1,r ) + v a (T 4,r − T 1,r )} −Nz b is negative for one curve and positive for the other.Proof.Figure 2 shows the fixed delay d as a function of clock drift θ D only which is in reality an intersection of 4D curves as explained above.The MLE ( θ D , d, θ S , θ O ) can be derived by the following observations.
) Since the boundary of the support region is formed by the curves in(12) in an order of decreasing slopes {z b }, the intersection where the sign of N r=1 {(T 2 4,r − T2  1,r ) + v a (T 4,r − T 1,r )} − N zb (and hence the sign of φ a,b ) changes from negative to positive occurs only once.Therefore, the MLE must be unique.(5)Let c = min a {v a } and s = {a}\c.Now comparing φ c,b and φ s,b on the boundary of the support region (see Figure 2) yields the following three options.(i) The signs of both φ s,b and φ c,b change at the same intersection of curves in (12).In this case, φ c,b < φ s,b since v c < v s .(ii) The sign change for φ s,b occurs at an intersection (say intersection 2 in Figure 2) of the curves in (12) to the right of the intersection (say intersection 1 in Figure 2) where the sign change for φ c,b occurs.This is not possible because for the same z b , φ s,b must have a sign change at or to the left of the intersection where the same occurs for φ c,b .(iii) The sign of φ s,b changes at an intersection of curves in (12) (say intersection 1 in Figure 2) which is to the left of the intersection where the sign change for φ c,b occurs (say intersection 2 in Figure 2).Now even on intersection 1, φ c,b < φ s,b since v c < v s .Due to the continuity of φ c,b (and hence the continuity of the likelihood function) on the support region, φ c,b can be further decreased by increasing θ D until it touches the intersection 2. Therefore, a = c should be used to find the index b corresponding to the minimum of the set φ c,b .

Finally, in the
light of above observations, by checking the sign of the expression N r=1 {(T 2 4,r − T 2 1,r ) + v c (T 4,r − T 1,r )} − Nz b for each b, we conclude that the MLE θ D can be expressed as

Figure 3 :
Figure 3: Mean square error (MSE) of θ D as a function of number of observations N.