CONTENTS

Acknowledgement

Most of all, I would like to thank Prof. Alain Lambert, Prof. Thomas Nowak, and Prof. Michel Kieffer who supported my work and helped me get results of better quality. I also like to thank the reviewers and the juries. With their help, I have better knowledge of the domain of wireless sensor network and related algorithms. They are the best models for my future career. I would like to thank my girl friend Dr. Min-I Lee. I could not finish this long journey without hers accompany. Studying in Paris-Sud University with her is the best decision in my life. I would like to thank my friends for the support and being my travel companions. Last but not the least, I would like to thank my family: my parents, my brother, aunts and uncles for supporting me spiritually throughout writing this thesis and my my life in general.

Chapter 1 Introduction

Wireless sensor networks (WSN) consist of many individual low-cost and low-power sensor nodes, which achieve goals such as environmental monitoring by running micro electrical chips and detection sensors. Due to wide applications in environmental monitoring, health care, and the IoT embed WSN, WSN has caught researcher's attention and been developed rapidly [START_REF] Han | A nonintrusive power supply design for self-powered sensor networks in the smart grid by scavenging energy from ac power line[END_REF].

In many applications, sensors in WSNs observe some physical phenomenon and gather the observations to an evaluate estimation of the original physical phenomenon. This process is known as data fusion [START_REF] Huang | Tree-Clustered Data Gathering Protocol. for wireless sensor networks[END_REF]. In order to perform this task, the ability of sensors to estimate their physical location in some well chosen frame and the knowledge of time for WSN to work under a universal synchronized clock is essential. Synchronization among sensors is usually more important than absolute synchronization with some reference high-accuracy clock. For instance, real-time controller systems usually assume that events occur at predetermined time instants, and that require sensors sampled at the assumed rate together, otherwise the calculations might be wrong. This necessity is reflected by the development of a significant number of algorithms for localization and time synchronization for sensor networks in the recent past.

For example, many sensors can only detect the proximity of an observed object. Higher-level information, such as speed, size, or shape of an object can then only be CHAPTER 1. INTRODUCTION obtained by correlating data from multiple sensor nodes whose locations are known. The shape or the velocity of a mobile object, for example, can be estimated by the ratio of the spatial and temporal distances between two consecutive object sightings by different sensor nodes. For example, a sensor node can detect measurements of a widespread object, and the size and shape can be approximated by the information of the sensor nodes.

In WSN, many applications are based on the ability of sensors nodes to estimate the physical location and to synchronize the time with each other. Thus, a significant number of algorithms for localization and clock synchronization for WSN were studied in the recent years. However, due to the complexity and dynamicity of the WSN, robust algorithms that can tolerate inaccurate measurements during the communication and computation of WSN still have a tremendous demand.

For the synchronization problem, sensor nodes synchronize the clock value to perform data fusion algorithms or speed computation [START_REF] Sundararaman | Clock synchronization for wireless sensor networks: a survey[END_REF][START_REF] Solis | A new distributed time synchronization protocol for multihop wireless networks[END_REF][START_REF] Su | Time-diffusion synchronization protocol for wireless sensor networks[END_REF]. The local physical clock is driven so that the clock values are synchronized among sensors. A straightforward solution is to use a source node to provide a coordinate universal time. This solution is not suitable for WSN because the location of the sensor nodes is uncertain, so a source node may not be able to send the message to all the nodes in the WSN and the propagation delay is uncertain. Moreover, this solution is not robust, because a failure of the source node can affect all the sensors in the network.

Considering the dynamicity, a distributed algorithm is proposed, which extend the previous work from Függer et al. [2015b] and reduce the clock skew. In the new pulsecoupled synchronization, time-wheel algorithm is introduced for sensors to exchange the round number and averaging the pulses of the same round. Moreover drift compensation is applied to reduce the clock skew. The convergence analysis shows that the drift compensation reduces the clock skew and improves the convergence rate. When considering a highly dynamic WSN, the algorithm converges the software clock faster then in previous works. The proposed algorithm is successfully applied to a vehicular network.

For the localization problem, sensors can try locating each other possibly with the help of some anchor sensors in order to obtain a position estimation. A simple method is to put a GPS receiver on each sensor [START_REF] Yick | Wireless sensor network survey[END_REF]. With GPS receivers, sensors can receive information from satellites and estimate the geographical position themselves. However, GPS receivers drain the battery very quickly, so usually, only a small portion of sensors in WSN have GPS receivers. For some applications, GPS signal may not be received by the sensors (e.g., indoor, underwater scenarios).

For the second part of the thesis, we estimate the location of each sensor using LSCR, which produces a non-asymptotic confidence region. The benefit of applying LSCR to localization problem is that the noise assumption on the noise of RSS signal is mild, so the confidence region can be estimated accurately. We analyze different correlation functions of LSCR and found that only some correlation functions can produce a well-shaped confidence region. We also implement alternative localization algorithms, involving, e.g., a maximum likelihood estimator (MLE), semi-definite programming (SDP), bounded error estimator (BE), robust bounded error estimator (RBE). We compared the alternative methods and found that LSCR provides the best localization results: its average error is smaller than other alternative methods.

In this thesis, we divide our work into two parts for viewers to understand the backgrounds and methodologies for both clock synchronization and localization. In the first part, the clock synchronization is discussed, and in the second part, the wireless sensor localization is discussed.

For clock synchronization part, Chapter 3 introduces a formal model of the wireless sensor network and specify the clock synchronization problem. In Chapter 4, the methodology used in this work is shown with the analysis of the clock skew and convergence rate with mathematical derivation, which is followed by some experimental results in Chapter 5, and finally, the conclusions of the survey and the direction of future work are given in Chapter 6.

For the localization part, in Chapter 8, tools from Interval analysis is described, which is used by some of the localization methods. With the knowledge of Interval analysis, Chapter 9 introduces RSS localization model and specify several localiza-tion approaches. In Chapter 10, the methodology used in this work is shown and compared with other localization approaches. The different correlation functions of LSCR are also simulated and analyzed. Finally, the conclusions of the survey is given in Chapter 12.

Part I Clock synchronization

Chapter 2

Introduction

This chapter presents the importance of the clock synchronization problem. Previous works and proposed techniques are explained briefly. The chapter ends with a description of the structure and contents of clock synchronization.

Motivation

Clock synchronization is a critical issue in sensor networks [START_REF] Ganeriwal | Timing-sync protocol for sensor networks[END_REF][START_REF] Saad Ahmed Munir | Mobile wireless sensor network: Architecture and enabling technologies for ubiquitous computing[END_REF][START_REF] Affan | Time synchronization for high latency acoustic networks[END_REF][START_REF] Sun | Fault-tolerant cluster-wise clock synchronization for wireless sensor networks[END_REF]. In particular, it is a cornerstone of data fusion, where algorithms work in a better way when data are precisely synchronized rather than only time stamped. Distributed data fusion of nonsynchronized data necessitates taking into account the shift between the time of the measurement. Such shifts increase the data fusion imprecision results [START_REF] Sundararaman | Clock synchronization for wireless sensor networks: a survey[END_REF][START_REF] Solis | A new distributed time synchronization protocol for multihop wireless networks[END_REF][START_REF] Su | Time-diffusion synchronization protocol for wireless sensor networks[END_REF]. A quite large number of clock synchronization protocols, mostly for wired networks, have been developed over the past few decades. However, most of the protocols are unsuitable for wireless sensor networks due to the dynamicity, low bandwidth, and often low-powered computing nodes [Elson andRömer, 2003, Ganesan et al., 2004].

On the other hand, the hardware of the sensors are usually very restricted (lim-CHAPTER 2. INTRODUCTION ited storage, slow computation) since the size of the sensors is small [START_REF] Bertrand | Applications and trends in wireless acoustic sensor networks: A signal processing perspective[END_REF][START_REF] Sundararaman | Clock synchronization for wireless sensor networks: a survey[END_REF]. Energy consumption of wireless sensor networks can be high when service constraints are high [START_REF] Ian F Akyildiz | Wireless sensor networks: a survey[END_REF][START_REF] Sundararaman | Clock synchronization for wireless sensor networks: a survey[END_REF]. The energy is usually provided by the battery with limited power for sensors, so synchronization must be achieved while preserving energy [START_REF] Ye | An energy-efficient MAC protocol for wireless sensor networks[END_REF]. For those sensors, the limit on bandwidth directly influences message exchanges among sensors [START_REF] Kamal | Routing techniques in wireless sensor networks: a survey[END_REF]. During the message delivery, the flux of intensity of electromagnetic radiation is reduced through the medium. Therefore, the message may not be delivered from the sender to the receiver if their distance is long [START_REF] Amulya | A model for the classification and survey of clock synchronization protocols in wsns[END_REF]. In addition, unlike wired networks, wireless sensor networks are dynamic. Sensors may move and communication links may change. In view of this, traditional clock synchronization methods for wired networks [START_REF] David | The network computer as precision timekeeper[END_REF][START_REF] Elson | Fine-grained network time synchronization using reference broadcasts[END_REF] are often not applicable for wireless sensor networks.

Background

Most clock synchronization algorithms for wireless sensor networks are packet-based methods, where sensors exchange digital data to agree on a common time base. Centralized approaches [START_REF] David | Internet time synchronization: the network time protocol[END_REF][START_REF] Cristian | Probabilistic clock synchronization[END_REF][START_REF] Kopetz | Clock synchronization in distributed real-time systems[END_REF] distribute the clock of a master node to other nodes and are not robust to any failure of the master node. In the light of this, an average consensus-based protocol is proposed in [START_REF] Schenato | Average timesynch: A consensus-based protocol for clock synchronization in wireless sensor networks[END_REF], where sensors exchange clock readings and update its clock drift and offset by averaging the value of its neighbors. A maximum consensus-based protocol based on [START_REF] Schenato | Average timesynch: A consensus-based protocol for clock synchronization in wireless sensor networks[END_REF] is proposed in [START_REF] He | Time synchronization in wsns: A maximum-value-based consensus approach[END_REF] to improve the convergence rate. A second-order consensus method is proposed in Carli and Zampieri [2014] to reduce the clock skew further. For all the above packet synchronization [START_REF] Schenato | Average timesynch: A consensus-based protocol for clock synchronization in wireless sensor networks[END_REF][START_REF] He | Time synchronization in wsns: A maximum-value-based consensus approach[END_REF][START_REF] Bibliography | Network clock synchronization based on the second-order linear consensus algorithm[END_REF], nodes exchange data information (global time, frequency, or system parameters) between its neighboring nodes and agree on some parameters of the model to have a common notion of time. In large sensor networks, 2.2. BACKGROUND message delivery might be unreliable due to both dynamicity and sparsity of these networks [START_REF] Wu | Clock synchronization of wireless sensor networks[END_REF]. Most wireless sensor networks usually have stationary sensor nodes. The network topology is usually fixed, so sensors can transmit packets without failure. However, for a sensor network, all the sensors move rapidly, which makes it difficult to apply a common clock synchronization protocols such as [START_REF] Schenato | A distributed consensus protocol for clock synchronization in wireless sensor network[END_REF], [START_REF] Elson | Fine-grained network time synchronization using reference broadcasts[END_REF]. Also, during intensive transmission, sensors suffer from the packet collision and delay [START_REF] Pešović | Effect of hidden nodes in ieee 802.15. 4/zigbee wireless sensor networks[END_REF]. When a sensor tries to send packets, it detects and waits until the medium is not occupied by the other sensors, so the information transmission might be postponed. Moreover, if sensors are far away, the sender cannot detect whether the medium is occupied near the receiver, so the packet collision might happen when two senders transmit packets to the receiver simultaneously.

An alternative to the packet-based synchronization approaches, pulse-coupled synchronization, is proposed in [START_REF] Li | Global clock synchronization in sensor networks[END_REF], where nodes only exchange pulses instead of packets. The convergence properties of [START_REF] Li | Global clock synchronization in sensor networks[END_REF] is studied in [START_REF] Simeone | Distributed time synchronization in wireless sensor networks with coupled discrete-time oscillators[END_REF] through linear algebra with different communication graphs. A different topology (scenario) of the network is investigated, and a practical implementation is discussed. Also, the dynamic and bi-directional properties of pulse-coupled synchronization are analyzed in Függer et al. [2015b]. The algorithm solves clock synchronization as long as the topology of the dynamic network contains a directed spanning tree. Idle listening is introduced in [START_REF] Wang | Energy-efficient pulse-coupled synchronization strategy design for wireless sensor networks through reduced idle listening[END_REF] to reduce total energy consumption.

Pulse-coupled synchronization [START_REF] Li | Global clock synchronization in sensor networks[END_REF][START_REF] Simeone | Distributed time synchronization in wireless sensor networks with coupled discrete-time oscillators[END_REF], Függer et al., 2015b[START_REF] Wang | Energy-efficient pulse-coupled synchronization strategy design for wireless sensor networks through reduced idle listening[END_REF] can be viewed as the first-order phase locked loops [START_REF] Simeone | Physical-layer distributed synchronization in wireless networks and applications[END_REF], which suffer from carrier frequency offset between sensors resulting in a constant phase lag. Moreover, pulse frequency synchronization is studied [START_REF] Bibliography | Distributed frequencylocked loops for wireless networks[END_REF].

Objectives

The primary objective of the part I is to improve the limitations of the current state-of-the-art techniques.

The specific objectives of this thesis are as follows.

• To implement a clock synchronization protocol is able to work well with any number of nodes in a dynamic network.

• To implement a decentralized clock synchronization protocol, which tolerates communication loss of the sensors in the network.

• To develop a new communication method to improve precision (clock skew) of the clock synchronization.

• To relax the assumption of pulse-coupled clock synchronization with more robust information exchange protocol.

• To implement an energy efficient clock synchronization protocol that only performs simple calculation instead of performing energy-consuming optimization.

This work extended and improved the algorithm of the previous time synchronization protocols [Függer et al., 2015b] based on theory and experiment results in clock skew and some other aspects.

Outline

Chapter 3 introduces a formal model of the wireless sensor network and specify the clock synchronization problem. In chapter 4, the methodology used in this work is shown with the analysis of the clock skew and convergence rate with mathematical derivation, which is followed by some experimental results in Chapter 5, and finally, the conclusions of the survey and the direction of future work are given in Chapter 6.

Chapter 3

Formal Model and Problem Formulation

Introduction

This chapter describes the formal model of the wireless sensor network and specifies the problem of clock synchronization. First, the local physical clock will be introduced. Second, the formal model of the wireless sensor network is defined. In the end, an error function is presented to evaluate the quality of the clock synchronization algorithm. Also, the existing state of the art techniques which deal with this problem are described.

Formal Model and Problem Specification

The goal of clock synchronization is to synchronize the time of all sensors because physical clocks of each sensor drift even if they have the same frequency, and these frequencies might even change over time. Figure 3.1 shows the physical clocks of sensor i and j. The frequency of the physical clocks for i and j are different due to clock drift or hardware setup, so the time period of pulses are different as well.

Figure 3.1: The physical clocks of the sensor i and j without any corrections. The frequency of the physical clock for both i and j can be different due to the clock drift or hardware setup, so the time period of pulses are different. The time difference of the ticks between sensors is evaluated using clock drift and clock skew. The clock drift is the difference of clock periods between sensors and the clock skew is the largest difference of time between pulses of sensors.

The quality of the clocks is evaluated using clock drift and clock skew. The clock drift is the difference in clock periods between sensors. The clock skew is the largest difference of time between pulses of sensors. The goal is to reduce the clock skew and the clock drift as much as possible.

Formal Model of Wireless Sensor Network

A finite set N = {1, 2, . . . , n} of sensors communicating by message passing is considered. They are endowed with imperfect physical local clocks evolving in the time-base

T = [0, 1) of nonnegative reals.
The sensor network is modeled by a set of graph G t = (N , V t ), t 2 T , where N is the set of sensors and V t is a set of directed edges. The relation (i, j) 2 V t means that sensor i can send a message to sensor j at time t. Mobility is modeled by having a (possibly different) communication graph of different time instants.

FORMAL MODEL AND PROBLEM SPECIFICATION

A local algorithm for a sensor comprises a set of local states S, a set of events E, and its state transition function f (s, e). The local state is the internal variables for each sensor. The local state is initialized at s 0 2 S. Local state changes only when an event e 2 E occurs.

An execution of an algorithm in scenarios G t is a sequence of events (e q ) q 0 , which are triggered by sensors i q 2 N with its present states s q 2 S at a specific time t q 2 T . The value of q is a sequence number of the event in the system. Thus, the first event is e 0 and the second event is e 1 and so on. The sequence of events is sorted by t q . At the beginning of the execution, every sensor triggers its first event with its initial state s (q) 0 at time t q = 0. There are three types of event e b , e m , e t ⇢ E. The event e b represents the beginning of the local algorithm. The event e m indicates that a sensor receives a message.

The event e t is triggered when the timer in a sensor reaches a certain value (sensors can schedule this event when executing state transition functions). A sensor i determines its next local state s 0 using the state transition function f (s, e). A sensor cannot access the global time, only the physical local clocks of the sensor are obtained.

We assume all the sensors have the following features: A sensor i can send a message to sensor j at time t if and only if (i, j) 2 V t . A sensor i produces and sends a message using its sending function: S ! M, where the S is the local state, and M is the message. Sending function maps the state s k to message content. A sensor j receives messages through a receive function when a message event e m occurs. The time difference between sending and receiving a message is assumed to be zero. The assumption is reasonable if the message contains only a few bits. Notice that, the noise caused the transmission delay leads to larger clock skew in the wireless sensor network. However, with the pulse synchronization, sensors only exchange few bits, so the transmission delay is neglectable.

A sensor i is equipped with a local physical clock with the clock period µ i (c), where µ i is a function, which subsumes manufacturing defects, hardware differences, and different temperatures for quartz oscillators. The variable c is a local physical counter, which is initiated at 0 and increases over time whenever the local physical clock ticks. For simplicity of notation, assume that the nominal real time duration of a local physical clock tick is equal to 1. With c, sensors read time t i as

t i (c) = c X m=0 µ i (m). (3.1) Assuming that (1 %)µ  µ i (c)  (1 + %)µ
, where % represents the bound of the clock period and µ is the manufacture specification of the clock period, one gets

(1 %)µc  t i (c)  (1 + %)µc, (3.2)
which indicates that the uncertainty becomes larger through time when c and t i increase.

The main problem for clock synchronization to solve is that c i for each sensor does not increase at the same pace. Therefore, in this framework, a sensor i is able to calculate and exchange the round number k i . The round number k i for sensor i can be viewed as a software clock for each sensor. The value of k i starts from zero and increases periodically for each sensor i. With state transition functions, a sensors are able to broadcast messages and try to increase k i in the same pace.

Evaluation Function

Sensors synchronize the value of k i , i 2 G t and let k i increase periodically at a similar pace using the state transition functions: f (s, e b ), f (s, e m ), and f (s, e t ). To evaluate the precision of an algorithm, the maximum skew of rounds between sensors is taken into account:

(k) = max i 1 ,i 2 2Gt |⌧ i 1 (k) ⌧ i 2 (k)| , (3.3)
where ⌧ i (k) represents the time instant when the value of k i changes to k. Maximum clock skew of each round cannot reveal the accuracy of the algorithm, which needs

STATE-OF-THE-ART

to consider the errors between the clock ticks and the real time.

Equation 3.3 evaluates the precision of the algorithm. Sensors converge the software clock. If the software clock is identical, then (k) is zero.

State-of-the-art

This section describes the state of the art of clock synchronization including Average TimeSynch Protocol (ATS) [START_REF] Schenato | Average timesynch: A consensus-based protocol for clock synchronization in wireless sensor networks[END_REF], Robust Average TimeSynch Protocol (RoATS) [START_REF] Garone | Clock synchronization protocol for wireless sensor networks with bounded communication delays[END_REF], and Average time synchronization (ATSP) [START_REF] Wu | Average time synchronization in wireless sensor networks by pairwise messages[END_REF].

Average TimeSynch Protocol (ATS)

In ATS [START_REF] Schenato | Average timesynch: A consensus-based protocol for clock synchronization in wireless sensor networks[END_REF], Sensor i is equipped with a local physical clock c i (t) 2 R, which is similar to the physical clock counter c but in real value.

c i (t) = ↵ i t + i , (3.4)
where ↵ i and i are the frequency and offset of the local physical clock, respectively. Both variables ↵ i and i might be different for each sensor due to the imperfection of the manufacturing, but the values are fixed during the synchronization. To correct this difference of the local physical clock, each sensor i also has a software clock

ĉi (t) = ↵i (t)c i (t) + ˆ i (t), (3.5)
where ↵i (t) and ˆ i are tunable variables. With those variables, sensor i adjusts the variables and converges its software clock with other sensors to achieve clock synchronization.

During the ATS time synchronization [START_REF] Schenato | Average timesynch: A consensus-based protocol for clock synchronization in wireless sensor networks[END_REF], each sensor sends a packet periodically. At the k-th updates, the sensor i receives a packet from its neighbors j, which is a tuple (ID j , ↵j (t k ), ˆ j (t k ), c j (t k )). The parameter ID j CHAPTER 3. FORMAL MODEL AND PROBLEM FORMULATION is denoted as sensor j 0 s identifier, which is a unique value for each sensor. The parameters ↵j , ˆ j and c j are the latest value of ↵j (t), ˆ j (t), and c j (t), respectively at the time instant when sensor j sends the packet. The packet is then saved to local storage with the current local physical clock reading c ij (t k ).

if k > 1, sensor i performs drift estimate to calculate a relative drift ↵ ij . The Drift and offset compensation is also applied to update the value of both ˆ i (t) and ↵i (t), respectively.

Drift estimation and compensation

During drift estimation, a relative drift ↵ ij (t k ) of sensor j to sensor i is estimated using the following equation:

↵ ij (t k ) = c j (t k ) c j (t k 1 ) c ij (t k ) c ij (t k 1 ) , (3.6)
where c ij (t k ) is the local clock time when sensor i receives the packet from sensor j.

With ↵ ij (t k ), the value of ↵i is updated as follows:

↵i (k) = ⇢ ↵ ↵(k-1) + (1-⇢ ↵ )↵ ij (t k )↵ i (k-1), (3.7)
where ⇢ ↵ 2 (0, 1) is a user-defined parameter that adjust the convergence rate.

Offset compensation

In the offset compensation, the offset i is updated using the difference of ĉj (t k ) and

ĉi (t k ). ˆ i (k) = ⇢ ˆ (k 1) + (1 ⇢ ) [ĉ j (t k ) ĉi (t k )] (3.8)
Assume that ↵ i and i remain constant during the synchronization process, the transmission delay between sensors is zero, and the scenario is fully connected, then the software clocks ĉi (t) in the wireless sensor network eventually converge.

ATS shows an elegant iterative method to converge the software clocks in the wireless sensor network. However, it does not take into account the transmission delay during the communication of the sensors. When considering the transmission delay, original ATS [START_REF] Schenato | Average timesynch: A consensus-based protocol for clock synchronization in wireless sensor networks[END_REF] cannot ensure the convergence of the software clocks as the incremental error is accumulated during the iteration process.

Robust Average TimeSynch Protocol (RoATS)

Robust Average TimeSynch Protocol from [START_REF] Garone | Clock synchronization protocol for wireless sensor networks with bounded communication delays[END_REF] is proposed to improve the ATS method by adding an extra information exchange between sensors. In RoATS, sensors take an extra step to exchange ↵ ij and ↵ ji , so that the clock drift ↵i (k) and ↵j (k) is updated symmetrically to each other for convergence.

