Embed Size (px)
Citation preview
Joint Localization of Multiple Vehicles Based on RangeDifference Measurements
Vasco Manuel Nunes dos Loios Ludovico
Dissertação para a obtenção do Grau de Mestre em
Engenharia Aeroespacial
Orientador: Prof. Doutor João Pedro Castilho Pereira Santos Gomes
Júri
Presidente: Prof. Doutor João Manuel Lage de Miranda LemosOrientador: Prof. Doutor João Pedro Castilho Pereira Santos GomesVogal: Prof. Doutor Paulo Jorge Coelho Ramalho Oliveira
Maio 2014
ii
Acknowledgments
I would like to thank to my thesis advisor, Professor João Pedro Gomes, for introducing these challenges
and for the numerous suggestions, guidance and help throughout this work.
I would also like to thank to Professor João Xavier for being available to help whenever needed and for
giving me the opportunity to assist to the lectures on Convex Optimization. I am also grateful to Thomas
Furfaro at NATO STO Centre for Maritime Research and Experimentation for taking the time to answer
to my questions regarding the MORPH’12 trials.
To all my friends for making me feel alive, to Lorina for giving color to my days, to Rita for the comforting
smiles and care.
Finally, my deepest thanks to all my family for the everyday love, especially to my Mother, my Father and
my Brother for making me feel part of something bigger than me.
iii
iv
Resumo
A localização de múltiplos veículos numa constelação depende frequentemente de sinais GNSS. Con-
tudo, esta solução pode não ser viável devido às características do ambiente de operação, requisitos
de elevada precisão ou limitações orçamentais. Nesse caso, uma alternativa válida é a localização
baseada em medidas de distância entre pares.
Quando o conjunto de medidas é incompleto e a característica da solução é imposta, o problema
torna-se NP-difícil. Uma estratégia usual é relaxar as imposições do problema de forma a obter uma
formulação semi-definida que um algoritmo de otimização convexa genérico consiga tratar.
As contribuições desta tese podem ser agrupadas em quatro tópicos. a) Refinamento dos resultados
experimentais dos ensaios MORPH’12. b) Proposta de um esquema de interrogação alternativo que
permita uma redução na quantidade de informação transmitida durante o processo de determinação de
distâncias entre pares. c) Avaliação do desempenho de três formulações convexas, uma delas original,
considerando medidas ruidosas e incompletas. d) Desenvolvimento de um algoritmo de otimização
anterior à relaxação Semi-Definida em que o objetivo é encontrar o conjunto de medidas em falta que
minimizam a dimensão do problema.
Os resultados indicam que a nova formulação convexa permite um menor erro de posicionamento
ainda que a formulação clássica se apresente como a mais robusta. Para além disso, as simulações
numéricas levadas a cabo mostram que o algoritmo de otimização desenvolvido pode reduzir significa-
tivamente o erro de posicionamento embora à custa de um maior esforço computacional.
Palavras-chave: Localização, Reconstrução de constelações, Preenchimento de EDM, Otimiza-
ção Convexa, Método de Powell, Método BFGS
v
vi
Abstract
The localization of multiple vehicles in a constellation often depends on Global Navigation Satellite Sys-
tem (GNSS) signals. However, this might not be a viable solution due to environmental characteristics,
high accuracy requirements or overall cost constraints. In that case, localization based on pairwise
distances might be a valid choice.
When the set of distances is incomplete and the rank of the solution is fixed, the problem becomes
NP-hard. A popular strategy is to relax the constrains of the problem in order to obtain a semidefinite
formulation capable of being handled by a generic convex optimization solver.
The contributions of this thesis can be summarized in four major topics. a) Refinement of the experi-
mental results of the MORPH’12 trials. b) Suggestion of an alternative interrogation scheme that allows
a reduction on the communications burden during the ranging process. c) Performance evaluation of
three convex formulations, one of them original, under noisy and incomplete measurements. d) Devel-
opment of an optimization algorithm prior to the convex relaxation where the goal is to find the set of
missing distances that minimizes the embedding dimension of the problem.
Results indicate that the novel convex formulation achieves a lower positioning error although the
classical formulation remains as the most robust method. Plus, numerical simulations show that the pro-
posed completion algorithm can achieve a significant reduction on the positioning error. Unfortunately,
this is gained at the expense of a larger computational effort.
Keywords: Localization, Framework reconstruction, EDM Completion, Convex Optimization,
Powell’s Method, BFGS Method
vii
viii
Contents
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
Resumo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii
Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv
Common symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii
1 Introduction 1
1.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Thesis Outline and Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2 Wireless Sensor Network Localization 7
2.1 Localization Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.1 Range-free Localization Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.2 Range-based Localization Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Mathematical Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.1 EDM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.2 Invariance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.3 Framework Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.4 Uniquely Realizable Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.5 Euclidean distance matrix completion problem . . . . . . . . . . . . . . . . . . . . 20
2.2.6 Convex Relaxations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Localization Schemes - State of the Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3 Collaborative localization based on differential range measurements 31
3.1 Morph Experimental Results Refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.2 Interrogation Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2.1 Master/Slaves Interrogation Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2.2 Symmetric Interrogation Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2.3 Communications Burden Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3 Performance Evaluation of Several EDM Completion Convex Formulations . . . . . . . . 40
ix
3.3.1 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.4 Pre-distance Matrix Rank Minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.4.1 Powell’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.4.2 BFGS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.4.3 Dealing with Local Minima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.4.4 Performance Evaluation of the Pre-distance Matrix Rank Minimization . . . . . . . 55
4 Conclusions & Future Work 61
4.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
Bibliography 68
A Matrix Differential Calculus 69
A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
A.2 The chain rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
A.3 The first identification theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
A.4 The differential of eigenvalues: symmetric case . . . . . . . . . . . . . . . . . . . . . . . . 70
A.5 Additional fundamental rules of differential calculus . . . . . . . . . . . . . . . . . . . . . . 70
A.5.1 Differential of the Hadamard product . . . . . . . . . . . . . . . . . . . . . . . . . . 70
A.5.2 Differential of the vectorization of a matrix function . . . . . . . . . . . . . . . . . . 70
x
List of Figures
1.1 MORPH 3D Habitat Mapping. Figure taken from [1]. . . . . . . . . . . . . . . . . . . . . . 2
2.1 Illustration of the Area Localization Scheme. Figure taken from [2]. . . . . . . . . . . . . . 8
2.2 Illustration of the time based methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Illustration of the 2D Trilateration Scheme. Figure taken from [3]. . . . . . . . . . . . . . . 11
2.4 Example of localization of an ordinary node by an hyperbola-based approach. Figure
taken from [3]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.5 Example of a framework in R2. Figure adapted from [4]. . . . . . . . . . . . . . . . . . . . 14
2.6 Importance of the affine independence. Figure taken from [5]. . . . . . . . . . . . . . . . . 16
2.7 Constrained Procrustes problem example. . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.8 Illustration of a reflection over axis t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.9 Example of an universally rigid instance in R2. Figure taken from [5]. . . . . . . . . . . . . 20
2.10 Probability distribution of signal strenght with distance. Figure taken from [6]. . . . . . . . 26
2.11 Effect of the sparsification heuristics suggested in [5]. Figure taken from [5]. . . . . . . . . 28
3.1 Collaborative localization based on differential range measurements. . . . . . . . . . . . . 32
3.2 Hardware configuration of a network test node. Figure taken from [7]. . . . . . . . . . . . 36
3.3 Regenerated node tracks following the procedure in [7]. Figures adapted from [7]. . . . . 36
3.4 Regenerated node tracks following an alternative procedure. . . . . . . . . . . . . . . . . 37
3.5 Master/Slaves interrogation scheme. Example for a 5 node network. . . . . . . . . . . . . 38
3.6 Symmetric Interrogation Scheme. Example for a 3 node network. . . . . . . . . . . . . . . 38
3.7 Illustration of the reasoning behind Lower Bound EDM with Plain Ranges (EDM-RLB). . . 42
3.8 Example of a constellation reconstruction using EDM with Squared Ranges (EDM-SR).
Average noise power of 60 decibels. Constellation of 10 nodes, 3 of them are anchors. . . 43
3.9 Root Mean Square Error (RMSE) over 100 Monte Carlo runs for an average noise power
from 80 to 0 decibels and complete sets of observed distances. . . . . . . . . . . . . . . . 43
3.10 Uncertanty ellipsoids for a 10 decibels noise level. . . . . . . . . . . . . . . . . . . . . . . 44
3.11 Uncertanty ellipsoids for a 20 decibels noise level. . . . . . . . . . . . . . . . . . . . . . . 44
3.12 Uncertanty ellipsoids for a 30 decibels noise level. . . . . . . . . . . . . . . . . . . . . . . 44
3.13 RMSE over 100 Monte Carlo runs for a constant noise power of 50 decibels, a number of
anchors from 3 to 8 out of 10 nodes and complete sets of observed distances. . . . . . . . 45
xi
3.14 RMSE over 100 Monte Carlo runs for an average noise power from 80 to 40 decibels with
5 randomly chosen missing ranges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.15 RMSE over 100 Monte Carlo runs for an average noise power from 80 to 40 decibels with
10 randomly chosen missing ranges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.16 RMSE over 100 Monte Carlo runs for an average noise power from 80 to 40 decibels with
20 randomly chosen missing ranges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.17 Pairwise distances mapping for a 5 node constellation. Unknown distances: 3− 4, 2− 5,
3− 5, 4− 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.18 Constellation reconstruction using EDM-SR and 80 decibels noise level. . . . . . . . . . . 48
3.19 Example of a plot of f (c) for a 5 node planar network with missing distances between
nodes 3 − 4, 2 − 5, 3 − 5 and 4 − 5. Null-space of A is two-dimensional. Average noise
level of 80 decibels. Minimum at c = [−1.0713,−0.8198]. . . . . . . . . . . . . . . . . . . . 50
3.20 Pairwise distances mapping for a 5 node constellation. Unknown distances: 3− 4, 3− 5,
4− 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.21 Example of a plot of f (c) for a 5 node planar network with missing distances between
nodes 3− 4, 3− 5 and 4− 5. Null-space of A is one-dimensional. Average noise level of
80 decibels. Minimum at c1 = [1.2440]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
xii
List of Tables
3.1 Ratio of solved instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2 Scenario 1: Complete set. Noise level: 60 decibels. . . . . . . . . . . . . . . . . . . . . . . 57
3.3 Scenario 1: Complete set. Noise level: 80 decibels. . . . . . . . . . . . . . . . . . . . . . . 57
3.4 Scenario 2: One dimensional null space. Noise level: 60 decibels. . . . . . . . . . . . . . 57
3.5 Scenario 2: One dimensional null space. Noise level: 80 decibels. . . . . . . . . . . . . . 59
3.6 Scenario 3: Two dimensional null space. Noise level: 60 decibels. . . . . . . . . . . . . . . 59
3.7 Scenario 3: Two dimensional null space. Noise level: 80 decibels. . . . . . . . . . . . . . . 59
xiii
xiv
Abbreviations
AoA Angle-of-Arrival
EDM Euclidean Distance Matrix
EDMCP Euclidean Distance Matrix Completion Problem
EDM-R EDM with Plain Ranges
EDM-RLB Lower Bound EDM with Plain Ranges
EDM-SR EDM with Squared Ranges
GNSS Global Navigation Satellite System
GPS Global Positioning System
LS Least-Squares Method
MPT Mean Processing Time
PMRM Pre-distance Matrix Rank Minimization
RF Radio Frequency
RMSE Root Mean Square Error
RSS Received Signal Strength
RSSI Received Signal Strength Indicator
RTT Round-Trip Times
SDP Semidefinite Programming
SNL Sensor Network Localization
SNR Signal-to-Noise Ratio
SVD Singular Value Decomposition
TDOA Time-Difference-of-Arrival
TOA Time-of-Arrival
UAV Unmanned Aerial Vehicle
UWSN Underwater Wireless Sensor Network
WSN Wireless Sensor Network
xv
xvi
Common symbols
diag (A) vector formed from the diagonal of A
A�B Hadamard product of A and B
A−1 inverse of A
A⊗B Kronecker product of A and B
A+ Moore-Penrose pseudo-inverse of A
rank (A) rank of A
σr (A) rth largest singular value of A
trace (A) trace of A
AT transpose of A
vec (A) vectorization of A, columns of A stacked from left to right
EDM cone of EDM matrices
S set of symmetric matrices
D pre-distance matrix of observed squared ranges
E euclidean distance matrix
G gram matrix
I identity matrix
J orthogonal projection matrix on the orthogonal complement of the vector of ones
X set of nodes coordinates
∇ gradient operator
J set of known edges
N set of nodes
r embedding dimension
L number of edges in a framework
M number of anchor nodes
N number of nodes
1 vector of all ones
dij euclidean distance between nodes i and j
xvii
xviii
Chapter 1
Introduction
In the past, the use of unmanned vehicles have been mostly related to military applications such as
reconnaissance, surveillance and location acquisition of enemy targets. Recently, interest for unmanned
vehicles systems has begun to grow also in the direction of civil applications. To properly georeference
the gathered data and also for precise navigation, accurate localization systems are essential. First of
all, a distinction must be made between navigation and localization. Navigation is the guidance between
one point to another, localization is the ability of the robot to localize itself in a reference system. The
issues related to navigation are out of the scope of this thesis and will be left for future work.
While Global Positioning System (GPS) is one of the most popular and widely accessible positioning
technologies, it is not an universal solution. There are several good reasons why GNSS enabled devices
may not be a valid alternative:
� The GNSS signals might not be available or may be highly attenuated. For example, in underwater
environment, in the presence of obstacles that block the line-of-sight from the satellites or in indoor
applications.
� The cost of implementing GNSS receivers in all the modules might not be acceptable.
� In the case of micro vehicles, the size of the GNSS receiver and its antenna might not be accept-
able.
� Formations of vehicles may demand relative position estimates more accurate than those provided
by GNSS receivers.
The alternative to a global positioning system usually falls into one of three main categories [8]:
� Inertial/dead reckoning: Inertial positioning and navigation uses accelerometers and gyroscopes
to propagate the current state. The main drawback is the unbounded position error growth.
� Communication localization: Because the distance between the beacon and the vehicle can be
estimated by communication, the location can then be derived by localization algorithms without
extra consumption. The ranging methods are usually based on the signal strength or time-of-flight.
Alternatively, the angle-of-arrival can be the parameter of interest.
1
Figure 1.1: MORPH 3D Habitat Mapping. Figure taken from [1].
� Geophysical: Techniques that use external environmental information as references for navigation.
Methods based on communications have gained increasing popularity since no additional equipment is
required (the vehicles are expected to have communications modems) and recent advances on localiza-
tion algorithms have achieved accurate and robust solutions (see [2, 9, 10]). Although these techniques
can be applied to all classes of vehicles (terrestrial, aerial or underwater vehicles), new developments
in this field are often related with underwater applications since the methods based on time-of-flight
have shown to be particularly well suited for the localization of robots in this environment. Localization
underwater is challenging, as autonomous systems cannot rely on radio communications and global po-
sitioning. Such signals propagate only short distances underwater and therefore communications must
depend on acoustic transmissions. An interesting opportunity to implement and test time-of-flight based
localization techniques appeared with the kick off of the MORPH Project, an European Project in which
IST/ISR plays an important role.
The MORPH Project
MORPH, short for Marine System of Self Organizing, Logically Linked Physical Nodes, advances the
novel concept of an underwater robotic system composed of a number of spatially separated mobile
robot-modules, carrying distinct and yet complementary resources. Instead of being physically coupled,
the modules are connected via virtual links that rely on the flow of information among them. Inter-module
interactions are enabled by underwater acoustic communication networks at far and close ranges and
vision at very close range (see Figure 1.1).
It will provide efficient methods to map the underwater environment with great accuracy in situations
that defy existing technology, namely, underwater surveys over rugged terrain and structures with full
3D complexity, including walls with negative slope. The MORPH networks are composed of surface and
underwater vehicles, with the reception of GNSS signals limited to the surface modules. Positioning of
the remaining elements relies on measured pairwise distances to the surface vehicle(s).
The techniques explored under this thesis can be applied to any kind of autonomous vehicle in
multiple environments as long as a communication modem is available. However, and due to the strong
2
connection between IST/ISR and the MORPH project, the challenges found during the design of the
localization systems for this supra-vehicle have acted as motivation for the development of solutions
presented later on this work. Whenever possible, the main benefits that some of the proposed methods
could bring to a MORPH formation will be highlighted. Additionally, the experimental data set collected
during the MORPH’12 experiments in Faial, Azores, was also of great utility, in particular, for testing
alternative methods for reconstructing constellations of vehicles.
1.1 Problem Statement
The challenge addressed in this thesis is to localize an entire constellation of vehicles given a set of mea-
sured pairwise distances and a subset of elements with known location called anchors. It is assumed
that the network has wireless communications capabilities, in particular, communications modems ca-
pable of recording the time-of-arrival of the messages received (a common feature available in some
commercial solutions). This is fundamental since the pairwise distances will be computed based on the
difference of times-of-arrival of messages received from different vehicles. Alternatively, a simpler pro-
cedure based on Round-Trip Times (RTT) could be employed. In this case, the measurement of each
distance only requires the intervention of the two vehicles separated by the distance in question. How-
ever, this process is less efficient, in the sense that implies the transmission of a larger amount of data
in order to perform the same task. This causes unnecessary cluttering of the bandwidth, particularly
harmful in low bandwidth channels such as the underwater acoustic channels.
The vehicles dynamics will be ignored, meaning that localization will only be based on pairwise
ranges (more exactly, range differences) obtained in each measurement cycle with no additional trajec-
tory prediction steps. In this case, the vehicles can be treated as generic nodes in a Wireless Sensor
Network (WSN) taking advantage of all the techniques developed in the field of Sensor Network Lo-
calization (SNL) (see [11] for a detailed description of the problem). The SNL problem perfectly fits our
problem since it consists in locating the positions of wireless sensors, given only the (squared) Euclidean
distances between sensors that are within a given radio range and given the positions of a subset of
the sensors. The solution must meet the dimensional constraints of the problem, this is, the estimated
constellation must be planar or tridimensional according to the known geometry of the network. The
assumptions and constraints considered for this work meet the needs of the MORPH project and will be
detailed next.
Dimensional Space
The networks can assume tridimensional or planar geometries. The proposed techniques can be applied
to both cases although simulations have considered only the two dimensional scenario.
3
Transmission Range
It is considered that all the nodes are within communication range of each other, which means that
the maximum transmission range can be considered infinite. However, this does not imply that all the
messages arrive at all the listeners. Due to interferences or signal blockage, some of the packets may
be lost or corrupted, in which case the messages will be discarded.
Noise
It is assumed that range measurements between nodes are affected by white noise. The ratio between
the noise covariance of each measurement and the respective (squared) distance is considered con-
stant, following the procedure in [12].
Clock Synchronization
All the techniques proposed in this work do not require clock synchronization since the time differences
of arrival consider timestamps recorded on the same node. This implies that although clock offsets may
assume any value, clock skews must be negligible.
Line-of-sight
Communications are assumed to travel in straight line and therefore the pairwise distances are a linear
function of time intervals and sound speed, considered to be constant. In practice, refraction due to
variations in the properties of the medium and ambiguities between primary and reflected signals might
introduce measurement errors (outliers) that are not accurately modeled by the white noise.
Small Scale Networks
The convex optimization algorithms extensively used throughout the work are computationally demand-
ing and therefore not suitable for networks with more than a few dozen nodes. This is not a serious
limitation for the scenarios envisaged in the MORPH project, and in most foreseeable deployments of
underwater vehicles in the near future.
Anchor nodes
The access to absolute location is restricted to a group of nodes generically called anchors. Whether it
is by the usage of GNSS signals, other kind of localization system, or simply by the precise deployment
of the nodes (in a scenario where the modules are fixed), the position of these elements in a reference
system is known without ambiguity.
Simulation Environment
The numerical simulations conducted during this work were performed using a 32-bit version of Matlab
7.12(R2011a) on a laptop running Windows 8, with a 2.40 GHz Intel Core i3-3110M processor and with 4
4
GB of RAM. The convex optimization problems were solved using CVX (see [13]). CVX (free software)
is a Matlab-based modeling system for convex optimization allowing constraints and objectives to be
specified using standard Matlab expression syntax. This tool works as a black box, in the sense that
the user does not need to care about the internal processes that lead to the solution but only about
posing the right problem. This means that the functions and sets must obey a set of rules that guarantee
that the problem can be solved under the convex optimization formulation. Only if the computational
burden of using a generic solver were excessive for the real-time constraints of the application would it
be necessary to consider specifically tailored solution methods. These methods are outside the scope
of this Thesis.
1.2 Thesis Outline and Contributions
Chapter 2 introduces the most common localization schemes applied in WSN, and reused in this Thesis
for vehicle formations. After a brief introduction to the range-free methods, that provide only coarse
estimates of the node positions, range-based techniques are presented in more detail. The survey
addresses both ranging and location estimation stages providing a general overview of the whole signal
processing chain. Simpler algorithms, such as the Least-Squares Method (LS) and hyperbolic methods,
are also introduced and will be useful as a benchmark against the more complex convex optimization
methods. Those methods will only be introduced after the required mathematical framework is set and
some of the fundamental results on SNL are properly explored. Chapter 2 ends with a survey on the
state-of-the-art on SNL addressing in particular the work on the Euclidean Distance Matrix Completion
Problem (EDMCP) and Semidefinite Programming (SDP) relaxations.
Chapter 3 presents the main contributions to localization of vehicles using a SNL approach. In par-
ticular, it highlights the most relevant developments for the MORPH project localization system design.
These achievements can be categorized in four different topics, each one of them trying to answer to a
different challenge. The first section is concerned with the framework reconstruction procedures used
to match the estimated constellation with the anchors location. At this stage, it is considered that the
localization algorithms have already provided the relative positioning of the nodes based only on the
available pairwise distances. It is easily understood, that up to this point the positioning cannot be global
since any distance preserving map (isometry) would return a different constellation while continuing to
respect the distances between nodes. To remove the remaining degrees of freedom, the anchor nodes
positions must be brought into play. A procedure of this kind is reported in [7], applied to the data
set collected during the MORPH’12 trials. Using the same data set, it was suggested an alternative
approach that resulted in a better agreement between the position estimates and the GPS data. The
following section describes the ranging scheme used during the MORPH’12 trials (see [7]) and proposes
an alternative one. The goal is to reduce the amount of data transmitted while maintaining the ability to
determine all the pairwise distances. This was accomplished by taking advantage of all the information
flowing around the network that was ignored in the initial ranging scheme. Still in Chapter 3, Section
3.3 addresses several localization methods to handle the SNL problem when the set of measurements
5
is noisy and maybe incomplete. In fact, only simulation scenarios are able to retrieve a complete and
noiseless set of measured distances. In real environments, inaccuracies in the range measurements
cannot be avoided and the loss of one or more packets can lead to undetermined pairwise distances.
In this case, before the framework reconstructions procedures can be applied, one must find the set of
coherent (squared) Euclidean distances closest to the available data. A class of popular methods uses
the so-called SDP relaxations to perform that task. These techniques describe as a convex formulation
a problem near to the original that can be later handled by a convex optimization solver (CVX in this
case). The ability of the formulation to approximate the original problem impacts on the accuracy of the
solution. Three different SDP formulations were analyzed in order to evaluate the performance of each
one of them under noisy and incomplete sets of data scenarios. The methods under test were the clas-
sical EDM-SR (described in [14]), the EDM with Plain Ranges (EDM-R) (proposed in [14]) and a novel
formulation denominated as EDM-RLB. Results have shown that the EDM-RLB achieves more accurate
estimates than the other two procedures although the well-known EDM-SR remains as the most robust,
in the sense that it conserves the best ratio of solved problems. The last section of Chapter 3 describes
a novel algorithm whose main goal is to lower the embedding dimension of the estimated solution. The
SDP relaxations techniques used to complete the set of (squared) Euclidean distances do not constraint
the dimension of the associated constellation. For that reason, the estimated set of distances often
corresponds to constellations that lie in a higher dimension. This algorithm intends to find a completion
of the set of distances, before the SDP problem, that leads to a solution with the correct embedding
dimension. Results are encouraging with reductions in the RMSE up to 84% at the expense of a severe
increase in the processing time.
Lastly, the main conclusions and an overview on future work are drawn in Chapter 4. The work
developed during this thesis is reported in a short paper submitted to the Underwater Communications
Conference and Workshop (UCOMMS’14) [15] with a planned extension as a submission to the Journal
of Oceanic Engineering [16].
6
Chapter 2
Wireless Sensor Network Localization
When the dynamics of each vehicle is ignored, the localization of vehicles in a formation can be treated
as a classical problem of localization of generic sensors in a WSN. In that case, only the pairwise
distances between nodes in each measurement cycle are relevant for the relative positioning of the
nodes, ignoring the past positioning of each one of them. This chapter gives an overview on the WSN
localization schemes introducing the most common methods and conducting a brief survey on the state
of the art in this field. Although the main interest in this thesis lies on positioning based on pairwise
distances, alternative and complementary methods will not be ignored, since in most cases remain as
valid solutions for similar problems.
2.1 Localization Methods
The node localization methods (or self-localization in some cases) are usually categorized in range-
based and range-free schemes. The former uses the measured distance/angle to estimate the location
while the latter uses the connectivity or pattern matching methods. Next, range-free techniques will be
succinctly characterized followed by range-based schemes in more detail.
2.1.1 Range-free Localization Methods
Range-free localization schemes do not use range or bearing information to estimate distances to other
nodes which, in some cases, might be a desirable feature. However, these techniques only provide
a coarse estimate of a node’s position. Range-free schemes can be broadly classified into hopcount-
based schemes and area-based schemes [10].
Hopcount based Schemes
Hopcount-based schemes are characterized by the placement of anchor nodes at the corners or along
the boundaries of a square grid. One of the most basic range-free schemes is the DV-Hop [10]. It
first employs a classical distance vector exchange so that all the nodes in the network get distances,
7
Figure 2.1: Illustration of the Area Localization Scheme. Figure taken from [2].
in number of hops, to the anchor nodes. The number of hops is the number of nodes involved in the
transmission of a message from one node to another. Each node maintains a table and exchanges
updates only with its neighbors. Once an anchor gets distances to other marks, it estimates an average
distance for one hop, which is then propagated as a correction to the entire network. Upon receiving
the correction, an arbitrary node then estimates its distances to the landmarks which can be used to
perform triangulation. The accuracy of this method is dependent upon the density and uniformity of the
network which in actual deployments might not be the best.
Centroid Scheme
In this scheme, anchor nodes are placed as a rectangular mesh. The anchor nodes send out beacon
signals at periodic intervals with their respective locations. From the beacon signals received, a receiver
node infers proximity to a collection of anchor nodes. The location of the node is then estimated to be
the centroid of the anchor nodes that it can receive beacon packets from. A high concentration of anchor
nodes is required for this scheme to work well.
Area-based Localization
Area-based methods reduce the possible locations of each node to a certain area based on the con-
nections established with the anchors. The most relevant example comes from the Area Localization
Scheme where anchor nodes send out beacon signals at varying power levels. The multiple power levels
divide the plane into many small sub-regions defined by intersecting circles. The sensors measure the
lowest power level that they receive from each anchor node and this information is forwarded with the
sensor data to the sink for processing. By analyzing the minimum power received from all the anchors
each ordinary node can be attributed to one sub-region (see Figure 2.1). The main drawbacks of this
scheme are related with its limited coverage which is dependent on the communication range of the
reference nodes and its poor accuracy when the number of anchors is low.
8
2.1.2 Range-based Localization Methods
Range-based localization1 typically comprises two steps: determining the distances (or other spatially-
related metrics) between nodes, also known as ranging, followed by the location estimation step. Meth-
ods for ranging (see [2]) are usually based on travel times (TOA, TDOA), angles of arrival (AOA), and
received signal strength (Received Signal Strength (RSS)).
Received Signal Strength (RSS)
Ranging based on received signal strength (RSS), also commonly known as received signal strength
indicator (RSSI), measures the power of the signal at the receiver and based on the known transmit
power, the effective propagation loss can be calculated. Next, by using theoretical and empirical models
this loss is translated into distance estimates. This method has been used mainly for Radio Frequency
(RF) signals. Since the path loss in underwater acoustic channels is usually time varying and may
significantly depart from simple spherical propagation models, the RSS method is not the primary choice
for underwater localization.
Angle of Arrival (Angle-of-Arrival (AoA))
AoA estimates the angle at which signals are received to work as input to the location estimation stage.
Usually, this is accomplished by measuring the difference in the signal’s phase at each element in a
antenna array. This way, the hardware required to perform such a task is usually more complex than the
equipment used in ranging by time-of-flight or RSS and may not be feasible to embed on each sensor
whether it is by cost or volume constraints.
Time based Methods (Time-of-Arrival (TOA),Time-Difference-of-Arrival (TDOA))
These methods record the time-of-arrival (TOA) or the time-difference-of-arrival (TDOA) of transmitted
signals. The propagation time can be directly translated into distance, based on the (approximately)
known signal propagation speed. TOA methods are based on RTT and may not be the most efficient
solution for ranging measurements. Traditional TDOA methods are based on the difference of times-of-
arrival at different receivers of a same signal, in which case clock synchronization is vital. An alternative
approach, adopted in this thesis, considers the differences between times-of-arrival of signals from
different sources at the same receiver. This way, clock offsets end up canceling each other avoiding the
need of clock synchronization. Figure 2.2 illustrates both range measurements schemes. In Figure 2.2a
(TOA) the interval between the interrogation signal from node 1 and the received response from node 2
(∆T ) is equal to 2 times the propagation delay (t) plus a known processing time (tp).
The spatial metrics obtained during the ranging stage constitute the input for the location estimation
step. Given a set of metrics and a set of nodes with known coordinates called anchors, the (single-
source) location estimation step aims to determine the position of the source in a reference system. The
1 This commonly used term is somewhat of a misnomer, as it encompasses other sensing modalities besides distance measure-ments.
9
(a) TOA range measurement method.(b) Non-classical TDOA range measurementmethod.
Figure 2.2: Illustration of the time based methods.
procedures detailed here will focus on distance-based methods, neglecting mixed distances/angles al-
gorithms and range-free schemes. Before discussing advanced localization estimation algorithms, some
basic methods which underlie various localization procedures will be introduced. A detailed explanation
of each one of these methods can be found in [3].
Trilateration
Fixed pairwise distances between reference nodes and anchors result in ambiguity regions given by
circles or spheres (2D or 3D respectively). Trilateration determines the location using the intersection
point of three circles or four spheres (2D and 3D) as exemplified in Figure 2.3. However, this method
is only applicable to noiseless range measurements, otherwise additional processes are required. A
simple and popular method to deal with data inconsistencies is the LS method which extends trilateration
to noisy range measurements. The position of an unknown node is computed given the available ranges
to anchor nodes. A solution exists if the number of available anchor nodes (M ) is at least equal to r + 1
where r is the embedding dimension of the problem.
Consider the planar problem (easily extended to the 3D case), where (x, y) are the coordinates of an
unknown node and (xi, yi) are the coordinates of the anchor nodes, the least-squares problem can be
formulated as:
√(x− x1)
2+ (y − y1)
2= d1√
(x− x2)2
+ (y − y2)2
= d2...√
(x− xm)2
+ (y − ym)2
= dm
(2.1)
After squaring both sides of the equation and applying some basic algebra, the problem can be
rewritten in the canonical least-squares form:
AX = b (2.2)
10
Figure 2.3: Illustration of the 2D Trilateration Scheme. Figure taken from [3].
where X =[x y
]Tand
A =
2 (x1 − xm) 2 (y1 − ym)
......
2 (xm−1 − xm) 2 (ym−1 − ym)
(2.3)
b =
x1
2 − xm2 + y1
2 − ym2 + dm
2 − d12
...
xm−12 − xm
2 + ym−12 − ym
2 + dm2 − dm−1
2
(2.4)
The optimal solution reaches the minimum of ‖AX− b‖2 and is given by the well-known expression (see
for example [17]):
X =(ATA
)−1AT b (2.5)
where X is the optimal location of the unknown node.
Hyperbola-based Approach
A hyperbola may be defined as the locus of points where the absolute value of the difference of the
distances to two points, known as foci, is constant. Therefore, by computing the difference in times-of-
arrival from different anchors to one reference node we may find its location by intersecting hyperbola
branches instead of circles, as in the Triangulation case. This way, the existence of the beacon intersec-
tion is ensured. Figure 2.4 exemplifies this method.
Optimization Methods
Localization of single sources or multiple sources can be formulated as optimization problems, with
the key idea being to minimize the discrepancies between measured variables and their predictions
for a hypothesized set of locations. Consider for example the single-source localization problem. The
11
Figure 2.4: Example of localization of an ordinary node by an hyperbola-based approach. Figure takenfrom [3].
optimization problem for a range-based approached can be formulated as,
minimizex
M∑m=1
|‖x− sm‖p − dpm|
q (2.6)
where x and sm denote the positions of the unknown source and them-th anchor, and dm is the observed
absolute distance between them. The choice of the parameters p and q leads to different formulations
with distinct speed and convergence properties. For example, p = 1 and q = 2 actually results in the
likelihood function of observations for a Gaussian noise model. Alternatively, if we consider localization
based on range-differences, the problem can be formulated as,
minimizex
M∑m=1
|‖x− sm‖p − (dm + ‖x‖)p|q (2.7)
An optimization-based approach provides an elegant and systematic framework to deal with inconsis-
tencies in measurements and outliers (see [14]). Furthermore, it works as a support tool for the devel-
opment of efficient algorithms, including distributed localization techniques. The methodical process of
developing, test and compare different formulations is one the main advantages that have supported the
choice of following this research direction during this thesis.
A variety of numerical optimization approaches were proposed, all of them with their own limited
scope of application (see [18] for a brief overview). One particularly relevant field is the domain of con-
vex optimization, which provides a simple and efficient framework for problems with noisy and incom-
plete measurements. In recent years, we have seen a growing interest in applying convex optimization
techniques due to their accuracy and ease of formulation. The major obstacle to the implementation of
these methods is the non-convexity of most of the standard formulations, a challenge usually tackled by
relaxing the original problem to a similar but convex formulation. By doing this, it is expected that the
solution of the relaxed problem is sufficiently near the solution of the original problem.
The next section presents background material (notation and mathematical framework) for the convex
12
formulation of the SNL problem. This forms the basis for the new methods developed in Chapter 3 for
the specific problem of jointly localizing multiple vehicles from differences of packet travel times.
2.2 Mathematical Framework
Consider a network of sensor nodes with wireless communication capabilities, the SNL problem consists
in locating the positions of N wireless sensors, xj ∈ Rr, j = 1, . . . , N , given only an incomplete and
noisy set of distances between sensors and given the positions of a subset of M sensors, xj , j =
N −M + 1, . . . , N , called anchors. The embedding dimension of the problem is denoted by r, which in
practical situations is 2 or 3.
The set of exact squared distances can be arranged in matrix form giving origin to a mathematical
object called Euclidean Distance Matrix (EDM). The properties of EDM matrices will be summarized
next, following the notation in [19]. A detailed overview on Euclidean distance geometry can also be
found in [4].
2.2.1 EDM
Let Rr be a finite-dimensional Euclidean space having an inner product defined on it and thus inducing
a metric [4]. An EDM, living in R+N×N , is a complete table of squared distances d2ij between pairs of
points taken from the list of N points (xj , j = 1, . . . , N ) in Rr. The entries d2ij are given by
d2ij = ‖xi − xj‖2
2(2.8)
The row or column index of an EDM, i or j = 1, . . . , N , individually addresses all the points in the list.
Due to symmetry, from the N2 entries, only N (N − 1) /2 of these may be distinct. As an example, the
EDM for the framework in Figure 2.5 is given by:
E =
0 1 5 2
1 0 4 1
5 4 0 1
2 1 1 0
(2.9)
Gram Matrix
Before defining the properties of an EDM matrix it is worthwhile to introduce some auxiliary notation.
Consider the matrix G given by,
G = XTX (2.10)
where X = [x1, . . . , xN ] so that Gij = xTi xj . This matrix is usually referred to as the Gram matrix of the
set of points. By construction, the Gram matrix is symmetric and positive semidefinite (see [19] for a
simple proof).
13
Figure 2.5: Example of a framework in R2. Figure adapted from [4].
The diagonals of G are given by the lengths of the vectors in X, denoted by lj , j = 1, . . . , N , and
therefore the pairwise distances can be expressed as
dij = ‖xi − xj‖2
=(l2i + l2j − 2xTi xj
)1/2=(l2i + l2j − 2Gij
)1/2 (2.11)
Conversely, Gij can be expressed in terms of dij as
Gij =l2i + l2j − d
2ij
2(2.12)
or in matrix form,
G =(z1T + 1zT − E
)/2 (2.13)
where z ∈ RN is the vector of squared lengths defined by zi = l2i , 1 is the vector of all ones and E ∈ SN
is defined by Eij = d2ij . Here S denotes the set of symmetric matrices.
EDM properties
The definition introduced in the beginning of Section 2.2.1 has already imposed the first two properties
of an EDM. A matrix E ∈ SN is an EDM if and only if all its elements are non-negative (Eij ≥ 0, i, j =
1, . . . , N ) and its diagonal is composed by zeros (Eii = 0, i = 1, . . . , N ). The third condition guarantees
that the associated Gram matrix is positive semidefinite and thus, the set of pairwise distances in E is
realizable. Consider that the centroid of the constellation is located at the origin, i.e., the average of
the vectors in X is zero. Since translation is a distance-preserving map, the corresponding EDM is not
affected by this assumption. This condition validates the following expression (see [19] for the complete
proof),
JGJ = − (1/2)JEJ (2.14)
14
where J is the orthogonal projection matrix on 1⊥ (the orthogonal complement of 1) given by J =
I − 1N 11T , where I denotes the identity matrix as usual. Recall that the Gram matrix must be positive
semidefinite, implying that,
JTEJ � 0 (2.15)
From now on, a matrix that respects the conditions above will be designated by E whereas the matrix of
squared observed ranges (noisy and possibly incomplete) will be referred to as D.
2.2.2 Invariance
An EDM represents an infinity of configurations related by isometric transformations also known as
distance-preserving maps, namely, rotations, reflections and translations.
Translation
Any translation α common among all the nodes in X (X ′ = X + α1) will be canceled in the formation of
each pairwise distance. The proof results directly from definition (2.8).
Rotation
The rotation or reflection of a node constellation can be achieved through the multiplication over an
orthogonal matrix Q. This operation also does not affect the original EDM as it will be shown next.
Recall (2.13) and rearrange the terms such that,
E = z1T + 1zT − 2G (2.16)
Recalling that z is the vector of the squared lengths of the elements in X, it follows that it will remain un-
changed since the multiplication over an orthogonal matrix, composed by orthonormal vectors, preserve
the original lengths. Furthermore, in an orthogonal matrix QT = Q−1 thus,
G (QX) = XTQTQX = XTX = G (X) (2.17)
2.2.3 Framework Reconstruction
Given a complete EDM, the reconstruction of the respective node configuration is unique up to an
isometric transformation since, as explained before, the information on translation, rotation and reflection
is lost in the construction of this matrix.
To eliminate the remaining degrees of freedom (which may be necessary or not, depending on the
application) anchor nodes are required. There must be at least M ≥ r + 1 anchors with affinely inde-
pendent coordinates. In particular, 3 non-collinear anchor nodes are needed for the planar case. Figure
2.6 illustrates the prevailing ambiguity when the 3 anchors are collinear, this is, its coordinates are not
affinely independent. The ordinary node x1 could be reflected over the anchor axis while continuing to
15
Figure 2.6: Importance of the affine independence. Figure taken from [5].
respect all distances constraints. In this case, node localization is ambiguous. The procedure to extract
the node coordinates from the EDM matrix follows [20] and will be detailed next.
The first step is to reconstruct the Gram matrix up to a translation mismatch (G). From (2.14),
− (1/2)JTEJ = JTGJ
= JTXTXJ
= AT A
= G
(2.18)
The extraction of A from G is achieved via a Singular Value Decomposition (SVD) decomposition. Since
G is positive semidefinite,
G = V ΣV T
= V√
Σ√
ΣV T
= AT A
(2.19)
and therefore A =√
ΣV T . In most cases the SVD will return a coordinate matrix whose rank is greater
than the embedding dimension of the problem, so valid coordinates are obtained by truncating the SVD
to the appropriate rank, in a procedure known as best rank r approximation [21]. Alternatively, the
positions in higher dimensions can then be projected into Rr as suggested in [22]. At this point it is
worthwhile to recall equation (2.17) that highlights the invariance of the Gram matrix to a multiplications
of the matrix of nodes coordinates by an orthogonal matrix Q. In fact, any constellation given by A =
Q√
ΣV T results in the same Gram matrix G. Since the distances between nodes are preserved, the
orthogonal matrix is often seen as a rotation/reflection matrix or simple rotation when its determinant is
+1. If the purpose of the location estimation is solely to find the positioning of the nodes relatively to
each other, or in the absence of anchors, the procedure ends here. However, to achieve an absolute
node location, it is necessary to eliminate the translation, rotation and reflection ambiguities. It should be
highlighted that, although the EDMCP has ignored the anchor positions, there may be some advantages
in incorporating the localization of the anchors in the optimization problem as described in [23] and
16
[22]. The matching procedure between the estimated constellation and the true anchor positions will be
described next.
Rotation Matching
The relative orientation of the framework is adjusted by solving a constrained Procrustes problem using
a set of procedures known as the Kabsch algorithm here described without proof (see [24] for a detailed
explanation). In a constrained Procrustes problem, the goal is to find the rotation matrix Q ∈ Rr×r
(instead of a general orthogonal matrix such as in the classical Procrustes formulation) which best
approximates, in the least-squares sense, a matrix B ∈ Rm×r to A ∈ Rm×r (see Figure 2.7). Formally,
the problem can be described as,
minimize ‖A−BQ‖2F
s.t. QTQ = I
det (Q) = 1
(2.20)
where ‖ · ‖F denotes the Frobenius norm and det (·) the determinant of a matrix as usual. To ensure that
Q is a rotation matrix, not only must it be orthogonal (QTQ = I) but also its determinant must be equal
to 1. This choice avoids unwanted reflections which allows the separation between the rotation and the
reflection matching procedures. On the other hand, if matrix Q were allowed to be a general orthogonal
matrix, the reflection ambiguity would be readily solved.
The Kabsch algorithm starts with a translation step, where the centroids of both set of points, A and
B, are relocated to match the origin of the reference system. In our case, the matrix A holds the anchors
known positions and the matrix B holds the estimated anchors positions, i.e., the subset of coordinates
retrieved by the framework reconstruction step corresponding to the anchor nodes. The second step
consists of calculating the SVD of the covariance matrix BTA,
BTA = USV T (2.21)
and a squared orthogonal matrix T given by,
Tij = 0, i, j = 1, . . . , r
Tii = 1, i = 1, . . . , r − 1
Trr = sgn(
det(UV T
)) (2.22)
where sgn (·) represents the sign function. The value in Trr (−1 or 1) ensures a right-handed coordinate
system. The optimum rotation matrix is finally calculated as,
Q = UTV T (2.23)
17
Figure 2.7: Constrained Procrustes problem example.
Figure 2.8: Illustration of a reflection over axis t.
Reflection Matching
The reflection ambiguity also results from the SVD decomposition of the Gram matrix. It can be simul-
taneously solved with the rotation ambiguity, as stated in Section 2.2.3, if no determinant constraints
are imposed to matrix Q. However, if for some reason it is desirable to separate both procedures, the
reflection ambiguity can be eliminated later by solving two constraint Procrustes problems in parallel.
One over the original set of coordinates (from the EDM decomposition), and the other over a reflected
(flipped) version of the set over an arbitrary axis of reflection (see Figure 2.8) . The chosen constellation
is the one that minimizes,
‖A−Bi‖F i = 1, 2 (2.24)
where the subscript i designates the original or flipped configurations. Finally, the estimated configura-
tion Bi is translated by a vector symmetric to the initial translation of A so that the centroid of Bi matches
the centroid of A.
18
2.2.4 Uniquely Realizable Networks
Until now, the problem of estimating the position of the sensors given a set of range measurements
was formulated and some components of the solution were discussed with the tacit assumption that
the problem is solvable. Yet, the SNL problem raises a related and fundamental question regarding the
unique localizability of a given network. Assuming that not all the pairwise distances may be available,
either due to packet losses or by a deliberated suppression (see [5]), under which conditions does the
incomplete set of data assure that nodes can be unambiguously localized?
Consider for now the noiseless case, suppose that we are given a partial EDM where only the
elements in the set J are known. In addition, suppose every fully specified principal submatrix has an
embedding dimension less or equal to r. Let N be the node set N := {1, . . . , N}. An unweighted graph
G = (N ,J ) (ignoring the lengths of the available distances) together with a mapping p : N → Rr is
called a framework in Rr and is denoted by G (p). Frameworks G (p) in Rr and G (q) in Rs are called
equivalent if
‖pi − pj‖ = ‖qi − qj‖, for all ij ∈ J (2.25)
which means that both frameworks respect the same set of available distances. Furthermore, G (p) and
G (q) are congruent if
‖pi − pj‖ = ‖qi − qj‖, for all i, j = 1, . . . , N (2.26)
that is, configuration q can be obtained from configuration p by applying a rigid motion such as a trans-
lation or a rotation in Rr [25].
A framework G (p) in Rr is called globally rigid in Rr if all equivalent frameworks G (q) in Rr are
congruent to it, i.e. they satisfy condition (2.26). Similarly, a framework G (p) in Rr is called universally
rigid in Rr if, for all s = 1, . . . , N − 1, all equivalent frameworks G (q) in Rs satisfy condition (2.26)
(example in Figure 2.9).
The nomenclature detailed above establishes the required background to finally introduce the notion
of uniquely localizable networks. A given framework is uniquely r-localizable if it has a unique localization
in Rs for all s ≥ r. This definition not only restricts the existence of alternative equivalent configurations
in r but also in any dimension s ≥ r. In [5] it was formally demonstrated that for a given instance to
be uniquely r-localizable, it needs to be universally rigid and have at least r + 1 anchors with affinely
independent coordinates. It was also showed that the converse is also true.
In most existing work however, the looser notion of global rigidity is employed to characterize a given
framework, ignoring possible realizations in higher dimensions. In [26], the unique localization of a planar
network is determined using the concept of rigidity matrix. Consider once more a graph G = (N ,J )
modeling a framework in R2 of N vertices and L edges. The rigidity matrix is defined with an arbitrary
ordering of the vertices and edges and has 2N columns and L rows. Each edge gives rise to a row, and
if the edge links vertices j and k, the nonzero entries of the row of the matrix are in columns 2j − 1, 2j,
2k − 1, and 2k and are, respectively, xj − xk, yj − yk, xk − xj , yk − yj . An alternative formulation for the
rigidity matrix can be found in [25].
Making use of the rigidity matrix as an auxiliary object, the global rigidity of a graph in R2 is guar-
19
Figure 2.9: Example of an universally rigid instance in R2. Figure taken from [5].
anteed if and only if the rigidity matrix has rank 2N − 3. In that case, the dimension of the kernel of
the matrix is 3, associated with the 3 degrees of freedom of a rigid network in two dimensions. When
the kernel dimension is greater than 3, independent motions in addition to translations and rotation are
possible, corresponding to some kind of flexing. The same kind of reasoning however, cannot be applied
for dimensions r ≥ 3. As it is explained in [27] , the necessary conditions to guarantee global rigidity in
R2 are not sufficient when in R3. In fact, the authors state that rigidity in a generic r-dimensional problem
is NP-hard, a class of problems which most people believe cannot be solved in polynomial time (see [28]
for more information on NP-hard problems). Furthermore, the notion of global rigidity, although widely
spread, is mentioned in [29] as being not entirely satisfactory. Alternatively, the authors proposed an
efficient algorithm using SDP that decides whether a given network localization instance is universally
rigid. For more on the topic of SDP see [19].
The work related to rigidity theory usually considers exact range measurements, this is, the distances
between vertices are free of noise. This assumption simplifies the analysis and intuitively seems valid.
In practice, however, distance measurements will never be exact, and the equations whose solutions
deliver sensor positions in the noiseless case will no longer have a solution in the general case. The
authors in [26] argue that, in the planar case at least, if the distance measurement errors are not too
great, the associated graph is globally rigid and if three or more noncollinear anchors are available,
the network will be approximately localizable, in the sense that estimates can be found for the sensor
positions which are near the correct values. It was also stated, that the error between the true and
estimated sensor positions goes to zero continuously as the amplitude of noise perturbations in the true
distances also go to zero. Nevertheless, the minimization problem is not guaranteed to have just one
local minimum, and the determination of an effective robust algorithm is a challenging task.
2.2.5 Euclidean distance matrix completion problem
The framework reconstruction procedure introduced in Section 2.2.3 assumes a complete and noiseless
EDM as an input. In practice however, range measurements give rise to noisy and possibly incomplete
sets of distances. The location estimation step must therefore be preceded by an EDM completion
algorithm. This problem is formally known as EDMCP and it has been extensively studied in the literature
20
(see for example [22, 30–32]).
Usually, it is of great interest to constrain the solution to have an embedding dimension r in order
to meet the physical arrangement of the network. This extension of the EDMCP is known as the low-
dimensional EDMCP. As an example, if it is known that the network is planar, the estimated EDM must
have an embedding dimension of 2.
Consider once more a partial pre-distance matrix D with zero diagonal entries and with certain
nonnegative elements equal to the square of observed ranges, Dij = d2ij , and the remaining elements
considered free. The nearest EDM problem is to find an EDM E that is nearest in the least-squares
sense to matrix D, when free variables are not considered. Formally,
minimizeE
‖W � (E −D) ‖
subject to E ∈ EDM
rank (JEJ) = r
(2.27)
where W is the 0− 1 mask matrix with zeros in the entries corresponding to the free elements of D and
ones elsewhere. As usual, (�) denotes the Hadamard product and EDM denotes the EDM cone, i.e.,
the cone of matrices that respect the EDM properties. The norm in (2.27) was not specified since there
are several valid options. Details on the formulation of the EDMCP optimization problem will be given
later in this section.
Recalling the properties of an EDM matrix detailed in Section 2.2.1, matrix E belongs to the EDM
cone (more details on convex sets in [19]) if it satisfies,
Eii = 0 Eij ≥ 0 −JEJ � 0 (2.28)
The rank constraint ensures that the solution is compatible with a constellation of sensor points in Rr.
This is easily seen by remembering (2.14),
rank (JEJ) = rank
(−1
2JEJ
)= rank (JGJ)
= rank(JXTXJ
) (2.29)
For a real matrix A, rank(ATA
)= rank (A). Since J = I− 1
N 11T ,
rank (JEJ) = rank(JXTXJT
)= rank
(XJT
)= rank (X)
(2.30)
Since X ∈ Rr×N , for generic configurations, rank(XJT
)= r. The term generic configurations excludes
some particular node arrangements that by linear dependence of the nodes vector coordinates reduce
21
the embedding dimension which translates into rank(XJT
)< r.
When the set of observed distances happens to be complete, a simple way to obtain the node
coordinates is to directly reconstruct the framework from D. However, since the distances are noisy, this
procedure does not guarantee that the matrix D will satisfy the EDM properties, and thus, solving (2.27)
as an intermediate step provides more accurate solutions.
It is important to note, that in this formulation there are no references to the anchor nodes. In fact,
during the EDM completion process this information is ignored and the method proceeds the same
way whether the anchor node coordinates are available or not. The anchors will only take part in the
constellation reconstruction step (see Section 2.2.3) by removing the remaining degrees of freedom.
During the EDMCP, anchors are treated as ordinary nodes. It should be noted however that, as stated
before, some methods directly incorporate the anchors location into the EDMCP itself ([23] and [22] for
example). The theoretical background that supports the dissociation between both problems is given in
[21]. Although this seems an intuitive approach, the formal proof results from an intricate process.
The low-dimensional EDMCP as formulated before, appears as a simple and elegant mathematical
framework fitting the needs of the SNL problems. In spite of that, the problem was shown to be NP-hard
for general r > 1 (see [33]). Due to these hardness results, it becomes necessary to find alternative
ways to tackle the problem. One popular strategy is to turn to convex relaxations which can be solved
efficiently, but may not solve the original problem.
2.2.6 Convex Relaxations
Convex relaxations are one of the most powerful techniques for designing polynomial time approximation
algorithms for NP-hard optimization problems [22]. The strategy consists in loosening the constraints
of the original problem in such way that the new formulation, although much more simpler to solve,
maintains its solution sufficiently near the solution of the original problem. Some of the formulations
most relevant to the SNL problem will be described next, while a comparative numerical analysis will be
performed later in Section 3.3.
The first two formulations, EDM-SR and EDM-R, try to solve the nearest EDM problem of (2.27).
In both methods, the rank constraint is dropped, leading to SDP formulations that can be handled by
standard convex optimization software (more details on the software used for the convex optimization
problems are given in Section 3.3). This semidefinite relaxation essentially allows the points to move
into Rl where l > r. That is, if a solution of the problem has rank (JEJ) = l > r, then we have found a
Euclidean distance matrix completion of D with embedding dimension l. This may occur even if D has
a completion with embedding dimension r.
EDM with Squared Ranges
EDM-SR is a classical formulation that operates on squared ranges (see [4] for example). The Frobenius
norm is employed in (2.27) leading to,
22
minimizeE
‖W � (E −D) ‖2F
subject to E ∈ EDM(2.31)
By squaring the observed ranges in the cost function of (2.31) the impact of noise will increase [7]. To
avoid this, one may modify the cost function to accept plain distances leading to the EDM-R formulation.
EDM with Plain Ranges
The EDM-R formulation was proposed in [14] and it is obtained from the relaxed nearest EDM problem
formulation,
minimizeE
∑i,j
(√Eij −
√Dij
)2subject to E ∈ EDM
(2.32)
Expanding the objective function in (2.32) results in,
minimizeE
∑i,j
(Eij − 2
√EijDij +Dij
)subject to E ∈ EDM
(2.33)
Finally, by introducing an epigraph variable T and dropping Dij (constant),
minimizeE,T
∑i,j
(Eij − 2Tij
√Dij
)subject to T 2
ij ≤ Eij
E ∈ EDM
(2.34)
The EDM-R formulation doubles the number of variables relative to EDM-SR. Even so, in [7] it was
found to have better numerical properties enabling the use of a faster convex optimization solver. A
comparison between the two methods is carried out in Section 3.3.
When distance measurements have errors, the distance constraints are usually inconsistent and
there is no solution in Rr satisfying them. However, since the SDP approach drops the rank constraint it
may still be possible to find a set of sensor coordinates in higher-dimensional space that exactly match
the observed distances and make the objective function value zero. The optimal solution in higher-
dimensional space always results in smaller objective function values than the one constrained to lie in
Rr. A good analogy is to think about a rigid structure consisting of a set of points in a plane with specified
distances from each other. If we perturb some of the specified distances, the configuration may need to
be readjusted by setting some of the points outside the plane.
In order to encourage having a solution of the semidefinite relaxation with low rank some heuristics
have been suggested and used with great success on the SNL problem [22]. Although they assume
noiseless range measurements it is still worthwhile to explore these approaches since the key ideas can
be adapted to the nearest EDM problem.
23
Flatten Graph Heuristic
This heuristics has been suggested by [34] under the well known manifold learning problem framework
and then used by [31] on the SNL problem. The idea is try to flatten the graph associated with a
partial Euclidean distance matrix by pushing the nodes of the graph away from each other as much as
possible. This flattening of the graph then corresponds to reducing the rank of the semidefinite solution
of the relaxation. This translates into trying to maximize the objective function,
N∑i,j=1
‖xj − xi‖2 = 1TE1 (2.35)
The optimization problem follows as,
maximizeE
1TE1
subject to W � E = W �D
E ∈ EDM
(2.36)
By fixing the known pairwise distances it is assumed that they are exact and therefore there is no need
for further refinement. In [31] however, the measurement errors are taken into account, reformulating
the idea in the form of a regularization term in the cost function of the nearest EDM problem. Assuming
that (2.31) is the original formulation, the flattening strategy would result in,
minimizeE
‖W � (E −D) ‖2F − λ1TE1
subject to E ∈ EDM(2.37)
where λ is a positive regularization parameter that is to be chosen. This choice is crucial for the effective-
ness of this method. If λ is too low, the regularization term will barely influence the solution. However, if
it is chosen to be too high, the regularization term will overwhelm the error term. The authors suggest an
heuristic to estimate the optimum value for λ. After solving a non-regularized instance of the optimization
problem, the choice of λ can be obtained from,
λ∗ =‖W �
(E −D
)‖2F
1T E1(2.38)
where E was the solution for the non-regularized instance. This heuristic however requires solving the
problem twice which seems to be a drastic solution.
Nuclear Norm Heuristic
The fundamental idea behind this heuristic is to try to reduce the rank of the solution of the semidefinite
relaxation by minimizing the nuclear norm of the Gram matrix G (see [22]). The nuclear norm of a
24
general matrix A is given by
‖A‖∗ =
k∑i=1
σi (A) (2.39)
where σi (A) is the ith largest singular value of A, and rank (A) = k. Since the Gram matrix G is positive
semidefinite,
‖G‖∗ =
k∑i=1
σi (G)
=
k∑i=1
|λi (G) |
=
k∑i=1
λi (G)
= trace (G)
(2.40)
It is interesting to compare this heuristic with the flatten graph strategy. Recalling (2.16),
1TE1 = 1T(z1T + 1zT − 2G
)1
= 2n1T z − 21TG1
(2.41)
If once again X1 = 0, and remembering that z holds the norms of the vectors in X,
1TE1 = 2n1T diag (G)− 21TXTX1
= 2n trace (G)(2.42)
where diag (G) denotes the vector formed from the diagonal of G. Thus, minimizing the nuclear norm of
G produces the opposite effect intended with the flatten graph heuristic. Geometrically, it is interpreted
as trying to bring all the nodes of the graph as close to the origin as possible, which is not an intuitive
approach. Even so, both heuristics have be implemented with success as stated in [22, pp. 15-16]
although a performance comparison is lacking.
The next section conducts a brief survey on the most relevant work developed under the SNL frame-
work. The main purpose is to introduce several different localization schemes, both for terrestrial and
underwater scenarios, with special focus on range-based localization methods, particularly, those which
make use of EDM completion approaches.
2.3 Localization Schemes - State of the Art
Localization in WSN is an active area of research with a vast amount of published work on several
topics related with this subject. An introductory theoretical background was provided in Section 2.1 and
the mathematical concepts required to fully understand the reasoning behind some of the strategies
presented in this section can be found in Sections 2.1.2 and 2.2.
The extension and diversity of localization schemes proposed so far raised the need to create some
25
Figure 2.10: Probability distribution of signal strenght with distance. Figure taken from [6].
form of common ground that would help in the standardization of taxonomy, evaluation criteria and rigor-
ous definition of concepts. Some surveys fulfill those demands and unify some of the methods under a
common framework, while serving as excellent introductory material to the SNL theme. In [18], an exten-
sive survey is conducted on localization schemes covering both range-based and range-free methods.
The advantages and drawbacks of each one of them are pointed out on a merely qualitative basis. The
underwater environments and the specific challenges associated with this medium are addressed in [2]
and [10]. Finally in [9], the survey is limited to energy based localization methods, addressing problems
like Non Line of Sight broadcasting and node resource management.
Although its performance is not as good as other ranging techniques, energy based methods (usu-
ally denoted as RSSI-based) are a relatively cheap solution that does not require extra devices, as all
sensor nodes are likely to have radios with some form of received power monitoring. In [6] a robust
and distributed RSSI-based position estimation algorithm is proposed for SNL in the presence of range
measurement inaccuracies. The authors consider a probabilistic approach that can be easily geared to
various environmental conditions. Figure 2.10 gives an example of a probability distribution of the signal
strength with distance. These data, obtained from experimental measurements, form the basis of the
proposed algorithm.
The SNL problem can also be restricted to a single source localization problem, where only one
source has unknown coordinates while all the other nodes act like anchors. The work in [35] introduces
a number of exact and approximate solutions for the single source localization problem based on least-
squares approaches. Numerical results show that the proposed algorithms outperform other existing
solutions. These ideas were adapted in [7] to obtain fast, robust and flexible methods that can be used in
underwater waveguides. In [12], the main concern lies on overcoming the need of clock synchronization
between anchors and a single source. Furthermore, the possibility of the source being compromised
and possibly providing erroneous timestamps is addressed. A round-trip time scheme followed by a
least-squares or an iterative least-squares fit is proposed, showing a good accuracy while ignoring the
timestamps from the unknown source. A closed-form solution with low complexity was preferred over
an alternative maximum likelihood estimator, leading to lower processing times and simpler hardware
requirements that better fit the low complexity and low cost of current nodes in (most) sensor networks.
Many range-based schemes for WSN require precise time synchronization among nodes for rang-
26
ing. In fact, the problem of time synchronization is itself not an easy one and may introduce additional
communication overhead. In [36], the authors suggest an energy efficient algorithm based on round-trip
times assuming that all the nodes are in the same broadcast domain and therefore can listen to all the
data exchange. The use of round-trip times eliminates the need for clock synchronization, converting
the computation of pairwise distances into a straightforward task. Then, by a simple least-squares pro-
cedure, node coordinates are extracted from the measured ranges. Although it may be computationally
simple, the method does not guarantee that the set of estimated distances lies on the EDM cone, some-
times resulting in poor position estimates. A comparison between a simple least-squares procedure and
an EDM completion problem is presented in Section 3.3. The strategy in [37] is radically different. A
two-stage localization scheme is proposed. The goal of the first stage is to estimate several parameter
related to range measurements such as node internal delays, clock frequency skew, clock offset and
non line-of-sight propagation. Those estimates are then used in an iterative refinement process. This
is one of the few works where these practical sources of error are estimated to provide more accurate
results.
Concerns with the scalability of solutions led to the development of an efficient localization method
specifically designed for large-scale underwater networks [38]. The problem is solved in a 3D scenario
through a recursive algorithm based on Euclidean distances that aims for low communication costs.
These ideas later evolved to a scalable localization scheme with mobility prediction (SLMP) based on
predictable mobility patterns of underwater objects in sea shore environments [39]. Taking advantage
of the fact that mobility patterns can be suitably modeled, a high accuracy range based scheme is pro-
posed. The predicted node mobility allows a decrease in data exchanges leading to low communication
cost. However, the mobility pattern only matches one specific environment and does not guarantee
good performance for other marine conditions. Simulation results for SLMP show an overall better be-
havior than a similar algorithm without mobility prediction, thus encouraging the pursuit of this research
direction.
The methods introduced so far keep the computational complexity low at the expense of relatively
coarse position estimates. In fact, the set of distances used in the positioning processes is subject to
error measurements and therefore there is no guarantee that the EDM properties are respected without
additional refinement. In addition, incomplete sets of distances are not properly exploited, although they
may contain all the information required to successfully localize all the nodes in the network. In the
last years, researchers have turned to alternative solutions based on EDMCP techniques. In particular,
SDP relaxations have been extremely successful in providing computationally efficient solutions for this
non-convex class of problems. In [22], the authors conduct a survey on the EDM completion problems
covering the most pertinent issues in this field. Some valuable information can also be found in [32].
Although that survey specifically focuses on the applications of EDMCP to the problem of molecular
conformation, it includes a comprehensive mathematical introduction and a section specifically dedi-
cated to the SNL problem.
The noiseless scenario is considered in [23], [30] and [5]. In the first one, a theoretical analysis
is conducted on uniquely localizable networks showing, for the first time, that the SDP method yields
27
an algorithm that is guaranteed to find the solution in polynomial time, if the input graph is uniquely
localizable. The proof employs SDP duality theory and properties of interior-point methods, a class of
algorithms intended to solve linear and non-linear convex optimization problems (see [19]). One early
attempt to solve the noiseless EDMCP via SDP through an interior-point method can be found in [30]. In
[21], the complexity of SNL problems for large-scale networks is kept low by eliminating the expendable
edges in a process called edge sparsification (see Figure 2.11). The algorithm identifies universally
rigid sub-graphs and eliminates all the edges that are not vital to maintain that property. By doing so,
the properties of the original framework remain unchanged while the computational complexities and
memory requirements of the SDP algorithms are significantly reduced.
(a) Before the sparsification process. (b) After the sparsification process.
Figure 2.11: Effect of the sparsification heuristics suggested in [5]. Figure taken from [5].
Noise affects range measurements in practical scenarios, leading to infeasible systems. The work
in [21] considers this case where the problem consists in finding the nearest EDM matrix to the given
range measurements. To estimate the sensor positions, an SDP relaxation together with rank reduction
techniques are applied in order to achieve a solution in low-dimensional space. Additionally, the authors
build a theoretical background that allows to treat anchors as ordinary nodes. While in [21], the SDP
was solved via a general convex optimization solver, in [11] it is showed that, in many cases, the solution
can be efficiently computed without using any SDP solver. The work in [31] explores a regularization
technique designed to achieve a low-rank solution for the SDP problem. Described in Section 2.2.6,
this heuristic is inspired in a geometric perspective and tries to stretch the graph forcing a r dimensional
solution while keeping the distances between nodes. This method can be compared to the stretching
of a towel over a table. While the towel remains with waves, its underlying dimension is 3, but when
stretched it lies on 2 dimensions only. The presence of noise further suggests the use of this method
since, as argued in the same paper, in the presence of noise the nodes estimated from the nearest EDM
problem tend to crowd together towards the center of the configuration.
SDP relaxation strategies are used in [20] and [14] solely as a initialization stage followed by a refine-
ment phase based on plain Maximum Likelihood estimation that proved to be efficient over a standard
EDMCP approach for large problems sizes. The authors propose an incremental scheme that allows
the integration of new nodes without the need for solving an entirely new EDM completion problem.
28
Completion algorithms for Gaussian and Laplacian noise densities were developed. For the initialization
stage, and under the Gaussian noise assumption, two cost-functions are suggested. The basic EDM
completion problem operating on squared ranges, and a modified formulation operating on plain ranges,
which ended up showing better accuracy under Gaussian noise and even in the presence of outliers.
Those formulations are detailed in Section 2.2.6 and constitute one of the main approaches to the SNL
problem pursued on this thesis. The Laplacian noise assumption led to the introduction of an l1 norm
on the cost-function that proved to be robust in the presence of outliers, as demonstrated by extensive
simulation results. An experimental evaluation of the algorithms proposed in [14] can be found in [40]
where a series of 3D indoor tests were realized using mixed ultrasound-RF ranging.
In [7], the plain and squared ranges EDM completion strategies of [14] were applied to an experi-
mental underwater dataset collected under the Morph project. Further refinement of these results can
be found in Section 3.1.
29
30
Chapter 3
Collaborative localization based on
differential range measurements
In single-source localization the transition from range measurements (using TOA) to range-difference
measurements (using TDOA) leads to significant changes in the geometry and properties of cost func-
tions, and also in the algorithms used to find their optimal points. As stated in Section 2.1.2, TDOA is a
natural option to consider when there is no synchronization between the transmitter and receiver nodes,
so that absolute travel times cannot be determined unless the endpoints cooperate in a two-step proce-
dure to measure round-trip times. This may or may not be feasible depending on whether the transmitter
is able or willing to collaborate with the receivers in determining those times. Even when the transmitter
collaborates, measuring RTTs to each of the anchors sequentially is a more cumbersome process than
measuring all TOAs simultaneously when clocks are synchronized, or all TDOAs to a common reference
node (which requires clock synchronization only at receiver nodes).
Multiple-source (i.e., collaborative) scenarios pose similar challenges regarding synchronization. In
our reconstruction framework the goal will be to build, at each sensor, a matrix with pairwise Euclidean
distances from which the full set of unknown positions may be extracted through factorization. In a fully
connected network configuration, envisaged for project MORPH and assumed in this Thesis, such a
matrix might be built as follows:
1. Each node sends a request packet to query its neighbours (all other nodes)
2. The neighbours sequentially reply, and distances to the originating node are determined through
RTTs
3. The node broadcasts the set ot RTTs (or a subset thereof) to make them available to other nodes
4. When all nodes have completed the query/broadcast phase the set of relevant pairwise distances
is known across the network
The number of transmissions required by this scheme is relatively large, which constitutes a problem
in underwater networks due to the scarce available bandwidth. Even for a small number of nodes in
31
!"
#" $"
%!#&%#$" %!$"
!ij(k)"'"%!#&%#$(%!$"
(a) Time differences of arrival between packets received at node j are proportional to the differences in path lengthsi→ k → j and i→ j.
!"
#" $"
%&'(")'*+"#"!"*!#",-../"
!"
#" $"
%&'(")'*+"$"!"*!$",-../"
!"
#" $"
*!#0*#$" *!$"
!"1"*!#0*#$2*!$"
34+)"*#$",5&'("6789'7*"!/"
(b) Determining all pairwise distances at node 1 in a network of 3 nodes.
Figure 3.1: Collaborative localization based on differential range measurements.
the network the shared channel may be mostly taken up by ranging-related transmissions, leaving few
resources for actual data transmission. Ranging data can be piggybacked on normal data packets
to improve the efficiency by reducing the time occupancy of the channel, but the underlying problem
remains.
An alternative TDOA-type of approach is pursued here where measured time differences between
packet arrival times provide information on the range differences between propagation paths. Unlike
TDOA, however, time differences are not measured between a sensor and a reference sensor for a given
transmission, which would require clock synchronization across nodes, but rather between different
(but related) packets received on a given node, making it immune to clock offset issues. Figure 3.1a
illustrates the concept: Node i transmits a packet that is received at several nodes, including nodes k
and j. Prompted by that event, node k sends a reply packet back to node i that is also received at node
j. Assuming that node j duly compensates for known packet durations and processing delays, the time
difference between the instants when the original transmission from node i and the reply packet from
node k are received is proportional to dik + dkj − dij , where dik denotes the pairwise distance between
nodes i and k, and similarly for the other terms.
Due to the nature of the sequential interrogation schemes considered in this Thesis and the full net-
work connectivity, any given node will overhear multiple interdependent packet transmissions from other
nodes, from which it may compute differences in arrival times that are obtained from sums/differences of
several unknown pairwise distances in the constellation. Sharing these measurements (and associated
32
linear equations) will eventually allow any node to gather enough linear equations and invert the system
to obtain the ranges. That is the essence of the differential approaches considered in the Thesis. The
approach, enabled by the collaborative nature of nodes communicating over a shared medium, is very
different from conventional TDOA techniques for single-source localization, where range differences
cannot be linearly inverted to get ranges.
Figure 3.1b shows a simple interrogation scheme for a network of 3 nodes, adopted in Section 3.1,
that allows node 1 to determine all pairwise distances. Initially, node 1 broadcasts a request that is
first answered by node 2. The distance d12 is determined from the RTT. When node 3 replies d13 is
similarly obtained from the RTT, but the packet additionally includes in the payload the locally measured
difference in arrival times between the original request from node 1 and the reply from node 2. As
argued in connection with Figure 3.1a, that interval is proportional to d12 + d23 − d13, from which node 1
determines d23 with the help of the 2 measured RTTs. This particular scheme can readily be extended
to networks with more nodes, but the precise chaining of operations to determine absolute ranges can
easily be disrupted if even a single transmitted delay is corrupted. The interrogation and processing
schemes proposed below are shown to be more resilient to data losses.
The remainder of this chapter reports contributions in different parts of the processing chain for col-
laborative localization, each addressing a different challenge. Section 3.1 examines a real data set
for a small network of 3 nodes, focusing on the reconstruction framework for determining the spatial
coordinates of nodes from (noisy) absolute pairwise ranges. These were obtained with the simple
RTT/differential scheme depicted in Figure 3.1b, followed by nearest-EDM completion and factoriza-
tion of the Gram matrix to estimate node positions as described in [7]. Up to this point the position-
ing cannot be global, since any distance preserving map (isometry) would return a different constel-
lation while continuing to respect the distances between nodes. To remove the remaining degrees of
freedom, the anchor positions must be brought into play to synchronize the estimated node positions
across different interrogation cycles. While a synchronization procedure of this kind (i.e., setting rota-
tion/reflection/translation) was already reported in [7] for the MORPH’12 data set under consideration,
this section proposes a modified approach that yields better agreement with the GPS data.
Section 3.2 describes the ranging scheme used during the MORPH’12 sea trial [7], extends it to more
than 3 nodes, and proposes an alternative one. The goal is to reduce the amount of data transmitted
while maintaining the ability to determine all pairwise distances, regenerating at any given node a so-
called Euclidean pre-distance matrix. This was accomplished by noting that the original master/slave
scheme is inefficient because it imposes an artificial hierarchy (nodes sequentially assume the role of
master to collect the required pairwise distances), whereas the natural arrangement is non-hierarchical.
In the resulting scheme nodes share the timing of events measured in their local timelines in a symmetric
way that captures all relevant data flowing on the network and avoids unnecessarily duplicating the
replies to queries.
Section 3.3 addresses several localization methods to handle the SNL problem based on pre-distance
matrices that may be noisy or incomplete. In fact, only perfect simulation scenarios are able to retrieve
a complete and noiseless set of measured distances. In real environments, inaccuracies in range mea-
33
surements cannot be avoided and the loss of one or more packets can lead to missing pairwise dis-
tances. In this case, before framework reconstructions procedures can be applied, one must find the set
of coherent (squared) Euclidean distances closest to the available data. The chosen approach adopts
semidefinite (SDP) relaxations to EDM completion problems, which can be very conveniently later han-
dled by convex optimization solvers such as CVX. Three different SDP formulations are benchmarked:
Classical EDM reconstruction with squared ranges (EDM-SR), reconstruction with plain ranges (EDM-R
[14]), and a novel formulation termed Lower bound EDM with plain ranges (EDM-RLB) that outperforms
the other two in some problem instances.
Finally, Section 3.4 addresses the topic of rank minimization of estimated EDMs, the main goal being
to lower the embedding dimension of the estimated solution. Factorization and rank truncation of the
Gram matrix attains best results if its rank is close to the dimension of the ambient space, but since rank
constraints cannot be handled in convex problems surrogate formulations are required. As an alternative
to popular approaches for approximate rank minimization using nuclear norms, a novel technique that
explicitly minimizes a key singular value of the Gram matrix is proposed. This formulation is non-convex
and iterative, thus requiring suitable initialization and leading to an increase in processing time, but it
may afford significant reductions in the RMSE of reconstructed constellations.
3.1 Morph Experimental Results Refinement
The work in [7] reports the experimental results obtained during the MORPH’12, conducted in Faial,
Azores, on July 2012. The goal was to evaluate a range-based positioning system in a real underwater
scenario comparing the position estimates with GPS data. In this section, an alternative method for
tracks synchronization is presented with improvements in the GPS/range-based positions matching.
The experimental setup available for data collection is described in detail in [7]. The system consisted
of three nodes deployed from two vessels. Each node included an acoustic modem, a depth cell and
an underwater audio recorder. The interface computer for each node was connected on the dry side
to a GPS receiver and a 2.4 GHz radio. The wireless links between the dry-side elements establish a
terrestrial network, used for oversight and control of the analogous underwater network composed of
the modems. Figure 3.2 presents the nodal hardware configuration. The two modems deployed from
a larger vessel (nodes 1 and 2) were suspended from outriggers on both sides, providing a baseline of
around 10 m perpendicular to the ship’s axis. The third modem (node 3) was in a second vessel. The
packet exchange followed the Master/Slave scheme that, as described later in Section 3.2.1, is immune
to internode clock offsets. The complete interrogation cycle lasted around 20 seconds. The sound speed
was considered to be constant and equal to 1500 m/s.
Figure 3.3 compares the range-based estimated positions following the reconstruction procedure in
[7] with the GPS location for the three nodes during part of the trajectory (Figure 3.3b). It is clear from
3.3b that the GPS from node 2 is unreliable, as the distance between 1 and 2 is variable along the path
when it should be constant (10 m) due to physical constraints. In fact, the authors reported frequent
GPS outages throughout the tests. The estimates based on range measurements are presented in
34
Figure 3.3a, merging all reconstructed constellations into a single timeline, i.e., they do not account for
which modem was in possession of the range data on a particular interrogation cycle. To estimate the
constellations at each measurement cycle, the authors started by solving an EDMCP based on plain
ranges (see Section 2.2.6) using the set of ranges measurements collected during the experimental
trials. Due to the reduced number of nodes, 3, only complete sets of distances were considered, and the
remaining ones discarded. Consider the framework reconstruction procedure described in Section 2.2.3.
Factorization of − 12J
TEJ yields the Gram matrix of node coordinates, recentered at the origin, up to an
unknown orthogonal matrix Q. To remove ambiguities, the authors combined the framework resulting
from the EDM completion step with some of the GPS information. First, for each measurement cycle,
the constellation was rotated until the pair of nodes 1 and 2 (sharing the same vessel, 10 m apart) lied
along the horizontal axis with node 2 ahead of node 1. The third node was forced have positive vertical
coordinates by allowing a reflection over the horizontal axis whenever needed. Next, from the GPS track
of node 1 the authors estimated the attitude of the larger vessel and used this information to translate
and rotate the reconstructed range-based constellation over time. The agreement with GPS data is
seen to be quite acceptable although some occasional large deviations are visible in the reconstructed
trajectories of node 3. The challenge addressed in this thesis was to develop an alternative framework
reconstruction procedure capable of eliminating that deviations most probably induced by errors in the
estimated attitude of the vessel.
The method proposed here starts with an EDM completion problem based on squared ranges, al-
though the results with plain ranges are similar. The main difference appears in the rotation and reflection
matching between the estimated constellation and GPS data. GPS data from node 2 was considered
unreliable and therefore not included in the rotation and translation matching. However, it was consid-
ered that its error was sufficiently low so that it could be used to remove the flip ambiguity. This method
starts by following the procedure described in Section 2.2.3. Two procrustes problems are solved based
on anchors 1 and 3, one with the original estimated constellation and the other with a reflected version
over an arbitrary axis. The choice between the original/reflected constellation takes into account the
GPS position of node 2. Finally, the translation matching uses the mean value between nodes 1 and 3.
Although conceptually simpler than the method in [7], it avoids the errors introduced by computing the
vessel’s attitude based on one node only. The results for the whole trajectory are presented in Figure
3.4.
3.2 Interrogation Schemes
Until now, the mechanisms to obtain the pairwise distances were not brought into discussion. This sec-
tion, introduces some interrogation schemes designed to collect the time intervals required to determine
the distances between nodes. As stated before, the signal speed is considered constant. Thus, once
the time intervals are completely determined, the computation of the distances becomes a straightfor-
ward process. The goal is to design a distributed interrogation scheme capable of collecting enough
information, so that all the nodes are able to determine all the pairwise distances with the minimum of
35
Figure 3.2: Hardware configuration of a network test node. Figure taken from [7].
(a) Range-based position estimates. (b) GPS position estimates.
Figure 3.3: Regenerated node tracks following the procedure in [7]. Figures adapted from [7].
transmissions and the maximum of redundancy of information. Throughout this work a complete net-
work scenario was assumed. In real applications, connectivity disruptions could be reduced by network
coding schemes.
The interrogation scheme must also take into account the goals and constraints of the network itself
such as latency, noise, vulnerability to interferences, among others. In high latency and very limited
bandwidth networks, such as underwater acoustic networks, the communications burden plays a critical
role and must be minimized, in particular, for large networks where the localization process requires a
huge amount of data. Next, it will be introduced two different interrogation schemes: the Master/Slaves
interrogation scheme [7] and an alternative symmetric scheme.
36
(a) Range-based position estimates. (b) GPS position estimates.
Figure 3.4: Regenerated node tracks following an alternative procedure.
3.2.1 Master/Slaves Interrogation Scheme
The Master/Slaves scheme is described in [7] and was evaluated in a real scenario during the MORPH’12
experiment. In each cycle, each one of the nodes has a position in the hierarchy, with rotating roles be-
tween cycles. The hierarchy (Master, Slave#1,...,Slave#N -1) defines the sequence of responses, with
increasing amount of data transmitted by each node. Figure 3.5 illustrates the interrogation process for
a 5 node network. The Master sends a Query at t0, Slave#1 replies after receiving the Query from the
Master (the processing times are ignored) which arrives back at the Master at t1. After receiving the
reply from Slave#1, Slave#2 replies including in the packet the time difference between the arrivals of
the reply from Slave#1 and the Query from the Master. The process goes on until all the nodes have
answered back. In the next cycle, the positions in the hierarchy shift until all the nodes have the chance
to be the Master since, in each cycle, only the Master is able to compute the distances between nodes.
It is relevant to note that time intervals always consider two instants recorded in the same local clock of
one of the nodes, which makes this process immune to clock offsets.
The time intervals ∆ij are assembled into a vector (∆) which is related with the vector of pairwise
distances (d) by a square matrix A of dimension 12
(N2 −N
)such that Ad = ∆, where N is the number
of nodes in the network. When all the time intervals are known, A is invertible, and the vector of pairwise
distances d can be determined from d = A−1∆. In the general case, d = A+∆ where A+ stands for the
pseudo-inverse of A. When some of the time intervals are missing, the solution for some of variables
d will not be unique since A will be rank deficient and its null space will not be empty. In that case,
the solutions can be expressed as the sum of a particular solution with an arbitrary element of the null
space. The variables of d with an unique solution correspond to the rows of the null space of A without
non-zero entries.
3.2.2 Symmetric Interrogation Scheme
The previous interrogation scheme does not take fully advantage of all the information flowing in the
network. Here, it is proposed an alternative scheme that shows that it is possible to reduce the amount
37
Figure 3.5: Master/Slaves interrogation scheme. Example for a 5 node network.
Figure 3.6: Symmetric Interrogation Scheme. Example for a 3 node network.
of data transmitted and still determine all the distances. This method does not establish an hierarchy in
each cycle, it is symmetric in the sense that all the nodes play the same role.
Each cycle starts with a round of queries followed by a second round where each one of the nodes
sends theN−1 time differences between arrivals. Figure 3.6 illustrates the process for a 3 node network.
As in the Master/Slaves scheme there is no need of clock synchronization.
Let ∆ be the vector grouping the time intervals and x the vector of unknown variables, containing not
only the pairwise distances divided by the signal speed but also the instants at which each one of the
nodes have sent its Query (tii). The vectors ∆ and x are related through a matrix A such that Ax = ∆.
A has dimensions(N2 −N, N
2+N2
)and therefore the system is overdetermined when all the ∆ entries
are known for networks with more than 3 nodes. In that case, x results from x = A+∆, where A+ stands
for the pseudo-inverse of A since the matrix is not invertible in the general case. Note that one of the
variables tii must be assumed as the origin of the time axis, causing x to lose a variable and A to lose
a column. Typically, each node will consider the moment at which have sent the Query (tii) as the origin
of the time axis. After this step, A becomes a full rank matrix when the sufficient ∆ entries are known.
The set of ∆ entries however, might not always be complete. Packets losses may lead to unknown
∆ entries, a scenario particularly relevant in harsh environments such as the underwater acoustic chan-
38
nels. Furthermore, some of the entries may not be transmitted on purpose, reducing the transmission
delays, while maintaining the vital data flowing around the network. A similar technique, called punc-
turing, is often employed in digital coding. In the case where some of the entries of ∆ are unknown,
the rows of A corresponding to the unknown ∆ entries are eliminated. If due to the row eliminations, A
becomes rank deficient, then for some of the variables in x the solution will not be unique since the null
space of A will not be empty as described before. Again, the solutions can be expressed as the sum of a
particular solution with an arbitrary element of the null space. The variables of x with an unique solution
correspond to the rows of the null space of A without non-zero entries.
3.2.3 Communications Burden Analysis
The minimum amount of data exchanged so that all the nodes are able to determine the complete set of
pairwise distances depends on the interrogation scheme in use. Here, the amount of data transmitted
is compared between the two methods analyzed before and the classical TOA scheme based on RTT
(see Section 2.1.2).
Classical Time-of-Arrival
The standard TOA method is based on RTT and imposes no hierarchy to the network. Therefore, all
the nodes determine the distances between themselves and the other elements. The RTT collected by
each node are then communicated to the network. In this case, there are N (N − 1) 2 = 2N2 − 2N
queries and N (N − 1) RTT transmitted per cycle. Therefore the amount of data transmitted depends
on O(N2)
.
Master/Slaves Interrogation Scheme
Only one node per cycle (Master) is able to determine the location of its peers which implies that the
amount of data transmitted for one cycle needs to be multiplied by the number of nodes (N ). There is
only one query andN∑
n=2(n− 2) time differences transmitted per cycle.
N∑n=2
(n− 2) =(N − 1)
2(N − 2) =
N2 − 3N + 2
2(3.1)
Therefore the amount of data transmitted depends on O(N3)
.
Symmetric Interrogation Scheme
In this scheme all the nodes determine the network configuration in just one cycle. The number of
queries needed is N and the number of time differences transmitted is N (N − 1) = N2 −N . Therefore
the amount of data transmitted depends on O(N2)
. This is a significant reduction when compared
with the Master/Slaves method. It is easy to see, that the Symmetric scheme is also preferable when
compared with the classical TOA since the required number of queries is lower although the number of
39
transmitted time-differences is similar. From now on, and whenever an interrogation scheme is neces-
sary, this will be the primary choice.
3.3 Performance Evaluation of Several EDM Completion Convex
Formulations
As described in Section 2.2.6 the EDMCP can be solved via convex relaxations by reducing the original
problem to a SDP formulation that can be solved in polynomial time. This section evaluates three
different formulations and justifies its usage over simpler methods such as a raw data approach or a LS
approach.
The convex optimization problems were solved using CVX with the standard SDPT3 solver. For the
RMSE only the solved runs were accounted for, this is, results from the CVX other than solved were
ignored (more on the CVX numerical results in [13]). All the simulations considered a planar constellation
with 10 nodes randomly placed on a unit square area with randomly chosen anchors. These scenarios
only contemplate the presence of Gaussian noise affecting the real distances. The level of noise in each
run is controlled by the average noise power defined as σ2 = 1/PP∑i=1
σ2i (see [12]), where σ2
i is the
noise covariance for the distance i and P is the number of pairs in the network P = CN2 = N(N−1)
2 . It
is assumed that σ2i is related to the distances according to the path loss law, so that σ2
i /d2i is a constant
ratio as suggested in [12].
The performance criterion is the RMSE of the estimated positions over a maximum of 100 Monte
Carlo runs which can be expressed as
RMSE =
√√√√ 1
MN
M∑j=1
‖X −X‖2F (3.2)
where M is the number of successful (solved) Monte Carlo runs, ‖ · ‖F denotes the Frobenius norm as
usual, X the estimated constellation and X the true constellation.
The formulations under analysis are EDM-SR and EDM-R described in Section 2.2.6 and also an
alternative formulation proposed next, under the name of EDM-RLB. To justify the usage of convex
optimization problem formulations over other simpler techniques, whenever feasible, the results using
a LS approach and a straightforward framework reconstruction will be presented as well. By this, it is
meant a direct framework reconstruction without a previous EDM completion procedure, in the cases of
course, where the set of observed distances is complete.
Finally, it should be noted that all the nodes in a constellation contribute to the computation of the
RMSE. Anchor nodes are not ignored, since for the convex optimization techniques, the estimated an-
chor positions are not forced to match the true anchor locations introducing errors that should be ac-
counted for.
40
Lower bound EDM with plain ranges
This formulation follows the same steps than EDM-R, starting from the nearest EDM problem,
minimizeE
∑i,j
(√Eij −
√Dij
)2subject to E ∈ EDM
rank (JEJ) = r
(3.3)
Again, by expanding the objective function, dropping the rank constraint, introducing an epigraph variable
T and ignoring constant terms,
minimizeE,T
∑i,j
(Eij − 2Tij
√Dij
)subject to T 2
ij = Eij
E ∈ EDM
(3.4)
The difference now lies in the relaxation of T . Instead of imposing an absolute value upper bound such
as in EDM-R, where T 2ij ≤ Eij , variable T is now restricted to assume positive values, upper bounded
once again by√Eij and lower bounded by the line between 0 and the estimated maximum value for√
Eij given by√Emax for Eij = Emax. Figure 3.7 illustrates this reasoning. The slope of the linear
lower bound is given by√
Emax
Emaxand therefore the constraints in T result in Eij√
Emax
≤ T ≤√Eij . The
EDM-RLB formulation is then given by,
minimizeE,T
∑i,j
(Eij − 2Tij
√Dij
)subject to Tij ≤
√Eij
Tij ≥Eij√Emax
E ∈ EDM
(3.5)
3.3.1 Simulation Results
The first experiment considered a scenario where the set of observed distances is complete but noisy,
with average noise powers ranging from 80 to 0 decibels. The goal is to evaluate the ability of each
formulation to deal with noisy measurements. The number of anchor nodes is fixed, 3 out of 10 nodes,
and the constellation changes randomly between runs. Figure 3.8 exemplifies the results for EDM-SR
with an average noise power of 60 decibels.
The results over 100 Monte Carlo runs are shown in Figure 3.9. The accuracy clearly degrades with
the increase in noise. Since the experiments consider a square of unit area, a RMSE near 1 means
the framework reconstruction has no practical value. As expected, the LS (see Section 2.1.2) and
the Raw data approaches exhibit lower accuracy than the convex optimization methods. And although
41
Figure 3.7: Illustration of the reasoning behind EDM-RLB.
the performance gap between the Raw method and the alternative convex approaches is small, it is
worthwhile to recall that a direct reconstruction from the set of observed ranges can only be performed
when all the pairwise ranges are available, thus making it an infeasible option for the general case.
EDM-SR and EDM-R exhibit similar results. In fact, for almost every noise level the points of the two
methods are coincident. For EDM-RLB however, the estimates are slightly better, which suggests the
use of this formulation over the others, at least, in the absence of missing ranges.
Similar conclusions can be drawn from Figures 3.10, 3.11 and 3.12 where the 95% uncertainty ellipses
were plotted for a randomly created constellation for average noise levels of 10, 20, and 30 decibels.
Again, the results are the outcome of 100 Monte Carlo runs with a complete but noisy set of distances.
Obviously, the constellation (including anchor nodes) is fixed.
An interesting phenomenon is related to the dependency with the scale verified in the formulations
based on plain ranges, namely EDM-R and EDM-RLB. In fact, once the distances between nodes start
growing the results achieved by the convex optimization solver may no longer be reliable even in the
absence of noise. This implies that, in practical situations, a preprocessing step must be carried out
before the EDMCP to ensure proper scaling of measured distances for the solver.
The effect of the number of anchors in the accuracy of the results is presented in Figure 3.13. As
expected, expanding the set of anchor nodes is advantageous only when the number of existing anchors
is significantly small, after which it has no significant effect. In this particular case, for LS, when the
number of anchors changes from 3 to 5 the RMSE decreases more than one order of magnitude. For
the other methods, however, only the first change, from 3 to 4 anchors, brings some significant reduction
in RMSE. The higher dependence of LS with the number of anchors can be explained by the role that
anchor nodes play in each method. LS is based on a triangulation procedure between anchors and
ordinary nodes, whereas the remaining techniques only make use of the anchors after the framework
reconstruction step to remove the translation and rotation/reflection ambiguities.
The convex optimization techniques were also tested assuming that some of the distances were un-
known, this is, when for some reason the set of observed ranges is incomplete. As in the previous case,
the simulations encompass a set of 100 Monte Carlo runs with randomly generated constellations and
42
Figure 3.8: Example of a constellation reconstruction using EDM-SR. Average noise power of 60 deci-bels. Constellation of 10 nodes, 3 of them are anchors.
Figure 3.9: RMSE over 100 Monte Carlo runs for an average noise power from 80 to 0 decibels andcomplete sets of observed distances.
an average noise level ranging from 80 to 40 decibels. Figures 3.14, 3.15 and 3.16 consider the cases
where 5, 10 and 20 randomly chosen distances are missing. The results for 5 missing distances (Figure
3.14) resemble those for a complete set of observed distances in Figure 3.9, showing a slight perfor-
mance advantage of plain ranges techniques (EDM-R and EDM-RLB) over EDM-SR. In subsequent
simulations, with 10 and 20 distances missing, the performance of the three methods can be considered
to be approximately equivalent. It is worthwhile to point out that for 20 missing distances (out of 45) the
RMSE is constant and relatively high throughout all the noise levels, indicating that these methods are
no longer reliable. A quick remark should also be made on the ratio of solved instances presented in
Table 3.1. EDM-SR stands out as the most robust formulation, while the plain ranges methods (EDM-R
and EDM-RLB) were unable to solve a significant share of the problems when the number of missing
distances is high (10 and 20).
43
Figure 3.10: Uncertanty ellipsoids for a 10 decibels noise level.
Figure 3.11: Uncertanty ellipsoids for a 20 decibels noise level.
Figure 3.12: Uncertanty ellipsoids for a 30 decibels noise level.
44
Figure 3.13: RMSE over 100 Monte Carlo runs for a constant noise power of 50 decibels, a number ofanchors from 3 to 8 out of 10 nodes and complete sets of observed distances.
Figure 3.14: RMSE over 100 Monte Carlo runs for an average noise power from 80 to 40 decibels with 5randomly chosen missing ranges.
45
Figure 3.15: RMSE over 100 Monte Carlo runs for an average noise power from 80 to 40 decibels with10 randomly chosen missing ranges.
Figure 3.16: RMSE over 100 Monte Carlo runs for an average noise power from 80 to 40 decibels with20 randomly chosen missing ranges.
46
Table 3.1: Ratio of solved instances
Missing distances 5 10 20Av. noise level (dB) 80 70 60 50 40 80 70 60 50 40 80 70 60 50 40
EDM-SR 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1EDM-R .98 .96 .98 1 1 .88 .85 .90 .90 .98 .49 .56 .52 .58 .65
EDM-RLB .99 1 1 1 1 .97 .95 .99 1 .99 .34 .39 .50 .66 .66
Figure 3.17: Pairwise distances mapping for a 5 node constellation. Unknown distances: 3 − 4, 2 − 5,3− 5, 4− 5.
3.4 Pre-distance Matrix Rank Minimization
The performance of the convex optimization techniques analyzed in Section 3.3 was found to be quite
satisfactory even when some distances are missing. However, since the rank constraint was dropped
in the nearest EDMCP, the solver is free to search for optimal solutions in higher dimensions, which
for some network geometries results in poorer localization performances. Figure 3.18 exemplifies this
effect for a 5 node network using the EDM-SR formulation. The set of observed distances in Figure
3.18a is complete leading to a low rank solution for the EDMCP and an accurate estimation of the
nodes positions. For the same constellation and set of observed distances, the ranges between nodes
3 − 4, 2 − 5, 3 − 5 and 4 − 5 were considered to be unknown with the results shown in Figure 3.18b.
The optimal solution for the EDMCP lies in a higher dimension and the estimates are now inaccurate
(RMSE = 0.0233). This configuration leaves node 5 with one connection only (see Figure 3.17) and,
as expected, the rigidity matrix for this constellation does not have maximum rank (6/7), which implies
that the constellation is not globally rigid and therefore not uniquely realizable in R2. In fact, node 5
can assume any position in the circumference centered in 1 with radius d15 without violating the known
pairwise distances. Even so, it is fair to suspect that better position estimates can be achieved by
reducing the rank of the solution towards the embedding dimension of the problem.
To evaluate the fit of an EDMCP solution to an embedding dimension r, we will analyze the magnitude
of the (r + 1)th, . . . , (N)
th largest singular values returned by a SVD of JEJ. Recall that from equation
(2.30), rank (X) = rank (JEJ), and since JEJ is a symmetric matrix, its singular values (σi) are given
by the magnitude of the eigenvalues (λi), σi = |λi|. The SVD seems a more appropriate choice than the
47
(a) Complete set of distances. RMSE = 7.1686e− 009.
(b) Missing distances: [2−5], [3−5] and [4−5]. RMSE =0.4962.
Figure 3.18: Constellation reconstruction using EDM-SR and 80 decibels noise level.
eigenvalue decomposition since it can be numerically computed with good robustness and in this case,
the interest lies uniquely on the magnitude of the values. As an example, the scenario of Figure 3.18b
led to the following set of singular values (by decreasing order),
σ (JEJ) = [1.6119 0.5232 0.0387 0.0000 0.0000] (3.6)
where the third value, 0.0387, indicates that this solution has an embedding dimension greater than
r = 2. As expected, the optimal EDM in the problem of Figure 3.18a leads to an SVD of JEJ with only
two significant values,
σ (JEJ) = [1.0706 0.5597 0.0000 0.0000 0.0000] (3.7)
It becomes clear that one reasonable approach is trying to reduce the embedding dimension of the
problem before solving a relaxed SDP instance. The proposed method, named Pre-distance Matrix
Rank Minimization (PMRM), proposes a pre-distance matrix completion that minimizes the rank of the
corresponding constellation.
48
Consider the Symmetric interrogation scheme introduced in Section 3.2.2, recall that the pairwise
distances (vector d) and the time intervals (vector ∆) share a linear relation through a matrix A such
that Ad = ∆. When some of the entries in the pre-distance matrix are missing so that matrix A has a
nontrivial null-space, the solution of the non-homogeneous system of equations Ad = ∆ is given by
d = d0 + c1v1 + . . .+ csvs (3.8)
where c1, . . . , cs ∈ R and {v1, . . . , vs} is the basis for the null space of A with dimension s. The idea
is now to find the optimal set c1, . . . , cs such that JDJ has minimum rank, where D is the pre-distance
matrix built from the squared ranges (d2ij). Although equation (2.30) holds only for an EDM, it is expected
that by minimizing the embedding dimension of the pre-distance matrix D as an initialization step, the
convex optimization problem that follows also achieves a low rank solution, leading to more accurate
position estimates. To achieve this, it is suggested to minimize the (r + 1)th singular value and neglect
the remaining ones by expecting that close to the optimal point their magnitudes are also negligible 1.
For simplicity, let c be the vector containing c1, . . . , cs, v the basis for the null space of A, and f (c) =
σr+1 (JDJ) where σr+1 (·) is the (r + 1)th largest singular value. The problem can be formulated as,
minimizec
f (c)
subject to d = d0 + vc
Dii = 0 Dij = d2ij
(3.9)
Figure 3.19 exemplifies the shape of f (c) for a 5 node planar network (r = 2). In this case, 5 specifically
chosen ∆ entries are missing (4 pairwise distances affected) so that the null space of A has two dimen-
sions (c = [c1, c2]). Besides the presence of a global minimum at c = [−1.0713,−0.8198] it is also evident
a local minimum in the vicinity. The existence of local minima will be the main drawback in this approach
since it will be highly sensitive to the starting point for the numerical search. The strategy to overcome
this problem is addressed in Section 3.4.3. The main idea is to try different initial points, based on a
randomized procedure, until a sufficiently low residue is attained by the numerical search.
The strategy is summarized in Algorithm 1. Steps 2 and 8 do not specify which SDP formulation
should be used. Any of the three methods described in Section 2.2.6 remains as a valid choice. The
simulations were carried out using EDM-SR since the performance evaluation of Section 3.3 indicates
that it is the most robust although it may not be the most accurate. The extraction of vector c from the
optimal EDM from step 2 follows a simple least-squares procedure formulated as,
minimizec
||d− dEDM ||2
subject to d = d0 + vc
(3.10)
where dEDM is the vector of pairwise distances from the optimal EDM at step 2. For the minimization
of function f (c) (step 7) we can no longer rely on the convex optimization techniques described before
1 see Section 2.2.6 for alternative rank minimization heuristics (Flatten Graph and Nuclear norm). As detailed in Section 4.1, thoseheuristics were tested but the results were not reported here, although remain as appealing research directions.
49
Figure 3.19: Example of a plot of f (c) for a 5 node planar network with missing distances between nodes3 − 4, 2 − 5, 3 − 5 and 4 − 5. Null-space of A is two-dimensional. Average noise level of 80 decibels.Minimum at c = [−1.0713,−0.8198].
since f (c) is clearly non-convex. Instead, two alternative numerical optimization algorithms will be
implemented. A zero-order procedure, Powell’s Method, and a Quasi-Newton method, BFGS.
Algorithm 1 PMRM algorithm.
1: procedure DCOMPLETION2: Solve a convex optimization problem with the available observed ranges.3: From the optimal EDM, extract the corresponding c vector.4: c0 ← c5: Find the vector c that minimizes f (c) with c0 as the initialization point.6: csol ← c7: Complete the pre-distance matrix using csol.8: Solve a convex optimization problem with the optimal pre-distance matrix from step 7.9: Reconstruct the framework based on the optimal EDM from step 8.
10: end procedure
3.4.1 Powell’s Method
Consider a problem in a n-dimensional design space. The objective is to minimize function f (x), where
the components of x are n design variables. To do so, Powell’s Method applies the basic strategy
common to all Conjugate Gradient Methods, as described in Algorithm 2:
50
Algorithm 2 Conjugate Gradient Methods basic algorithm. Adapted from [41].
1: procedure CONJUGATE
2: choose a starting point x0
3: for i = 1, 2, 3, . . . do
4: choose a vector vi
5: minimize f (x) in the direction of vi
6: let the minimum point be xi
7: if |xi − xi−1| < ε then
8: exit loop
9: end if
10: end for
11: end procedure
The choice of the vectors vi has given rise to several different numerical optimization methods, each
one of them suited for a particular class of problems. In the case of Conjugate Gradient Methods the
directions vi are mutually conjugate (non-interfering). The implication is that, once we have minimized
f (x) in one direction, we can move along the next direction without ruining the previous minimization.
In order to keep the complexity of the problem low, the first choice relied on a zero-order procedure,
this is, a method that only requires the evaluation of the function at each point without the need for com-
puting first derivatives. According to [41], the overwhelming choice on zero-order methods is Powell’s
Method. The standard procedure is described in Algorithm 3.
Algorithm 3 Powell’s Method basic algorithm. Adapted from [41].
1: procedure POWELL2: choose a starting point x03: choose the starting vectors vi, i = 1, 1, . . . , n . the usual choice are the unit vectors4: for i = 1, 2, . . . , n do5: Minimize f (x) along the line through xi−1 in the direction of vi6: Let the minimum point be xi7: end for8: vn+1 ← xn − x09: Minimize f (x) along the line through x0 in the direction of vn+1
10: Let the minimum point be xn+1
11: for i = 1, 2, . . . , n do12: vi ← vi+1 . v1 is discarded, the other vectors are reused13: end for14: if ||xn+1 − x0||2 > ε then15: return to 416: end if17: end procedure
Powell has demonstrated that the vectors vi are mutually conjugate, which implies that for a quadratic
problem with n variables, convergence will be achieved in less than or equal to n cycles [41]. Since
design optimization problems are seldom described by a quadratic polynomial, the solution should not
be expected to converge quadratically.
Alternatively, when the first derivatives are available, the Quasi-Newton methods provide dramatic
51
improvements with a moderate increase in complexity [42]. Unlike Newton Methods, second derivatives
are not required as the model of the objective function is constructed by measuring the changes in
gradients. One of the most popular Quasi-Newton procedures is BFGS, described later in Section 3.4.2.
Finally, the minimization along a line can be accomplished with any one-dimensional optimization
algorithm. One of the most commonly used is the Golden Section Search which is a technique for
finding the extremum of a strictly unimodal function by successively narrowing the range of values inside
which the extremum is known to exist. Again, function f (c) exhibits local minima that can lead to a line
search over a function which is not unimodal. Even so, this method is expected to perform well in the
cases where the starting point is sufficiently close to the solution.
Golden Section Search
The Golden Section Search conducts the search by evaluating the function at a triplet of points (a, b, c),
where a and c are the limits of the interval where the solution is known to exist and b is a middle point
placed at an optimal distance from a and c. The details of this method can be found in [43]. The
main idea is to maintain the same proportion of spacing throughout the procedure in such a way that
in each step, the two sub-intervals where the extremum might be are equal. By doing this, not only
it is guaranteed that the interval width shrinks by the same constant proportion in each step, but also
avoids slowing down the rate of convergence by setting uneven sub-intervals. The optimal ratio between
the distances (c − b) and (b − a) can be easily derived and is equal to the well known Golden Ratio
ϕ = 1+√5
2 ≈ 1.618 (or its inverse if (b − a) > (c − b)). The algorithm terminates the search when the
interval between the extreme points of the triplet is narrow enough for the problem at hand. Algorithm 4
summarizes the procedure for a minimum search.
3.4.2 BFGS
BFGS is one of the most successful Quasi-Newton methods. It starts by adjusting a quadratic model of
the objective function at each iterate k:
mk (p) = fk +∇fTk p+1
2pTBkp (3.11)
where Bk is an n × n symmetric positive definite matrix that approximates the Hessian matrix without
the need for explicit second derivatives.
The search direction is given by the minimizer of (3.11),
pk = −B−1k ∇fk (3.12)
which will be used in a one-dimensional search routine. More details on the derivation of BFGS, in
particular, of matrix B, can be found in [42].
It is advised to impose the Wolfe or strong Wolfe conditions when performing the line search, as
it is the only way to guarantee the self-correcting properties of BFGS [42]. The Wolfe conditions force
52
Algorithm 4 Golden Section Search algorithm. Adapted from [43].
1: procedure GSS2: consider the triplet (a, b, c) with a < b < c where b = c+ϕa
ϕ+1 or b = ϕ+1c+ϕa
3: while |c− a| > ε do4: if c− b > b− a then5: x = b+ (2− ϕ)(c− b)6: else7: x = b− (2− ϕ)(b− a)8: end if9: if f(x) < f(b) then
10: if c− b > b− a then11: (a, b, c)← (b, x, c)12: else13: (a, b, c)← (b, x, c)14: end if15: else16: if c− b > b− a then17: (a, b, c)← (a, b, x)18: else19: (a, b, c)← (x, b, c)20: end if21: end if22: end while23: return (c+ a)/224: end procedure
function f (x) to be evaluated at a point that verifies a sufficient decrease in f (x) with a sufficiently sharp
curvature. However, during the development of the PMRM, it was verified that respect for the Wolfe
conditions frequently would lead to an infinite loop where a solution could not be attained. Therefore,
a simple Golden Section Search was implemented as the line search procedure even at the risk of
compromising some of the properties of BFGS. The basic procedure for BFGS is described in Algorithm
5.
Gradient of the Objective Function
To fully specify the BFGS method an expression for the gradient of f (c), ∇f , must be derived. Recall
f (c) as f (c) = σr+1 (JDJ) where σr+1 (·) is the (r + 1)th singular value, J is the orthogonal projection
matrix on 1⊥ and D is the pre-distance matrix. Let for simplicity Z = JDJ and φ (Z) = λr+1 (Z), such
that f (c) = |φ (Z) |. The real absolute value function g (x) = |x| has a derivative for every x 6= 0, but is
not differentiable at x = 0. Its derivative for x 6= 0 is given by the step function,
d|x|dx
=
1 if x > 0
−1 if x < 0(3.13)
By the chain rule,
df (c) =
dφ (Z) if φ (Z) > 0
−dφ (Z) if φ (Z) < 0(3.14)
53
Algorithm 5 Adapted BFGS algorithm. Adapted from [43].
1: procedure BFGS2: Given a starting point x0, convergence tolerance parameters ε1 > 0 and ε2 > 0 and an inverse
Hessian approximation H0
3: k ← 04: while ||∇fk|| > ε do5: Compute search direction pk = −Hk∇fk6: xk+1 = xk + αkpk where αk is computed from a line search procedure7: Define sk = xk+1 − xk8: if sk < ε2 then break9: end if
10: Define yk = ∇fk+1 −∇fk11: Define ρk = 1
yTk sk
12: Hk+1 =(I− ρksky
Tk
)Hk
(I− ρkyks
Tk
)+ ρksks
Tk
13: k ← k + 114: end while15: end procedure
The mathematical results used throughout this section are almost exclusively based on [44]. The proofs
will be omitted here for the sake of brevity. For an introductory background see Appendix A. For a real
symmetric matrix Z,
dφ = dλr+1 (Z)
= uTr+1 (dZ)ur+1
=(uTr+1 ⊗ u
Tr+1
)vec (dZ)
= (ur+1 ⊗ ur+1)T
vec (dZ)
(3.15)
where (⊗) denotes the Kronecker product and ur+1 is the normalized eigenvector associated with λr+1
such that
Zur+1 = λr+1ur+1 uTr+1ur+1 = 1 (3.16)
Let D′ be the pointwise square root of D such that D = D′ � D′ and recall that the vector of pairwise
distances can be written as d = d0 + c1v1 + . . .+ csvs where c1, . . . , cs ∈ R and {v1, . . . , vs} is the basis
for the null space of A with dimension s. The generalization to the matrix form yields D′ = D′0 +s∑
i=1
ciD′i.
Hence,
dZ = 2J(D′ � dD′
)J
= 2J
(D′ �
s∑i=1
D′idci
)J
=
s∑i=1
(2J(D′ �D′i
)J)dci
(3.17)
54
Let Vi =(2J(D′ �D′i
)J), then,
vec (dZ) =
s∑i=1
vec (Vi) dci = V dc (3.18)
Therefore dφ = (ur+1 ⊗ ur+1)TV dc. From the first identification theorem in [44, pp. 199], the derivative
can be obtained from the differential such that,
Dφ = (ur+1 ⊗ ur+1)TV (3.19)
where in this case D stands for the derivative operator and should not be mistaken with the pre-distance
matrix. Finally, the gradient results from the derivative [44, pp. 99],
∇φ = DφT (3.20)
and the gradient of f (c) can be written as,
∇f (c) =
(
(ur+1 ⊗ ur+1)TV)T
if φ (Z) > 0
−(
(ur+1 ⊗ ur+1)TV)T
if φ (Z) < 0(3.21)
3.4.3 Dealing with Local Minima
Function f (c) suffers from the presence of local minimums in the vicinity of the global minimum point.
Therefore, the performance of search procedures for minimization of f (c) is highly dependent on the
initial point. To deal with this problem, a randomized initialization procedure was developed. The main
idea is to randomize the initialization point for the numerical search until the optimal value for f (c)
returned by the search is low enough. An extended version of Algorithm 1 is presented as Algorithm 6.
The randomized procedure considers the vicinity of the first estimate for c based on the EDMCP solved
with the available set of measures, denominated as coriginal (step 3). Simulations have shown that a
range of[0, 2coriginal
]produces good results although there were not conducted systematic studies to
support this choice.
3.4.4 Performance Evaluation of the Pre-distance Matrix Rank Minimization
Numerical simulations were conducted in order to evaluate the improvements in node positioning brought
by PMRM when compared against a simple EDM-SR technique. All the simulations considered a 5 node
network on a unit square area. To fully reproduce the results, the search parameters used during the
simulations are detailed next.
55
Algorithm 6 Pre-distance matrix completion algorithm with randomized initialization.
1: procedure DCOMPLETION2: Solve a convex optimization problem with the available observed ranges.3: From the optimal EDM extract the corresponding c vector.4: coriginal ← c5: c0 ← coriginal6: k ← 17: while k < Nc do8: Find the vector c that minimizes f (c) with c0 as the initialization point.9: if f (c) < ε then break
10: else11: c0 ← crand where crand is a vector given by random values in the range
[0, 2coriginal
]12: end if13: k ← k + 114: end while15: csol ← c16: Complete the pre-distance matrix using csol.17: Solve a convex optimization problem with the optimal pre-distance matrix from step 16.18: Reconstruct the framework based on the optimal EDM from step 17.19: end procedure
Powell Parameters
Powell’s Method was directly applied as suggested in [41] (see Algorithm 3) with an additional stopping
criterion that limits the number of iterations to 100s where s is the number of variables to optimize.
This avoids infinite loops and keeps the processing time under control. The main stopping criterion,
||xn+1 − x0||2F < ε, uses ε = 10−4.
BFGS Parameters
The implementation of BFGS adopted some modifications with respect to the algorithm suggested in
[42]. The main stopping criterion, ||∇fk|| < ε, uses ε = 10−4. However, due to the non-smooth shape
of f (c) near the true minimum (see Figure 3.19) this condition is seldom achieved. Therefore, two
additional stopping criteria have been set: sk < ε2 with ε2 = 10−4, this is, when the optimal point, xk,
changes less than ε2 between iterations, and once again, a limit on the number of iterations to 100s due
to the same reasons as described before. It is relevant to remember that the Wolfe conditions were not
verified and a simple Golden Section Search was used as the line search procedure.
Golden Section Search Parameters
The stopping criterion for the Golden Section Search was |c− a| < ε with ε = 10−3.
Scenario Testing
Three different scenarios were considered, each with two noise levels (60 and 80 decibels), totaling six
different test cases. As usual, the results for each test case consider 100 Monte Carlo runs with different
randomly created constellations for each run. A description and results for each scenario are given next.
56
Scenario 1: Complete Set
The first scenario considers a complete set of observed distances. In this case, PMRM is dispensable
since matrix D is already complete. The results are presented in Tables 3.2 and 3.3 for 60 and 80
decibels, respectively. The evaluation comprises the RMSE, the ratio of solved instances, the Mean
Processing Time (MPT) and the ratio between the RMSE of both methods. As expected the RMSE is
similar for the three methods while the MPT approximately doubles for the Powell and BFGS methods
with respect to EDM-SR. This is easily understood since the PMRM requires two convex optimization
problems to be solved.
Table 3.2: Scenario 1: Complete set. Noise level: 60 decibels.
RMSE Solved MPT (s) RMSE ratioEDM-SR 1.89× 10−3 1.00 0.33 -Powell 1.89× 10−3 1.00 0.68 1.00
BFGS 1.89× 10−3 1.00 0.67 1.00
Table 3.3: Scenario 1: Complete set. Noise level: 80 decibels.
RMSE Solved MPT (s) RMSE ratioEDM-SR 5.33× 10−4 1.00 0.33 -Powell 5.33× 10−4 1.00 0.68 1.00
BFGS 5.33× 10−4 1.00 0.68 1.00
Scenario 2: One Dimensional Null Space
The second scenario considers the elimination of 4 time intervals (affecting distances 3 − 4, 3 − 5 and
4−5) so that the null space of A is one dimensional. Figure 3.20 illustrates the known pairwise distances
for a random geometry and Figure 3.21 exemplifies the shape of function f (c) for the one dimensional
case. The results are presented in Tables 3.4 and 3.5 for 60 and 80 decibels, respectively. It is clear that
the accuracy improves when PMRM is applied, both with Powell’s and BFGS search. The difference
between RMSE values for the two procedures is not significant and therefore does not allow to draw a
definite conclusion regarding which search method is best suited for this application.
Table 3.4: Scenario 2: One dimensional null space. Noise level: 60 decibels.
RMSE Solved MPT (s) RMSE ratioEDM-SR 0.269 1.00 0.32 -Powell 0.043 1.00 0.89 0.16BFGS 0.077 1.00 1.24 0.29
Scenario 3: Two Dimensional Null Space
The third scenario considers the elimination of 5 time intervals (affecting distances 3− 4, 2− 5, 3− 5 and
4 − 5) in such a way that the null space of A is two dimensional. This scenario was already illustrated
in Figure 3.17 with an example of the shape of f (c) on the two dimensional case in Figure 3.19. The
57
Figure 3.20: Pairwise distances mapping for a 5 node constellation. Unknown distances: 3 − 4, 3 − 5,4− 5.
Figure 3.21: Example of a plot of f (c) for a 5 node planar network with missing distances betweennodes 3 − 4, 3 − 5 and 4 − 5. Null-space of A is one-dimensional. Average noise level of 80 decibels.Minimum at c1 = [1.2440].
58
Table 3.5: Scenario 2: One dimensional null space. Noise level: 80 decibels.
RMSE Solved MPT (s) RMSE ratioEDM-SR 0.266 1.00 0.36 -Powell 0.072 1.00 0.99 0.27BFGS 0.102 1.00 1.38 0.38
results are presented in Tables 3.6 and 3.7 for 60 and 80 decibels respectively. Again, the completion of
the pre-distances matrix before the convex optimization problem allowed a reduction on the positioning
error with both search methods at the expense of a considerable increase in the MPT. In this case,
BFGS seems to perform slightly better, although the difference is not significant enough to justify any
conclusions. More test cases would be needed to definitely prefer one method over the other.
Finally, one remark should be made about the error metric. The use of RMSE is very common and
it makes an excellent general purpose error metric for numerical predictions. But we must keep in mind
that compared to the similar Mean Absolute Error, the RMSE amplifies and punishes large errors. In fact,
a closer look at the results shows that the PMRM is extremely accurate for a great number of runs. The
value of the RMSE is almost exclusively due to a number of particular runs where the search procedure
was not capable of finding the global minimum of f (c).
Table 3.6: Scenario 3: Two dimensional null space. Noise level: 60 decibels.
RMSE Solved MPT (s) RMSE ratioEDM-SR 0.310 1.00 0.36 -Powell 0.202 1.00 0.93 0.65BFGS 0.143 1.00 1.22 0.46
Table 3.7: Scenario 3: Two dimensional null space. Noise level: 80 decibels.
RMSE Solved MPT (s) RMSE ratioEDM-SR 0.302 1.00 0.34 -Powell 0.176 1.00 0.87 0.58BFGS 0.156 1.00 1.06 0.52
59
60
Chapter 4
Conclusions & Future Work
This thesis addresses the problem of localization of multiple vehicles based on distance measures.
Although GNSS-based solutions effectively answer this need for a variety of scenarios, their usage is
not universal. Environments where the GNSS signals are not available, low-cost networks that do not
allow expensive equipment in every single node or even formations of vehicles that demand a high
accuracy relative positioning system are examples of scenarios that exclude GNSS localization as an
option. A viable alternative is positioning based on pairwise distance measures that locates the nodes
relative to each other when sufficient pairwise distances are collected. To achieve global localization,
however, it is required that some of the vehicles in the constellation have known positions (anchors) so
that the remaining degrees of freedom are resolved.
The framework reconstruction can assume very simple algorithms such as the LS or the hyperbola
methods [3] that provide straightforward and fast solutions, or on the other hand, complex and more
accurate algorithms such as the optimization-based techniques. This thesis is mainly focused on the
latter where the accuracy of the position estimates is the most important driving force.
In real applications the ranging process is not ideal and the sets of range measurements are often
noisy and incomplete. This situation falls into the well-known category of problems, EDMCP, where
despite the enormous amount of work published in this area and the variety of algorithms proposed
(see the state-of-the-art in Section 2.3) it is still an open problem. In fact, the completion of an EDM
subject to a rank constraint, in order to guarantee that the solutions lies in the embedding dimension of
the problem, was proven to be NP-hard [33]. Several alternative methods were proposed to overcome
this challenge. One of the most popular strategies are the SDP relaxations, that neglect some of the
constraints such that the problem is capable of being handled by a generic convex optimization solver.
The SDP formulations were the main approach to the node localization throughout this work.
The main contributions of this thesis can be summarized in four different topics, each one of them
trying to answer to a different challenge. First, an alternative constellation reconstruction procedure is
proposed for the MORPH’12 experimental results. In [7], GPS data is compared with the position esti-
mates obtained through a framework reconstruction based on distance measurements obtained under
the MORPH’12 trials. Using the same data set, a different approach for the anchors position matching
61
has been proposed, achieving a better agreement between the estimated trajectory and the GPS data.
Second, it has also been proposed a different ranging scheme (or interrogation scheme) that super-
sedes the Master/Slaves used on the MORPH’12 experiments. The greatest advantage is the reduction
in the amount of data transmitted, avoiding the bandwidth cluttering with ranging related transmissions.
Third, three different convex formulations have been tested under noisy and incomplete measurements
conditions. The simulations covered several noise powers and incomplete data scenarios. The for-
mulations under analysis were the classical EDM-SR (see [4] for example), the EDM-R (proposed in
[14]) and the novel EDM-RLB, that showed higher accuracy than the remaining ones. The conventional
EDM-SR however, appeared to be the most robust with a higher ratio of instances solved. Finally, the
most relevant achievement in this dissertation is the development of the pre-distance matrix completion
algorithm, the PMRM. Prior to the SDP problem, this algorithm searches for the pre-distance matrix
completion that minimizes the embedding dimension of the problem. This way, we prevent the solver
to search for a solution in higher dimensions. In practice, the search translates into a minimization of
a non-convex function where all the collected time difference of arrivals contribute to narrow the search
space. The minimization followed two alternative methods, the Powell’s method and the BFGS both
using a randomized initialization procedure. The performance of each one of them is analyzed in the
end of Chapter 3. The simulations considered three scenarios: complete pre-distance matrix, a one
dimensional search and a two-dimensional search. The results indicate that the PMRM achieves a sig-
nificant improvement in the positioning accuracy, with reduction in the RMSE up to 84% and 71% for the
Powell’s method and the BFGS respectively, at the expense of a higher processing time (between two
and three times the processing time for the EDM-SR).
4.1 Future Work
The figures presented for the PMRM algorithm are encouraging but are far from being closed. Not only
the procedure is not completely optimized, through a fine tuning of the search parameters for example,
but also the simulation results are not complete since the spectrum of the scenarios analyzed must be
wider. Besides, Powell and BFGS methods, as well as the randomized initialization procedure might
not be the most appropriate techniques, and therefore other options should also be evaluated. It should
also be noted that the inclusion of a pre-distance completion algorithm is only required because it does
not exist yet a convex formulation that effectively achieves low-rank solutions on a simple and elegant
mathematical framework. That accomplishment would represent a tremendous step in the range-based
localization research and a valuable tool in the development of all kinds of autonomous vehicles. The
research for this thesis highlighted two promising directions on the path to accomplish that goal.
The first is the development of an effective regularization term forcing a low-rank solution. This kind
of approach was described in Section 2.2.6 reporting the work in [31]. As stated before, the authors
add a regularization term that forces the flattening of the graph by pushing the nodes away from each
other. In order to be effective, the weight of the regularization parameter must lie in an admissible range
of values. This reasoning looks appealing and some experiments were actually conducted in order to
62
understand the potential of this formulation. However, during the course of this thesis, it was impossible
to achieve encouraging results since the optimum value for the weight of the regularization parameter
was shown to be highly dependent of the constellation itself. Plus, the heuristic suggested in [31] has
not produced any satisfactory results either. Alternatively to the Flatten Graph formulation, the same
kind of reasoning should be tried for a regularization term based on the nuclear norm. Formally,
minimizeE
‖W � (E −D) ‖+ λ1TE1
subject to E ∈ EDM(4.1)
As detailed in Section 2.2.6, the reasoning behind this approach is the opposite of the strategy applied
in the Flatten Graph heuristic and may not be intuitive. However, has received much attention lately with
some proofs of success (see [22]) which encourages a closer look at this technique.
The second appealing research direction is a natural development of the EDM-RLB proposed here.
Using a similar formulation, the idea is now to include the time intervals collected during the ranging
process in the convex formulation. A first draft of the this formulation can be written as,
minimizeE,T
∑i,j
(Tij − Lij (d)
)2subject to Tij ≤
√Eij
Tij ≥Eij√Emax
E ∈ EDM
Ad = ∆
(4.2)
This way, the relation between the time intervals and pairwise distances can now be used to narrow the
search space. This method has not been tested yet and at this point it only constitutes a proposal for
future work.
63
64
Bibliography
[1] Morph, January 2014. URL http://morph-project.eu/.
[2] Hwee-Pink Tan, Roee Diamant, Winston K.G. Seah, and Marc Waldmeyer. A survey of techniques
and challenges in underwater localization. Ocean Engineering, 38(14,15):1663 – 1676, 2011.
[3] Sen Wang and Huosheng Hu. Wireless sensor networks for underwater localization: A survey.
Technical report, School of Computer Science and Electronic Engineering, University of Essex,
2012.
[4] Jon Dattorro. Convex Optimization & Euclidean Distance Geometry. Meboo Publishing USA, 2011.
[5] Zhisu Zhu, Anthony Man-Cho So, and Yinyu Ye. Universal rigidity and edge sparsification for sensor
network localization. SIAM J. on Optimization, 20(6):3059–3081, October 2010.
[6] Vaidyanathan Ramadurai and Mihail L. Sichitiu. Localization in wireless sensor networks: A prob-
abilistic approach. In Weihua Zhuang, Chi-Hsiang Yeh, Olaf Droegehorn, C.-T. Toh, and Hamid R.
Arabnia, editors, International Conference on Wireless Networks, pages 275–281. CSREA Press,
2003.
[7] João Gomes, Ehsan Zamanizadeh, José M. Bioucas-Dias, João Alves, and Thomas C. Furfaro.
Building location awareness into acoustic communication links and networks through channel delay
estimation. In Proceedings of the Seventh ACM International Conference on Underwater Networks
and Systems, WUWNet ’12, pages 24:1–24:8. ACM, 2012.
[8] Liam Paull, Sajad Saeedi, Mae Seto, and Howard Li. AUV Navigation and Localization: A Review.
IEEE Journal of Oceanic Engineering, 39(1):131–149, January 2014.
[9] Long Cheng, Chengdong Wu, Yunzhou Zhang, Hao Wu, Mengxin Li, and Carsten Maple. A Survey
of Localization in Wireless Sensor Network. International Journal of Distributed Sensor Networks,
2012:1–12, 2012.
[10] Vijay Chandrasekhar, Winston KG Seah, Yoo Sang Choo, and How Voon Ee. Localization in un-
derwater sensor networks: survey and challenges. In Proceedings of the 1st ACM international
workshop on Underwater networks, WUWNet ’06, pages 33–40. ACM, 2006.
[11] Nathan Krislock and Henry Wolkowicz. Explicit sensor network localization using semidefinite rep-
resentations and facial reductions. SIAM J. on Optimization, 20(5):2679–2708, July 2010.
65
[12] Yiyin Wang, Xiaoli Ma, and G. Leus. Robust time-based localization for asynchronous networks.
Trans. Sig. Proc., 59(9):4397–4410, September 2011.
[13] Cvx research, January 2014. URL http://cvxr.com/.
[14] P. Oguz-Ekim, J. P. Gomes, J. Xavier, and P. Oliveira. Robust localization of nodes and time-
recursive tracking in sensor networks using noisy range measurements. Trans. Sig. Proc., 59(8):
3930–3942, August 2011.
[15] Vasco Ludovico, João Gomes, João Alves, and Thomas Furfaro. Joint Localization of Underwater
Vehicle Formations Based on Range Difference Measurements. Submitted to Proceedings of IEEE
Underwater Communications Conference and Workshop (UComms), 2014.
[16] João Gomes, Vasco Ludovico, João Alves, and Thomas Furfaro. Collaborative Localization based
on Differential Range Measurements. In preparation to be submitted to IEEE Journal of Oceanic
Engineering, 2014.
[17] Karl Johan Åström and Björn Wittenmark. Computer-controlled systems: theory and design; 3rd
ed. Prentice-Hall information and system sciences series. Prentice-Hall, 1997.
[18] Amitangshu Pal. Localization Algorithms in Wireless Sensor Networks : Current Approaches and
Future Challenges. Journal of Network Protocols and Algorithms, 2(1):45–74, 2010.
[19] Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
[20] Pinar Oguz-ekim, João Gomes, João Xavier, and Paulo Oliveira. ML-BASED SENSOR NETWORK
LOCALIZATION AND TRACKING : BATCH AND TIME-RECURSIVE APPROACHES. In Proceed-
ings of the 17th European Signal Processing Conference, number Eusipco in EUSIPCO’09, pages
80–84, 2009.
[21] Yichuan Ding, Nathan Krislock, Jiawei Qian, and Henry Wolkowicz. Sensor network localization,
euclidean distance matrix completions, and graph realization. In Proceedings of the first ACM
international workshop on Mobile entity localization and tracking in GPS-less environments, MELT
’08, pages 129–134. ACM, 2008.
[22] Nathan Krislock and Henry Wolkowicz. Euclidean distance matrices and applications. In Handbook
of Semidefinite, Cone and Polynomial Optimization, 2010.
[23] Anthony Man-Cho So and Yinyu Ye. Theory of semidefinite programming for sensor network lo-
calization. In Proceedings of the sixteenth annual ACM-SIAM symposium on Discrete algorithms,
SODA ’05, pages 405–414. Society for Industrial and Applied Mathematics, 2005.
[24] W. Kabsch. A discussion of the solution for the best rotation to relate two sets of vectors. Acta
Crystallographica Section A, 34(5):827–828, September 1978.
[25] A.Y. Alfakih. On the dual rigidity matrix. Linear Algebra and its Applications, 428(4):962 – 972,
2008.
66
[26] Brian D. O. Anderson, Iman Shames, Guoqiang Mao, and Baris Fidan. Formal theory of noisy
sensor network localization. SIAM J. Discret. Math., 24(2):684–698, June 2010.
[27] Bill Jackson and Tibor Jordán. On the rank function of the 3-dimensional rigidity matroid. Int. J.
Comput. Geometry Appl., 16(5-6):415–430, 2006.
[28] Jeff Erickson. Algorithms. Departement of Computer Science, University of Illinois, 2004.
[29] Zhisu Zhu, Anthony Man-Cho So, and Yinyu Ye. Universal rigidity: Towards accurate and efficient
localization of wireless networks. In Proceedings of the 29th Conference on Information Communi-
cations, INFOCOM’10, pages 2312–2320. IEEE Press, 2010.
[30] Abdo Y. Alfakih, Amir Khandani, and Henry Wolkowicz. Solving euclidean distance matrix com-
pletion problems via semidefinite programming. Comput. Optim. Appl., 12(1-3):13–30, January
1999.
[31] Pratik Biswas, Tzu-Chen Liang, Kim-Chuan Toh, Yinyu Ye, and Ta-Chung Wang. Semidefinite
programming approaches for sensor network localization with noisy distance measurements. IEEE
T. Automation Science and Engineering, 3(4):360–371, 2006.
[32] Leo Liberti, Carlile Lavor, Nelson Maculan, and Antonio Mucherino. Euclidean distance geometry
and applications. SIAM Review, 56(1):3–69, 2014.
[33] James Saxe. Embeddability of weighted graphs in k-space is strongly nphard, 1980.
[34] K.Q. Weinberger, F. Sha, and L. K. Saul. Learning a kernel matrix for nonlinear dimensionality
reduction. In Proceedings of the Twenty First International Conference on Machine Learning (ICML-
04), pages 839–846, 2004.
[35] A. Beck, P. Stoica, and Jian Li. Exact and approximate solutions of source localization problems.
Trans. Sig. Proc., 56(5):1770–1778, May 2008.
[36] Bin Liu, Hongyang Chen, Ziguo Zhong, and H. Vincent Poor. Asymmetrical round trip based
synchronization-free localization in large-scale underwater sensor networks. Trans. Wireless.
Comm., 9(11):3532–3542, November 2010.
[37] K. Yu, Y. J. Guo, and M. Hedley. Toa-based distributed localisation with unknown internal delays
and clock frequency offsets in wireless sensor networks. Signal Processing, IET, 3(2):106–118,
2009. ISSN 1751-9675.
[38] Zhong Zhou, Jun-Hong Cui, and Shengli Zhou. Efficient localization for large-scale underwater
sensor networks. Ad Hoc Networks, 8(3):267 – 279, 2010.
[39] Zhong Zhou, Zheng Peng, Jun-Hong Cui, Zhijie Shi, and Amvrossios Bagtzoglou. Scalable lo-
calization with mobility prediction for underwater sensor networks. IEEE Transactions on Mobile
Computing, 10(3):335–348, March 2011.
67
[40] P.O. Ekim, J. Gomes, J. Xavier, and P. Oliveira. Experimental evaluation of simultaneous 3d local-
ization of sensor nodes and tracking moving targets. In Signal Processing Conference (EUSIPCO),
2012 Proceedings of the 20th European, pages 180–184, 2012.
[41] P Venkataraman. Applied Optimization with MATLAB Programming. Wiley, 2002.
[42] J. Nocedal and S. J. Wright. Numerical Optimization. Springer, 2nd edition, 2006.
[43] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes
3rd Edition: The Art of Scientific Computing. Cambridge University Press, 3 edition, 2007.
[44] J.R. Magnus and H. Neudecker. Matrix Differential Calculus with Applications in Statistics and
Econometrics. Wiley Series in Probability and Statistics: Texts and References Section. Wiley,
1999.
[45] Abdo Y. Alfakih. On the uniqueness of euclidean distance matrix completions: the case of points in
general position. Linear Algebra and its Applications, 397(0):265 – 277, 2005.
[46] I.F. Akyildiz, W. Su, Y. Sankarasubramaniam, and E. Cayirci. Wireless sensor networks: a survey.
Computer Networks, 38(4):393 – 422, 2002.
[47] G. Laman. On graphs and rigidity of plane skeletal structures. Journal of Engineering Mathematics,
4(4):331–340, October 1970.
[48] Isy linkoping university, January 2014. URL http://www.isy.liu.se/en/.
[49] K. B. Petersen and M. S. Pedersen. The matrix cookbook, nov 2012. Version 20121115.
68
Appendix A
Matrix Differential Calculus
A.1 Introduction
This Appendix introduces the most relevant mathematical results used in Section 3.4.2. The notation
follows [44] except for the transpose operator which is here defined as (.)T . Interested readers can also
refer to [44] for detailed proofs of the results.
A.2 The chain rule
Let S be a subset of Rn, and assume that f : S → Rm is differentiable at an interior point c of S. Let T
be a subset of Rm such that f (x) ∈ T for all x ∈ S, and assume that g : T → Rp is differentiable at an
interior point b = f (c) of T . Then the composite function h : S → Rp defined by
h (x) = g (f (x)) (A.1)
is differentiable at c, and
Dh (c) = (Dg (b)) (Df (c)) (A.2)
A.3 The first identification theorem
Let f : S → Rm be a function defined on a set S in Rn, and differentiable at an interior point c of S. Let
u be a real n× 1 vector. Then,
df (c;u) = (Df (c))u (A.3)
where df (c;u) is called the (first) differential of f at c and (Df (c)) is an m × n matrix whose elements
Djfi (c) are the partial derivatives of f evaluated at c, known as the Jacobian matrix of f at c. The
transpose of the m × n Jacobian matrix is an n × m matrix called the gradient of f at c, denoted by
∇f (c). Thus,
∇f (c) = (Df (c))T (A.4)
69
The first identification theorem has great practical value due to the fact that if f is differentiable at c and
we have found a differential df at c, then the value of the partials at c can be immediately determined. In
particular, for a scalar function φ (x) of a vector variable x. Let the differential be written as,
dφ = aT dx (A.5)
then the derivative comes as
Dφ (x) = aT (A.6)
A.4 The differential of eigenvalues: symmetric case
Let X0 be a real symmetric n × n matrix, and let u0 be a (normalized) eigenvector associated with an
eigenvalue λ0 of X0. Since X0 is real and symmetric, its eigenvalues are real and its eigenvectors can
always be taken to be real. This assumption fits our problem and simplifies the derivation. A real-valued
function λ and a vector function u are defined for all X in some neighborhood N (X0) ⊂ Rn×n of X0,
such that
λ (X0) = λ0, u (X0) = u0 (A.7)
and
Xu = λu, uTu = 1 (X ∈ N (X0)) (A.8)
The function λ is∞ times differentiable on N (X0) as long as λ0 is simple, and the differentials at X0 are
dλ = uT0 (dX)u0 (A.9)
A.5 Additional fundamental rules of differential calculus
A.5.1 Differential of the Hadamard product
Let U and V be two matrix functions. Then, the differential of the Hadamard product between U and V
can be written as,
d (U � V ) = (dU)� V + U � dV (A.10)
A.5.2 Differential of the vectorization of a matrix function
Let U be a matrix function. Then, the differential of the vectorization of U can be written as,
d vecU = vec dU (A.11)
70
