Maxime Cohen
email: maxcohen@tx.technion.ac.il

Nahum Shimkin
email: shimkin@ee.technion.ac.il

Decentralized algorithms for sequential network time synchronization

Keywords: network clock synchronization, decentralized algorithms, Kalman filtering, recursive estimation. I

Accurate clock synchronization is important in many distributed applications. Standard algorithms, such as the Network Time Protocol (NTP), essentially rely on pairwise offset estimation between adjacent nodes. Some recent work introduced more elaborate algorithms for clock offset estimation, which take into account the algebraic constraints imposed on the sum of offsets over network cycles, using a least-squares framework. These algorithms are iterative and decentralized in nature, requiring several cycles of local communication among neighbors for convergence. In this paper, we extend this approach towards a sequential estimation framework, which allows to incorporate initial time estimates along with their uncertainty, as well as multiple rounds of pairwise measurements. We propose a decentralized implementation of the estimation algorithm that employs only local broadcasts and establish its convergence to the optimal centralized solution. We also present some simulation results to illustrate the performance benefits of the suggested algorithms.

INTRODUCTION

Accurate clock synchronization has been extensively studied and applied in the context of communication networks, from the Internet [START_REF] Mills | Internet time synchronization: the network time protocol[END_REF] to sensor networks [START_REF] Sivrikaya | Time synchonization in sensor networks: a survey[END_REF]. The task of synchronizing clocks in distributed systems is usually accomplished via the exchange of time-stamped messages (probe packets) between the distributed entities in order to coordinate their time. There is a large literature on how to synchronize clocks in traditional networked systems; among these, the "Network Time Protocol" (NTP) is the most widely accepted standard for synchronizing clocks over the Internet [START_REF] Mills | Internet time synchronization: the network time protocol[END_REF], [START_REF] Mills | Network Time Protocol (version 3) specification, implementation and analysis[END_REF]. This protocol essentially uses a so-called hierarchical approach by sending probe messages along a layered spanning tree of the network.

More recently, a novel approach for time synchronization termed CTP -Classless Time Protocol [START_REF] Gurewitz | Network time synchronization using clock offset optimization[END_REF] was proposed. This non-hierarchical approach exploits convex optimization theory in order to minimize a quadratic objective function of clock offsets. It was shown that CTP substantially outperforms hierarchical schemes such as NTP in terms of clock accuracy without increasing complexity. A related approach that relies on Least-Squares fit was proposed in [START_REF] Giridhar | Distributed clock synchronization over wireless networks: algorithms and analysis[END_REF], [START_REF] Solis | A new distributed time synchronization protocol for multihop wireless networks[END_REF]. The accuracy of clock synchronization was improved by exploiting global network-wide constraints (i.e., the relative offsets must sum up to zero over network loops). The central characteristic of these methods is the use of a distributed algorithm that requires only local broadcasts among neighboring nodes. The algorithms are iterative, and typically converge within a small number or rounds. The work in [START_REF] Barooah | Distributed estimation from relative measurements in sensor networks[END_REF] extends the same LS approach to a Weighted Least-Squares (WLS) framework, where each measurement may be assigned a different weight.

It is interesting to note, following [START_REF] Barooah | Distributed estimation from relative measurements in sensor networks[END_REF], that the time synchronization problem is mathematically equivalent to a general class of distributed estimation problems of additive quantities over a network, a class that includes sensor localization over a sensor network (in Cartesian coordinates).

The basic Least-Squares framework of [START_REF] Giridhar | Distributed clock synchronization over wireless networks: algorithms and analysis[END_REF] considered the offset estimation problem using only a single set of measurements. Our goal here is to extend this framework into a sequential estimation one, which handles prior estimates of clock offsets as well as multiple measurement sets. These goals can be cast in the framework of Kalman Filtering, which indeed readily provides an optimal centralized solution. However, as the obtained equations are not readily amenable to a decentralized implementation, we resort to the equivalent least-squares formulation, and employ localized least squares iterations to obtain a decentralized algorithm. This algorithm employs only local broadcasts between neighbors when the initial covariance matrix is diagonal.

A rich literature exists on distributed implementation of the Kalman Filter (KF), dating back to [START_REF] Hashemipour | Decentralized structures for parallel Kalman filtering[END_REF], [START_REF] Rao | Fully decentralized algorithm for multisensor Kalman filtering[END_REF]. More recent work, such as [START_REF] Olfati-Saberi | Distributed Kalman filtering for sensor networks[END_REF], focused on network structures and developed consensus-based algorithms that employ only local communication. This line of work allows measurement to be processed locally, however each node keeps a copy of the entire state vector. In our framework, each node keeps only data related to its own offset. Interesting work in this line is presented in [START_REF] Khan | Distributing the Kalman Filter for large-scale systems[END_REF], where decentralized approximation schemes are presented. Our goal here is to develop exact algorithms, for the more limited problem that we consider.

Finally, we present simulation results over several network topologies for evaluating and comparing the accuracy of the proposed time synchronization schemes. We provide several interesting comparisons, where the Kalman Filter approach outperforms the existing algorithms.