The main advance of RoATS is the drift compensation. To describe the modification of RoATS to ATS, let ↵ ⇤ ij (t k ) be the relative drift which considers noise. Equation 3.7 is then reformulated as

↵i (k) = ↵(k-1) + (1-⇢ ↵ ) ⇥ ↵ ⇤ ij (t k )↵ i (k-1) ↵(k-1) ⇤ , (3.9) with ↵ ⇤ ij (t k ) = ✏ ij (t k )↵ ij (t k ), (3.10)
where

✏ min = 1 1 + ij,max T min  ✏ ij (t k )  1 1 ij,max T min = ✏ max . (3.11)
The variable ij is the transmission delay, and T is the update period. To consider noise during transmission, the bound of ↵ ij (t k ) and ↵ ji (t k ) can be derived using Equation 3.11 and Equation 3.10.

↵ ⇤ ij (t k ) ✏ max  ↵ j (t k ) ↵ i (t k )  ↵ ⇤ ij (t k ) ✏ min , (3.12) CHAPTER 3. FORMAL MODEL AND PROBLEM FORMULATION and ✏ min ↵ ⇤ ji (t k )  ↵ i (t k ) ↵ j (t k )  ✏ max ↵ ⇤ ji , (3.13)
from which it follows

↵ ⇤ ij (t k )  ↵ ij (t k ) = ↵ i (t k ) ↵ j (t k )  ↵ ⇤ ij (t k ), (3.14)
where

↵ ⇤ ij (t k ) = max ⇢ ↵ ij ✏ max , ✏ min ↵ ji , (3.15)
and

↵ ⇤ ij (t k ) = min ⇢ ↵ ij ✏ min , ✏ max ↵ ji . (3.16)
With the error bound calculated, the correction direction and magnitude are computed locally. The direction sign s(k) is used to determine the convergence direction. The range of direction sign s(k) is ( 1, 0, 1), which is designed as follows:

s(k) = 1 2 sign ⇥ ↵ ⇤ ij (t k )↵ i (t k 1 ) ↵(t k-1 ) ⇤ + 1 2 sign ⇥ ↵ ⇤ ij (t k )↵ i (t k-1 ) ↵(t k-1 ) ⇤ . Note that if the correction term ↵ ⇤ ij (t k )↵ i (t k-1 ) ↵(t k-1
) is close to zero, then the update direction cannot be decided due to the noise, so the two sign functions become opposite and s(k) = 0.

The magnitude is computed using the smallest correction for both sensors i and j to guarantee the value of update steps are the same.

STATE-OF-THE-ART

m(t k ) = min { ↵ ⇤ ij (t k )↵ i (t k-1 ) ↵(t k-1 ) , ↵ ⇤ ij (t k )↵ i (t k-1 ) ↵(t k-1 ) , 1 ↵ ⇤ ij (t k ) ↵i (t k-1 ) ↵(t k-1 ) , 1 ↵ ⇤ ij (t k ) ↵i (t k-1 ) ↵(t k-1 ) } .
Finally, the drift compensation is reformulated as

↵i (t k ) = ↵(t k-1 ) + (1-⇢ ↵ ) [s(t k )m(t k ) ↵(t k-1 )] (3.17)
With the term s(k)m(k), RoATS converges ĉi (t) considering the transmission delay. RoATS inherit the benefit of the ATS, which is a decentralized and one hop recursion algorithm. However, RoATS uses additional communication to perform the synchronization, which is not possible when the communication graph is direct. With the direct edge, when a sensor receives a message, it might not able to send back the calculations. In RoATS, sensors need to perform back-and-forth communication to achieve their clocks synchronized.

Average time synchronization (ATSP)

In ATSP [START_REF] Wu | Average time synchronization in wireless sensor networks by pairwise messages[END_REF], the two neighboring nodes exchange their clock values and calculate the average value. In the k-th update round, three packets are transmitted between sensor i and sensor j. The three packet protocol is as follows:

1. For the first packet, sensor i sends a message at its local physical clock c s1 i (k 1) to sensor j, and sensor j receives the message at its local physical clock c r1 j (k 1). Then c r1 j (k 1) = c s1 i (k 1) + d + , where d is the transmission delay and is the clock offset.

2. For the second packet, sensor j sends the second packet back to sensor i at c s2 j (k 1) as soon as the first packet is received. The second packet contains the information of both c r1 j (k 1) , and c s2 j (k 1). When sensor i receiving the second packet at c r1 i (k 1), it estimates the average delay as

d = 1 2 c r1 j (k 1) c s1 i (k 1) + 1 2 c r2 i (k 1) c s2 j (k 1) (3.18)
and updates its clock to (3.19) 3. In order to obtain the time average synchronization, the third packet is required to send the information back from sensor i to sensor j. Similar to the second packet, when sensor i receives the second packet, it sends the third packet at

c i (k) = 1 2 c r2 i (k 1) + 1 2 c s2 j (k 1) + d ,
c s3 i (k 1) immediately.
The third packet contains the information of c r2 i (k 1), c s3 i (k 1), and the updated clock c i (k). When receiving the packet, sensor j updates its clock as follows:

c j (k) = 1 2 c r3 j (k 1) + 1 2 c s3 i (k 1) + d . (3.20)
ATSP uses pairwise messages exchanges to reduce the clock offset. ATSP apply the third packet transmission to update the clocks of both the paired sensor. The disadvantage is that the communication method is relatively complicated. In the real environment, communication between sensors in the wireless network might have many communication faults, which reduce the convergence rate significantly.

Conclusion

This chapter states the problem formulation of clock synchronization and reviews three state-of-the-art methods widely used. We also highlight their main advantages and drawbacks. Clock synchronization is still an open problem due to its complexity, 3.4. CONCLUSION limitation of the energy of sensors, and high wireless communication fault rate compared to the static wired networks. Research on energy efficient and fault tolerant protocol that can be applied to dynamic sensor network are the objectives of many research works. Next chapter presents our pulse synchronization protocol, which intends to overcome the limitation of the state-of-the-art algorithms.

Chapter 4

Time Wheel Algorithm with Drift Compensation

In this chapter, we introduce the pulse-coupled algorithm. First, the synchronization function for different events is detailed. The synchronization function is executed by a local sensor, which predicts the pulse of next round. The synchronization function saves the timing information of pulses in the time wheel and use it to calculate the next pulse. The synchronization function calculates the offset and the time period of pulses, and then synchronizes the pulse of the next round with the neighboring sensor nodes. To synchronize the pulses, both offset and drift compensation are applied. Finally, the convergence rate of the algorithm is analyzed.

Introduction

In the previous chapter, the events e b , e m , e t ⇢ E are introduced, which occur during the clock synchronization. A sensor uses synchronization functions to handle the event. When executing synchronization functions, a sensor can broadcast messages to its neighbors and to setup a timer. The timer is used to estimate the timing for the next pulse to broadcast. The timer is set by the timer set function to correct the

Synchronization functions

We start by describing the general round-structure that our algorithm operates in. At the beginning of each round, i.e., when the round number k i changes, sensor i broadcasts an electromagnetic pulse to its neighbors, where the pulse is a short burst of electromagnetic energy. By recording the local time of pulses in the present round, sensors are able to converge the timing and switch to the next round. Sensors start in round k i = 0 and broadcast a pulse in round k i 1. A sensor i waits for a period T i (k), which is subject to its local clock drift. Because the durations T i (k) and the round start times may not be the same for different sensors, pulses are not perfectly synchronized with each other. Therefore, corrections must be computed using information of previous pulses from neighboring sensors.

In previous work [Függer et al., 2015b], round numbers are assumed to be known for all sensors, which is not always realistic, because rounds estimated by sensors can be faster or slower. In our work, sensor i broadcasts a value p i 2 Z, 0  p  p max 1. The value p max is usually small so p can be represented by different frequency of 4.2. SYNCHRONIZATION FUNCTIONS the pulse. The value p i serves a similar role as k i , which increases every round, except when p i reaches p max , it returns to zero for the next round. When p max = 1, it is not possible for sensors to distinguish the different rounds. With higher p max , sensors are able to converge the p i . However, if p max is higher, it takes more time to converge due to the algorithm converge process. In our experiment, we discover that p max = 2 is sufficient for the sensors to synchronize efficiently. If p max is large, then sensors need to transmit more bits, which leads to transmission delay. If p max = 0 or p max = 1, then it is difficult for sensors to distinguish the round numbers. In a view of implementation, broadcasting p max can be represented using different frequency or wave pattern of the pulses.

A sensor i uses a circular buffer called time wheel w i to record p j from neighboring sensor j. Figure 4.1 shows the structure of the time wheel w i . A time wheel of sensor i has p max slots. The pointer in the center of time wheel is directed to the value of p j just broadcasted by sensor j. When sensor i receives a pulse with a tag p j from sensor j, it appends the local clock time c i to the list of p j -th slot of w i . The timing information in w i will be used to estimate the pulse for the next round.

Each sensor runs the pseudo-code given in Algorithm 4.1. The code is divided into three parts: initial function (line 1), receive function (line 5), and state transition function (line 9). These three parts are not executed in sequence. Instead, they are triggered by different events. The initial function is executed once when the sensors start to synchronize the clock. The receive function is executed whenever a pulse with a pointer value p j is received. The state transition function is executed when a new round starts.

Initial function

At the beginning of the algorithm, sensor i initializes an empty time wheel w i with p max slots (line 2). The sensor reset p i (line 3) to zero. The variable c i is the counter of the physical clock started at initialization, which increases subject to the clock drift. The sensor also broadcasts p i to its neighbors (line 4). w i = TimeWheel(s max )

3:

p i = 0 4:
sending_function(p i ) 5: else if e m then 6:

Obtain p j from m 7:

Append c i to the list in p j -th slot 8:

⌧ (k + 1) = timer_set_function(w i , ⌧ (k)) and set the count down timer 9: else if e t then 10:

Clear the p i -th slot 11:

p i = (p i + 1) mod (p max ) 12:
sending_function(p i ) 13: end if Receive function When sensor i receives a message from sensor j, with a tag p j (line 6). Sensor i appends the current local time c i of the received message to the list of the p j -th slot (line 7). Recall that the time between broadcasting and receiving a pulse is assumed to be negligible, since sensors only exchange zero or few bit pulses. A sensor also schedules its next state transition by setting its count down timer to ⌧ i (k + 1) c i , which is the time period for the next pulse to broadcast. (line 8). Details of the timer set function can be found in section 4.2.

State transition function

When the timer reaches zero, the sensor executes its state transition function and start the next round. First, the sensor clears the current slot from the current slot p i (line 10), because the information in the slot is used. Second, it updates p i (line 11). In the end, it broadcasts p i to its neighboring sensors (line 12).

In previous work from Függer et al. [2015b], the rounds are assumed known by sensors and the pulses in the same round is assumed to have small differences, in this work, the time wheel is applied to improve the robustness of the algorithm, so the assumption is no longer needed. In the second round, R is equal to 5 for all the sensors. However, the constant clock ticks R do not represent the constant time period due to the manufacturing defects and clock drifts. The difference of clock period causes divergence of the timing of pulses. In the third round, sensors adjust its R 0 = R + corr 1 i + corr 2 i , so that the clock skew becomes smaller. Notice that we setup k max = 3, so the round number can be 3.

The timer set function in Algorithm 4.1 (line 8) estimates the counter value for next pulse ⌧ (k + 1) according to the following formula:

⌧ i (k + 1) = ⌧ i (k) + R + corr 1 i (k) + corr 2 i (k), (4.1)
where R is a constant representing the time period between the next and present rounds. All sensors use the same value of R, but the time period between pulses T i (k + 1) may be different due to frequency variation of local clocks described in Equation 3.2. The corrections corr 1 i (k) and corr 2 i (k) are the adjustment of the timing to reduce the clock skew. These corrections adjust the pulse period using the information saved in the time wheel w. Figure 4.2 illustrates the pulse synchronization method. In the first round, sensors broadcast pulses at a similar time. In the second round, R is equal to 5 for all the sensors. However, the constant clock ticks R do not represent the constant time period due to the manufacturing defects and clock drifts. The difference of clock period causes divergence of timing of pulses. At the third round, sensors adjust its R 0 = R + corr 1 i (k) + corr 2 i (k), so that the clock skew becomes smaller. Notice that both corr 1 i (k) and corr 2 i (k) depend on the previous round k.

Each sensor only communicates through pulses and additional round information c i . The overall benefit of the method is that sensors do not require to identify their neighboring sensor and their pulses, so the message is short and robust.

Offset compensation

The term corr 1 i is used for correcting the offset of the pulses for all the sensors in previous work [Függer et al., 2015b]. corr 1 i is a time difference between the average pulse of the neighboring sensors and the pulse of sensor i itself in round k:

corr 1 i (k) = 1 |In i (k + 1)| X j2In i (k+1) (⌧ j (k) ⌧ i (k)) (4.2)
where |In i (k + 1)| is the number of incoming messages at the k + 1-th round and

P j2In i (k+1) (⌧ j (k) ⌧ i (k)
) is calculated by averaging the values in the current slot of wheel w i . Compared to the original next pulse ⌧ i (k + 1) = ⌧ i (k) + R, the equation considering corr 1 i becomes

⌧ i (k + 1) = ⌧ i (k) + corr 1 i + R = 1 |In i (k + 1)| X j2In i (k+1) ⌧ j (k) + R.
The equation shows that instead of using its previous pulse, sensor i uses the average of the previous pulses of its neighbors to estimate the timing of the next pulse. Thus, the precision of the next pulse tends to increase. However, the constant clock ticks R do not represent the constant time period due to the manufacturing defects and 4.3. CONVERGENCE ANALYSIS clock drifts. In order to further reduce the clock skew of the next pulse, R should also adjust to adapt the variation of local physical clocks.

Drift compensation

The goal of drift compensation is to let sensors agree on a global clock cycle. Since sensor i converges the previous pulses

⌧ i (k) with (⌧ i (k) + corr 1 i (k))
, where corr 1 i (k) is the correction for the previous pulse ⌧ i (k). We obtain an approximation of time period (in terms of the number of clock ticks) (R + corr 2 i (k)) in the current round as follows:

R + corr 2 i (k) = ⌧ i (k) + corr 1 i (k) 1 k , (4.3) or corr 2 i (k) = ⌧ i (k) + corr 1 i 1 k R. (4.4)
In our experiment, if drift compensation is enabled at the beginning, the consensus may not converge under some extreme cases, so the round filter is applied to Equation 4.4:

corr 2 i (k) = 8 < : 0 k  (⌧ i (k) + corr 1 i (k)) 1 k R else . (4.5)
By using Equation 4.5, the drift compensation only enabled after the `-th round.

The reason to add this filter will be derived and discussed mathematically in the next section.

Convergence Analysis

In this section, we analyze first the effect of OC algorithm (the correction with only offset correction), then analyze the effect of OCDC algorithm (the correction for both offset Correction and drift Correction).

Convergence Analysis for OC algorithm

In OC algorithm, only corr 1 i is enabled. Recall that the relation between time and the counter of the local clock of sensors is as follows:

t i (k) = ⌧ i (k) X c=1 µ i (c), (4.6)
where t i (k) represents the time of k-th pulse broadcasted by sensor i. ⌧ i (k) is the local physical clock counter of sensor i when it broadcasts k-th pulse. µ i (c) is the time period of a clock cycle, which is bounded by

(1 %)µ  µ i (c)  (1 + %)µ. (4.7)
The incremental physical clock counter in each sensor is assumed to be discrete. Each value of the physical clock counter represents a period of time.

Definition 1. Denote by ⌧ i (k) the local tick number of sensor i at which sensor i starts its k th round and by t i (k) the corresponding real time. Denote by ⌧ i j (k) the local tick number of sensor i at which it receives the pulse of sensor j for round k.

A sensor i reads the clock timing of the pulse from sensor j within the bound:

⌧ i j (k) 1 X c=1 µ i (c)  t j (k)  ⌧ i j (k) X c=1 µ i (c), (4.8)
with the assumption that the transmission delay between sensors is zero. The variable t j (k) is the real world time of k-th pulse broadcasted by sensor j, and ⌧ i j (k) is the physical local clock counter of sensor i when it receive the k-th pulse from sensor j. The variable t j (k) and ⌧ i j (k) indicate that Equation 4.8 only holds if sensor i receives the message from sensor j in round k.

CONVERGENCE ANALYSIS

Lemma 1. The relationship of the local clock difference ⌧ i j (k) ⌧ i (k) and the time difference

t j (k) t i (k) is ⌧ i j (k) ⌧ i (k) = (t j (k) t i (k)) / i j (k) + ↵ i j (k), (4.9)
where

i j (k) 2 [µ(1 %), µ(1 + %)] is caused by clock skew and ↵ i j (k) 2 [ 1, 1
] is caused by imprecision of the discrete time.

Proof. According to Equation 3.1 and Equation 4.8, the timing of the pulse is assumed to be within the bound:

⌧ i j (k) 1 X c=1 µ i (c)  t j (k)  ⌧ i j (k) X c=1 µ i (c), (4.10)
and

⌧ i (k) 1 X c=1 µ i (c)  t i (k)  ⌧ i (k) X c=1 µ i (c). (4.11)
Combining Equation 4.11 and Equation 4.10, we have

⌧ i j (k) 1 X c=1 µ i (c) ⌧ i (k) X c=1 µ i (c)  t j (k) t i (k)  ⌧ i j (k) X c=1 µ i (c) ⌧ i (k) 1 X c=1 µ i (c), (4.12) If ⌧ i j (k) 1 ⌧ i (k)
, the left-hand side of Equation 4.12 is reformulated using Equation 4.7:

⌧ i j (k) 1 X c=1 µ i (c) ⌧ i (k) X c=1 µ i (c) = ⌧ i j (k) 1 X c=⌧ i (k)+1 µ i (c) (4.13) ⌧ i j (k) ⌧ i (k) 1 (1 %) . (4.14) If ⌧ i j (k) 1  ⌧ i (k)
, with the similar derivation as Equation 4.14, the left-hand side of Equation 4.12 is

⌧ i j (k) 1 X c=1 µ i (c) ⌧ i (k) X c=1 µ i (c) = ⌧ i (k) X c=⌧ i j (k) µ i (c) (4.15) ⌧ i (k) ⌧ i j (k) + 1 (1 %) µ (4.16) = ⌧ i j (k) ⌧ i (k) 1 (1 %) µ, (4.17)
With Equation 4.14 and Equation 4.17, we have

⌧ i j (k) 1 X c=1 µ i (c) ⌧ i (k) X c=1 µ i (c) ⌧ i j (k) ⌧ i (k) 1 (1 %) µ, (4.18)
With the similar derivation as Equation 4.18, the right-hand-side of Equation 4.12 is

⌧ i j (k) X c=1 µ i (c) ⌧ i (k) 1 X c=1 µ i (c)  ⌧ i j (k) ⌧ i (k) + 1 (1 + %) µ. (4.19)
Combining Equation 4.18, Equation 4.19 with Equation 4.12, we have

⌧ i j (k) ⌧ i (k) 1 (1 %) µ  t j (k) t i (k)  ⌧ i j (k) ⌧ i (k) + 1 (1 + %) µ (4.20)
From Equation 4.20, we get

t j (k) t i (k) = ⌧ i j (k) ⌧ i (k) ↵ i j (k) i j (k) (4.21) with i j (k) 2 [µ(1 %), µ(1 + %)] (4.22) 4.3. CONVERGENCE ANALYSIS and ↵ i j (k) 2 [ 1, 1] (4.23)
to build Equation 4.9.

Lemma 2. The relationship of the local clock difference ⌧ i j (k + 1) ⌧ i (k) and the time difference t j (k + 1) t j (k) is

⌧ i (k + 1) ⌧ i (k) = (t i (k + 1) t i (k)) / i (k) + ↵ i (k) (4.24) where i (k) 2 [µ(1 %), µ(1 + %)] and ↵ i (k) 2 [ 1, 1] .
Proof. The timing of the pulse from sensor j to i is received as a discrete time event.

According to Equation 3.1 and Equation 4.8, we get:

⌧ i (k+1) 1 X c=1 µ i (c)  t i (k + 1)  ⌧ i (k+1) X c=1 µ i (c). (4.25)
Combining with Equation 4.7, we have

⌧ i (k+1) 1 X c=1 µ i (c) ⌧ i (k) X c=1 µ i (c)  t i (k + 1) t i (k)  ⌧ i (k+1) X c=1 µ i (c) ⌧ i (k) 1 X c=1 µ i (c), (4.26) If ⌧ i (k + 1) ⌧ i (k) + 1, the left-hand side of Equation 4.26 is ⌧ i (k+1) 1 X c=1 µ i (c) ⌧ i (k) X c=1 µ i (c) = ⌧ i (k+1) 1 X c=⌧ i (k)+1 µ i (c) (4.27) (⌧ i (k + 1) ⌧ i (k) 1) (1 ⇢) µ. (4.28) If ⌧ i (k + 1) 1  ⌧ i (k)
, the left-hand side of Equation 4.26 is

⌧ i (k+1) 1 X c=1 µ i (c) ⌧ i (k) X c=1 µ i (c) = ⌧ i (k) X c=⌧ i (k+1) µ i (c) (4.29) (⌧ i (k) ⌧ i (k + 1) + 1) (1 ⇢)µ (4.30) (⌧ i (k + 1) ⌧ i (k) 1) (1 ⇢)µ, (4.31)
With Equation 4.28 and Equation 4.31 we have

⌧ i (k+1) 1 X c=1 µ i (c) ⌧ i (k) X c=1 µ i (c) (⌧ i (k + 1) ⌧ i (k) 1) (1 ⇢)µ, (4.32)
With the similar derivation, the right-hand-side of Equation 4.26 is

⌧ i (k+1) X c=1 µ i (c) ⌧ i (k) 1 X c=1 µ i (c)  (⌧ i (k + 1) ⌧ i (k) + 1) (1 + ⇢)µ. (4.33)
Combining Equation 4.32, Equation 4.33 with Equation 4.26, we have 4.34, we get

(⌧ i (k + 1) ⌧ i (k) 1) (1 %) µ  t j (k) t i (k)  (⌧ i (k + 1) ⌧ i (k) + 1) (1 + %) µ (4.34) From Equation
t i (k + 1) t i (k) = (⌧ i (k + 1) ⌧ i (k) ↵ i (k)) i (k) (4.35) where i (k) 2 [µ(1 %), µ(1 + %)] (4.36)
and

↵ i (k) 2 [ 1, 1] (4.37)
to build Equation 4.24.

CONVERGENCE ANALYSIS

For the duration of round k + 1, the algorithm sets its timer to the following number of local clock ticks:

⌧ i (k + 1) ⌧ i (k) = 1 |In i (k)| X j2In i (k) ⌧ i j (k) ⌧ i (k) + R, (4.38)
where R is the original number of clock ticks between pulses.

To calculate the maximum pulse skew differences of sensors in round k, the Dobrushin semi-norm is introduced as follows defined on real vectors x by

(x) = max i,j |x i x j |. (4.39) By setting x = t(k), t(k) = [t 1 (k), ..., t n (k)]
T , we calculate the maximum pulse skew in the round k. Using properties of the Dobrushin seminorm, we are able to calculate the upper bound of the maximum pulse skew. The following conditions hold for all vector x and y 2 R n , and any scalar c.

• 1. (x) 0,

• 2. (cx) = |c| (x), and

• 3. (x + y)  (x) + (y)
To estimate the averaging operation of sensors, the network G t of WSN can be interpreted as A(k):

A(k) i,j = 8 < : 1 |In i (k)| (i, j) 2 V 0 else . (4.40)
Note that with (i, j) 2 V, the matrix A(k) is stochastic because the sum of row P n j=1 A(k) i,j is equal to one for 1  j  n. The Dobrushin semi-norm can be apply to a real matrix A as follows.

(A) = sup x2R n , (x)6 =0 (Ax) (x) . (4.41)
The properties of the stochastic matrix and the Dobrushin semi-norm described as follows:

Lemma 3. (Charron-Bost [2013]). Let A, B 2 R n⇥n . If A is stochastic, we have (AB)  (A) (B) .
When A is stochastic, we have (A)  1. We call stochastic matrices as a scrambling matrices if the strict inequality holds.

Lemma 4. [START_REF] Hajnal | Weak ergodicity in non-homogeneous markov chains[END_REF]) Let A be a stochastic matrix with and minimal positive entry . Then,

(A) = 1-min i 1 ,i 2 X j min {A i 1 ,j , A i 2 ,j } . (4.42)
Since the network is a non-split graph, the matrix A has the scrambling property.

With scrambling property and Lemma Equation 4.42, we have (A)  1 , where is the minimal positive entry of all matrices A(k) during the execution.

Since 0  A(k) i,j  1 according to Equation 4.40, it is a scrambling matrix, then k Y `=1 A(k) ! = k Y `=1 (A(k))  (1 ) k , (4.43)
In the worse case scenario, sensors are disconnect from each other. Performing averaging method will not help since

A = I, (A) = 1 and ⇣ Q k `=1 A(k) ⌘ = 1.
In a better scenario, a connected network (connected graph), each sensor has neighbors, i.e., |In i (k)| > 1, so < 1 and the value can be converge through the averaging process.

. Lemma 5. Equation 4.38 can be transformed to real time in matrix form

t(k + 1) = A(k)t(k) + b(k) + c(k), (4.44)
where A(k) is a weighted adjacency matrix, which satisfies

A i,j (k) = 8 < : 1/|In i (k)| if j 2 In i (k) 0 else . (4.45) and b i (k) = i (k) i j (k) 1 ! ⇣ A(k)t(k) i t i (k) ⌘ (4.46)
and

c i (k) = ⇣ ↵ i (k) + ↵ i j (k) + R ⌘ (k). (4.47)
with the existence of some

i j (k) 2 [µ(1 %), µ(1 + %)] and ↵ i j (k) 2 [ 1, 1]. Proof. Substitute ⌧ i j (k) ⌧ i (k) in Definition 1 using Lemma 1, we get ⌧ i (k + 1) ⌧ i (k) = 1 |In i (k)| X j2In i (k) ⌧ i j (k) ⌧ i (k) + R = 1 |In i (k)| X j2In i (k) (t j (k) t i (k)) / i j (k) + ↵ i j (k) + R. = 1 |In i (k)| X j2In i (k) (t j (k) t i (k)) / i j (k) + 1 |In i (k)| X j2In i (k) ↵ i j (k) + R.
The lower bound can be calculated by taking the largest i j (k) and smallest ↵ i j (k).

⌧ i (k + 1) ⌧ i (k) 1 |In i (k)| X j2In i (k) (t j (k) t i (k)) /µ (1 + %) 1 |In i (k)| X j2In i (k) 1 + R = 1 |In i (k)| 0 @ X j2In i (k) (t j (k) t i (k)) 1 A |In i (k)| µ (1 + %) 1 |In i (k)| |In i (k)| + R = 1 |In i (k)| 1 µ (1 + %) 0 @ X j2In i (k) (t j (k) t i (k)) 1 A 1 + R. (4.48)
With the similar derivation, the upper bound can be calculated as well.

⌧ i (k + 1) ⌧ i (k)  1 |In i (k)| 1 µ (1 %) 0 @ X j2In i (k) (t j (k) t i (k)) 1 A + 1 + R (4.49)
Combining Equation 4.48 and Equation 4.49, we get the existence of some

i j (k) 2 [µ(1 %), µ(1 + %)] and ↵ i j (k) 2 [ 1, 1] such that ⌧ i (k + 1) ⌧ i (k) = 1 |In i (k)| i j (k) X j2In i (k) (t j (k) t i (k)) + ↵ i j (k) + R (4.50)
Combining Lemma 2 and Equation 4.50 , we get

(t i (k + 1) t i (k)) = i (k) |In i (k)| i j (k) X j2In i (k) (t j (k) t i (k)) + ⇣ ↵ i (k) + ↵ i j (k) + R ⌘ i (k) 4.3. CONVERGENCE ANALYSIS Putting i (k) = i (k) i j (k)
1 and adding t i (k) on both sides of the equation, we get

t i (k + 1) = 1 |In i (k)| X j2In i (k) t j (k) + i (k) |In i (k)| X j2In i (k) t j (k) t i (k) + ⇣ ↵ i (k) + ↵ i j (k) + R ⌘ i (k) = A(k)t(k) i + i (k) ⇣ A(k)t(k) i t i (k) ⌘ + ⇣ ↵ i (k) + ↵ i j (k) + R ⌘ • i (k) (4.51)
or in vector form as Equation 4.44.

Lemma 6. Semi-norm of b(k) is b(k) = i (k) i j (k) 1 ! ⇣ A(k)t(k) t(k) ⌘ !  4% 1 % (k), (4.52)
where (k) is a simplified symbol represents the clock skew (t(k)) in round k.

Proof. According to Equation 4.39, suppose we have t i 1 and t i 2 such that (k) = (t i 1 (k) t i 2 (k)), and since A(k) is a stochastic matrix, i.e., P j [A(k)] i,j = 1, we have

t i 2 (k)  A(k)t(k)  t i 1 (k), which implies A(k)t(k) t(k) i  (k). (4.53)
The bound of

✓ i (k) i j (k) 1 ◆ is 2% 1 + %  i (k) i j (k) 1 !  2% 1 % , (4.54)
Considering Equation 4.53 and Equation 4.54 with Equation 4.46, we get

2% 1 % (k)  b i (k)  2% 1 % (k) , (4.55)
which means that

(b(k))  max i 1 b i 1 (k) min i 2 b i 2 (k)  4% 1 % (k). (4.56)
Lemma 7. By assuming original number of clock ticks between pulses R > 2, the maximum difference of the vector c(k) is

(c(k))  (4 + 2%R) µ (4.57) Proof. If R > 2, we have c i (k) 0 in Equation 4.47, so c i (k) is within a bound: ( 2 + R)µ(1 %)  c i (k)  (2 + R)µ(1 + %), (4.58)
which means that

(c(k))  max i 1 ,i 2 |c i 1 (k) c i 2 (k)| (4.59)  (2 + R)µ(1 + %) ( 2 + R)µ(1 %) (4.60) = (4 + 2%R) µ (4.61)
Theorem 1. The clock skew (k) is within a bound:

(k)  ✓ 1 + 4% 1 % ◆ k (0) + (4 + 2%R) µ 1 (1 ) k (4.62)
Proof. Transforming Equation 4.44 to precision domain using Lemma 6 and Lemma 7, we get

(k + 1)  ✓ (A(k)) + 4% 1 % ◆ (k) + (c(k)) , (4.63) 50 4.3. CONVERGENCE ANALYSIS
Expanding the right-hand side of the equation gives:

(k)  k 1 Y m=0 ✓ (A(k)) + 4% 1 % ◆ (0) + k 1 X `=0 k 1 Y m=` (A(m)) (c(k)) ! = ✓ 1 + 4% 1 % ◆ k (0) + (4 + 2%R) µ k 1 X `=0
(1

) k  ✓ 1 + 4% 1 % ◆ k (0) + (4 + 2%R) µ 1 (1 ) k , (4.64)
As a sanity check, by assuming > 4% 1 % , we can let k ! 1 and find

lim k!1 (k)  (4 + 2%R) µ . (4.65)

Convergence Analysis for OCDC Algorithm

In previous section the converge property is analyzed for corr 1 i . In this section the converge property of the OCDC algorithm with both corr 1 i and corr 2 i is also analyzed. We start by a slightly different relation between ⌧ i (k + 1) and ⌧ i (k) due to the additional correction corr 2 i .

Lemma 8. With corr 2 (k), Equation 4.38 becomes

⌧ i (k + 1) ⌧ i (k) = 1 |In i (k)| X j2In i (k) ⌧ i j (k) ⌧ i (k) k + 1 k . (4.66)
Lemma 9. Lemma 8 can be reformulated with matrix form

t(k + 1) = A(k)t(k) + b(k) + c(k), (4.67)
where

A i,j (k) = 8 < : 1/|In i (k)| if j 2 In i (k) 0 else . (4.68) and b i (k) = (k + 1) i (k) k i j (k) 1 ! ⇣ A(k)t(k) i t i (k) ⌘ (4.69)
and

c i (k) = ✓ ↵ i (k) + k + 1 k ↵ i j (k) ◆ (k) (4.70)
with the existence of some

i j (k) 2 [µ(1 %), µ(1 + %)] and ↵ i j (k) 2 [ 1, 1]. Proof. Substitute ⌧ i j (k) ⌧ i (k)
in Definition 8 using Lemma 1, we get

⌧ i (k + 1) ⌧ i (k) = 1 |In i (k)| X j2In i (k) ⌧ i j (k) ⌧ i (k) k + 1 k = 1 |In i (k)| X j2In i (k) (t j (k) t i (k)) / i j (k) + ↵ i j (k) k + 1 k
For the average, we get the existence of some

i j (k) 2 [µ(1 %), µ(1 + %)] and ↵ i j (k) 2 [ 1, 1] such that ⌧ i (k + 1) ⌧ i (k) = 1 |In i (k)| i j X j2In i (k) (t j (k) t i (k)) k + 1 k + k + 1 k ↵ i j (k) (4.71)
Combining Lemma 2 and Equation 4.71 , we get

(t i (k + 1) t i (k)) = i (k) |In i (k)| i j (k) X j2In i (k) (t j (k) t i (k)) k + 1 k + ✓ ↵ i (k) + k + 1 k ↵ i j (k) ◆ i (k) 4.3. CONVERGENCE ANALYSIS Putting i (k) = (k+1) i (k) k i j (k)
1 and adding t i (k) on both sides of the equation, we get

t i (k + 1) = 1 |In i (k)| X j2In i (k) t j (k) + i (k) |In i (k)| X j2In i (k) t j (k) t i (k) + ✓ ↵ i (k) + k + 1 k ↵ i j (k) ◆ i (k) = A(k)t(k) i + i (k) ⇣ A(k)t(k) i t i (k) ⌘ + ✓ ↵ i (k) + k + 1 k ↵ i j (k) ◆ i (k) (4.72)
or in vector form as Equation 4.67.

Lemma 10. The variable b(k) has the following constraint:

b(k)  4% 1 % k + 1 k (k). (4.73)
Proof. With (k) is the maximum difference of pulses in the same round, we have

A(k)t(k) i t i (k)  (k). (4.74)
The bound of i = (k+1) i (k)

k i j (k) 1 is 2% 1 + % k + 1 k  i (k)  2% 1 % k + 1 k , (4.75)
Considering Equation 4.74 and Equation 4.75 with Equation 4.46, we get

2% 1 % k + 1 k (k)  b i (k)  2% 1 % k + 1 k (k) , (4.76) which means that (b(k))  max i 1 b i 1 (k) min i 2 b i 2 (k)  4% 1 % k + 1 k (k). (4.77)
Lemma 11. By assuming R > 2, the maximum difference of the vector c(k) is 4.70, so c i (k) is within a bound:

(c(k))  ✓ 4k + 2 k ◆ µ (4.78) Proof. If R > 2, we get c i (k) 0 in Equation
( 1 k + 1 k )µ(1 %)  c i (k)  (1 + k + 1 k )µ(1 + %), (4.79)
which means that

(c(k))  max i 1 ,i 2 |c i 1 (k) c i 2 (k)| (4.80)  (1 + k + 1 k )µ(1 + %) ( 1 k + 1 k )µ(1 %) (4.81) = ✓ 4k + 2 k ◆ µ (4.82)
Theorem 2. The clock skew (k) is within a bound:

(k)  ✓ 1 + 4% 1 % k + 1 k ◆ k (0) + (4 + 2%R) µ 1 (1 ) k (4.83)
Proof. Transforming Equation 4.67 to precision domain using Lemma 10 and Lemma 11, we get

(k + 1)  ✓ (A(k)) + 4% 1 % k + 1 k ◆ (k) + (c(k)) , (4.84)
Expanding the right-hand side of the equation gives:

4.4. CONCLUSION (k)  k 1 Y m=0 ✓ (A(k)) + 4% 1 % k + 1 k ◆ (0) + k 1 X `=0 k 1 Y m=` (A(m)) (c(k)) ! = ✓ 1 + 4% 1 % k + 1 k ◆ k (0) + 4k + 2 k µ k 1 X `=0
(1

) k  ✓ 1 + 4% 1 % k + 1 k ◆ k (0) + 4k + 2 k µ 1 (1 ) k , (4.85)
As a sanity check, by assuming > 4% 1 % k+1 k , we can let k ! 1 and find

lim k!1 (k)  4 µ , (4.86)
which is much smaller than Equation 4.65.

Conclusion

In this chapter, a new pulse-coupled clock synchronization method is proposed (OCDC algorithm), which improves the pulse-coupled synchronization (OC algorithm) by introducing a time-wheel algorithm and drift compensation. The time-wheel algorithm enables sensors to exchange the information of the rounds during clock synchronization. Drift compensation applies an additional correction to clock skew. We also 1 To obtain the Equation 4.86, the first term in Equation 4.85 is removed because

lim k!1 ⇣ 1 + 4% 1 % k+1 k ⌘ k = 0, where the term 1 + 4% 1 % k+1 k
< 0 since we assume that

> 4% 1 % k+1 k . It is a reasonable assumption because % is much smaller than 1, so the value of 4% 1 % k+1 k
is small. Moreover, the term µ 1 (1 ) k in Equation 4.85 becomes µ when k approach to infinity because

lim k!1 (1 ) k = 0. The term 4k+2 k in Equation 4.85 becomes 4 since lim k!1 4k k = 4 and lim k!1
derive the convergence rate of both algorithms to show the improvements of the algorithm. In order the view the properties of our proposed algorithm, in the next chapter, simulations are performed to compare our proposed algorithm with the previous work Függer et al. [2015a].

Chapter 5

Simulation Results

Introduction

In this chapter, we perform several simulations to compare the proposed approach to its previous work [Függer et al., 2015a]. The simulation environment is written in Python3 and run on the MacOS operating system with 2 GHz Intel dual Core i7 processor and 8 gigabyte RAM. A different scenario is build to test the limit of the algorithm. At first, a simple environment is setup and simulated. We setup different environment to test the robustness of the algorithm. In the end, a dynamic real-world vehicular network is simulated to compare the algorithms.

Dense Scene

Figure 5.1 shows the layout of the simple environment, n = 10 sensors are randomly scattered. The average connectivity of scene kVk kN k is 5.2 per node. Each node has different sensing range to detect the pulses from neighboring sensors. Sensors can only sense the neighbors within the distance r i 2 R, which is bounded within a range:

r(1 ✏)  r i  r(1 + ✏).
The average of sensing range for all the sensors r = 30 and ✏ = 0.5. The size of the environment is 50 ⇥ 50 m 2 . We set

T i (k) = R(1 + i + k ) Figure 5.1: Scene in the experiment with | i |  30% and | k |  1% and R = 0.
03 second for different sensors. i is a sensor varying variable and k is a sensor and time varying variable. The variation i is the different frequency for different sensor, which only be randomly chosen once when sensor is initialized. The variation k is the different frequency in each round for each sensor, which affects the frequency of the clock every round. A realistic value of k is 10 2 % to 10 4 % according to Mills [1998]. We setup all the sensors waking up at time 0 and starting to broadcast at time T i (0) for every sensor.

We compare three different algorithms in the same environment:

1. an algorithm without any correction 2. clock synchronization using only corr 1 i (OC algorithm)

CLOCK DRIFT ANALYSIS

3. clock synchronization using corr

1 i + corr 2 i (OCDC algorithm)
The second algorithm is the previous work from Függer et al. [2015a], and the third is our proposed technique. We take corr 2 i into account at the sixth round. One can observe that if no clock synchronization technique is applied, the pulses for each round is diverging (the first figure in Figure 5.3). Pulses in the same round for algorithm corr 1 i are successfully trapped within an interval, but the interval will not reduce as time goes on (the second figure in Figure 5.3). In contrast, pulses from the OCDC algorithm are broadcasted more orderly than OC algorithm as time goes on (the third figure in Figure 5.3).

Figure 5.4 shows the size of clock skew for each round. We calculate clock skew by considering the nearest pulse from all sensors in the same round k to the target sensor. Clock skew for first five rounds are the same, because OCDC algorithm does not enable corr 2 yet. As soon as corr 2 is enable, the clock skew reduces dramatically. The clock skew of OCDC algorithm is approximately 7 times smaller than the clock skew of OC algorithm after they are stable.

Clock Drift Analysis

In this experiment, we analyze the methods: OC and OCDC algorithm in situations that i and k is different. We set i and k at different values to verify the robustness i . The last figure shows the pulses with much smaller clock skew thanks to both corr 1 i and corr 2 i .

Figure 5.3: Timing of pulses at each round. This figure shows the timing of pulses for each round. Dots linked by lines with different colors indicate that the pulses is broadcasted by a same sensor. One can observe that if no clock synchronization technique is applied, the pulses for each round is diverging (column 1). Pulses in the same round for OC algorithm are successfully trapped within an interval, but the interval will not reduce as time goes on (column 2). In contrast, pulses from the OCDC algorithm are broadcasted more orderly than OC as time goes on (column 3). of OC and OCDC algorithm. In order to test the edge cases, the total variation of i + k is ±20% of the constant R. Figure 5.5 on the facing page shows the five different situations. The environment is the same except i and k are changed. In these five simulations, the error for OC algorithm becomes small when i decrease, which shows that it can handle random uncertainty of clock drift k better than OCDC algorithm. In contrast, the error for OCDC algorithm increases when k increases. It shows that large amount of additional random clock drift has chance to let OCDC algorithm become worse than OC algorithm. Notice that we test the edge cases here, in normal case, k is smaller than 10 4 %.

Sudden Change of Clock Frequency

Clock drift can suddenly change due to a sensor turning into power-saving mode. There we setup an extreme experiment where all sensors change their frequency in the environment with i = 30%, and k = 1% during the 13th and 14th round in the simulation. Figure 5.6 shows the result of this experiment. Both of the methods suffer from the instant change of clock frequency at the 14-th round but converge afterward.

CLOCK DRIFT

ANALYSIS case1 case2 case3 case4 case5 i,max  20%  15%  10%  5%  0% k,max  0%  5%  10%  15%  20% Figure 5
.5: This figure shows the five different situations. The environment is exactly the same except i and k is changed. In these five simulations, the error for method corr 1 i becomes small when the i decreases, which shows that it can handle random uncertainty of clock drift k better than corr 1 i + corr 2 i . In contrast, the error for method corr 1 i + corr 2 i increase when k increases. It shows that large amount of additional random clock drift has chance to let corr 1 i + corr 2 i become worse than corr There is one line diverge at the fifth second of simulation, which means that there is a sensor that does not follow the pace of the network. The figure on the right shows that with OCDC algorithm, all sensors remain stable even if the connectivity is zero because corr 2 i modifies the local clock period according to local history data. It does not need to use the pulses from neighboring sensors to compute corr 2 i , so the round remains regular even if the network is changed. If the connectivity remains zero, the pulses in the same round diverge through time if the clock skew is not zero.

Sudden Change of a communication graph

In the wireless sensor network, sensors will not remain at one location, so the neighboring sensors will not be the same in different rounds. In this experiment, we want to discover the stability of the two methods OC and OCDC algorithm considering the changes of communication graph G(k). In this environment, all the parameters are the same as in the first experiment. In order to test our algorithm under the extreme case, all the sensors are relocated suddenly, and communication graph G(k) becomes different at the fifth second in the simulation.

Figure 5.7 shows the comparison of the two methods. Figure 5.7(left) shows the round and the pulses using OC algorithm. There is one line diverge at the fifth second of simulation, which means that there is a sensor that does not follow the pace of the network. Figure 5.7 (right) shows that with OCDC algorithm, all sensors remain stable even if the connectivity is 0, because corr 2 i modify the local clock period shows the round and the pulses using OC algorithm. We discover that some sensors may have few or 0 neighboring node. In that case, corr 1 is difficult to correct the timing of pulses, because corr 1 only depend on the pulses from neighboring sensors. The right figure shows the round and the pulses using OCDC algorithm. It is obvious that the pulses are aligned.

according to local history data. It does not need to use the pulses from neighboring sensors to compute corr 2 i , so the round remains regular even if the network is changed.

Random walk communication graph

We further analysis a rational environment to approximate the movement of the animals. In the environment, sensors are initially setup as Figure 5.1 but start to random walk with a speed s 2 [0, 2] m/s. Figure 5.8 shows the comparison of the two methods. Figure 5.8 (left) shows the round and the pulses using only corr 1 . We discover that some sensors may have few or zero neighboring node. In that case, corr 1 is difficult to correct the timing of pulses, because corr 1 only depends on the pulses from neighboring sensors. Figure 5.8 (right) shows the round and the pulses using both OC and OCDC algorithm. It is obvious that the pulses are aligned. That is because the calculation of corr 2 only involve the pulses of the sensor itself, which is accessible even if the sensor has no neighbor.

EXTREME ENVIRONMENT

Figure 5.9: Experiment in harsh environment. A comparison between pulses broadcasted using corr 1 i (upper curve) and using corr 1 i + corr 2 i (lower curve). The y-axis represents the value of the maximum clock skew during the simulation.

Extreme Environment

In the simulation, to drive our algorithm to the extreme, 50 sensors are randomly spread in a 217 ⇥ 217 m 2 environment. The connection for every sensor is distance r = 15 ± 7.5 m 2 . Moreover, the sensors are dynamic and randomly relocated every 0.05 second. We set k max = 3 and i = 20% and k = 1%. Figure 5.9 shows the convergence result, where corr 1 + corr 2 not only converge faster but also have much smaller clock skew then corr 1 .

Vehicular Network

Vehicular network is a categories of wireless sensor network, where vehicles communicate to each other to perform data fusion. With vehicular network, the vehicles are able to corporate to increase the road safety or improve the traffic. We simulate

Static Vehicle Network

Figure 5.10 shows the simulation environment, where n = 75 vehicles are randomly placed with an average number of 2.84 neighbors per vehicle. The average sensing range for all the vehicles is r = 300 (m) and sensing range variation ✏ = 0.2. Sensing range for each sensor node i is bounded

(1 ✏)r  r i  (1 + ✏)r.
(5.1)

The number of slots of time wheel is s max = 2. of the vehicles. The parameters in this simulation are listed as follows:

• n = 75 vehicles

• Low connectivity c = 2.84 (nodes in average)

• sensing range r = 300 (m)

• clock drift % = 0.2

• maximum round number s max is set to 2

• sensing range variation ✏ = 0.2 Figure 5.12: Clock skew comparison between OC (red) and OCDC (blue) algorithm.

Clock skew for the first five rounds are the same because OCDC algorithm includes the term corr 2 = 0. As soon as corr 2 is enabled, the clock skew reduces significantly.

The clock skew of OCDC algorithm is approximately 7 times smaller than the clock skew of OC algorithm after both clock skew are stabilized.

We compare the two different algorithms in the same environment:

1. OC algorithm

OCDC algorithm