The paper is organized as follows. In section II, we describe the model and formulate the estimation problem in state-space form. In section III, we present an iterative, decentralized algorithm that combines prior offset estimates with a single set of measurements, and establish its convergence. While the development is carried out for a general initial state covariance matrix, the local communication structure is maintained only if this matrix is diagonal. Hence, a recursive extension of this algorithm to the case of multiple measurement sets is not apparent. The case of multiple measurement sets is taken up in Section IV, where an exact algorithm is developed. We also consider for comparison purposes a simple sub-optimal algorithm that neglects the off-diagonal terms of the inverse covariance matrix. Section V presents some simulation results that evaluate and compare the performance of the proposed algorithms. Finally, the conclusions and some notes on future directions are reported in section VI.

II. MODEL AND PROBLEM DEFINITIONS

We model the network as a directed graph ( )

, G V = Ε with N V = nodes denoted { } 1 2
, ,..., N Λ Λ Λ , and m = Ε edges. Each edge represents the ability to transmit and receive packets between the corresponding pair of nodes. The edge connecting nodes i Λ and j Λ is denoted by ij e . We assume that all the edges are bidirectional, so that ij e ∈ Ε implies ji e ∈ Ε , and that the network graph is connected, namely there exists a path between any pair of nodes in the network. Denote by i N the set of nodes which are the neighbors of i Λ , i.e., one edge away from node i Λ , and let | | i N be the number of such neighbors. We assume that each node in the network keeps a local clock, and our goal is to estimate the time offset of each clock with respect to some global reference. Without loss of generality, we may assign node 1 Λ as the reference time node, so that all clock offsets are estimated with respect to this node's clock.

A. Clock model

A standard model for the clock drift at a node follows the linear form: ( )

i i i T t t α τ = +
, where i α and i τ are the skew (rate deviation) and the offset parameters respectively, t is the real time (or reference time) and ( ) i T t is the local time at node i Λ . This model is known as the two parameters linear model (see [START_REF] Solis | A new distributed time synchronization protocol for multihop wireless networks[END_REF] and the references therein). As node 1 serves as the reference, we have by default 1 0 τ = and 1 1 α = .

We focus here on the simplified model where all clocks run at the same speed, so that there is no skew ( 1 i α = for all i ). This assumption is appropriate when the measurement time span is short so that the rate deviation is small relative to the skew.

B. Measurement model

We consider a two-way offset measurement scheme. Each network node ( , 1, 2,...,

i i N Λ =
) sends probe packets to each of its neighbors. Upon sending a packet m k the sender i Λ stamps the packet with its local time ( )

i m
T k , and the receiver j Λ stamps the packet upon receiving it with its local time ( ) j m R k . Then, node j Λ retransmits the packet back to the source with time stamp ( )

j m
T k , and the source stamps its local time ( ) i m R k when receiving the packet back. We thus obtain

( ) ( ) ( ) . ij j m i m ij m i j ij T R k T k D k τ τ ε ∆ - = -+ + Here ( ) ij m D k is the propagation delay over link ij e , ij
ε is an additive noise that represents the random queuing delay (and the other unknown influences) and j i τ τ

is the difference between the two clock offsets. Assuming that ( ) ( )

ij m ji m D k D k =
(symmetric propagation delay) we obtain:

( )

1 ˆ2 ij ij ji j i ij O T T v τ τ ∆ -∆ = -+ (1)
where

( ) / 2 ij ij ji v ε ε = - is the effective measurement noise.
We shall assume that ij v is zero mean, with covariance 0 ij r > .

Remark:

In practice, the effective measurement ˆij O may be obtained by processing several subsequent packet exchanges.

One way to combine these successive measurements is to consider their average. More refined methods such as a minimum filter may favor measurements that reflect smaller delays, to account for changing congestion conditions (cf. [START_REF] Mills | Internet time synchronization: the network time protocol[END_REF]). In any case, the combined measurement is taken to be of the above form.

C. State-space equations

Our objective is to synchronize all the clocks in the network with the reference time. This is equivalent to estimating the offset i τ at each network node. We start by formulating the statistical estimation problem in state space form. Define the state as the column vector ( ) 2 ,...

T N

x τ τ (recall that 1 0 τ = by definition). As x is assumed to be constant for the time frame of interest, we can write ( 1) ( ) x n x n + =

, where 0 n ≥ is the time (or step) index. The initial state (0)

x is assumed to have known first and second order statistics:

[ ] 0 (0) E x x = , [ ] 0 cov (0) x P = .
This allows to take into account prior information about the accuracy of the initial clock offsets in different nodes. For example, the network may contain several reference nodes which keep accurate time, which translates to small covariance entries in 0 P .

The effective measurement for each pair of neighboring nodes is given by equation ( 1). Thus, ˆij O measures the offset difference for these nodes, plus an additive noise ij v . The noise ij v is assumed to have zero mean, with finite covariance 0 ij r > , and is independent across node pairs. We collect the individual measurements in a column vector ( ,

) ij y O ij = ∈ Ε ,
which represents a single measurement set. We assume that a new measurement set ( ) y n is obtained at each step 1 n ≥ . We note that n need not refer to actual time, but rather corresponds to the epoch when the n-th measurement set ( ) y n becomes available.

To express the measurement equation in vector form, define the graph incidence matrix A % whose dimensions are N (nodes) × m (edges), and where in the row corresponding to node i Λ , we have an entry +1 for all edges ij e originating in i Λ , an entry -1 for all edges ji e terminating in i Λ , and 0 otherwise. For a connected graph, the rank of the incidence matrix is 1 N -. Thus, deleting any row from the incidence matrix yields a full row rank matrix, which is called the reduced incidence matrix. Here, we will work with the ( )

1 N m -×
matrix A obtained by deleting the row corresponding to the reference node 1 Λ .

The model may now be summarized in state space form:

( ) ( 1) ; 1 ( ) ( ) ( ) T x n x n n y n A x n v n = -  ≥  = +  with initial conditions [ ] 0 (0) E x x = , [ ] 0 cov (0) x P = . The measurement noise sequence { } ( )
v n is assumed to be a white noise sequence with zero mean and covariance ( ) 0 R n R = > . We further assume that there is no correlation between the measurement noise over different links, so that diag{ ,

} ij R r ij E =
∈ , as a diagonal matrix with elements 0 ij r > along the diagonal. Finally, the measurement noise and initial state (0)

x are uncorrelated.

Our goal is to estimate the offset vector ( )

x n x = based on the measurement sets (1), , ( ) y y n K and the prior information 0 0 , x P . This is of course a classical problem in sequential estimation theory. As is well known, the Kalman filter provides the optimal linear solution the MMSE (Minimal Mean Square Error) sense, which further coincides with the optimal (conditional expectation) solution under the Gaussian assumption. However, while the centralized KF equations are easily written, it is not readily seen how they may be distributed. We thus proceed to develop distributed algorithms that converge to the optimal KF solution.

III. SINGLE MEASUREMENT SET

We start by considering a single measurement update, namely the problem of estimating [START_REF] Barooah | Distributed estimation from relative measurements in sensor networks[END_REF] x x = based on the measurement (1) y y = and the prior information 0 0 , x P . As is well known (e.g., [START_REF] Jazwinski | Stochastic processes and filtering theory[END_REF], Section 5.3), the KF equations are equivalently obtained as a solution to a Least-Squares deterministic problem, which in our case reduces to the minimum of the following objective function:

1 1 0 0 0 ( ) ( ) ( ) ( ) ( ) T T T T J x x x P x x y A x R y A x - - = - - + - - (2)
The first term is related to the initial information regarding the clock offsets, whereas the second term is associated with the single set of measurements and its corresponding covariance matrix R .

In the development of a distributed algorithm, we will find it more convenient to manipulate the above deterministic LS problem rather than starting with the KF equations.

A. Baseline algorithm

We first review the existing algorithm introduced in [3], [START_REF] Gurewitz | Network time synchronization using clock offset optimization[END_REF], [START_REF] Solis | A new distributed time synchronization protocol for multihop wireless networks[END_REF]. In this case, the objective function is given by:

( ) 2 , ( ) ( ) ( ) i T T T ji i j i j j N J x y A x y A x O τ τ ∈ = - - = -+

∑

The first order optimality conditions are ( ) ( ) ( )

2 0 i T ji i j i i j N i J AA x A y O τ τ τ ∈ ∂ = - = - -+ = ∂

∑

From this, we obtain:

( )

1 î i ji j j N i O N τ τ ∈ = ⋅ + ∑ (3)
The above equation must be satisfied by the optimal solution of the offset estimation problem. While this is a set of linear equations, a direct solution cannot be carried out in a decentralized manner. Instead, a decentralized iterative algorithm was suggested and shown to convergee to the optimal centralized solution. This algorithm simply iterates the above equation, which can be simply interpreted as follows. Each node obtains from its neighbors their current estimates for their own clock offsets, and computes its offset estimate as the average of all its neighbors' estimates plus the corresponding relative measurements. This procedure is the same as in [START_REF] Giridhar | Distributed clock synchronization over wireless networks: algorithms and analysis[END_REF], [START_REF] Solis | A new distributed time synchronization protocol for multihop wireless networks[END_REF] and one can easily show that this is equivalent to the algorithm in [START_REF] Gurewitz | Network time synchronization using clock offset optimization[END_REF].

Our objective is to extend the previous result to a wider framework and we will obtain this procedure as a special case of a more general algorithm.

Next, we will consider the more general framework that includes the initial covariance matrix in the objective function in addition to a weighting matrix 1 R -. The analysis is divided in two cases: non-diagonal and diagonal initial covariance matrix.

B. Adding initial conditions and measurement weights

We now consider the distributed solution of equation ( 2), that incorporates the initial estimate 0

x with covariance 0 P , as well as possibly different weights for the measurement as expressed by the covariance matrix R . We develop an iterative algorithm, and establish its convergence to the optimal solution of (2). However, it will be seen that the algorithm is not truly distributed, in the sense that each node must receive data from non-neighboring nodes, unless the initial covariance matrix is diagonal. The latter case will be considered in the next subsection.

Proceeding similarly to above, the first-order optimality conditions for (2) are

( ) ( ) ( ) 1 0 1 1 1 ˆ(0) 0 i i N i ji j k k ik j N j N k i ji ji J O P r r τ τ τ τ τ - ∈ ∈ = ∂ = - + + - = ∂ ∑ ∑ ∑ or, equivalently, ( ) ( ) ( ) (
)

1 1 0 0 1 1 1 ˆ(0) (0) i N i ji j i m m ii im j N m i ji m i O P P I r τ τ τ τ τ - - ∈ = ≠     = + + - -       ∑ ∑ (4)
where ( )

1 0 1 . i i ii j N ji I P r - ∈ = +

∑

These equations motivate the following synchronous iterative algorithm for their solution:

( ) ( ) ( ) ( ) ( 1) ( ) 1 0 1 ( ) 0 1 1 1 ˆˆ(0) ˆ(0) [ ] i k k i ji j i ii j N i ji N k m m im m m i O P I r P τ τ τ τ τ + - ∈ - = ≠ = + + - - - ∑ ∑ K (5)
with initialization (0) ˆ(0), 2, 3,...

i i i N τ τ = = .
Here, 0 k ≥ is the iteration number. This algorithm can be interpreted as the classical Jacobi iteration for the solution of the linear equations ( 5), or, equivalently, as a local least-squares algorithm where each node minimizes at each iteration the objective function (2) over its own offset, given the current offset estimates of the other nodes.

The following convergence result is next established.

Theorem 1. Suppose that:

(1) The matrix 1 Ris diagonal and Positive Semi-Definite, that is:

( ) 1 0 , ji r i j - ≤ < ∞ ∀ .
(2) The initial covariance matrix 0 P is an M-matrix, namely:

( ) ( ) ( ) ( ) 1 0 1 1 0 0 0 0 0 ij j ii ij P P and P i j - - -  ≥    ≥ ≤ ≠  ∑ (3)
The clock adjustment operation in (5) is applied synchronously by all nodes ( 2,3,...

i N =
) in all iterations.

Then, the iterated estimators ( ) ˆ( ) 2,3,...

k i n i N τ = converge (as k → ∞
) to the optimal offsets that minimize the objective function in [START_REF] Cohen | Network Time Synchronization Using Decentralized Kalman Filtering[END_REF]. The proof is provided in the Appendix.

The main problem with the last iterative algorithm is that, in general, each node needs to communicate with all the other nodes and not only with its neighbors. Thus, each node is required to be aware of the global topology of the network, and the algorithm does not satisfy the requirement of local information exchanges only. Fortunately, this problem does not exist when initial covariance matrix 0 P is diagonal.

C. Diagonal 0 P

We henceforth specialize the discussion to the case where the initial covariance matrix 0 P (hence its inverse 1 0 P -) is diagonal, with diagonal elements ( ) i p . This will be the case when the initial estimates of clock skews are obtained by the different nodes independently. For example, some nodes may have a GPS receivers which allows then to obtain an accurate estimate of the time. Or, starting from an initially accurate estimate, each node has after some time an added uncertainty due to the estimated drift of its clock.

The decentralized iterative procedure in (5) now reduces to:

( )

( 1) ( ) (0) 1 1 ˆ1 1 i i k k i i ji j j N ji i j N ji i O r p r p τ τ τ + ∈ ∈   = ⋅ + +       + ∑ ∑ (6)
In words, i τ is computed at each iteration as a weighted average between the modified estimates obtained from adjacent nodes, and its prior estimate. Observe that if the matrix 1 0 Pis equal to zero and 1 R I -= , we obtain the equation ( 3) as in the basic LS case described above. This algorithm requires only local broadcasts between adjacent nodes. Evidently, the convergence result in Theorem 1 applies here as a special case.

IV. MULTIPLE MEASUREMENT SETS

In the previous section we have developed a decentralized algorithm for the case of a single measurement set, under the assumption that 0 P is diagonal. We next consider the case when multiple measurement sets become available sequentially, with the goal of presenting a recursive version of the previous decentralized algorithms. As its turns out, a simple recursive extension of the single-measurement case, in the style of the Kalman Filter, will not work here. The reason is that even if the initial covariance matrix 0 P is diagonal, which we assume, after the first iteration the state covariance matrix 1 P will not be diagonal any more (as can be easily verified, the estimates will become correlated). Thus, although the iterative procedure in (5) can be formally repeated, it will not be distributed any more, as each node has to communicate and store information related to all the other nodes over the network and not only with the one-hop neighbors.

Our objective is therefore to derive an alternative recursion that is suitable for decentralized implementation. This is done in the next subsection. For comparison purposes we also mention, in the subsequent subsection, an approximate algorithm that simply neglects the off-diagonal elements of the covariance matrix after each measurement update.

A. Optimal decentralized algorithm

For the multiple measurement update case, the equivalent Least Squares objective function is given by:

1 1 0 0 0 1 ( ) ( ) ( ) ( ( ) ) ( ( ) ) n T T T T k J n x x P x x y k A x R y k A x - - = = - - + - - ∑ (7)
For notational simplicity assume that the matrix 1 R -is identical for each set of measurements. We propose the following iterative algorithm: (

1 0 1 ˆˆ( ) ) ( ) ( ) ( ) ( 1) 1 ( ) ( ) ( )
} i i i i k n k i i i ji i j j N ji k j j j N ji i i ii j N j N ji ji n n I n O n n r n n n r I n I n P n n r r τ τ τ τ τ τ + - ∈ ∈ - ∈ ∈ = -+ ⋅ ⋅ - -- + -⋅ ⋅ - - = -+ = + ⋅ ≥ ∑ ∑ ∑ ∑ (8) for 2, 3,..., , i N = initialized with (0) ˆ(0) i i τ τ = , 1 0 (0) ( ) i ii I P - = . ( 1) ( ) [ ( 1) ( ) ] 1 ˆ1 [ ( ) ( 1) ] , 1 1 ( ) ( 1) , 1 {
The above set of equations is a decentralized, synchronous and recursive algorithm that computes at each step the estimated offsets and their corresponding error variances. The main advantage of this algorithm is its local nature; each network node needs to communicate only with its neighbors.

We now describe in words the iterative procedure in [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF]. At time n , we assume that the estimate of ˆ [START_REF] Barooah | Distributed estimation from relative measurements in sensor networks[END_REF] i n τis given.

Then, ( ) ˆ( ) 1, 2,...

k i n k τ = is computed based on ˆ( 1) i n
τand the last measurement set ( ) y n . We assume that a sufficient number of iterations is performed at each time n , so that the estimate ˆ( ) i n τ is accurate.

Remarks:

1. We point out that the suggested recursion slightly deviates from a standard recursive estimation scheme due to the presence of n (time or measurement count) as a factor in the recursion.

2. It may be shown that the elements of ( ) i I n are the diagonal entries of the inverse covariance matrix in the Kalman filter equations, namely 1 ( ) ( )

i n ii I n P - =
. In fact, the above iteration for ( ) i I n is the same as the iteration over the diagonal elements in the information form of the Kalman covariance update:

( ) ( ) 1 1 1 1 T n n P P AR A - - - + = + (9)
Observe however that we do not compute the non-diagonal elements of the inverse covariance matrix.

The proposed algorithm may be derived by differentiating ( 1) J n -and ( )

J n with respect to the offsets vector x and setting the partial derivatives to zero. The algebraic details (which can be found in [START_REF] Cohen | Network Time Synchronization Using Decentralized Kalman Filtering[END_REF]) are omitted since the procedure is similar to the previous case.

An alternative derivation, also presented in [START_REF] Cohen | Network Time Synchronization Using Decentralized Kalman Filtering[END_REF], may start with the KF equations in information form. However, the recursion in [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF] is not equivalent to the KF recursion, and is not readily seen from these equations.

We next address the convergence the set of equations in [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF] to the optimal centralized solution.

Theorem 2. Suppose that:

(1) Assumptions (a) and (b) from Theorem 1 hold. [START_REF] Cohen | Network Time Synchronization Using Decentralized Kalman Filtering[END_REF] The clock adjustment operation in [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF] is applied synchronously by all nodes ( 2, 3,...

i N =
) in all iterations, recursively for n sets of measurements.

(3) A sufficient number of iterations is performed after each measurement set n , so that ( ) ˆ( )

k i n τ converges to ˆ( ) i n τ .
Then, for each 1 n ≥ , the iterated estimators

( ) ˆ( ) 2,3,... k i n i N τ = converge (as k → ∞
) to the optimal offsets that minimize the objective function in [START_REF] Jazwinski | Stochastic processes and filtering theory[END_REF].

The proof is provided in the Appendix.

Next, we propose a simple sub-optimal algorithm for the case where multiple sets of measurements are available.

B. A sub-optimal decentralized algorithm

For the case where 0 P is non-diagonal, we obtained in ( 5)) that the estimated offset of node i Λ depends on all the other offsets and not only on those of its neighbors. One can consider the naïve sub-optimal algorithm that neglects the off-diagonal terms of the inverse covariance matrix:

( ) ( ) ( ) ( 1) ( ) 1 0 1 0 1 1 ˆˆ(0) 1 i i k k i ji j i ii j N ji ii j N ji O P r P r τ τ τ + - - ∈ ∈   = + +       + ∑ ∑
The decentralized sub-optimal recursive algorithm for the multiple measurement scenario is given by:

( ) ( ) ( 1)
( ) ( )

1 0 ˆ( ) ( 1)
1 1 ˆˆ( 1) ( ) 1 i i k i i n k ji i j j N ji ii j N ji n n O n n r P n r τ τ τ τ + - ∈ ∈ = -       + ⋅ ⋅ - --         + ⋅ ∑ ∑
Here, we may neglect the off-diagonal terms before inverting the information matrix ( )

1 1 n P - -
, in order to reduce the algorithm complexity. In this case, we will invert a diagonal matrix and hence the time computation will significantly decrease.

V. NUMERICAL RESULTS

In this section, we implement some of the algorithms that we previously developed for typical problems and we compare the results with the existing algorithms. More extensive comparisons can be found in [START_REF] Cohen | Network Time Synchronization Using Decentralized Kalman Filtering[END_REF]. The convergence rates of the decentralized algorithm are not presented here as convergence is achieved after a relatively small number of iterations and the results are very similar than [START_REF] Gurewitz | Network time synchronization using clock offset optimization[END_REF]. Consider two different network topologies:

-Network 1: a 400 node network with 997 edges.

-Network 2: a 170 node network with 1200 edges.

The first case we analyze is the one where 10% of the nodes are perfectly synchronized to the global time (through a GPS satellite receiver for example), and the remainder are not synchronized at all. Namely, for these arbitrary 40 nodes we take the initial variances to be very small (0.01) and the offsets equal to zero, and for the rest of the nodes, the variances tend to infinity and the offsets are randomly chosen according to a uniform distribution. The graphical comparison between the decentralized CTP algorithm (equation ( 3)) and the Decentralized Kalman Filter (DKF) (equation ( 6)) is presented in Fig. 2. As expected, the DKF algorithm outperforms the decentralized CTP method in terms of clock accuracy. Fig. 2 shows the fraction of nodes with clock offset with respect to the reference time node that is not grater than t for the different algorithms. In other words, the y-axis represents the fraction of nodes with clock offset, relating to the global time, not greater than the value described by the x-axis. The second part of this section is devoted to the comparison of the recursive Centralized Kalman Filter (CKF) algorithm to the Sub-Optimal Algorithm (SOA) that neglects the offdiagonal terms of the inverse covariance matrix (see section V. B). We consider the topology of Network 2 and we check several values of n (the number of measurements). The queuing delay is randomized in accordance with the Kalman Filter assumptions, namely normally distributed with zero mean and covariance matrix R :

[ ] ( ) 0.01,12 0, delay R U Q N R
In addition, we consider that 10% of the nodes are perfectly synchronized to the global time and the remainder are not synchronized at all (similar to the case in Fig. 2). In this analysis, we also compare the results to the Centralized Least-Squares (CLS) algorithm. Fig. 3 presents the results for the offsets obtained by applying the optimal CKF method, the SOA and the CLS algorithms for two different values of n . As expected, the optimal algorithm gives the best results. The suboptimal algorithm gives relatively poor results but reduces the complexity and is not diverging. Moreover, we obtained that the sub-optimal algorithm is even worse (in terms of clock accuracy) than the basic centralized Least-Squares method (that does not take into account the initial covariance matrix).

VI. CONCLUSION

In this paper, we have developed several decentralized algorithms for estimating the offset at each network node with respect to the reference time, utilizing a sequential estimation framework. The essential characteristic of these algorithms is their decentralized nature; each node can estimate its clock offset by exchanging packets with its one-hop neighbors only. We extend the existing Least-Squares based algorithms so that we may assign different weights to the measurements according to their accuracy, include a-priori information, and provide a recursive estimation scheme. The main algorithm is both decentralized (requires only local broadcasts), recursive (works in on-line applications) and converges to the optimal centralized solution. Finally, some numerical results were presented to show that, as expected, the proposed algorithm outperforms the existing methods.

We close the paper by mentioning several extensions of interest. A discount factor is easily incorporated into the objective function [START_REF] Jazwinski | Stochastic processes and filtering theory[END_REF] in order to give a higher weight to the more recent measurements, and leads to similar algorithms. This will be useful when the offsets are time-varying. The proposed algorithms may also be extended to handle dynamic changes in the communication topology by considering temporary link failures, following the treatment in [START_REF] Barooah | Distributed estimation from relative measurements in sensor networks[END_REF]. Further details related to these two issues may be found in [START_REF] Cohen | Network Time Synchronization Using Decentralized Kalman Filtering[END_REF]. One may also consider more elaborate state dynamics to model possible time variations in the clock offsets. The simplest is adding a white system noise in the state space model. Interestingly, the results of this paper are not easily extendable to this model. Another major issue is the incorporation of the clock skew parameter into the clock model (see section II. A). These scenarios were partially investigated in [START_REF] Cohen | Network Time Synchronization Using Decentralized Kalman Filtering[END_REF] and may be considered as directions for future research.

Figure 1 .

 1 Figure 1. Comparison between the decentralized CTP and DKF algorithms (with 10% of nodes synchronized via GPS) in Network 1.

Figure 2 .

 2 Figure 2. Comparison between CKF, SOA and CLS (with

Proof of Theorem 1: Let us recall that the general objective function is given by:

Let us analyze the convergence properties of the general case, where 0 P is not necessarily assumed to be a diagonal matrix. We recall that iteration (5) cannot be easily decentralized when 0 P is not diagonal as we previously explained. However, the iteration is still well defined mathematically. The synchronous iteration can be written in vector form:

(

) ( )

The optimal solution (equivalent to performing the centralized protocol) is given by: ( ) ( )

Then we obtain after some manipulations:

Thus, the convergence of the sequence ( ) ˆk τ to * τ is equivalent to the convergence of ( ) k τ to the zero vector, which is determined by the matrix M . The necessary and sufficient condition for this convergence is that the spectral radius of M is strictly smaller than 1. The following result is well known (see, e.g., [START_REF] Horn | Matrix analysis[END_REF], chapter 6).

Proposition 1.

Consider a non-negative square matrix A with the following properties: a) All the row sums of A are smaller or equal than 1. b) At least in one row this sum is strictly smaller than 1. c) The matrix A is irreducible (i.e., there exists a path from any node to any other node in the network). Then, ( ) 1

According to Proposition 1, ( ) 0 k τ → if the sufficient conditions apply to the matrix M . In order to show that the spectral radius of M is strictly smaller than 1, we will require that the matrix M is both non-negative and sub-stochastic (the row sums are smaller than one). The elements of the matrix M can be determined by inspection as the following:

Let us find the conditions for the row sums of the matrix M to be smaller than 1:

if and only if:

( )

In other words, we obtained that the necessary condition is that for each node i Λ , the row sum of the matrix 1 0 P -has to be non-negative.

Requiring that all the entries of the matrix M are nonnegative leads to:

Hence, we can write: ( ) ( )

The above requirement can be seen as a diagonal dominance condition over the matrix 1 0 P -. In the case that the node i Λ is adjacent to the reference node, the corresponding row sum of the M matrix is given by:

In the case that i Λ is not adjacent to the reference node, the corresponding row sum of the M matrix is given by: ( ) ( ) ( ) ( ) ( )

Hence, we have shown that at least in one row, the row sum of M is strictly smaller than 1. Actually, we proved that the iteration matrix verifies all the sufficient conditions for convergence. Namely, the row sums of the matrix M are less or equal than 1 (and at least in one row this sum is strictly smaller than 1), the matrix M is irreducible and all its entries are non-negative.

As a result, we proved the convergence of the decentralized algorithm to the optimal solution performed by the centralized Kalman Filter for the most general case.

To sum up, the convergence conditions are given by:

The case of 1 n = was treated in Theorem 1. Our proof relies on the following lemma.

Lemma 1.

Suppose that 0 P satisfies the convergence conditions of Theorem 1, namely 0 P is an M-matrix. Let n P be computed using [START_REF] Khan | Distributing the Kalman Filter for large-scale systems[END_REF], then n P is an M-matrix for all 1 n ≥ . Proof Equation ( 9) corresponds to the measurement update equation of the inverse covariance matrix of the KF. Recalling that the matrix 1 R -is assumed to be diagonal, let us analyze the properties of the matrix 1 T AR A -. For the reduced incidence matrix, we have:

M

Here, v is a vector with non-negative components. The structure of the matrix

The row sums are ( )

for each node i Λ that is not adjacent to the reference node. Moreover, if i Λ is adjacent to the reference, this sum is a strictly positive number. We conclude that if the a-priori inverse covariance matrix ( ) This lemma immediately implies the convergence of the recursive extension (for several measurement sets) of equation ( 5) to the optimal solution, where at each step, the new covariance matrix is computed according to [START_REF] Khan | Distributing the Kalman Filter for large-scale systems[END_REF]. Since the iterations in [START_REF] Kelley | Iterative methods for linear and nonlinear equations[END_REF] are equivalent to the procedure in [START_REF] Hashemipour | Decentralized structures for parallel Kalman filtering[END_REF], we obtain the claimed convergence in Theorem 2. ■


Y Léguillon

P Besnard

L Provino

A Monteville

D Méchin

D Trégoat

M Doisy

F.-X Launay

Phase sensitivity to axial strain of microstrustured optical silica fibers

Keywords: Microstructured silica fiber, optical fiber sensors, phase sensitivity, acoustic sensors

We compare, thanks to a Sagnac interferometer, the phase sensitivity to strain of different microstructured optical silica fibers (MSF) that we design and fabricate. Our results show that when a same elongation is applied to different MSF, the induced phase change is equal or lower than the one obtained for a standard fiber, showing no advantage on this parameter for sensing applications.

INTRODUCTION

Fiber optic hydrophone based on interferometric sensors have advantages over conventional piezoelectric ones. [START_REF] Kirkendall | Overview of high performance fibre-optic sensing[END_REF] However, this interferometric optical technology is complex and has a relatively big size. Distributed feedback fiber laser could be a solution to these inconvenients. [START_REF] Cranch | Fiber laser strain sensors: enabling a new generation of miniaturized high performance sensors[END_REF] These fiber laser sensors have already been used in multiplexed hydrophone arrays, [START_REF] Hill | Fiber laser hydrophone array[END_REF] even if improvements on the size of the mechanical system have to be made. Preliminary work [START_REF] Pang | Phase sensitivity of hollow-core photonic bandgap fiber to acoustic pressure[END_REF] seems to indicate that the phase sensitivity of an hollow-core photonic bandgap silica fiber to an axial strain due to acoustic pressure is improved because of the presence of air holes. The goal of this paper is to experimentally check whether MSF, on which we can potentially inscribe a Bragg grating, are much more sensitive to an axial strain or not than standard silica fibers. The experimental set-up consists in a Sagnac interferometer shown on figure 1, commonly used as a fiber optic gyroscope or acoustic sensor. [START_REF] Blin | Pickup suppression in sagnacbased fiber-optic acoustic sensor array[END_REF] Light from a broadband source is split into clockwise (CW) and counterclockwise (CCW) beams by a 3×3 coupler. The CW and CCW have asymmetric path as the CW beam crosses first the delay coil and then the fiber under test (FUT) while the CCW propagates through these different elements in the reverse order. Each FUT is wrapped around the same cylindrical piezoelectric transducer (PZT), cleaved and directly spliced to the standard fiber in the interferometer. The PZT is driven by an alternative voltage to periodicaly pull and release the FUT, thus allowing the intrinsic sensitivity of the fiber to an axial strain to be measured. The fiber-length modulation causes a change in the propagation phase. It follows, at a given time, a phase difference between the two counter-propagating beams, because they do not meet the PZT at the same time. The phase modulation is converted into a power modulation when the beams recombine at the 3×3 coupler. If an alternative voltage, modulated at the frequency Ω and applied to the PZT induces a phase modulation φ h cos(Ω), the output power of the interferometer is then given by 5

EXPERIMENTAL SET-UP

P (t) = P 0 1 + v • cos 2 • φ h • η • sin Ωt + Ωτ 2 + π 3 , ( 1
)
where P 0 is a power scaling factor, v the interference visibility, φ h the phase modulation amplitude, τ the fiber delay time, η = sin(Ωτ/2) the sensitivity of the Sagnac interferometer.

The 3×3 coupler is used in order to passively bias the interferometer close to the quadrature at the phase bias π/3. To avoid acoustic perturbations, the test part and the interferometer are acoustically isolated. A polarization beam splitter (PBS) and a Faraday rotator mirror (FRM) reduce polarization induced signal fading, which simplified measurements by stabilizing the interferometer. Polarization controllers were used to maximize the amplitude and visibility of the detected signal. Note that as the Sagnac is a common-path interferometer, it is insensitive to the source phase noise.

By interchanging the FUT, it is possible to compare the results obtained with a standard silica fiber and with different MSF.

THE MICROSTRUCTURED OPTICAL SILICA FIBERS

Microstructured fibers are nowadays well-known speciality fibers. Their design could bring several advantages such as a smaller effective area in order to reach higher intensities in the waveguide commonly to lower nonlinear threshold, or on the contrary large effective area while maintaining single mode operation, or the realization of a desired optical bandwidth, or lower losses to bends. Three MSF have been tested. Their scanning electron microscopy images are given on figure 2 and their design parameters in table 1 . They have not been manufactured especially for this study. They have been chosen among MSF singlemode at 1550 nm available at PERFOS. Nevertheless, their different design could give us an idea of their sensitivity to an axial strain.

EXPERIMENTAL RESULTS

The fiber delay coil was about 5 km long. Because light passes twice through it, the effective loop length was twice that of an unfolded coil, which is 10 km. It corresponds to a loop proper frequency of ≈ 10.3 kHz, which is the frequency to set the sensitivity η at 1.

The interferometer output power is converted into a voltage by the photodiode and is read onto an oscilloscope. A typical response observed for a fiber SMF 28 is given on figure 3. When the PZT is not modulated, we obtain a continuous voltage called V dc . When it is modulated by a voltage V P ZT cos(2πf p t)/2, we obtain "pseudo"-sines (see equation ( 1)), whose amplitude increases when V P ZT increases. The power received by the photodiode depends on the FUT because of splice losses. So, to compare the different fibers, we measure, for a given amplitude V P ZT , the peak-to-peak voltage V p-p of interferences normalised by the voltage V dc . Measurements for each fiber are presented on figure 4 where V p-p /V dc is plotted versus the applied voltage on PZT. We can observe that MSF are not more sensitive to strain than the standard fiber SMF 28. For small phase shift amplitude, which corresponds to the linear part of the curve 4(a), we have

V p-p V dc ≈ 4 √ 3 • v v + 2 • η • φ h , ( 2
)
where v is the visibility, η the sensitivity and φ h the phase modulation amplitude. On figure 4(b), we linearly fitted curves in order to obtain the slope for each fiber. Since we fixed both visibility v and sensitivity η to be equal to 1, we can determine the sensitivity φ h /V P ZT (in rad/V) to an axial strain for each fiber from equation (2). These sensivities are 0.125 rad/V for both SMF 28 and HF 148, 0.123 rad/V for HF 202 and 0.098 for HF 174. For a standard fiber, the relationship between the relative phase shift and the applied axial strain z is given by 6

∂φ φ = [1 -P e ] • ε z ( 3
)
where P e = (n 2 /2)[p 12 -ν(p 11 -p 12 )] is the effective strain-optic constant, n is the effective index, p ij are the Pockel's coefficients of the stress-optic tensor and ν is Poisson's ratio. With the PZT used here, and considering P e for standard silica fibers 6 equal to 0.22, we should measure a sensitivity φ h /V P ZT of ≈0.117 rad/V. Our results are in agreements with theory. However, the sensitivity of the HF 174 is lower than others. As a matter of fact, the high birefringence of the fiber implies that not all the optical power is on the same birefringent axis even if a polarization controller is used. As in our set-up, the normalization is set over the total intensity, while the measurement concerns the axial strain along with an unique birefringent axis, the phase sensitivity is under estimated. We think that if we could align properly the polarization axis, the results would have been similar to the other fibers. It shows the limitation of our bench when highly birefringent fibers are under test.

Note that the use of all polarization maintaining fibers and components will not improve anything as far as the broadband source is not linear polarized. However, we made up the choice of showing the obtained results for such birefringent MSF in order to give an order of magnitude of phase sensivity. For MSF, if we consider more in details the different parameters of P e , this study confirms a posteriori that P e cannot signicantly be changed even with the presence of air-holes.

CONCLUSION

In this communication, we have experimentally shown that microstructured silica optical fibers have the same sensitivity than standard fibers when an axial strain is applied. These results can be summed up by the fact that the induced phase change for a given strain is, to the first order, determined by the silica material elasto-optic properties. However, future prospects with MSF are the use of radial stress or the use of multimaterial fiber taking benefits of the aire holes. Another idea is to expose MSF directly to acoustic pressure. The induced axial strain to MSF (thus induced phase change) will be as large as the air-filling ratio.

Figure 1 .

 1 Figure 1. Sagnac interferometer. PC: Polarization controller; PD: Photodiode; PBS: Polarization beam splitter; FRM: Faraday rotator mirror; PZT: Piezoelectric transducer; CW: clockwise; CCW: counterclockwise.

Figure 2 .

 2 Figure 2. Scanning electron microscopy images of microstructured silica fibers.

Figure 3 .

 3 Figure 3. Example of results with standard fiber SMF 28.

  (a) Set of measures (b) Zoom on linear part

Figure 4 .

 4 Figure 4. Sensitivity to strain; comparison of the 3 fibers.

Table 1 .

 1 Design parameters of the microstructured silica fibers. d is the air-holes diameter. Λ is the pitch between air-holes. Nr is the number of holes rings.

	Fiber d (µm) Λ (µm) N r
	HF 148	≈ 1.8	≈ 6.8	5
	HF 174	-	≈ 2	5
	HF 202	≈ 2	≈ 6	2

ACKNOWLEDGMENTS

The authors want to acknowledge S. Blin for all his advices and discussions. The authors thank the French territorial and governmental organizations that partially funded this work through the "Sea Innovation & Business cluster (Pôle Mer) ATOS project.