The first algorithm is a pulse synchronization algorithm considering a highly dynamic network published by Függer et al. [2015b] and the second is our proposed technique. We take corr 2 i into account after the sixth round. Figure 5.12 shows the comparison of these algorithms. The OC algorithm (Figure 5.12 red) successfully converges all the pulses in the same round within an interval, but the interval will not reduce as time goes on (Figure 5.12 blue). In contrast, pulses from the OCDC algorithm are broadcasted more orderly than corr 1 i as time goes on. Figure 5.12 shows the clock skew for each round. Clock skew for first five rounds are the same, because OCDC algorithm includes the term corr 2 = 0. As soon as corr 2 is enabled, the clock skew reduces significantly. The clock skew of OCDC algorithm is approximately seven times smaller than the clock skew of OC algorithm after both clock skews are stabilized.

Pulse Frequency and Clock Skew

In this section, different frequencies of the pulse are considered. Same parameters as the previous section are used. If the vehicles exchange pulses more frequently, then the clock skew temp to converge faster. Fig Figure 5.13 illustrates four simulations where the scene is setup the same as Fig Figure 5.10 except T (k), the default time period between two pulses, which is different in the four simulations. If T (k) is smaller, then the vehicle network converges faster then when T (k) is larger. However, higher pulses frequently cost higher energy consumption considering the same time period. The parameter T (k) can be chosen to balance the convergence rate and the energy consumption. In the vehicle network, all pulses in the same round converge regardless T (k). There is also a possibility to slow down the pulse exchange after the clock synchronization has converged.

Moving Vehicle Network

Vehicles are initially setup as Figure 5.10. The acceleration and the deceleration of each vehicle is 2.6 (m/s 2 ) and 4.5 (m/s 2 ), respectively, which is the default parameters in SUMO simulator. Vehicles slow down when they turn and stop when traffic lights are red. There are 214 traffic lights on the map. The realist vehicle traces are generated by SUMO.

Figure 5.14 shows the clock skew of OC and OCDC algorithm in the simulation. The clock skew is initialized at 0.07 seconds and drops efficiently due to the correction term corr 1 i . The communication graph is fully connected at the beginning, so the clock skew is smaller after few seconds. During the simulation, some mobile vehicle has no or few neighbors. In that case, the clock offset and drift cannot converge, because there is no communication between the networks, so the clock skew starts to increase. For OC algorithm, the frequency of the pulses reverts to its original frequency, causing relatively large clock skew. For OCDC algorithm, the clock skew increase slowly due to unconnected topology of the communication graph. When the nodes in the network are connected, pulses are re-synchronized, so the clock skew drop. Overall OCDC algorithm has better performance than OC algorithm due to the additional correction of corr 2 i , which is more robust in the dynamic vehicle network.

Conclusion

In this chapter, the proposed clock synchronization is simulated and compared with the previous work of Függer et al. [2015a]. From the simulation, it is clear to observe that the clock skew is reduced dramatically compared to the previous work. A random walk environment is simulated, which shows that with corr 2 i , sensors are able to remember the corrected frequency even if they do not have any neighboring Figure 5.14: Time of each round of a dynamic vehicle network for 2 seconds. The clock skew is initialized at 0.07 seconds and drop efficiently due to the correction term corr 1 i . The communication graph is fully connected at the beginning, so the clock skew is smaller after few seconds. During the simulation, some mobile vehicle has no or few neighbors. In that case, the clock offset and drift cannot converge, because there is no communication between the networks, so the clock skew starts to increase. For OC algorithm, the frequency of the pulses reverts to its original frequency, causing relatively large clock skew. For the OCDC algorithm, the clock skew increase slowly due to the unconnected topology of the communication graph. sensors. Two different types of clock drift are also analyzed to show the performance of the corrections. A vehicle network is also simulated to show the benefit of adding corr 2 i .

Chapter 6

Conclusion

This chapter summarizes the principal contributions of the first part of this thesis.

Suggestions for further work in some of the subjects encompassed by the research contained in this thesis are also included. At the end of the chapter, a list of related publications is provided.

Conclusion and Contributions

We proposed a new pulse-coupled synchronization by introducing a time-wheel algorithm and drift compensation. The time-wheel algorithm enables sensors to exchange the information of the rounds during clock synchronization. Drift compensation applies an additional correction to clock skew. The simulation results show that the proposed algorithm performs better than the compared algorithm, especially when considering a highly dynamic wireless sensor network. The proposed algorithm is successfully applied to the vehicular network with a real geometric location and dynamic vehicles.

This thesis represents a step forward in research into pulse-coupled clock synchronization. The second correction reduces the clock skew of the sensors in a natural and elegant way. Moreover, using a comparison with the state-of-the-art technique, CHAPTER 6. CONCLUSION it was proved that the proposed algorithm was more precise in most of the cases. This work is published on the 29th IEEE Intelligent Vehicle Symposium Han et al. [2018b].

Future Works

For future work, we intend to evaluate the algorithm using a practical wireless sensor network. The next step will be to select an appropriate communication hardware and to closely simulate it in a large scale network.

Chapter 7 Introduction

Wireless sensor networks (WSNs) play an important role in applications such as environmental monitoring, source tracking, etc. In WSNs, the locations of sensing nodes are often required to make the collected information useful. Node localization algorithms have thus been developed, usually relying on a number of anchor nodes, whose locations are known. Localization is mainly performed considering the received signal strength (RSS), the time of arrival (TOA), the time difference of arrival (TDOA) of electromagnetic signals transmitted by a node to the anchor nodes or by the anchor nodes to the node to localize, see [START_REF] Gezici | A survey on wireless position estimation[END_REF][START_REF] Li | Contributed review: Source-localization algorithms and applications using time of arrival and time difference of arrival measurements[END_REF][START_REF] Iliev | Review and comparison of spatial localization methods for low-power wireless sensor networks[END_REF].

In this chapter, one considers guaranteed confidence region characterization for the location estimate of nodes using RSS measurements. In most of the papers, location estimates are evaluated first, and approximate characterization of confidence regions are then provided. Here, one evaluates first accurate confidence regions with prescribed confidence levels. Estimates may then be obtained using classical techniques, or may be evaluated as the barycenter of the confidence region. For that purpose, the Leave-out Sign-dominant Correlated Regions (LSCR) [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF] approach is employed to define a region in the parameter space guaranteed to contain the true value of the parameter with an arbitrary confidence level. Confidence regions as defined by LSCR are valid even in the non-asymptotic regime, i.e., when only a small set of measurements is available.

Motivation

In most of the papers, the RSS measurement noise is assumed to be normal or lognormal [START_REF] Gholami | Wireless network positioning as a convex feasibility problem[END_REF][START_REF] Vaghefi | Cooperative received signal strength-based sensor localization with unknown transmit powers[END_REF]. This allows one to apply classical least-squares (LS), maximum likelihood (ML), or maximum a posteriori (MAP) estimation techniques. The main difficulty is that the cost function associated to the LS, ML, or MAP estimation problem may be non-convex, especially when the path-loss exponent as well as the reference transmission power have to be estimated jointly with the location. The differential RSS measurement localization approach proposed in [START_REF] Wang | Robust device-free wireless localization based on differential RSS measurements[END_REF]. With the differential RSS, one avoids estimating the reference power and facilitates estimation, e.g., via the total least-squares approach considered in [START_REF] Rahimi | A generalized total least squares algorithm for RSS-based localization with unknown path-loss model parameters[END_REF]. Convex [START_REF] Tomic | RSS-based localization in wireless sensor networks using convex relaxation: Noncooperative and cooperative schemes[END_REF] or semi-definite [START_REF] Kim | Robust localization with unknown transmission power for cognitive radio[END_REF][START_REF] Lee | Location estimation using differential RSS with spatially correlated shadowing[END_REF][START_REF] Wang | Efficient semidefinite relaxation for energy-based source localization in sensor networks[END_REF], [START_REF] Ouyang | Received signal strength-based wireless localization via semidefinite programming: Noncooperative and cooperative schemes[END_REF][START_REF] Vaghefi | Cooperative received signal strength-based sensor localization with unknown transmit powers[END_REF] relaxation techniques have been proposed, to facilitate the optimization. The non-linear localization problem is relaxed to a convex problem and solved by the semi-definite programming. Using a semi-definite relaxation, [START_REF] Lohrasbipeydeh | A minimax SDP method for energy based source localization with unknown transmit power[END_REF] presents a minimax approach involving again the differential RSS measurements. Minmax approach iteratively minimize the possible loss for a worst case senario. In most of these approaches, the estimator confidence region is then defined using the Cramér-Rao lower bound (CRLB) [START_REF] Torrieri | Statistical theory of passive location systems[END_REF]. Nevertheless, confidence regions obtained from CRLB are usually optimistic and provide a relatively inaccurate confidence region when there are few measurements. Alternative bounded-error estimation techniques have been proposed in [START_REF] Mourad | Anchor-based localization via interval analysis for mobile ad-hoc sensor networks[END_REF], [START_REF] Léger | Guaranteed robust distributed estimation in a network of sensors[END_REF], Mourad et al. [2011a] to get guaranteed set estimates. In these latter approaches, the RSS measurement errors are assumed to be bounded with known bounds. One tries to evaluate the set of all source locations that are consistent with the measurement model and the noise bounds. The resulting set is guaranteed to contain the true location of the source node, provided that the hy-7.2. LEAVE-OUT SIGN-DOMINANT CORRELATED REGIONS potheses on the noise bounds are satisfied. Nevertheless, in practice, noise samples do not really follow a normal or lognormal distribution, and the bounds considered in bounded-error estimation are either often violated because they are too small and result in an empty set estimate, or lead to huge solution sets when the bounds are too large.

Leave-out Sign-dominant Correlated Regions

LSCR belongs to a family of recently-proposed techniques to define non-asymptotic confidence regions (NACR) [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF], 2006, 2010[START_REF] Csáji | Non-asymptotic confidence regions for the least-squares estimate[END_REF]. The main assumptions of LSCR are (i) that the system output corresponds to the model output for a given true value of the parameter vector (ii) that the measurement noise samples are independent and symmetrically distributed. Considering a given confidence level, and a finite set of measurements, LSCR defines a subset of the parameter space that is guaranteed to contain the true value of the parameter vector with a probability equal to the confidence level. Considering relatively mild conditions on the measurement noise, one is thus able to obtain guaranteed confidence regions. Nevertheless, the confidence region as defined by LSCR may be non-convex or even may consist of several disconnected components. As shown in [START_REF] Kieffer | Guaranteed characterization of exact non-asymptotic confidence regions as defined by LSCR and SPS[END_REF], using tools from interval analysis [START_REF] Jaulin | Applied interval analysis: with examples in parameter and state estimation, robust control and robotics[END_REF], one is able to obtain inner and outer approximations of the confidence regions as defined by LSCR using subpavings, i.e., set of non-overlapping interval vectors.

Outline

In Chapter 8, tools from Interval analysis is described, which is used by some of the localization methods. With the knowledge of Interval analysis, chapter 9 introduces RSS localization model and specify several localization approaches. In chapter 10, the methodology used in this work is shown and compared with other localization approaches. The different correlation functions of LSCR are also simulated and analyzed. Finally, the conclusions of the survey is given in Chapter 12.

Interval

The interval [x] is a connected subset of R. Even when the interval is not closed, we shall keep the notation [x]. Whether the empty set ; should also be considered as an interval is still a subject of discussion. The lower bound lb ([x]) of an interval [x], also denoted by x, is defined as

x = lb[x] , sup{a 2 [x]|8x 2 [x], a  x} (8.1)
Its upper bound ub ([x]), also denoted by x, is defined as

x = ub[x] , inf{b 2 [x]|8x 2 [x], x  b} (8.2)
Thus, x is the largest number on the left of [x] and x is the smallest number on its right. The width of an interval

[x] is w([x]) , x x (8.3)
and the centre of any bounded and non-empty interval [x] is defined as

mid([x]) , x x 2 . (8.4)
Also, the set-theoretic operations can be applied to intervals. The intersection of two intervals [x] and [y] is defined by

[x] \ [y] , {z 2 R|z 2 [x] and z 2 [y]}, (8.5)
and the union of two intervals [x] and [y] is defined by

[x] [ [y] , {z 2 R|z 2 [x] or z 2 [y]}. (8.6)
There are four classical operations of real arithmetic, namely addition (+), subtraction ( ), multiplication (⇤) and division (/). All of them can be extended to 8.2. INTERVAL intervals as follows:

[x] + [y] = [x + y, x + ȳ] [x] [y] = [x ȳ, x y] [x] ⇤ [y] = [min(x ⇤ y, x ⇤ ȳ, x ⇤ y, x ⇤ ȳ)] [x]/[y] = [x] ⇤ (1/[y]).
(8.7)

An interval vector [x] is a subset of R n that can be defined as the Cartesian product of n closed intervals. An interval vector [x] is also called a box, which is written as

[x] = [x 1 ]⇥[x 2 ]⇥...⇥[x n ], with [x i ] = [x, x] for i = 1, ..., n. Its ith interval components [x i ] is the projection of [x]
onto the ith axis. Many of the operations extend without difficulty to boxes.

x = lb(x) , (lb(x 1 ), lb(x 2 ), ..., lb(x 3 )) T x = ub(x) , (ub(x 1 ), ub(x 2 ), ..., ub(x 3 ))

T [x] + [y] = ([x 1 ] + [y 1 ]) ⇥ ([x 2 ] + [y 2 ]) ⇥ ... ⇥ ([x n ] + [y n ]) [x] [y] = ([x 1 ] [y 1 ]) ⇥ ([x 2 ] [y 2 ]) ⇥ ... ⇥ ([x n ] [y n ]) [x] T ⇤ [y] = [x 1 ] ⇤ [y 1 ] + [x 2 ] ⇤ [y 2 ] + ... + [x n ] ⇤ [y n ] mid([x]) , (mid(x 1 ), mid(x 2 ), ..., mid(x n )) T w([x]) , max 1in w([x i ]) (8.8)
Let R m⇥n be the set of all matrix with real coefficients, m rows and n columns. A m ⇥ n dimensional interval matrix is a subset of R m⇥n that can be defined as the Cartesian product of m ⇥ n closed interval matrix. The interval matrix [A] is written in the following form:

[A] = 0 B B @ [a 11 ] . . . [a 1n ] . . . . . . [a m1 ] . . . [a mn ] 1 C C A , (8.9)
Many of the operations can be extended.

A = lb([A]) , 0 B B @ [a 11 ] . . . [a 1n ] . . . . . . [a m1 ] . . . [a mn ] 1 C C A Ā = ub([A]) , 0 B B @ [ā 11 ] . . . [ā 1n ] . . . . . . [ā m1 ] . . . [ā mn ] 1 C C A [A] + [B] = 0 B B @ [a 11 + b 11 ] . . . [a 1n + b 1n ] . . . . . . [a m1 + b m1 ] . . . [a mn + b nm ] 1 C C A [A] [B] = 0 B B @ [a 11 b 11 ] . . . [a 1n b 1n ] . . . . . . [a m1 b m1 ] . . . [a mn b nm ] 1 C C A [A] ⇤ [x] = n X j=1 [a ij ] ⇤ [x j ] ! 1im mid([A]) , (mid([a ij ])) 1im,1jn w([A]) , max 1im,1jn
w([a ij ]) (8.10)

Inclusion Function

Natural Inclusion Function

Natural inclusion function is similar to normal functions, except that variables and functions are interval. One can easily convert a normal mathematical equation to Natural inclusion function by replacing the variables to intervals. Natural inclusion function are not minimum in general cases, because of the dependency effects. The accuracy of the natural inclusion function depend on the expression of f .

Example Consider the function f = (x + 1) 2 , g = x 2 + 2x + 1. It is obvious that f = g. However, if one applies natural inclusion function to evaluate [f ] and [g] with [x] = [ 1, 1] then we get:

[f ]([x]) = ([ 1, 1] + 1) 2 = [0, 4] [g]([x]) = [ 1, 1] 2 + 2 ⇤ [ 1, 1] + 1 = [0, 1] + [ 2, 2] + 1 = [ 2, 4] (8.12)
The results are different due to different expression of f and g. For same function but different expression, if x occurs more times, the result tend to be more inaccuracy.

Centred Inclusion Function

Centred inclusion function can be applied to any function. Centred inclusion function linearize the function respect to the input parameter [x]. Given any function f , we can obtain its centred inclusion function as follow:

[f ]([x]) = f (m) + [J f ]([x]) ⇤ ([x] m), (8.13)
where m can be any point in [x] and [J f ] is the Jacobian matrix of [f ].

Mix-Centred Inclusion Function

The centred inclusion function can be improved by a more complicated formulation.

The main idea for mix-centred inclusion function is to apply Equation 8.13 for n times, where n is the number of the variables. The case n = 3 will be treated first, to simplify the exposition. Next, we show an example for the mix-centred inclusion function for any

[f ]([x 1 ], [x 2 ], [x 3 ]).
To derive the mix-centred form, the function

[f ]([x 1 ], [x 2 ], [x 3 ]) is modified to 1 variable form, for example f [x 1 ],m 3 ([x 2 ]
) , where the [x 2 ] is the only variable for this function and the subscripts are the fixed variables.

The interval [x 1 ] is fixed and the value of [x 3 ] is fixed at an point m 3 . According to this formation, we enumerate the following 3 valid functions using Equation 8.13:

[f m 2 ,m 3 ]([x 1 ]) = f m 2 ,m 3 (m 1 ) + [g 1 m 2 ,m 3 ]([x 1 ])([x 1 ] m 1 ) (8.14) [f [x 1 ],m 3 ]([x 2 ]) = [f [x 1 ],m 3 ](m 2 ) + [g 2 [x 1 ],m 3 ]([x 2 ])([x 2 ] m 2 ) (8.15) [f [x 1 ],[x 2 ] ]([x 3 ]) = [f [x 1 ],[x 2 ] ](m 3 ) + [g 3 [x 1 ],[x 2 ] ]([x 3 ])([x 3 ] m 3 ) (8.16) Since [f [x 1 ],[x 2 ] ](m 3 ) in Equation 8.16 is equal to [f [x 1 ],m 3 ]([x 2 ]
) in Equation 8.15, it can be substituted using Equation 8.15. Moreover,

[f [x 1 ],m 3 ](m 2 ) in Equation 8.15 is equal to [f m 2 ,m 3 ]([x 1 ]
) in Equation 8.14, so it can be substituted to Equation 8.14.

After the substitutions, the equation become:

[f [x 1 ],[x 2 ] ]([x 3 ]) = f m 2 ,m 3 (m 1 ) +[g 1 m 2 ,m 3 ]([x 1 ])([x 1 ] m 1 ) +[g 2 [x 1 ],m 3 ]([x 2 ])([x 2 ] m 2 ) +[g 3 [x 1 ],[x 2 ] ]([x 3 ])([x 3 ] m 3 ) (8.17)
It is the same equation as:

[f ]([x 1 ], [x 2 ], [x 3 ]) = f (m 1 , m 2 , m 3 ) + [g 1 ]([x 1 ], m 2 , m 3 ) ⇤ ([x 1 ] m 1 ) + [g 2 ]([x 1 ], [x 2 ], m 3 ) ⇤ ([x 2 ] m 2 ) + [g 3 ]([x 1 ], [x 2 ], [x 3 ]) ⇤ ([x 3 ] m 3 ). (8.18)
The expression can be generalize as follow:

f ([x]) ⇢ f (m) + n X i=1 [g i ]([x 1 ], ..., [x i ], m i+1 , ..., m n ) ⇤ ([x i ] m i ) (8.19)

Monotonic Inclusion Function

Monotonic inclusion function is a type of inclusion function that reduce the pessimism efficiently. Monotonic inclusion function performs well when the function is monotonic, because it reduce the interval pessimism by checking whether

f ([x]) is monotonic (lb( @[f ] @[x i ] ) > 0 or ub( @[f ] @[x i ] ) < 0) for each interval [x i ]. If f ([x]
) is monotonic for variable x i , we can further find out the x + i , the desired value of x i which maximize f ([x]). We can also find out the x i the desired value of

x i which minimize f ([x]). If f ([x]
) is monotonic increase, then we let x + i = xi and

x i = x i . If f ([x]
) is monotonic decrease, then we let x + i = x i and x i = xi . Let f be a function defined on variables V of domains [V ]. Let X ✓ V be a subset of monotonic variables. Consider the values x + i and x i such that: if x i 2X is an increasing variable, then x + i = xi and x i = x i , and if x i 2 X is an decreasing variable, then x + i = x i and x i = xi . Consider W = V \X the set of variables not detected monotonic. Then, f min and f max are functions defined by:

f min (W ) = f (x 1 , x 2 , ..., x n , W ) f max (W ) = f (x + 1 , x + 2 , ..., x + n , W ) (8.20)
Finally, the monotonicity-based extension [f ] of f in the box [V ] produces the following interval vector:

[f ] = [[f min ](W ), [f max ](W )] (8.21) Example Consider for example f (x 1 , x 2 , w) = x 2 1 + x 1 x 2 + x 2 w 3w in the box [V ] = [6, 8] ⇥ [2, 4] ⇥ [7, 15]. 8.4. SIVIA ALGORITHM [f ]([x 1 ], [x 2 ], [w]) = [6, 8] 2 + [6, 8] ⇥ [2, 4] + [2, 4] ⇥ [7, 15] 3 ⇥ [7, 15] = [-83, 35] @f @x 1 (x 1 , x 2 ) = 2x 1 + x 2 (8.22) [ @f @x 1 ]([6, 8], [2, 4]) = [ 14, 8]. (8.23)
Since [ 14, 8] < 0, we deduce that f is decreasing w.r.t. x 1 . Similarly, we can get

@f @x 2 (w) = 1 + w (8.24) [ @f @x 2 ]([w]) = [8, 16] (8.25)
Since [8, 16] > 0, we deduce that x 2 is increasing w.r.t x 2 . Finally, 0 2

[ @f @w ]([x 1 ], [x 2 ], [w]) = [ 1, 1]
, so that w is not deduced monotonic. The monotonicity-based evaluation yields: (6, 4, [7, 15])] = [-79, 27]

[f ]([V ]) = [[f ](x 1 , x 2 , [w]), [f ](x 1 , x 2 , [w])] = [[f ](8, 2, [7, 15]), [f ]

SIVIA algorithm

A set-inversion problem can be solved by SIVIA algorithm proposed by Jaulin and Walter [1993b]. Let f be a possibly non-linear function from R n to R m and let Y be a subset of R m . Set inversion is the characterization of

X = {x 2 R n |f (x) 2 Y} = f 1 (Y). (8.26) Algorithm 8.1 SIVIA algorithm Require: [f ], Y, [x] Ensure: X, X 1: if [f ]([x]) \ Y = ; return; 2: [f ]([x]) ⇢ Y then X := X [ [x]; X := X [ [x]; return; 3: w([x] < ") then X := X [ [x]; return; 4: SIVIA([f ], Y, L[x], ", X, X) 5: SIVIA([f ], Y, R[x], ", X, X)
For any Y ⇢ R m and for any function f admitting a convergent inclusion function

[f ], two regular subpavings X and X such that X ⇢ X ⇢ X. can be obtained with the algorithm SIVIA (Set Inverter Via Interval Analysis).

SIVIA require a search box [x] to which X is guaranteed to belong. The general cases may be encountered:

1. If [f ]([x]
) has a non-empty intersection with Y, but is not entirely in Y, then [x] may contain a part of the solution set. [x] is said to be undetermined. If [x] has a width greater than a pre-specified percision parameter ", then it should be bisected and the test should be recursively applied to these newly generated boxes.

If [f ]([x]

) has an empty intersection with Y, then [x] does not belong to X and can be cut off from the solution tree.

If [f ]([x]

) is entirely in Y, then [x] belongs to the solution subpaving X, and is stored in X and X.

4. If the box considered is undetermined, but its width is lower than ", then it is deemed small enough to be considered in the outer approximation X of X.

Algorithm 8.1 shows the SIVIA based on an inclusion function, where X and X have been initialized as empty.

8.5. CONTRACTORS

Contractors

Forward-backward contractor is a simple technique to contract interval using constraints. It can be used in wireless sensor localization problem. However, this contractor only consider one equation at a time and solve equation independently. On the other hand, Newton contractor and Krawczyk contractor are powerful since they can consider multiple equations at the same time. In this thesis, we aim to compare the forward-backward contractor between Newton contractor and Krawczyk contractor.

Constraint Satisfaction problem (CSP)

Consider n x variables x i 2 R, i 2 {1, ..., n x } and n f constraint f j (x 1 , x 2 , ..., x nx ) = 0, j 2 {1, ..., n f } (8.27)

We assume that the prior domain for

x is [x] = [x 1 ] ⇥ [x 2 ] ⇥ ... ⇥ [x nx ].
Let f be the function whose coordinate functions are the f j s. Equation 8.27 can be formulated as a constraint satisfaction problem

H : (f (x) = 0, x 2 [x]) (8.28)
The solution set of H is defined as

S = {x 2 [x]|f (x) = 0} (8.29)
Contracting S means replacing [x] by a smaller domain which still contain the solution set. The optimal contraction is to replace [x] by the smallest box that contains S.

Forward-Backward Contractor

The forward-backward contractor is based on constraint propagation. It decomposes the function into forward constraints and backward constraints to reduce the size of the potential search space. It tries to find a consistent solution for each equation in ICSP. The forward-backward contractor only considers the equations in ICSP one by one, so it can be applied on the ICSP where the number of functions and the number of variables are different. Consider a simple constraint x 3 = x 1 x 2 . This constraint can be rewritten as forward phase and backward phase: Forward phase: 8,40] is the intervals we want to contract. The contractor can be build as Forward phase: In Figure 8.2, two anchors a and b are located at (4, 3) and ( 4, 3), respectively. Both of the anchors can measure the distance d + e to agent and send this information to the agent, where e is the measurement noise with zero mean and variance 2 . There is one agent located at x = (0, 0), whose location need to be determined. The observation location of the agent for anchors a and b are (u a1 , u a2 ) and (u b1 , u b2 ), re-8.5. CONTRACTORS spectively. We can model scenario of Figure 8.2 to a constraint satisfaction problem:

x 3 = x 1 x 2 Backward phase: x 1 = x 3 /x 2 x 2 = x 3 /x 1 If the box [x] = [x 1 ] ⇥ [x 2 ] ⇥ [x 3 ] = [1, 4] ⇥ [1, 4] ⇥ [
[x 3 ] = [x 3 ] \ ([x 1 ][x 2 ]) Backward phase: [x 2 ] = [x 2 ] \ ([x 3 ]/[x 1 ]) [x 1 ] = [x 1 ] \ ([x 3 ]/[x 2 ]) (
H : 8 > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > : (u a1 4) 2 + (u a2 3) 2 (d 2 a + N (0, 2 )) = 0 (u b1 + 4) 2 + (u b2 3) 2 (d 2 b + N (0, 2 )) = 0 u a1 x 1 = 0 u b1 x 1 = 0 u a2 x 2 = 0 u b2 x 2 = 0 d a 5 = 0 d b 5 = 0 (8.31)
Since N (0, 2 ) is uncertain, it is difficult to calculate Equation 8.31. Instead, we can model the Equation 8.31 to a interval constraint satisfaction problem (ICSP) given

3 = 0.2 H : 8 > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > : ([u a1 ] 4) 2 + ([u a2 ] 3) 2 [d 2 a ] = 0 ([u b1 ] + 4) 2 + ([u b2 ] 3) 2 [d 2 b ] = 0 [u a1 ] [x 1 ] = 0 [u b1 ] [x 1 ] = 0 [u a2 ] [x 2 ] = 0 [u b2 ] [x 2 ] = 0 [d a ] = [4.8, 5.2], [d b ] = [4.8, 5.2] (8.32)
We assume the agent's location is already known within the box

[x] = [[ 2, 2], [ 2, 2]]
T , but still have high uncertainty. By applying the forward-backward contractor, we can reduce the x down to [ 1.1, 1.1][ 1.3, 2.0]. It is apparent that the result is not optimized. Forward-backward contractor fails to reduce [x] to the minimum size of the solution, so it is not an optimal contractor for this ICSP. Forward-backward can-Figure 8.2: A scenario of localization CSP problem. The two anchors a and b are located at (4, 3) and ( 4, 3), respectively. Both of the anchors can measure the distance d + e to agent and send this information to the agent, where e is the measurement noise with zero mean and variance 2 . There is one agent located at x = (0, 0), whose location need to be determined. not consider multiple equations at the same time. Krawczyk and Newton contractor are numerical contractors that consider ICSP as a matrix. They use the Jacobian matrix to find the solution of ICSP. Since the Jacobian matrix can evaluate the slop using multiple constraints instead of one. However, Krawczyk and Newton contractor cannot contract successfully when the slop of the constraints are dramatically variant.

Newton Iteration Function

Here we introduce the newton iteration function because it is a preliminary used for building Krawczyk and Newton contractors. Figure 8.3 shows the idea of Newton iteration: one starts with an initial guess (x n ) which is reasonably close to the true root, then the function (f ) is approximated by its tangent line ( f ), and one computes the x-intercept of this tangent line (a point xn+1 such that f (x n+1 ) = 0 ). This Figure 8.3: Newton Iteration in one dimensional. One starts with an initial guess (x n ) which is reasonably close to the true root, then the function (f ) is approximated by its tangent line ( f ), and one computes the x-intercept of this tangent line (a point xn+1 such that f (x n+1 ) = 0 ). This x-intercept is a better approximation to the function's root than the original guess, and the method can be iterated.

Figure 8.4: A scenario that newton iteration fails to converge the solution to root, where x 1 , x 2 , and x 3 are the results of each iteration. We can observe that from the second to the third iteration, x 3 is farther away from the root than x 2 , so the solution is not converged.

x-intercept is a better approximation to the function's root than the original guess, and the method can be iterated. If the initial guess is far away from the true root and the function is not monotonic, this iterative procedure may not converge to the solution point due to the bad approximation of the tangent lines. Figure 8.4 shows a scenario that newton iteration fails to converge the solution to root, where x 1 , x 2 , and x 3 are the results of each iteration. We can observe that from the second to the third iteration, x 3 is farther away from the root than x 2 , so the solution is not converged.

To introduce the newton iteration in a more formal way, we first assume that f : R k ! R k is differential, and we have the first guess of the solution xn . Since f (x n ) 6 = 0, we want to find a better solution xn+1 . First, we approximated f by its tangent line at xn

f (x) ⇡ f (x) = f (x n ) + J f (x n )(x xn ) (8.33)
Where J f is the Jacobian matrix of f . Then, we solve f (x n+1 ) = 0 to get the x-intercept of this tangent line. Equation 8.33 is reformulated as

xn+1 = xn J 1 f (x n )f (x n ) (8.34)
where xn+1 is the better approximation of the xn . We compute Equation 8.34 iteratively until an sufficient accurate value is obtained.

In Equation 8.34, the one has multiply with inverse of the k-by-k Jacobian matrix J 1 f (x). This is the Newton iteration method solving a system of k equations. Assume that the system has m variables and k equations with m > k. The Newton's method can still be used if the algorithms use the generalized inverse of the non-square Jacobian matrix J + = ((J T J) 1 )J T (from Wikipedia) instead of the inverse of J. If the non-linear system has no solution, the method attempts to find a solution in the non-linear least squares sense. For subsection 8.5.4 on the next page and subsection 8.5.5 on page 100, we will demonstrate Krawczyk contractor and Newton contractor using a system with k variables and k equations (invertible matrix) for simplification.

Krawczyk Contractor

Theoretically, we can compute the interval version of newton iteration by the following equation, which is the inclusion function of Equation 8.34 where xn and xn+1 are is replace by the interval [x] and [x 0 ], respectively:

[x 0 ] = [x] [J 1 f ]([x])[f ]([x]) (8.35)
where [x 0 ] is a fixed-point solver for CSP problem.

[x] is the interval which contain the solution set S, [f ] is the inclusion function of f and [J f ] is the inclusion function for the Jacobian matrix J of function f . However, Equation 8.35 is suffered from too many occurrence of [x]. To reduce the occurrence, Krawczyk contractor replace 8.35 to reduce the occurrence of x:

J 1 f ([x]) by J 1 f (x 0 ) in Equation
[ ]([x]) = [x] J 1 f (x 0 )[f ]([x]) (8.36)
where x 0 is a point near x. This equation is still a fixed-point solver because

[f ]([x]) = 0 imply [ ]([x]) = [x] J 1 f (x 0 ) ⇤ 0 = [x].
Then, it linearize Equation 8.36 using centred inclusion function. Combine Equation 8.36 and Equation 8.13 to get

[ ]([x]) = (x 0 ) + [J ]([x]) ⇤ ([x] x 0 ) (8.37)
where [J ] is an inclusion function for the Jacobian matrix of . We can derive [J ] using Equation 8.36. Notice that the J f (x 0 ) is only a scalar matrix, the differential operation of Equation 8.36 become easier than Equation 8.34 and the occurrence of

[x] is reduced. J (x) = Inclusion_Jacobian(x J 1 f (x 0 )f (x)) (8.38) = I J 1 f (x 0 )J f (x) (8.39)
From the centred inclusion fix-point function is obtained, the Krawczyk contractor can be formulate as follow.

CHAPTER 8. INTERVAL ANALYSIS

C K : [x] 7 ! [x] \ ( (x 0 ) + [J ]([x]) ⇤ ([x] x 0 )) (8.40)
Equation 8.40 can be further modified by substitute with Equation 8.36 and by substitute [J ]([x]) with Equation 8.38.

C K : [x] 7 ! [x] \ (x 0 J 1 f (x 0 )f (x 0 ) + (I J 1 f (x 0 )[J f ]([x])) ⇤ ([x] x 0 )) (8.41)

Newton Contractor

Newton contractor is not a direct interval version of the newton iteration function.

Still, it approximates the function by its tangent lines. Given any x and x 0 2 R and a continuous function f , we can find a point ⇠ that satisfy the following equation:

f (x) = f (x 0 ) + J f (⇠)(x x 0 ). (8.42)
Since ⇠ is guaranteed within the interval [x], we turn Equation 8.42 into the interval version as follow.

[f ](

[x]) = f (x 0 ) + [J f ]([x])([x] x 0 ) (8.43)
where x 0 can be any point in the interval [x]. According to mean-value theorem, we can find an approximated slop of function within [x] such that the approximate function has the same root as the original function, so J f (⇠) is replaced by

[J f ]([x]). Since the Jacobian matrix [J f ]([x]
) is a interval instead of a point. It enumerate all the possibility of the approximated tangent lines. In order to find all the x-intercepts of those tangent lines, we solve [ f ]([x]) = 0 using Equation 8.43 to get

[x] = x 0 [J 1 f ]([x])f (x 0 ) (8.44)
Figure 8.5: A simple example of newton contractor of a k = 1 (1 variable and 1 function) system. We use colors to distinguish different lines. The curve black line represents the function f . The orange lines represent the interval [x] before (after) the contraction. The green lines are the tangent lines with maximum and minimum slop. The red lines are maximum and minimum slop shifted to the point (x 0 , f (x 0 )).

There are two gray triangle area represent all possibility of the tangent lines that pass through the point (x 0 , f (x 0 )). Since the tangent lines are the approximation of function f , we can find the better solution by finding all the x-intercepts of the tangent lines.

By applying Equation 8.30, the Newton Contractor can be derived as follows.

C N : [x] 7 ! [x] \ (x 0 [J 1 f ]([x])f (x 0 )), x 2 [x] (8.45)
Figure 8.5 is a simple example of newton contractor of a k = 1 (1 variable and 1 function) system. We use colors to distinguish different lines. The curve black line represents the function f . The orange lines represent the interval [x] before and after the contraction. The green lines are the tangent lines with maximum and minimum slop. The red lines are maximum and minimum slop shifted to the point (x 0 , f (x 0 )).

There are two gray triangle area represent all possibility of the tangent lines that pass through the point (x 0 , f (x 0 )). Since the tangent lines are the approximation of function f , we can find the better solution by intersecting x-axis and the triangle area.

Gauss-Seidel Contractor

Gauss-Seidel contractor can be applied to calculate [J 1 f ]([x])f (x 0 ) in Equation 8.45 efficiently. Gauss-Seidel contractor is a iterative interval matrix solver. Given a problem:

Ap = b, where A 2 [A], b 2 [b], p 2 [p].
Assume A is invertible. We can decompose A as the sum of a diagonal matrix and a matrix with zeros on its diagonal:

A = diag(A) + extdiag(A) (8.46) Solving Ap = b is equivalent to diag(A)p + extdiag(A)p = b (8.47)
The equation can be rewritten as

p = diag(A) 1 (b extdiag(A)p) (8.48)
8.6. CONCLUSION An inclusion function for Equation 8.48 is

[p] = diag([A]) 1 ([b] extdiag([A])[p]) (8.49)
We can contract p by applying Equation 8.49 iteratively until p is contracted.

To calculate

[J 1 f ]([x])f (x 0 ) using Gauss-Seidel contractor, one can take [A] = [J f ]([x]), [b] = f (x 0 ), and [p] = x 0 [x]. The reason for [p] = x 0 [x] is because Equation 8.43 imply that [J 1 f ]([x])f (x 0 ) = x 0 [x].

Conclusion

This chapter introduces the tools of interval analysis. Interval analysis allow us to consider the uncertainty as bounds. With the inclusion function, the bounded uncertainty can propagate to different domains. A contractor is able to reduce the solution space using multiple inclusion functions. For the next section several localization algorithms, which are based on interval analysis will be discussed.

Chapter 9

Problem Formulation and Related Approaches

Introduction

Localization in a network of sensors is usually performed considering the received signal strength (RSS), the time of arrival (TOA), the time difference of arrival (TDOA) of electromagnetic or acoustic signals transmitted by some anchor nodes, see [START_REF] Gezici | A survey on wireless position estimation[END_REF] for an overview. In this thesis, one will consider localization from RSS, which has received significant attention for more than a decade. By measuring the received power of agent node from a specify location, we are able to estimate the location of an agent.

Problem Formulation and Assumptions

Consider a set of n a anchor nodes situated on a plane, with fixed and known locations ✓ i , i = 1, . . . , n a , in some reference frame attached to the plane. This paper aims at providing a guaranteed NACR for the unknown location ✓ 0 of some static node in the x-y plane, i.e., for a given confidence level, to provide a subset of locations to which the node belongs with the considered confidence level. For that purpose, the node periodically emits some electromagnetic or acoustic signal, which is received by the anchors. The k-th RSS measurement performed by anchor node i is denoted y (i, k). The anchor nodes transmit the RSS measurements to some central processing unit in charge of estimating ✓ 0 and the associated confidence region from y (i, k), i = 1, . . . , n a , k = 1, . . . , n.

We assume that the k-th measurement obtained by anchor node i is described by the Okumura-Hata model [START_REF] Okumura | Field strength and its variability in VHF and UHF land-mobile radio service[END_REF]

y (i, k) = P 0 10 P log 10 k✓ 0 ✓ i k d 0 + " (i, k) , (9.1)
where P 0 is the signal power measured at some known reference distance d 0 , P is the path-loss exponent, k•k is the Euclidian norm, and " (i, k) is the measurement noise.

One assumes that P is the same for all anchors. Usually, the " (i, k)s are assumed to be realizations of some normal or log-normal random variables. In this paper, one only assumes that the " (i, k)s are realizations of independently, not necessarily identically distributed random variables with distributions symmetric around zero.

Since P 0 and P are not necessarily known by the anchors nodes, one considers the parameter vector p = ⇥ ✓ T 0 , P 0 , P ⇤ T for which a confidence region has to be estimated and from which a confidence region for ✓ 0 may then be deduced. Notice that d 0 = 1 is a fixed parameter linked with P 0 , thus it does not appear in p.

Assuming that there exists some true value p ⇤ = ⇥ ✓ ⇤T 0 , P ⇤ 0 , ⇤ P ⇤ T of the parameter vector, Equation 9.1 can be rewritten as

y(i, k) = y m (i, p ⇤ ) + " (i, k) (9.2) with y m (i, p ⇤ ) = P ⇤ 0 10 ⇤ P log 10 k✓ ⇤ 0 ✓ i k d 0 . (9.3)
One assumes further that the search for the confidence region for p ⇤ has to be done in some known compact set P 0 , which may be chosen arbitrarily large.

STATE-OF-THE-ART APPROACHES

State-of-the-art Approaches

Range-based approaches consider the localization problem in mobile sensor networks using interval analysis and other tools. Range-based approaches can take outliers into account, i.e., data not matching the observation model. Some range-based approaches perform a set-membership estimation, where only the maximal number of outliers is required to be known. Mourad et al. [2011b], [START_REF] Léger | Guaranteed robust distributed estimation in a network of sensors[END_REF] Using these algorithms, estimates consist of sets of boxes whose union surely contains the correct location of the sensor, provided that the considered hypotheses are satisfied.

Bounded-Error Approach

BE estimation techniques [START_REF] Halbwachs | Bounded-error estimation for mobile vehicule localization[END_REF], Jaulin and Walter [1993a] assume that the noise samples " (i, k) in Equation 9.1 are bounded with known bounds, i.e., that " (i, k) 2 [e (i, k) , e (i, k)], with w (i, k) and w (i, k) known for all i = 1 . . . n a and k = 1 . . . n. The set of all parameter vectors p belonging to some initial search set P 0 that are consistent with the measurement model, the measurements, and the noise bounds is then characterized. This set may be defined as

P BE = {p 2 P 0 |y m (i, p) 2 y (i, k) [w (i, k) , w (i, k)] , i = 1 . . . n a , k = 1 . . . n} .
(9.4) Introducing where

w (i) = ( w (i, 1) , . . . , w (i, n)) T w (i) = (w (i, 1) , . . . , w (i, n)) T
⌧ (p) = na X i=1 n X k=1 ⌧ i,k (p) , (9.8) ⌧ i,k (p) = 8 < : 1 if y m (i, p) 2 y (i, k) [w (i, k) , w (i, k)] ,
0 else (9.9)

and

Y ⇠ = [n a n ⇠, n a n].
Again, (Equation 9.4) can be viewed as a set-inversion problem

P RBE,⇠ = P 0 \ ⌧ 1 (Y ⇠ ) , (9.10)
which may be solved using SIVIA.

This approach does not need to specify in advance which measurements have to be considered as outliers. Considering a fixed value of ⇠, one obtains an estimator robust against ⇠ outliers. If ⇠ is not fixed a priori, a possible strategy is to start with ⇠ = 0, and to increase ⇠ as long as an empty solution is provided by SIVIA. When a non-empty solution set has been obtained for some ⇠ 0 > 0, one may consider ⇠ = ⇠ 0 + 1 or ⇠ = ⇠ 0 + 2 and evaluate P RBE,⇠ to further increase the robustness of the estimated set P RBE,⇠ .

The resulting set P RBE,⇠ is guaranteed to contain the true value p ⇤ of the parameter vector provided that there are actually less than ⇠ outliers.

Nonlinear ML estimate and Cramér-Rao Lower Bound

Maximum Likelihood Estimation (MLE)

In MLE [START_REF] Ziskind | Maximum likelihood localization of multiple sources by alternating projection[END_REF], Assuming that the noise samples w(i, k) in Equation 9.1 are independently and identically distributed (iid) as zero-mean log-normal 109 CHAPTER 9. PROBLEM FORMULATION AND RELATED APPROACHES with variance 2 , one may obtain the MLE of p ⇤ as pMLE = arg max ⇡(y|p), (9.11) where ⇡(y|p) is the likelihood function, which requires a precise knowledge of the noise probability distribution. The cost function to minimize is usually non-convex and has several local minimizers. Nevertheless, one may try to minimize it using, e.g., lsqnonlin of Matlab. The function applies Levenberg-Marquardt algorithm (LMA) to solve the optimization problem. First, considering the logarithm of ⇡(y|p), we get

l(p) = log(⇡(y|p)) = log na Y i=1 n Y j=1 1 p 2⇡ exp (y(i, j) y m (i, p)) 2 2 2 !! ) (9.12) = c na X i=1 n X j=1
(y(i, j) y m (i, p)) 2 , (9.13)

where c = log ⇣ 1 p 2⇡ ⌘ 2 2 .
In order to find pMLE , we need to maximum ⇡(y|p), which implies minimizing

l 0 (p) = na X i=1 n X j=1 (y(i, j) y m (i, p)) 2 . (9.14)
Given a start point of p, LMA replaces p by a new estimation p+ . To determine , the functions y m (i, p + ) are approximated by their linearization. y m (i, p + ) ⇡ y m (i, p) + @y m (i, p) @p (9.15)

STATE-OF-THE-ART APPROACHES

Combining Equation 9.14 and Equation 9.15, we get

l 0 (p) ⇡ c na X i=1 n X j=1 ✓ y(i, j) y m (i, p) @y m (i, p) @p ◆ 2 . (9.16)
Reformulate Equation 9.16 in vector notation, we get l 0 (p) ⇡ ky y m (p) J k 2 , (9.17)

= (y y m (p)) T (y y m (p)) 2(y y m (p))J + T J T j (9.18)
where J is the Jacobian matrix of y m (i, p). The elements in the matrix is calculated as

[J] ij = [y m (i, p)] i [p] j . (9.19)
The operator [•] k in Equation 9.19 represents the k-th element in the vector or matrix. Since we want to minimize the function Equation 9.17, we can take the derivative with respect to and setting the result to zero:

@l 0 (p) @ = 0 (9.20) = 2(y y m (p))J + J T J (9.21) = 2(J T J) 1 [y y m (p)]J. (9.22)
We apply Equation 9.20 for every iteration to replace p to p+ , until p has converged.

Cramér-Rao Lower Bound (CRLB)

Once an estimate pMLE has found, the CRLB [START_REF] Torrieri | Statistical theory of passive location systems[END_REF] may be used to derive an ellipsoidal confidence region in the parameter space. Assuming that pMLE is unbiased, and introducing the estimator covariance matrix ⌃, one has 111 CHAPTER 9. PROBLEM FORMULATION AND RELATED APPROACHES (9.23) where J(p) is the dim(p) ⇥ dim(p) Fisher information matrix, with components

⌃ < J 1 (p),
[J(p)] i 1 ,i 2 = E  @ 2 ln ⇡(y|p) @p i 1 @p i 2 . (9.24)
The variable y demotes the observation vector in Equation 10.16, and ⇡(y|p) is the join conditional probability density function of the observation given p. The variable p i 1 and p i 2 are the i 1 -th and i 2 -th elements in p, respectively. When applying CRLB, we need to determine the probability function of the noise. We assume the noises are independent and identically distributed Gaussian random variables. The probability density function of an observation for a RSS measurement in Equation 9.1 is

⇡ i,j (y(i, j)|p) = 1 p 2⇡ e (y(i,j) y m (i,p)) 2 2 2 , (9.25)
where i represents the anchor node who measures the RSS value of a agent node. Since all probability density functions for different RSS value measured by anchor nodes are independent to each other. We can multiply all probability density function to get the joint conditional probability density function.

⇡(y|p) = na Y i=1 n Y j=1 ⇡ i,j (y(i, j)|p) (9.26)
Since each agent is localized individually. For a specific agent. By combining Equation 9.25, Equation 9.26 and Equation 9.1, the joint conditional probability density function of the observations is

⇡(y|p) = Na Y a=1 1 p 2⇡ exp ( (y(a) P 0 + 10 p log 10 k✓ 0 ✓ i k d 0 ) 2 2 2
) .

(9.27)

After the joint conditional probability density function of the observations is build,

STATE-OF-THE-ART APPROACHES

the fisher information matrix for localization is derived. We demonstrate fisher information matrix using three unknown variables for p.

J(p = [✓ 0,1 , ✓ 0,2 , P 0 ] T ) = E 0 B B @ 2 6 6 4
@ 2 ln ⇡(y|p) @✓ 0,1 @✓ 0,1 @ 2 ln ⇡(y|p) @✓ 0,2 @✓ 0,1 @ 2 ln ⇡(y|p) @P 0 @✓ 0,1 @ 2 ln ⇡(y|p) @✓ 0,1 @✓ 0,2 @ 2 ln ⇡(y|p) @✓ 0,2 @✓ 0,2 @ 2 ln ⇡(y|p) @P 0 @✓ 0,2 @ 2 ln ⇡(y|p) @✓ 0,1 @P 0 @ 2 ln ⇡(y|p) @✓ 0,2 @P 0 @ 2 ln ⇡(y|p) @P 0 @P 0

3 7 7 5 1 C C A (9.28)
and

J(p = [✓ 0,1 , ✓ 0,2 , p ] T ) = E 0 B B @ 2 6 6 4
@ 2 ln ⇡(y|p) @✓ 0,1 @✓ 0,1 @ 2 ln ⇡(y|p) @✓ 0,2 @✓ 0,1 @ 2 ln ⇡(y|p) @ p@✓0,1 @ 2 ln ⇡(y|p) @✓ 0,1 @✓ 0,2 @ 2 ln ⇡(y|p) @✓ 2 @✓ 0,2 @ 2 ln ⇡(y|p) @ p@✓0,2 @ 2 ln ⇡(y|p) @✓ 0,1 @ p @ ln ⇡(y|p) @✓ 0,2 @ p @ 2 ln ⇡(y|p) @ p@ p 3

7 7 5 1 C C A . (9.29)
The derivative elements in

J(p = [✓ 0,1 , ✓ 0,2 , P 0 ] T ) and J(p = [✓ 0,1 , ✓ 0,2 , p ] T ) are @ 2 ln ⇡(y|p) @✓ 1 @✓ 1 = (10 p ) 2 ln 10 Na X a=1 (✓ 0,1 ✓ a,1 ) 2 k✓ 0 ✓ a k 4 , (9.30) @ 2 ln ⇡(y|p) @✓ 2 @✓ 2 = (10 p ) 2 ln 10 Na X a=1 (✓ 0,2 ✓ a,2 ) 2 k✓ 0 ✓ a k 4 , (9.31) @ 2 ln ⇡(y|p) @✓ 0,1 @✓ 0,2 = @ 2 ln ⇡(y|p) @✓ 0,2 @✓ 0,1 = (10 p ) 2 ln 10 Na X a=1 (✓ 0,1 ✓ a,1 )(✓ 0,2 ✓ a,2 ) k✓ 0 ✓ a k 4 , (9.32) @ 2 ln ⇡(y|p) @✓ 0,1 @P 0 = @ 2 ln ⇡(y|p) @P 0 @✓ 0,1 = 10 p 2 ln 10 Na X a=1 (✓ 0,1 ✓ a,1 ) k✓ 0 ✓ a k 2 , (9.33) @ 2 ln ⇡(y|p) @✓ 0,2 @P 0 = @ 2 ln ⇡(y|p) @P 0 @✓ 0,2 = 10 p 2 ln 10 Na X i=1 (✓ 0,2 ✓ a,2 ) k✓ 0 ✓ a k 2 , (9.34)
@ 2 ln ⇡(y|p) @P 0 @P 0 = Na X a=1 2 , (9.35) @ 2 ln ⇡(y|p) @ p @✓ 0,1 = @ 2 ln ⇡(y|p) @✓ 0,1 @ p = 100 p 2 ln(10)

Na X a=1 (✓ 0,1 ✓ a,1 ) ln(k✓ 0 ✓ a k) k✓ 0 ✓ a k 2 , (9.36)
@ 2 ln ⇡(y|p) @ p @✓ 0,2 = @ 2 ln ⇡(y|p) @✓ 0,2 @ p = 100 p 2 ln(10)

Na X a=1 (✓ 0,2 ✓ a,2 ) ln(k✓ 0 ✓ a k) k✓ 0 ✓ a k 2
, (9.37) @ 2 ln ⇡(y|p) @P 0 @ p = @ 2 ln p(y|p) @ p @P 0 = 0, (9.38) and @ 2 ln ⇡(y|p) @ p @ p = ( 10 ln( 10)

) 2 Na X a=1 ln(k✓ 0 ✓ a k 2 ) k✓ 0 ✓ a k 2 .
(9.39)

As long as we have J(p = [✓ 0,1 , ✓ 0,2 , P 0 ] T ) and J(p = [✓ 0,1 , ✓ 0,2 , p ] T ), we can obtain the minimum covariance matrix of ⌃ by calculating the inversion of J(p = [✓ 0,1 , ✓ 0,2 , P 0 ] T ) and J(p = [✓ 0,1 , ✓ 0,2 , p ] T ), then we can illustrate CRLB as an error ellipse. An error ellipse is often used as a visual aid or performance measure to depict the accuracy or performance of an estimator, or other stochastic system. It implies that there are two (or more) variables involved, and that they are normally (Gaussian) distributed with known covariance.

Let us start with an 2-D zero mean random variable x = "

x 1 x 2

#

. We assume the covariance matrix of

x is C x = " 2 1 1 2 ⇢ 1 2 ⇢ 2 2

#

. It has the probability density function:

p(x) = 1 (2⇡) 2 p |C x | exp( 0.5x T C 1 x x), (9.40)
where C 1 x is:

C 1 x = " 2 1 ⇢ 1 2 ⇢ 1 2 2 2 # .
(9.41)

Combining Equation 9.40 and Equation 9.41, we get

9.3. STATE-OF-THE-ART APPROACHES pdf(x 1 , x 2 ) = 1 2⇡ 1 2 p 1 ⇢ exp ✓ 1 2(1 ⇢ 2 )  1 2 1 x 2 1 2⇢ 1 2 x 1 x 2 + 1 2 1 x 2 2 ◆ . (9.42)
We fix the value of the probability, i.e., p(x 1 , x 2 ) = p ⇤ , so we can obtain

pdf ⇤ = 1 2⇡ 1 2 p 1 ⇢ exp ✓ 1 2(1 ⇢ 2 )  1 2 1 x 2 1 2⇢ 1 2 x 1 x 2 + 1 2 1 x 2 2 ◆ , (9.43)
which can be rewritten as

k 2 = 1 2(1 ⇢ 2 )  1 2 1 x 2 1 2⇢ 1 2 x 1 x 2 + 1 2 1 x 2 2 , (9.44)
where

k = 2 ln(p ⇤ 2⇡ 1 2 p 1 ⇢). (9.45)
We can observe that Equation 9.44 is an ellipse with x 1 and x 2 being coordinate axes.

A more practical way to describe Equation 9.44 for n random variables is in terms of matrix. To draw the ellipse to unknown parameters of the localization problem, we set x = p p ⇤ , so that x is a vector of random variables with zero means.

k 2 = x T C x x. (9.46)
Given a fixed value pdf ⇤ and a covariance matrix C x , we can find a set of point b such that {b|k 2 = b T C x b} to illustrate the lower bound of the confidence level (n = 1 for an interval, n = 2 for an error ellipse, n = 3 for an ellipsoid).

To find the set b , [START_REF] Torrieri | Statistical theory of passive location systems[END_REF] apply a transformation matrix A to let z = Ax, where z's components are independent standard normal. i.e., no correlation between any of elements of z, and the variations of the elements in z are 1. After finding A, the desired b can easily be found using b = A 1 a, where a is a set of point of a circle with unit radius.

To find the transformation matrix A, we first apply the covariance of cov(z) = I. 8: Compute the error ellipse, b = A 1 a. 9: Add the offset by the mean of the random variable x.

Recall the properties of multivariate transformation. If x is a random vector of n variables with covariance C x , and A is an n ⇥ n matrix, then z = Ax is a random vector of n variables with covariance AC x A T . With this properties, we obtain

cov(z) = AC x A T = I. (9.47)
Then, we turn to Eigendecomposition, which allows any square matrix C x to be expressed as C x = VDV 1 , where V is an orthonormal square matrix of eigenvectors, and D is a diagonal matrix of eigenvectors. Then we have

cov(z) = AVDV 1 A T = AVD 1 2 D 1 2 V 1 A T = I, (9.48) which satisfied if D 1 2 VA T = I. Thus, setting A = D 1 2 V (9.49)
allows z = Ax to have a standard normal distribution. Working with z is much easier than working with x because z has a standard normal distribution, and the error ellipse is a circle. on the current page summarizes the procedure for finding the error ellipse.

We also want to use CRLB to calculate the lower bound of the location error (e).

STATE-OF-THE-ART APPROACHES

For a 2-D dimension, we can obtain the lower bound of the variance of ✓ 0,1 and ✓ 0,2 , which are

[J 1 ] 1,1 , [J 1 ] 2,2
, respectively. (The [M] 1,1 means the element in first row and first column; [M] 2,2 means the element in the second row and second column)

p E(e 2 ) = q E[(✓ ⇤ 0,1 ✓ 0,1 ) 2 ] + E[(✓ ⇤ 0,2 ✓ 0,2 ) 2 ] = q [J 1 ] 1,1 + [J 1 ] 2,2 . (9.50)
CRLB evaluates the lower bound of the minimum estimator.

Semidefinite programming (SDP) with Cramér-Rao Lower Bound

The approach considered in [START_REF] Vaghefi | Cooperative received signal strength-based sensor localization with unknown transmit powers[END_REF] consists in solving the localization problem via semidefinite programming. The confidence region is then evaluated considering confidence ellipsoid associated to the CRLB. This approach is briefly recalled here.

Only ✓ 0 and P 0 are assumed to be unknown. First, Equation 9.1 is reformulated as

log 10 d 2 i i = P 0 5 P + w(i, k) 5 P , (9.51)
where d i = k✓ 0 ✓ i k and i = 10 y(i,k)/5 P . Introducing ↵ = 10 P 0 /5 P , one deduces that

d 2 i i = ↵10 w(i,k)/5 P . (9.52)
Assuming that the noise is sufficiently small, the right-hand side of Equation 9.52 is expanded as

d 2 i i = ↵ + ✏(i, k), (9.53)
where ✏(i, k) = ln 10 5 P w(i, k). An estimator for p = P 0 , ✓ T 0 T is then expressed as b p = arg min

P 0 ,✓ 0 X i=1...na,k=1...n d 2 i i ↵ 2 , (9.54)
Finding b p is difficult in general, since the cost function is non-linear and non-convex. To address this issue, [START_REF] Vaghefi | Cooperative received signal strength-based sensor localization with unknown transmit powers[END_REF] introduce the auxiliary variable z = ✓ T 0 ✓ 0 to get the SDP problem

⇣ b ↵, b z, b ✓ 0 ⌘ = min ↵,z,✓ 0 X i=1...na,k=1...n d 2 i i ↵ 2 (9.55) subject to d 2 i = " ✓ i 1 # T " I 2 ✓ 0 ✓ T 0 z # " ✓ i 1 # , (9.56)
" I 2 ✓ 0 ✓ T 0 z # ⌫ 0. (9.57)
The solution of this problem may be obtained with standard SDP solvers such as SDPT3 [START_REF] Toh | SDPT3-a matlab software package for semidefinite programming, version 1.3[END_REF]. An estimate b p SDP is then easily deduced from b ↵, b z, and b ✓ 0 . The confidence region in this case is again obtained from the CRLB.

Conclusion

In this chapter, state-of-the-art algorithms are described. Some of the methods [START_REF] Vaghefi | Cooperative received signal strength-based sensor localization with unknown transmit powers[END_REF][START_REF] Torrieri | Statistical theory of passive location systems[END_REF][START_REF] Ziskind | Maximum likelihood localization of multiple sources by alternating projection[END_REF] allows one to apply maximum likelihood (ML), or maximum a posteriori (MAP) estimation techniques to the localization problem. The main difficulty is that the cost function associated to the ML, or MAP estimation problem may be non-convex, especially when the path-loss exponent, as well as the reference transmission power have to be estimated jointly with the location.

In the light of this, bounded error methods [START_REF] Mourad | Anchor-based localization via interval analysis for mobile ad-hoc sensor networks[END_REF], [START_REF] Léger | Guaranteed robust distributed estimation in a network of sensors[END_REF], Mourad et al. [2011a] are proposed, where the RSS measurement errors are assumed to be bounded with known bounds. One tries to evaluate the set of all source locations that are consistent with the measurement model and the noise bounds. The resulting set is guaranteed to contain the true location of the source node, provided that the hypotheses on the noise bounds are satisfied. Nevertheless, in practice, noise samples do not really follow a normal or lognormal distribution, and the bounds considered in bounded-error estimation are either often violated because they are too small and result in an empty set estimate, or lead to huge solution sets when the bounds are too large. In the next chapter, we introduce our proposed technique, which only uses mild assumptions on the noise, so the localization accuracy and the robustness are improved.

• Bounded-error approach (BE)

• Robust Bounded-error approach (RBE)

The goal of our experiments is to show the advantages of LSCR for localizing agents using the measurements from anchors.

Confidence Regions as Defined by LSCR

First consider the vectors

y (i) = (y (i, 1) , . . . , y (i, n)) T (10.1) y m (i, p) = (y m (i, p) , . . . , y m (i, p)) T (10.2)
gathering the RSS measurements and the model output for the i-th anchor node. Then, consider the vectors When p = p ⇤ , w (p ⇤ ) contains the n a ⇥ n realizations of the measurement noise for all anchors at all measurement time instants.

To define a NACR, LSCR introduces three integers s > 0, m > 0, and q > 0, and a vector e = (e 0 , ..., e s ) T of non-negative integers such that at least one e j , 0  j  s is odd. Considering some p 2 P 0 , for t = 1, ..., ny s, one can evaluate o of these subsets forms a group under symmetric difference, i.e., 8I j 2 G ny s and 8I j 0 2 G ny s , one has (I j [ I j 0 ) (I j \ I j 0 ) 2 G ny s . A procedure to build such G ny s is described in section 10.3. In Equation 10.7, I ny s j is the cardinal number of I ny s j .

c e t,e (p)
Then, the set

P e,q,m = ( p 2 P 0 | m X i=1 ⌧ i,e (p) > q and m X i=1 ⌧ + i,e (p) > q ) , (10.8)
where

⌧ i,e (p) = 8 > > > < > > > : 0 if s i,e (p) < 0 B 0.5 if s i,e (p) = 0 1 if s i,e (p)>0
(10.9) and

⌧ + i,e (p) = 8 > > > < > > > : 0 if s i,e (p) > 0 0 if s i,e (p) = 0 1 if s i,e (
p) < 0 (10.10) has been shown in [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF] to be such that .11) In ( 12) and ( 13), B 0.5 is a Bernoulli random variable with Pr (B 0.5 = 0) = Pr (B 0.5 = 1) = 0.5 to solve ties.

Pr (p ⇤ 2 P e,q,m ) = 1 2q m . ( 10

CHAPTER 10. PROPOSED LSCR APPROACH

The set P e,q,m contains all values of p 2 P 0 such that at least q out of the m s i,e (p) are positive and q are negative. It is a NACR containing p ⇤ with a probability that is exactly 1 2q m .

Construction of the groups G N

Original LSCR uses G N to construct correlation functions s j,e (p) in Equation 10.7.

With the properties of G N , the probability of the confidence region inEquation 10.11 holds. The analysis of the behavior of the LSCR approach in the considered case where different measurement models have to be considered requires a slight adaptation of the construction of the groups G N used by the LSCR approach as presented in Dalai et al. [2007b].

Consider the set I N = 1, ..., N and an integer m = 2 µ .To build a group G N for the symmetric difference operator, as in Dalai et al. [2007b], one considers first Gordon's iterative construction [START_REF] Gordon | Completely separating groups in subsampling[END_REF]. One starts with R(1) = [1] and evaluate

2 6 4 R(k 1) R(k 1) 0 R(k 1) J R(k 1) 1 0 T 1 T 1 3 7 5 , k = 1, . . . , µ, (10.12)
where J is the all-one matrix of appropriate size, and 0 and 1 are respectively the all-zero and all-one vectors of appropriate size. Then one considers

Q m 1 = " R(k) 0 T # = [q 1 , . . . , q m 1 ], (10.13)
where the q i s are the column vectors is repeated dN/(m 1)e times so that Q has at least N columns. The matrix G obtained selecting the first N columns of Q is the incidence matrix of the group

G N . Q is such that Q i,j = 1 if j 2 I N i .
10.4. CHARACTERIZING CONFIDENCE REGIONS

Characterizing Confidence Regions

The characterization of confidence regions as defined by LSCR is performed in [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF] by sampling points in P 0 and determining whether they belong to P e,q,m . One thus gets an inner approximation of P e,q,m , consisting of a union of points in P 0 which have been proved to belong to P e,q,m . When the confidence region is small, consists of several disconnected subsets, and when the number of sampling points is not enough, there is a risk to be unable to find any point in P e,q,m , except with a fine gridding of P 0 .

An alternative approach based on interval analysis has been proposed in [START_REF] Kieffer | Guaranteed characterization of exact non-asymptotic confidence regions as defined by LSCR and SPS[END_REF]. Its main property is to yield inner and outer approximations of P e,q,m consisting of unions of non-overlapping boxes. This approach is briefly recalled in what follows.

SIVIA

Characterizing P e,q,m may be formulated as a set-inversion Jaulin and Walter [1993b] problem:

P e,q,m = {p 2 P 0 |f r,m (p) 2 Y q } = P 0 \ f 1 e,m (Y q ) , (10.14)
where

f e,m (p) = " f + e,m (p) f e,m (p) # = " P m i=1 ⌧ + i,e (p) P m i=1 ⌧ i,e (p) # , (10.15)
and

Y q = " Y + q Y q # = " [q, 1[ [q, 1[ # , (10.16)
which may be efficiently solved via interval analysis [START_REF] Jaulin | Applied interval analysis: with examples in parameter and state estimation, robust control and robotics[END_REF][START_REF] Moore | Interval analysis. englewood cliff[END_REF] using the SIVIA algorithm [START_REF] Jaulin | Applied interval analysis: with examples in parameter and state estimation, robust control and robotics[END_REF]. For that purpose, inclusion functions for the f e,m are required. SIVIA in Algorithm 10.1 recursively partitions P 0 into 3 different kinds of boxes: boxes proved to belong to P e,q,m , boxes proved to have no intersection with P e,q,m , and undetermined boxes for which no conclusion can be CHAPTER 10. PROPOSED LSCR APPROACH Algorithm 10.1 LSCR with SIVIA Require: [f e,m ]([p]), [p], q, " Ensure: P e,q,m , P e,q,m 1:

if [f e,m ]([p]) \ Y q = ; return; 2: if [f e,m ]([p]
) ⇢ Y q then P e,q,m := P e,q,m [ [p]; P e,q,m := P e,q,m [ [p]; return; 3: if width([p]) < " then P e,q,m := P e,q,m [ [p]; return; 4: SIVIA([f e,m ]([p]), L[p], q, ", P e,q,m , P e,q,m ) 5: SIVIA([f e,m ]([p]), R[p], q, ", P e,q,m , P e,q,m ) obtained. Undetermined boxes [p] are bisected into two subboxes L[p] and R[p] on which SIVIA is applied again, until their width is less than some precision parameter ".

Contractor

Indetermination often results from range overestimation by inclusion functions. As a consequence, boxes have to be bisected by SIVIA many times to allow one to conclude on the position of the resulting boxes with respect to P e,q,m . This may entail a significant computational complexity, even for a moderate dimension of p.

Contractors may be useful in this context.

One may rewrite Equation 10.14 as follows P e,q,m = P + \ P (10.17)

P + = {p 2 P 0 |f + e,m (p) 2 Y + q } = P 0 \ f + r,m 1 Y + q (10.18) P = {p 2 P 0 |f e,m (p) 2 Y q } = P 0 \ f r,m 1 Y q . (10.19) A contractor C f + e,m ,Y +
q associated, e.g., with the generic set-inversion problem Equa-126 10.5. SIMULATION SETUP tion 10.18 is a function taking a box [p] as input and returning a box

C f + e,m ,Y + q ([p]) ⇢ [p] such that [p] \ P + = C f + e,m ,Y + q ([
p]) \ P + , so no part of P + in [p] is lost after contraction. Contractors eliminate parts of the candidate box [p] that do not belong to P + , without bisection.

Since f +

e,m and f e,m are not differentiable, most classical contractors, such as the forward-backward contractor, Newton's contractor, etc. [START_REF] Jaulin | Applied interval analysis: with examples in parameter and state estimation, robust control and robotics[END_REF] and [START_REF] Chabert | Contractor programming[END_REF]] cannot be used directly. To use contractors for the characterization of P + , considering a generic box [p], one first builds a set of m subboxes of [p]. The i-th subbox is obtained as the output of some contractor C s i,e ,[0,1[ , trying to eliminate all p 2 [p] such that s i,e (p) < 0, i = 1, ..., m. Similarly, for P , m subboxes of [p] are built, each obtained as the output of some contractor C s i,e ,] 1,0] , trying to eliminate all p 2 [p] such that s i,e (p) > 0, i = 1, ..., m. Assuming that y m (p) is differentiable, the forward-backward contractor [START_REF] Jaulin | Applied interval analysis: with examples in parameter and state estimation, robust control and robotics[END_REF], the contractor based on the centred form [START_REF] Kieffer | Guaranteed characterization of exact non-asymptotic confidence regions as defined by LSCR and SPS[END_REF], and the 3BCid contractor [START_REF] Chabert | Contractor programming[END_REF] may then be put at work. In a second step, the union of all intersections of q out of the m subboxes is evaluated and has been shown in [START_REF] Kieffer | Guaranteed characterization of exact non-asymptotic confidence regions as defined by LSCR and SPS[END_REF] to provide a contractor C f + e,m ,Y + q associated with the generic set-inversion problem Equation 10.18. A similar procedure may be considered for C f e,m ,Y q associated with the generic set-inversion problem Equation 10.19.

Simulation Setup

The simulation setup considered here is the same as that in [START_REF] Vaghefi | Cooperative received signal strength-based sensor localization with unknown transmit powers[END_REF]. Five anchor nodes (n a = 5) are placed in the corners and in the center of a square of 20 m⇥ 20 m; they are represented by black circles in the following Table 10.4 and Table 10.5. n a = 32 agents are regularly placed in the square; their location is indicated by crosses. At the beginning of the localization process, each agent broadcasts n = data are then transmitted to a central processing unit to determine an estimate of the position of the agents and its corresponding confidence region.

One has thus n y = n a ⇥ n = 50 measurements available for each agent. In the simulations, the reference power P 0 = 30 dBm at d 0 = 1 m is the same for all agents. The path loss exponent is taken as P = 4. Data are corrupted by two types of noise samples, assumed independent and identically distributed. In a first set of experiments, the noise samples are realizations of iid zero-mean log-normal variables with standard deviation 0 = 2 dBm. In a second set of experiments, the noise samples are realizations of iid Gaussian-Bernouilli-Gaussian (GBG) variables. In the latter case, with a probability p 0 = 0.9, the noise is zero-mean log-normal with standard deviation 0 = 2 dBm and with a probability p 1 = 0.1, it is zero-mean log-normal with standard deviation 1 = 5 dBm.

Three estimation problems are considered:

1. Only the location ✓ 0,i , i = 1, . . . , n a of each agent has to be estimated, P and P 0 are assumed to be known. The search box for ✓ 0,i is [0, 25] ⇥ [0, 25] m 2 . 2. ✓ 0,i and P 0,i , i = 1, . . . , n a have to be determined for each agent, the unknown reference power is not assumed to be the same for all agents, contrary to P which is assumed known. The search interval for the P 0 , i = 1, . . . , n a is [0 dBm, 40 dBm].

3. ✓ 0,i and P , i = 1, . . . , n a have to be determined for each agent, the reference power P 0,i is assumed to be known for all agents (e.g., transmitted in the packet broadcast by each agent). The search interval for P is [2, 8].

The goal is to get an estimate of the parameters in these three cases, as well as an estimate of the associated confidence region. The proposed LSCR technique is run for each agent taking q = 3, r = 1, and m = 63. To build G used in [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF], Gordon's procedure [START_REF] Gordon | Completely separating groups in subsampling[END_REF] is used, but each I j 2 G consists only of the n y r elements provided by the original procedure. The choice m = 63 10.6. SELECTION OF THE PARAMETERS OF THE LSCR APPROACH allows an exploitation of all data, and avoids having duplicate elements in G. As a consequence, (m 2q) /m ' 90 % confidence regions are provided by LSCR.

Selection of the parameters of the LSCR approach

To apply the proposed LSCR technique, several parameters have to be chosen. The most difficult to choose is e, which determines the considered statistic of the prediction residuals. Considering s = 1 and e = (1, 1) T corresponds to the correlation function considered in Equation 9.1, which provides satisfying results in the considered application, as will be seen in what follows.

Other choices are possible, but two aspects have to be taken into account. The first is related to the characterization by SIVIA of NACRs. In the inclusion function [f e,m ] of f e,m used by SIVIA, all occurrences of the parameters are considered as independent. The pessimism with which the range of f e,m is evaluated depends on the number of occurrences of each parameter to be estimated in the formal expression of f e,m , see [START_REF] Kieffer | Guaranteed characterization of exact non-asymptotic confidence regions as defined by LSCR and SPS[END_REF]. Thus, statistics with large values of P s i=0 e i should be avoided to improve the efficiency of SIVIA. The second is related to the LSCR method itself: the correlation functions should lead to small and regularlyshaped confidence regions. Some preliminary experiments have been conducted to evaluate the impact of e. One has considered s = 0 and s = 1 as well as q = 3, and m = 63. To build G used in Dalai et al. [2007b], [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF], Gordon's procedure [START_REF] Gordon | Completely separating groups in subsampling[END_REF] is used, but each I j 2 G consists only of the ny s elements provided by the original procedure. The choice m = 63 allows an exploitation of all data, and avoids having duplicate elements in G. As a consequence, (m 2q)/m ' 90% confidence regions are provided by LSCR.

Table 10.1 represents the NACR for the location ✓ 0 of an agent represented by a cross for different statistics Equation 10.6 used by LSCR. The correlation functions are defined by the vector e considering the organization of the measurement vector introduced in Equation 10.3. The shape of the confidence regions is satisfying only when e = (1, 1) T or when e = (3, 1) T . The sets in blue and red represent respectively the projections on the (✓ 1 , ✓ 2 )-plane of the inner (blue boxes) and outer (blue and red boxes) approximations, respectively of 90% NACRs.

Table 10.2 represents similar NACR but considering the following alternative organization of the measurement vector y(j) = (y(1, j), . . . , y(n a , j)) T (10.20) (10.21) In that case, the shapes of the confidence regions are satisfying for none of the vectors e. To interpret the previous results, one may consider the asymptotic analysis of the behavior of the LSCR approach provided in [7]. A condition to get P e,q,m shrinking to p ⇤ when n y = n a ⇥n tends to infinity is that for any p 6 = p ⇤ , all s j,e (p) should have the same sign as n y tends to infinity. In subsection 10.10.1 and subsection 10.10.2, one shows that this property is satisfied when considering e = (1, 1) T or e = (3, 1) T in Equation 10.6. Moreover, it is not possible to consider an organization of the measurement vector as described in Equation 10.20.

y = y T (1), . . . , y T (n) T .

Evaluation of the LSCR approach

First, one evaluates NACR as defined by LSCR for p = ✓ 0 only, using the plain SIVIA algorithm and compares its performance to a characterization via gridding, see Table 10.3. Both schemes provide sets of comparable shapes. When " decreases, the set of boxes provided by SIVIA and the set of points obtained by gridding become similar. Nevertheless, with " = 0.2 (resp. " = 0.05), characterization via SIVIA is 10 (resp. 100) times faster than via gridding. Gridding is faster than SIVIA only for large values of ", which confirms results already obtained in [START_REF] Kieffer | Guaranteed characterization of exact non-asymptotic confidence regions as defined by LSCR and SPS[END_REF]. Table 10.3: NACR as defined by LSCR: comparison of the sets obtained by gridding and using SIVIA; computing time as a function of e (left) and set of points or subpavings (right) for e = 0.2 and e = 0.05. Notice that this agent has five anchors. The number of the agent can affect the size of the confidence region too. With more anchors, one is able to obtain a smaller confidence region. Now, one aims at comparing SIVIA algorithm without contractor (NC), and variant of SIVIA involving the forward-backward contractor (FB), the 3BCid contractor (3B), and the contractor based on the centred form (CF). In all cases, these contractors are used to build, for some box [p], the m subboxes of [p] such that s e i,r (p) > 0, i = 1, ..., m and the m subboxes of [p] such that s i,e (p) 6 0, i = 1, ..., m, see Equation 10.9 and Equation 10.10. The q-intersection contractor [START_REF] Kieffer | Guaranteed characterization of exact non-asymptotic confidence regions as defined by LSCR and SPS[END_REF] is then applied on both subsets of m boxes.

Figure 10.1 shows the size of the confidence region for p = ✓ 0 when using the plain SIVIA, SIVIA with FB, 3B, and CF contractors. For each technique, two lines are provided. The upper lines represent the size of the outer approximation and the lower lines represent the size of the inner approximation of the NACR. The outer and inner approximation converge when the parameter " becomes smaller. For a given value of ", the confidence region obtained with SIVIA and using contractors has a smaller size than without contractors. The 3B and CF contractors provide the best results.

Figure 10.2 shows the evolution of size of the outer and inner approximations of the NACR as a function of the time required to obtain it with the different variants of SIVIA without and with contractors. Clearly, NC provides the best compromise between accuracy and computing time.

To understand the relative inefficiency of contractors in the context of NACR characterization, Figure 10.3 shows the contraction obtained on each box processed by SIVIA using the CF contractor. The x-axis represents the original surface of each processed box and the y-axis represents the surface of the box after contraction. One sees that large boxes are left unchanged. The CF contractor is only able to contract relatively small boxes, which helps obtaining a more accurate description of the NACR, but with a significantly increased computational complexity. The computing time penalty to apply contractors is not compensated by the gain in terms of reduction of the size of the search space provided by contractors.

Comparison with alternative methods

Table 10.4 summarizes an example of estimation results obtained with the approaches described in chapter 10 and section 9.3 considering Gaussian noise only. The columns correspond to different estimation problems. The rows correspond to the results provided, respectively by LSCR, the BE, and the RBE estimators as well as the estimate obtained using a MLE approach with confidence regions evaluated from the CRLB approach presented in subsection 9.3.3 on page 109.

The intervals for the measurement noise are taken as [ 3 0 , 3 0 ] for the BE approach and as [ 2 0 , 2 0 ] for the RBE approach. The plain SIVIA algorithm without contractors is used with the LSCR, BE, and RBE approaches with " = 0.05. The sets in green and yellow represent respectively the projections on the (✓ 1 , ✓ 2 )plane of the inner (blue boxes) and outer (blue and red boxes) approximations, respectively, of 90 % NACRs (for LSCR), or of the set estimates (for the BE or RBE approaches). For some cases, yellow part is too thin to verify. A triangle close to a node location indicates that the set provided for the node is either empty or does not contain the actual node location.

The initialization for the nonlinear MLE is taken at random within the convex hull of the anchors. Optimization is performed using Matlab lsqnonlin non-linear least-squares optimization function, which is detailed in subsection 9.3.3. With some initialization, lsqnonlin provides estimates which projection on the (✓ 1 , ✓ 2 )-plane are outside the [ 5, 25] 2 box. In that case, a new search is performed with an other random initialization. This reinitialization is performed up to three times. The green ellipses represent the 90 % confidence region derived from the CRLB. The lines link the estimated positions ✓ and true positions ✓ ⇤ of the agents.

Table 10.5 is similar to Table 10.4 but with results obtained from measurements corrupted by GBG noise. The intervals for the measurement noise are again taken as [ 3 0 , 3 0 ] for the BE approach and as [ 2 0 , 2 0 ] for the RBE approach.

Table 10.6 provides the estimates of ✓ 0 obtained with the SDP approach assuming that both ✓ 0 and P 0 have to be estimated and the projection of the 90 % 10.8. COMPARISON WITH ALTERNATIVE METHODS confidence ellipsoid derived from the CRLB in the case of Gaussian noise only (left) and Gaussian-Bernoulli-Gaussian noise (right).

Table 10.4 shows that when the number of the unknown variables increases, the size of the sets provided by the LSCR, BE, and RBE approaches increases. This is due to the fact that an uncertainty in P 0 or P translates directly in an uncertainty in terms of node localization. One also see that an uncertainty in P 0 or in P have comparable impacts on the size and shape of the obtained sets. This is consistent with the identifiability issue mentioned in section 10.5. Sometimes, the sets provided by LSCR are empty or do not contain the actual node location. This is consistent with the definition of a 90% confidence region, which should contain the actual location only in 90 % of the cases.

Considering GBG noise also increases the size of the sets obtained: confidence regions in Table 10.5 tend to be larger then those in Table 10.4. This is not surprising, since data corrupted by GBG noise contain less information about the node location than data corrupted by log-normal noise.

Nonlinear MLE provides good results when only ✓ 0 has to be estimated. Nevertheless, when both ✓ 0 and P 0 or ✓ 0 and P have to be estimated, the estimation errors increases significantly. SDP is more robust than nonlinear MLE when several parameters have to be estimated. In general, confidence regions derived from the CRLB are by far too optimistic: agents seldom belong to the 90 % confidence ellipsoid.

In Table 10.4, the sets evaluated using BE are in general larger than those obtained with RBE or than the confidence regions defined by LSCR. In the case of GBG noise, most of the sets provided by the BE approach is empty. This is not surprising, since the assumption on the noise bounds is violated in many cases. The size of the sets evaluated with the RBE approach is not much affected by the nature of the noise, showing the robustness of the approach to variations of the noise characteristics. Nevertheless, in some cases, in presence of GBG noise, the sets provided by the RBE approach do not contain the true value of the location of some nodes. This is due to the fact that not enough outliers where tolerated. Table 10.6: Projection on the (✓ 1 , ✓ 2 )-plane of estimates and confidence ellipsoids obtained by the SDP method in Gaussian noise (left) and GBG noise (right)

10.9 Effect of the noise variance

In this chapter, one studies the effect of the noise variance on the size of the confidence region and on the location estimation error. The same simulation setup as in section 10.5 is considered, but now, N = 100 agents to be localized are randomly placed in the square. Noise samples are realizations of iid zero-mean Gaussian random variables with standard deviation starting from 0 = 0.5 to 0 = 4 dBm for each simulation. Both ✓ 0,i and P 0,i , i = 1, . . . , N have to be determined for each agent.

Localization error

To evaluate the localization error of the LSCR, BE, and RBE approaches, one considers the outer approximations P r,q,m , P BE , and P RBE of P e,q,m , P BE , and P RBE and evaluate the center of mass of these sets defined, e.g., for P BE as 10.23) of the difference between the estimated location b ✓ 0,i and the actual location ✓ 0,i as a function of the measurement noise standard deviation. For the LSCR and BE techniques, the average is only evaluated when the obtained sets are not empty.

b p BE = 1 vol P BE X [p]2P BE vol ([p]) mid ([p]) , ( 10
= 1 N N X i=1 b ✓ 0,i ✓ 0,i (
For all estimation techniques, decreases with . The LSCR and RBE methods provide the best results. The performance of SDP is due to the relaxation of the constraints. MLE fails to provide reasonable estimates for some of the sensors due to local minima of the cost function.

Characteristics of the confidence region

Figure 10.5 shows the proportion of agents actually located in the set estimates or in the 90 % confidence region provided by the LSCR or derived from the CRLB. For BE, RBE, and LSCR, lower and upper bounds are provided, corresponding to the inclusion in the inner and outer approximation of the sets P BE , P RBE , and P e,q,m , respectively. Figure 10.5 shows that the confidence regions associated to the estimates provided by the MLE and SDP approaches are often too optimistic. LSCR, BE and RBE can produces confidence regions or set estimates which are reliable.

Figure 10.6 shows the average surface of the projection on the (✓ 1 , ✓ 2 )-plane of the NACR, of the BE and RBE set estimates, and of the confidence region derived from the CRLB. The size of the confidence regions provided by the CRLB (MLE and SDP approaches) are the smallest. Nevertheless, the regions are not realistic as shown in Figure 10.6. The LSCR and RBE approaches provide sets of similar size, Figure 10.5: Proportions of realizations for which the true value of the agent location is contained in the 90 % NACR (left), in the BE set estimates (up), or in the 90 % confidence regions derived from the CRLB (down) Figure 10.6: Evolution of the average surface of the projection on the (✓ 1 , ✓ 2 )-plane of the NACRs, the set estimates, or the confidence region derived from the CRLB while the size of the set provided by the BE approach is larger due to the larger noise bounds considered.

Non-Asymptotic Confidence Region Analysis

In this section, one analyzes the behavior of s j,e (p), j = 1, ..., m when n y tends to infinity, assuming that n a anchor nodes provide each the same number n of measurements tending to infinity. The parameters m and q are fixed. Our aim is to show that for properly chosen vectors e, when p = p ⇤ the functions s j,e (p), j = 1, ..., m tend to have all the same sign as n y tends to infinity. Moreover, to ensure that the NACR shrinks toward p ⇤ , one has to determine in which cases, p = p ⇤ is the unique solution of the equations s j,e (p) = 0, j = 1, ..., m. To perform this analysis, one assumes that the w j are independent and identically distributed.

First organization of the vector of measurements

One considers first an organization of the measurement vector as described in Equation 10.3 on page 122 and e = (1, 1) T . Moreover, the groups G ny s are assumed to be constructed as described in section 10.3 on page 124. For a given `, one considers the elements I ny s of G ny s with increasing n y and focus on values of n y s = kn a (m 1), with k = 1, .... In that case, for a given j, one gets

s j,e (p) = 1 |I kna(m 1) | X `2I kna(m 1)
c `,e (p).

(10.24)

Considering e = (1, 1) T , one get

c t,(1,1) (p) = w t (p)w t+1 (p), (10.25)
which, with the measurement vector (6), may take two forms, depending on whether w t (p) and w t+1 (p) involve measurements from the same anchor node (most of the cases), or from different anchor nodes, which is the case when t 2 {n, 2n, ..., (n a 1)n}. In the first case 10.26) where i t is the index of the anchor node which has produced measurement t. Similarly, in the second case, 10.29 on the preceding page may be rewritten as

c t,(1,1) (p) = ✓✓ P ⇤ 0 10 ⇤ P log 10 k✓ ⇤ 0 ✓ it k d 0 + w t ◆ ✓ P ⇤ 0 10 ⇤ P log 10 k✓ 0 ✓ it k d 0 + w t ◆◆ ✓✓ P ⇤ 0 10 ⇤ P log 10 k✓ ⇤ 0 ✓ it k d 0 + w t+1 ◆ ✓ P ⇤ 0 10 ⇤ P log 10 k✓ 0 ✓ it k d 0 + w t ◆◆ = ✓ 10 ⇤ P log 10 k✓ 0 ✓ it k k✓ ⇤ 0 ✓ it k + w t ◆ ✓ 10 ⇤ P log 10 k✓ 0 ✓ it k k✓ ⇤ 0 ✓ it k + w t+1 ◆ , (
c t,(1,1) (p) = ✓ 10 ⇤ P log 10 k✓ 0 ✓ it k k✓ ⇤ 0 ✓ it k + w t ◆ ✓ 10 ⇤ P log 10 k✓ 0 ✓ it+1 k k✓ ⇤ 0 ✓ it+1 k + w t+1 ◆ , (
s (1) j,e (p) = 1 |I kna(m 1) | na X i=1 X `2I kna(m 1) j \[(i 1)n+1in[ c `,e (p) = 1 |I kna(m 1) | na X i=1 k X =1 X `2I kna(m 1) j ,(i 1)n+1(` 1)k+<in c (` 1)k+,e (p) = 1 k |I kna(m 1) | na X i=1 X `2I kna(m 1) j X =1,(i 1)n+1(` 1)k+<in c (` 1)k+,e (p) = 1 |I kna(m 1) | na X i=1 1 k X `2I
kna(m 1) j X =1,(i 1)n+1(` 1)k+<in c (` 1)k+,e (p).

Then, using similar arguments as those considered in [Dalai et al. [2007a] Appendix A.2], for all p such that the associated ✓ 0 6 = ✓ i , i = 1, ..., n a , in the case e = (1, 1) T , one has

lim k!1 1 k X =1,(i 1)n+1(` 1)k+<in c (` 1)k+,(1,1) (p) = ✓ 10 ⇤ P log 10 k✓ 0 ✓ i k k✓ ⇤ 0 ✓ i k + w t ◆ 2 ,
(10.30) which is the expected value of c t,(1,1) (p) for some t such that i t = i.

Moreover, since s

(2) j,e (p) consists only in n a terms, for all p such that the associated

✓ 0 6 = ✓ i , i = 1, ..., n a , one has lim k!1 1 k |I na(m 1) | na X i=1 X `2I
kna(m 1) j \{n,...,(na 1)n} c `,e (p) = 0.

(10.31)

Let ↵ `,i > 0 be the proportion of indices in I ny s corresponding to c t,(1,1) (p) involving measurements of anchor node i only. When n y s = kn a (m 1), these proportions are independent of k. As a consequence, for all p such that the associated .32) which vanishes for all ✓ 0 such that .33) Provided that for each `there are at least three non-aligned anchor nodes such that ↵ `,i > 0, ✓ 0 = ✓ i is the unique solution of the system of Equation 10.33. Any ✓ 0 6 = ✓ i leads to strictly positive s 1 j,(1,1) (p).

✓ 0 6 = ✓ i , i = 1, ..., n a one has s 1 j,(1,1) (p) = lim k!1 s j,(1,1) (p) = na X i=1 ↵ `,i ✓ 10 ⇤ P log 10 k✓ 0 ✓ i k k✓ ⇤ 0 ✓ i k ◆ 2 , ( 10
k✓ 0 ✓ i k = k✓ ⇤ 0 ✓ i k , i 2 {1, . . . , n a |↵ `,i > 0}. ( 10
In the case e = (3, 1) T , using similar derivations, one obtains

s 1 j,(3,1) (p) = na X i=1 ↵ `,i ✓ 10 ⇤ P log 10 k✓ 0 ✓ i k k✓ ⇤ 0 ✓ i k ◆ 4 ,
and again, ✓ 0 = ✓ i is the unique solution of the system of Equation 10.33.

In the case e = (2, 1) T , one obtains now

s 1 j,(3,1) (p) = na X i=1 ↵ `,i ✓ 10 ⇤ P log 10 k✓ 0 ✓ i k k✓ ⇤ 0 ✓ i k ◆ 3 .
where each term in the function can be positive or negative causing multiple solution points. For this choice of e, the NACR will not shrink to p ⇤ when n y tends to infinity.

When both ✓ 0 and P 0 have to be estimated, in the case e = (1, 1)

✓ it k d 0 + w t+1 ◆ ✓ P 0 10 ⇤ P log 10 k✓ 0 ✓ it k d 0 + w t ◆◆ = ✓ P ⇤ 0 P 0 + 10 ⇤ P log 10 k✓ 0 ✓ it k k✓ ⇤ 0 ✓ it k + w t ◆ ✓ P ⇤ 0 P 0 + 10 ⇤ P log 10 k✓ 0 ✓ it k k✓ ⇤ 0 ✓ it k + w t+1
◆ .

(10.34) when measurements from the same sensors are combined. Then one gets

s 1 j,(1,1) (p) = na X i=1 ↵ `,i ✓ P ⇤ 0 P 0 + 10 ⇤ P log 10 k✓ 0 ✓ i k k✓ ⇤ 0 ✓ i k ◆ 2 , (10.35)
which vanishes for all (P 0 , ✓ 0 ) such that (10.36) This condition may be rewritten as

P ⇤ 0 P 0 + 10 ⇤ P log 10 k✓ 0 ✓ i k k✓ ⇤ 0 ✓ i k = 0, i = 1, . . . , n a .
k✓ 0 ✓ i k = k✓ ⇤ 0 ✓ i k ⇥ 10 P 0 P ⇤ 0 10 ⇤ P , i = 1, . . . , n a .
(10.37)

If P 0 is known to belong to some a priori search interval [P 0 , P 0 ], then ✓ 0 has to belong to the intersection of annulus centered in ✓ i with radius between k✓ ⇤ 0 ✓ i k ⇥ 10 P 0 P ⇤ 0 10 ⇤ P and k✓ ⇤ 0 ✓ i k ⇥ 10 P 0 P ⇤ 0 10 ⇤ P . The intersection of these annuli defines the NACR for ✓ 0 , and may also lead to a reduced confidence region for P 0 . Nevertheless, in this case, the confidence region does not shrink to a single point in the search space. Similar conclusions can be obtained when both ✓ 0 and P or ✓ 0 , P 0 , and P have to be estimated. 10.10.2 Alternative organization of the vector of measurements Consider now an organization of the measurement vector as described in Equation 10.20 and e = (1, 1) T . The expression of c t,e (p) is still given by Equation 10.25, but now, all c t,e (p) will involve measurements from two different anchor nodes. After some derivations, one obtains now:

s 1 j,(1,1) (p) = na 1 X i=1 `,i ✓ 10 ⇤ P log 10 k✓ 0 ✓ i k k✓ ⇤ 0 ✓ i k ◆ ✓ 10 ⇤ P log 10 k✓ 0 ✓ i+1 k k✓ ⇤ 0 ✓ i+1 k ◆ , (10.38) + `,na ✓ 10 ⇤ P log 10 k✓ 0 ✓ na k k✓ ⇤ 0 ✓ na k ◆ ✓ 10 ⇤ P log 10 k✓ 0 ✓ 1 k k✓ ⇤ 0 ✓ 1 k ◆ , (10.39)
where `,i is the proportion in s 1 j,(1,1) (p) of terms c t,e (p) involving a measurement provided by the i-th and the i + 1 mod n a anchor nodes. In that case, s 1 j,(1,1) (p) is positive for values of ✓ 0 far away from the anchor nodes, i.e., such that k✓ 0 ✓ i k > k✓ ⇤ 0 ✓ i k for all i = 1, . . . , n a . When ✓ 0 is close to ✓ ⇤ 0 the sign of s 1 j,(1,1) (p) is difficult to predict and unlikely to be the same for all s 1 j,(1,1) (p), j = 1, . . . m.

Conclusion

In this chapter, the localization algorithm using LSCR is discussed. With a user defined probability, the non-asymptotic confidence region NACR is defined. Originally, NACR is characterized using gridding method, where the points in the solution space are generated and verify whether it is in the NACR one by one. With the help of Interval analysis, the characterization is much faster when the required precision is high. Moreover, by using SIVIA, the outer and inner approximation of NACR is characterized, so the result is guaranteed.

This chapter also discusses the implementation and the result of LSCR localization. Considering the system parameters of LSCR, only a subset of parameters 10.11. CONCLUSION are able to generate quality results. We also show the runtime of both SIVIA and gridding under different precision. SIVIA is preferred than gridding if high precision and dimension of the solution space are considered. We also compare LSCR with some state-of-the-art approaches to show the robustness of LSCR. Finally, the error of the center of the confidence region is compared, LSCR estimates the locations of the sensors with higher accuracy than the alternative approaches. Different kinds of LSCR arrangement are analyzed and show that only a subset of parameters e lead to a robust and precise NACR.

Chapter 11

Conclusions

The second part of the thesis showed that NACR may be obtained using the LSCR approach applied to the problem of source localization from RSS measurements. LSCR provides confidence locations of the sensors under mild assumptions on the noise distribution. The noise is assumed to be symmetrically distributed around zero. Inner and outer-approximations of the NACR are provided using tools from interval analysis. Contrary to confidence regions provided by CRLB associated to MLE or SDP techniques, NACR provided by LSCR contains the true value of the parameter vector to be estimated with the prescribed confidence level. Furthermore, LSCR provides the best localization result: its average error is smaller than other compared methods (BE, RBE, MLE, SDP). This work is published in Han et al. [2018a].

Chapter 12 Conclusion for this thesis

This thesis focus on two important problems related to WSN: localization and time synchronization. For clock synchronization, robustness of the algorithm is essential. In WSN, mobile and low-cost sensors cause communication errors. Considering the dynamicity, a distributed algorithm is proposed, which extends the previous work from Függer et al. [2015b] and reduce the clock skew. In the new pulse-coupled synchronization, time-wheel algorithm is introduced for sensors to exchange the round number and average the pulses of the same round. Moreover drift compensation is applied to reduce the clock skew. The convergence analysis shows that the drift compensation reduces the clock skew and improves the convergence rate. When considering a highly dynamic WSN, the algorithm converges faster than in previous works. The proposed algorithm is successfully applied to a vehicular network.

For the second part of the thesis, we estimated the location of each sensor using LSCR, which produces a non-asymptotic confidence region. The benefit of applying LSCR to localization problem is that the assumptions on the noise of RSS signal is mild, so the confidence region can be estimated accurately. We analyze different correlation functions of LSCR and found that only some correlation function can produce a satisfying confidence region. We also implement alternative localization algorithms, such as maximum likelihood estimator (MLE), semi-definite programming (SDP), bounded error estimator (BE), robust bounded error estimator (RBE). We compare the alternative methods and found that LSCR provides the best localization result: its average error is smaller than other compared methods.

Summary and Contribution

WSNs have been widely studied in the past years, but research still needs to find solutions to some of the challenges posed by power limitation, dynamicity, and low computing clock rate. With the aim of contributing to the research on WSN, this thesis proposes new algorithms for both clock synchronization and localization. The detailed contributions of this thesis are summarized as follows.

• Chapter 3 described the formal model of the wireless sensor network and specified the problem of clock synchronization. The local physical clocks and the formal model of the wireless sensor network is defined. An error function is presented to evaluate the quality of the algorithm. Also, we summarize the existing state-of-the-art techniques such as ATS, RoATS, ATSP.

• Chapter 4 introduced time-wheel algorithm and drift compensation to extend the existing pulse-coupled synchronization [Függer et al., 2015a]. With our time-wheel algorithm, sensors can synchronize round number during consensus. With drift compensation, sensors estimate a better timing to broadcast pulse, which improves the accuracy of clock synchronization.

• Chapter 5 gave the simulation result that compare this work to the previous work [Függer et al., 2015a]. The result shows that the clock skew is reduced dramatically compared to the previous work. A random walk environment is simulated, which shows that with corr 2 i , sensors are able to remember the corrected frequency even if they do not have any neighboring sensors. A vehicle network is also simulated to show the benefit of adding corr 2 i .

• Chapter 8 introduced the tools of interval analysis. Interval analysis allow us to consider the uncertainty as bounds. With the inclusion function, the bounded 12.2. FUTURE WORK uncertainty can propagate to different domains. A contractor is able to reduce the solution space considering multiple inclusion functions. SIVIA algorithm is introduced to characterize the solution set.

• Chapter 9 described the localization algorithms. Numerical methods such as MLE and SDP with CRLB are discussed. Interval-based methods such as BE and RBE are also described. The main benefit and limitation of these algorithm is described.

• Chapter 10 proposed the localization algorithm using LSCR. With a user defined probability, the non-asymptotic confidence region NACR is defined under mild assumptions. With the help of Interval Analysis, the characterization is much faster when the required precision is high. Moreover, by using SIVIA, the outer and inner approximation of NACR is characterized, so the result is guaranteed. This chapter also discusses the implementation and the result of LSCR localization. Applying LSCR to solve localization problem achieve higher accuracy than the alternative approaches. This chapter also analyzed different kind of LSCR arrangement and find the best correlation function arrangement.

Future Work

This thesis represents a small contribution to the field of wireless sensor network. Therefore, future research is certainly needed. Some subjects where the research discussed in this thesis can continue are described in the following subsections.

Improve characterization speed of LSCR

The characterization of NACR using interval analysis is for the time being relatively time consuming. More efficient contractors are currently searched to speed-up the characterization process. Applications of the LSCR to recursive estimation is also considered. In this thesis, we used many contractors and tried to reduce the runtime.

However, the correlation functions are difficult to contractor, because there are too many occurrence of the same variable in these functions. With many occurrence, it is difficult for contractor to contract the boxes.

12.2.2 Applying LSCR to TOA or TDOA This thesis applies LSCR to RSS signal measurements, one can also apply the technique to time of arrival (TOA) or time-difference of arrival (TDOA) localization.

With TOA, the distance of sensors is calculated from the time of arrival as radiation signal travel with a known velocity. The transmission delay estimated from TOA or TDOA can also be used in clock synchronization algorithms.

Implement the Approaches in Real World

Simulation is not enough to micmic the noise samples in the real world environment.

For future work, we intend to evaluate the algorithm using a practical wireless sensor network. The first step will be to select an appropriate communication hardware and to closely simulate it in a large scale network. With the real world environment, we can start to consider the measurements that is affected by different obstacles such as buildings in the vehicular network or walls in the indoor localization problem.

Dynamic Pulse Frequency for Clock Synchronization

For the wireless sensor network, the power that drive the sensors is usually limited.

With limited battery, it is difficult for sensors to explore the real world. Our proposed clock synchronization method handles the problem of clock skew with the fixed pulse frequency. If the frequency is dynamic, then sensors are able to reduce the pulse frequency after the clock is synchronized, so the energy is preserved. 160 12.3. PUBLICATIONS

Considering Transmission Delay

In the proposed method for clock synchronization, we assume the transmission delay is zero because sensors only communicate with pulses. However, in order to communicate with pulses, the new protocol for wireless sensor network has to be setup.

If the proposed clock synchronization algorithm is applied on the regular WIFI or ZEEBEE, then transmission delay is unavoidable. In order to increase the robustness of the algorithm, we can also consider transmission delay in the future work.

Combining Approach

In the thesis, the clock synchronization and localization are two separated methods. However, many research shows that the clock synchronization and localization face the same difficulties and can be solved together. The both of the algorithms need to consider transmit delay. It is possible for sensors to exchange messages for both clock synchronization and localization at once to reduce the required bandwidth. For moving sensors, better clock synchronization reduce the clock skew, which leads to better localization result. It is interesting to merge the two different techniques to solve the problems at once, so the algorithm becomes more efficient and robust.

Publications

• Guaranteed confidence region characterization for source localization using LSCR.

-Han, Cheng-Yu, Michel Kieffer, and Alain Lambert.

-Summer Workshop on Interval Methods. 2016.

• Pulse synchronization for vehicular networks.

-Cheng-Yu Han, Thomas Nowak, and Alain Lambert.

-In Intelligent Vehicle Symposium, 2018b.

• Guaranteed confidence region characterization for source localization using RSS measurements.

-Cheng-Yu Han, Michel Kieffer, and Alain Lambert.

-Signal Processing, 152:104-117, 2018a. Pour effectuer ces tâches, la localisation et la synchronisation sont essentiels.

Figure 4

 4 Figure 4.1: When sensor i receives a pulse with a tag p j from sensor j, it appends the physical local clock time c i to the list of p j -th slot of w i .

CHAPTER 4 .

 4 TIME WHEEL ALGORITHM WITH DRIFT COMPENSATION Algorithm 4.1 State transitions for sensor i 1: if e b then 2:

Figure 4

 4 Figure 4.2: Illustration of a pulse synchronization.At the first round, sensors broadcast pulses at a similar time. In the second round, R is equal to 5 for all the sensors. However, the constant clock ticks R do not represent the constant time period due to the manufacturing defects and clock drifts. The difference of clock period causes divergence of the timing of pulses. In the third round, sensors adjust its R 0 = R + corr 1 i + corr 2 i , so that the clock skew becomes smaller. Notice that we setup k max = 3, so the round number can be 3.

  1

  Figure 5.2 and Figure 5.3 show the comparison of these three algorithms. Figure 5.2 shows the waveform of the pulses for the three algorithms. The x-axis indicates time in seconds, and the y-axis shows a value one whenever any pulse is broadcasted, otherwise it shows 0. The blue and red colors indicate the different round: k = 0 and k = 1. The thiner the width of a colored rectangle (blue or red) is, the more synchronized the pulses (of the various sensors) are, which means low clock skews. The clock skew of OCDC (the third figure of Figure 5.2) is much smaller than OC (the second figure of Figure 5.2). Figure 5.3 shows the timing of pulses for each round. Dots linked by lines with different colors indicate that the pulses are broadcasted by the same sensor.

Figure 5

 5 Figure 5.2: The figure of waveform of the pulses for all sensors. The x-axis indicates time in seconds and y-axis shows a value one whenever any pulse is broadcasted, otherwise it shows zero. The blue and red colors indicate the different round: k = 0 and k = 1. The first figure shows the waveform without any correction. The second figure shows the pulses are bounded in an interval thanks to corr 1 i . The last figure shows the pulses with much smaller clock skew thanks to both corr 1 i and corr 2 i .

  Figure 5.4: Clock skew comparison between OC algorithm (red) and OCDC algorithm (blue). The x-axis represents the round number and the y-axis represent the maximum clock skew of each round

Figure 5

 5 Figure5.5: This figure shows the five different situations. The environment is exactly the same except i and k is changed. In these five simulations, the error for method corr 1 i becomes small when the i decreases, which shows that it can handle random uncertainty of clock drift k better than corr 1 i + corr 2 i . In contrast, the error for method corr 1 i + corr 2 i increase when k increases. It shows that large amount of additional random clock drift has chance to let corr 1 i + corr 2 i become worse than corr 1 i .

Figure 5

 5 Figure 5.7: Experimental results of OC algorithm (left) and using OCDC algorithm (right). The figure on the left shows the round and the pulses using OC algorithm.There is one line diverge at the fifth second of simulation, which means that there is a sensor that does not follow the pace of the network. The figure on the right shows that with OCDC algorithm, all sensors remain stable even if the connectivity is zero because corr 2 i modifies the local clock period according to local history data. It does not need to use the pulses from neighboring sensors to compute corr 2 i , so the round remains regular even if the network is changed. If the connectivity remains zero, the pulses in the same round diverge through time if the clock skew is not zero.

Figure 5

 5 Figure 5.8: Time of each round of a random walk sensor network. The left figureshows the round and the pulses using OC algorithm. We discover that some sensors may have few or 0 neighboring node. In that case, corr 1 is difficult to correct the timing of pulses, because corr 1 only depend on the pulses from neighboring sensors. The right figure shows the round and the pulses using OCDC algorithm. It is obvious that the pulses are aligned.

Figure 5 .

 5 Figure 5.10: Static vehicle network simulation scenario. The width and height of the scenario is 1500 and 1000 meters. The pink squares represent the blocks. Black lines represent the road, and the blue circles are the position of the vehicles. The green area is the Jardin du Palais Royal.

  Figure 5.11: Static vehicle network, where the x-axis and the y-axis is the location represented in meters.

Figure 5 .

 5 Figure 5.13: Simulation result for different frequency of pulse. If T (k) is smaller, then the vehicle network converges faster then when T (k) is larger.

  Figure 8.1: Inclusion functions

  STATE-OF-THE-ART APPROACHESP RBE,⇠ = {p 2 P 0 |⌧ (p) 2 Y ⇠ } ,(9.7)

Algorithm 9. 1

 1 Error Ellipse 1: Procedure for finding error ellipse 2: Input: desired confidence probability pdf ⇤ , C x 3: Output: b 4: Compute the eigenvalues D and eigenvectors V. 5: Compute the transform matrix, A = D 1 2 V, and its inverse, A 1 . 6: Compute k according to Equation 9.45. 7: Compute some points around a circle a = k 2 [ cos(✓) sin(✓) ], ✓ = [0, ..., 2⇡] T .

y

  = y T (1) , . . . , y T (n a ) T (10.3) y m (p) = y T m (1, p) , . . . , y T m (n a , p) T (10.4) gathering in n y entries the n a ⇥ n measurements and n a ⇥ n model outputs. Considering some value p of the parameter vector, one may introduce the prediction residual w (p) = y y m (p) . (10.5)

  CONFIDENCE REGIONS AS DEFINED BY LSCR where w t (p) is the t-th component of w (p). Then, using Equation10.6, one may compute several estimates of the correlation function of the prediction residuals, for j = 1, . . of the set I ny s = {1, . . . , n y s} of indexes such that the collection

  NACR as defined by LSCR for the location of a given agent represented by a cross; the NACR are obtained for different statistics defined by e, with an organization of the measurement vector as described by Equation10.

Figure 10 . 1 :Figure 10 Figure 10

 1011010 Figure 10.1: Size of the outer (upper lines) and inner (lower lines) approximation of the confidence region as defined by LSCR provided by different variants of SIVIA, without and with contractors for different values of the precision parameter "

  considering measurements disturbed by GBG noise: Projections on the (✓ , ✓ 2 )-plane of inner and outer-approximations of set estimates obtained by SIVIA in the case NACR as defined by LSCR, and in the case of Bounded-Error (BE), and Robust Bounded-Error (RBE) estimation, as well as MLE with confidence regions derived from the CRLB; nodes for which the set estimate is empty or does not contain its actual location are marked by a triangle.

  .22) where [p] is a box contained in the list of non-overlapping boxes forming P BE , vol ([p]) and mid ([p]) are the volume and midpoint of [p]. The estimated location b✓ 0,BE is the deduced from b p BE . When an empty set is provided by the LSCR or the BE estimation techniques, no point estimate is evaluated.

Figure 10 .

 10 Figure 10.4 shows the average value of the norm

  kna(m 1) \{n,...,(na 1)n} c `,e (p) = 1 k |I na(m 1) | X `2I kna(m 1) \{n,...,(na 1)n} c `,e (p). In s (1) j,e (p), the second sum with indices `2 I kna(m 1) j T [(i 1)n + 1, in[ considers only the c `,e (p) involving measurements from the same anchor node i. In s (2) j,e (p), the sum with indices `2 I kna(m 1) j T {n, ..., (n a 1)n} considers only the c `,e (p) involving measurements from different anchor nodes. Since I kna(m 1) j contains sequences of kn a consecutive indices and I na(m 1) j contains sequences of n a consecutive indices, using I

Titre:

  Synchronisation d'horloge et localisation pour réseau de capteurs sans fil Mots clés : synchronisation d'horloge, localisation, couplage d'impulsions synchronisation, LSCR, réseau de capteurs sans fil Résumé : Les réseaux de capteurs sans fil (WSN) jouent un rôle important dans des applications telles que la surveillance de l'environnement, le suivi de sources et la surveillance médicale. Dans les WSN, les capteurs ont la capacité d'acquérir des mesures, de réaliser des traitements simples et de communiquer les résultats de ces traitements.

  Cette thèse se concentre sur deux problèmes importants liés à WSN: la localisation et la synchronisation de l'heure. Pour la synchronisation d'horloge, la robustesse de l'algorithme est essentielle. Dans WSN, les capteurs mobiles et à faible coût provoquent des erreurs de communication. Compte tenu de la dynamicité, un algorithme distribué est proposé, qui prolonge les travaux précédentsde Függer et al. (2015) et réduire le décalage d'horloge. Dans la nouvelle synchronisation à couplage d'impulsions, un algorithme de roue temporelle est introduit pour que les capteurs échangent le nombre de tours et calculent la moyenne des impulsions du même tour. De plus, une compensation de dérive est appliquée pour réduire le décalage d'horloge. L'analyse de convergence montre que la compensation de dérive réduit le décalage d'horloge et améliore le taux de convergence. Quand on considère un WSN hautement dynamique, l'algorithme converge plus rapidement que dans les travaux précédents. L'algorithme proposé est appliqué avec succès à un réseau de véhicules. Pour la deuxième partie de la thèse, nous estimons l'emplacement de chaque capteur à l'aide de LSCR, ce qui produit une région de confiance non asymptotique. L'application de la correction LSCR au problème de localisation présente l'avantage que les hypothèses sur le bruit du signal RSS sont faibles, de sorte que la région de confiance peut être estimée avec précision. Nous analysons différentes fonctions de corrélation de LSCR et avons constaté que seule une fonction de corrélation peut produire une région de confiance satisfaisante. Nous mettons également en oeuvre d'autres algorithmes de localisation, tels que l'estimateur de maximum de vraisemblance (MLE), la programmation semi-définie (SDP), l'estimateur d'erreur bornée (BE), l'estimateur robuste d'erreur bornée (RBE). Nous comparons les méthodes alternatives et avons constaté que LSCR fournit le meilleur résultat de localisation: son erreur moyenne est inférieure à celle des autres méthodes comparées.

Table 10

 10

	.4: Simulations considering measurements disturbed by Gaussian noise: Pro-
	jections on the (✓ , ✓ 2 )-plane of inner and outer-approximations of set estimates ob-
	tained by SIVIA in the case NACR as defined by LSCR, and in the case of Bounded-
	Error (BE), and Robust Bounded-Error (RBE) estimation, as well as MLE with con-
	fidence regions derived from the CRLB; nodes for which the set estimate is empty
	or does not contain its actual location are marked by a triangle.

  Les WSN ont été largement étudiés ces dernières années et la littérature scientifique rapporte de nombreux résultats qui les rendent applicables pour de nombreuses applications. Pour d'autres, la recherche doit encore trouver des solutions à certains des défis posés par la limitation énergétique, la dynamicité et la faible puissance de calcul. Dans le but de contribuer à la recherche sur les WSN, cette thèse propose de nouveaux algorithmes pour la synchronisation d'horloge et la localisation. La synchronisation d'horloge est utile pour résoudre les problèmes de localisation, de fusion de données et de consensus. En appliquant l'algorithme de synchronisation d'horloge, les capteurs etablissent un consensus temporel et travaillent donc au même rythme. Compte tenu de la dynamicité, des faibles capacités de calcul et de la parcimonie des WSN, un nouvel algorithme de synchronisation décentralisée à impulsions couplées est proposé pour réduire le décalage d'horloge de la synchronisation. L'avantage de ce type d'algorithme est que les capteurs échangent des impulsions au lieu de paquets, de sorte que non seulement la communication est efficace, mais aussi robuste à toute défaillance des capteurs dans le réseau. La localisation de capteurs a été largement étudiée dans la littérature scientifique. Cependant, la qualité et la précision de la localisation peuvent encore étre améliore. Cette thèse applique l'algorithme LSCR (Régression de régions corrélées à signes dominants) au problème de localisation. Avec LSCR, on évalue des régions de confiance avec des niveaux de confiance prescrits, qui fournissent non seulement on emplacement mais aussi la confiance en cet emplacement. Dans cette thèse, plusieurs approches de localisation sont implémentées et comparées. Le résultat de la simulation montre que, sous hypothèses modérées, LSCR obtient des résultats compétitifs par rapport à d'autres méthodes.

k = 0. With these calculations, Equation4.86 is obtained.

times a message containing its identifier. We assume that, using appropriate signaling, there is no collision between messages from different agents. Each anchor measures the RSS (in dBm) of the messages as well as the agent identifier. These

Part II

Guaranteed confidence region characterization for source localization using LSCR Chapter 8

Interval Analysis

To describe the problem formulation and the state-of-the-art solutions, it is necessary to understand the powerful tool, interval analysis, which is used in several localization approaches. Interval analysis help speed up the process for finding the confidence region. This chapter describes the fundamental concept and the tools in interval analysis.

Introduction

Interval analysis is a branch of mathematics considering intervals instead of numbers. It allows us to put bounds on rounding and measurements errors, thus produces reliable computed results. The guaranteed numerical methods are powerful tools for approximating solution sets. Guaranteed means outer and inner approximations of the sets of interest are obtained. It is possible for interval methods to characterize a set of solutions for a problem described by non-linear equations and inequalities.

The main idea to be used for interval analysis is to enclose real numbers in intervals and vectors in boxes. One of the benefits is that it is possible to obtain guaranteed results directly using interval variables.

w = wT (1) , . . . , wT (n a )

T w = w T (1) , . . . , w T (n a )

T and

[y] = [y w, y w] , (9.5) one may formulate (Equation 9.4) as a set-inversion problem

which may be characterized using SIVIA in section 8.4.

The resulting set P BE is guaranteed to contain the true value p ⇤ of the parameter vector provided that all noise samples are within the considered bounds.

Robust Bounded-Error Approach

BE approaches may provide an empty set P BE = ; as a result. Such result may be due to noise bounds which were violated by some measurements, to the initial search set which does not contain the solution set, or to the model which is a too coarse approximation of reality. RBE estimation methods [START_REF] Jaulin | Guaranteed robust nonlinear parameter bounding[END_REF][START_REF] Pronzato | Robustness to outliers of bounded-error estimators and consequences on experiment design[END_REF][START_REF] Léger | Guaranteed robust distributed estimation in a network of sensors[END_REF] have been proposed to improve the robustness to outliers of BE estimation techniques. In (Equation 9.4), the solution set is defined as the set of parameter vectors that are consistent with all measurements. With RBE estimation, one tries to find the set of parameter vectors that are consistent with all but ⇠ measurements and related noise bounds. In this context, the solution set may be formulated as

Chapter 10

Proposed LSCR Approach

This chapter proposes the localization method based on LSCR [Campi andWeyer, 2005, 2006]. The confidence region is defined and characterized using the tools from Interval analysis.

Introduction

LSCR [Campi andWeyer, 2005, 2006] defines a set P ⇢ P 0 that can be proved to contain p ⇤ with a prescribed probability. The main idea of LSCR is to consider various estimates of the correlations of the prediction residuals, i.e., the differences between the measurements and model outputs for some value p 2 P 0 . Then regions of P 0 are eliminated when too many correlation estimates have the same sign, which contradicts the hypothesis that the noise samples are independent and symmetrically distributed. The regions that were not eliminated constitute the non-asymptotic confidence region P. In this chapter, LSCR for localization is analyzed. Different configuration of the correlation functions in LSCR are discussed. Moreover, the LSCR is compared to different localization approaches: