Range-based Coordinate Alignment for Cooperative Mobile Sensor Network Localization

TL;DR

This paper introduces a novel range-based coordinate alignment method for cooperative mobile sensor network localization, transforming a non-convex problem into a convex one and proposing algorithms validated through simulations.

Contribution

It presents a new convex reformulation of the coordinate alignment problem and develops iterative and recursive algorithms for localizing target nodes in sensor networks.

Findings

Algorithms outperform existing methods in simulations

Effective in localizing multiple target nodes

Convex reformulation simplifies the non-convex optimization problem

Abstract

This paper studies a coordinate alignment problem for cooperative mobile sensor network localization with range-based measurements. The network consists of target nodes, each of which has only access position information in a local fixed coordinate frame, and anchor nodes with GPS position information. To localize target nodes, we aim to align their coordinate frames, which leads to a non-convex optimization problem over a rotation group $\text{SO}(3)$. Then, we reformulate it as an optimization problem with a convex objective function over spherical surfaces. We explicitly design both iterative and recursive algorithms for localizing a target node with an anchor node, and extend to the case with multiple target nodes. Finally, the advantages of our algorithms against the literature are validated via simulations.

Tables

Table 1. TABLE I: Relative errors of translation (%).

Algorithms	SDP	PPA	SDP+PPA	SDP+GD
SNR=20	7.85	5.87	5.87	5.99
SNR=30	6.20	4.48	4.48	5.12
SNR=80	0.09	0.06	0.06	0.06

Table 2. TABLE II: Distribution of target-anchor measurements over target nodes

# of target-anchor range measurements	$0$	$\le 5$	$\le 10$	$\le 15$	$\le 20$	$>20$
# of target nodes	56	23	17	10	3	1
percentages (%)	50.9	20.9	15.5	9.1	2.7	0.9

Equations

$$m_i(t)=\{d_{ij}(t),r_{ia}(t)|j\in {\mathcal{T}}_i(t),a\in {\mathcal{A}}_i(t)\}.$$

$$\text{SO}(3)=\{R\in {\mathbb{R}}^{3\times 3}|R{R}^{\prime }=I,\text{det}(R)=1\}$$

$$d_{ij}(t)$$

$$r_{ia}(t)$$

$$\begin{array}{rl}{f}_{ij}^{\mathcal{T}}(t,R,T)& =(d_{ij}(t)-\Vert R_i{p}_i^l(t)+T_i-R_j{p}_j^l(t)-T_j\Vert {)}^2,\\ f_{ia}^{\mathcal{A}}(t,R,T)& =(r_{ia}(t)-\Vert R_i{p}_i^l(t)+T_i-p_a^g(t)\Vert {)}^2.\end{array}$$

$$\begin{array}{rl}& \underset{R_i,T_i}{\mathrm{minimize}}f(R,T):=\sum _{i=1}^n{f}_i(R,T)\\ & \text{subject to}{R}_i\in \text{SO}(3),T_i\in {\mathbb{R}}^3,\forall i\in \mathcal{T}\end{array}$$

$$f_{i}(R,T)=\sum_{t=1}^{\overline{t}}\big{(}\sum_{j\in\mathcal{T}_{i}(t)}{f_{ij}^{\mathcal{T}}(t,R,T)}+\sum_{a\in\mathcal{A}_{i}(t)}{f_{ia}^{\mathcal{A}}(t,R,T)}\big{)}.$$

$$\mathcal{I}(\overline{t})={\bigcup }_{t=1}^{\overline{t}}\{p^l(t),p^g(t),r(t)\},$$

$$\underset{R,T}{\mathrm{minimize}}{\sum }_{t=1}^{\overline{t}}{f}_t(R,T)$$

$$\text{subject to}R\in \text{SO}(3),T\in {\mathbb{R}}^3$$

$$f_t(R,T)=(r(t)-\Vert R{p}^l(t)+T-p^g(t)\Vert {)}^2.$$

$$f_t(R,T)=\underset{y\in \mathcal{S}(t)}{\min }\Vert R{p}^l(t)+T-y{\Vert }^2,$$

$$\mathcal{S}(t)=\left\{y\in {\mathbb{R}}^3|\Vert y-p^g(t)\Vert =r(t)\right\}.$$

$$\begin{array}{rl}& \underset{R,T,y_{1:\overline{t}}}{\mathrm{minimize}}{\sum }_{t=1}^{\overline{t}}\Vert R{p}^l(t)+T-y(t){\Vert }^2\\ & \text{subject to}R\in \text{SO}(3),T\in {\mathbb{R}}^3,y(t)\in \mathcal{S}(t),t\in [1:\overline{t}].\end{array}$$

$$P_{\mathcal{S}(t)}(y)=p^g(t)+\frac{r(t)}{\Vert y-p^g(t)\Vert }(y-p^g(t))\text{if}y\notin \mathcal{S}(t).$$

$$\begin{array}{rl}{P}_{\text{SO}(3)}(\Omega )& =\underset{R\in \text{SO}(3)}{\mathrm{argmin}}\Vert R-\Omega {\Vert }_F^2\end{array}$$

$$P_{\text{SO}(3)}(\Omega )=UD{V}^{\ast },$$

$$D=\{\begin{array}{l}\text{diag}(1,1,+1),\text{ if }\text{det}(U{V}^{\ast })=1,\\ \text{diag}(1,1,-1),\text{ if }\text{det}(U{V}^{\ast })=-1.\end{array}$$

$$y^k(t)$$

$$=P_{\mathcal{S}(t)}(R^k{p}^l(t)+T^k),t\in [1:\overline{t}],$$

$$\underset{R\in \text{SO}(3),T}{\mathrm{minimize}}{\sum }_{t=1}^{\overline{t}}\Vert R{p}^l(t)+T-y^k(t){\Vert }^2.$$

$$\begin{array}{rl}{R}^{k+1}& =P_{\text{SO}(3)}(P^k),\\ T^{k+1}& ={\overline{y}}^k-R^{k+1}{\overline{p}}^l,\end{array}$$

$$P^k={\sum }_{t=1}^{\overline{t}}(y^k(t)-{\overline{y}}^k)(p^l(t)-{\overline{p}}^l{)}^{\prime }.$$

$$y^k(t)=P_{\mathcal{S}(t)}(R^k{p}^l(t)+T^k),t\in [1:\overline{t}]$$

$$\begin{array}{rl}{R}^{k+1}& =P_{\text{SO}(3)}(P^k),\\ T^{k+1}& ={\overline{y}}^k-R^{k+1}{\overline{p}}^l.\end{array}$$

$$M_T(R)=\frac{1}{2}(M-R{M}^{\prime }R).$$

$$y(t)=P_{\mathcal{S}(t)}(R(t-1)p^l(t)+T(t-1)),$$

$$(R(t),T(t))=\underset{R\in \text{SO}(3),T}{\mathrm{argmin}}\sum _{k=1}^t\Vert R{p}^l(k)+T-y(k){\Vert }^2,$$

$$\begin{array}{rl}\bar{y}(t)& =\bar{y}(t-1)+\frac{1}{t}(y(t)-\bar{y}(t-1)),\\ {\bar{p}}^l(t)& ={\bar{p}}^l(t-1)+\frac{1}{t}(p^l(t)-{\bar{p}}^l(t-1)),\\ P(t)& =P(t-1)+(1-\frac{1}{t})\\ & \times (y(t)-\bar{y}(t-1))(p^l(t)-{\bar{p}}^l(t-1){)}^{\prime }\end{array}$$

$$\begin{array}{rl}R(t)& =P_{\text{SO}(3)}(P(t)),\\ T(t)& =\bar{y}(t)-R(t){\bar{p}}^l(t).\end{array}$$

$$y(t):=P_{\mathcal{S}(t)}(R(t-1)p^l(t)+T(t-1)),$$

$$\begin{array}{rl}R(t)& =P_{\text{SO}(3)}(P(t)),\\ T(t)& =\bar{y}(t)-R(t){\bar{p}}^l(t).\end{array}$$

Full text

I Introduction

Cooperative localization is an important positioning technology [1, 2, 3, 4]. In the past decades, there are many methods for cooperative localization, such as semidefinite programming (SDP) [5], second-order cone programming [6], sum of squares [7], multidimensional scaling (MDS) [8], convex relaxation [9] and parallel projection algorithms (PPA) [10, 11]. Among them, PPA is reported to yield comparable accuracy to SDP and MDS with much shorter running time [11], and is an attractive localization approach.

By using both target-anchor and target-target range measurements, this work is concerned with cooperative localization problems over mobile sensor networks where anchor nodes are encoded with GPS positions and each target node is only aware of its position information in a local fixed coordinate frame, whose orientation and position relative to the global frame of the GPS are unknown. This framework is of great importance in both the underwater [12] and aerial localization [13]. For example, in case of multiple autonomous underwater vehicles (AUVs) the GPS information is often available to a very limited number of AUVs. Then, it is sensible to use cooperative methods to localize other AUVs with the inter-AUV range measurements [14, 15, 16, 17]. For unmanned aerial vehicles (UAVs), the target UAV in [13] is assumed to access to the Inertial Navigation System (INS), but the INS may continuously drift after initiation and lose the connection with the global coordinate system. That is, the GPS position of the target UAV is unavailable and requires to use inter-UAV range measurements for localization.

If a series of consistent positions in a local fixed coordinate frame can be obtained for a target node, its GPS position can be localized by aligning its local frame with the global frame of the GPS by using target-anchor measurements. To this end, a natural way is to parameterize the local coordinate frame by a rotation matrix $R\in\text{SO}(3)$ and a translation vector $T\in\mathbb{R}^{3}$. Then, the alignment problem reduces to the estimation of $(R,T)$, which is the key idea of [12, 13] and is also closely related to the idea of estimating the deviation of the local coordinate from the global coordinate in [14, 18, 19].

This work starts from investigating the problem of localizing a target node with an anchor node. The least squares estimate of $(R,T)$ can be obtained by solving an optimization problem with a non-convex objective function and non-convex constraints. Such a non-convex optimization problem in [13] is firstly relaxed as a SDP problem with $11$ equality constraints and the decision vector is a $17\times 17$ positive semi-definite matrix, hoping that the solution to the SDP problem can provide a good suboptimal solution. To further refine the SDP solution, they design a gradient descent algorithm over the rotation group $\text{SO}(3)$. Differently from [13], we exploit the geometric relations between nodes and reformulate the non-convex problem as a well-structured optimization problem with a convex cost over spherical constraints. This idea was presented in our conference paper [12], the major results of which are all contained in Section III(A)-(B) of this work.

The striking feature of our approach is that we are able to simultaneously solve the coordinate alignment problem for multiple target nodes in a general sensor network by using both target-anchor and target-target range measurements. Note that the authors in [13] only consider the case with only a target node, and is unclear how to extend to the general case with multiple target nodes.

With the aid of the block coordinate descent method [20], we propose a parallel projection algorithm (PPA) to solve the above reformulated problem. The projection is with respect to the spherical surfaces and can be explicitly written in a simple form, after which the constraint $R\in\text{SO}(3)$ can also be easily resolved. Overall, the iteration of the PPA is given in a simple form and can be implemented with a low computational cost, which is important to the sensor network. In comparison with [13], the PPA requires a much lower computational cost with comparable localization accuracy, both of which have been validated via numerical experiments.

Interestingly, the PPA can easily incorporate new measurements to update our estimate of $(R,T)$. Specifically, we propose a recursive version of the PPA, which is termed as recursive projection algorithm (RPA), to approximately solve the optimization problem for coordinate alignment. More importantly, we are able to extend our method to the case of multiple target nodes in a mobile sensor network. For a time-varying network, we further use the block coordinate descent idea to design the PPA to reduce the computational load. For a time-invariant network, we jointly use the Jacobi iterative method to run the PPA and obtain a distributed PPA, which only requires each target node to exchange information with its neighboring target nodes.

The rest of this paper is organized as follows. In Section II, we formulate the coordinate alignment over a time-varying network as a non-convex optimization. In Section III, focusing on two-node coordinate alignment problem, we propose the PPA and RPA. In Section IV, we extend them to the multi-node setting and propose a PPA algorithm by using the block coordinate descent idea. For a fixed communication graph, a distributed method with the Jacobi iteration is designed. The numerical experiments are conducted in Section V. Finally, some concluding remarks are drawn in Section VI.

II Problem Statement

II-A The mobile sensor network

The mobile sensor network is represented by a sequence of time-varying graphs $\mathcal{G}(t)=(\mathcal{V},\mathcal{E}(t))$ where $\mathcal{V}$ is the set of a fixed number of mobile nodes and $\mathcal{E}(t)$ is the set of edges between nodes at discrete time $t\in\mathbb{N}$. Specifically, $\mathcal{V}$ is the union of a target set $\mathcal{T}=\{1,\ldots,n\}$ and an anchor set $\mathcal{A}=\{n+1,\ldots,n+r\}$ where an anchor node can access its position information in the GPS while a target node does not and is only aware of its position information in a local fixed coordinate frame whose orientation and position relative to the global frame of the GPS are unknown. See an example of collaborative UAVs in Section I. For brevity, the former is called the global position and the later is called local position. Our objective is to localize the global positions of target nodes under information flow constraints, which are modeled by the graph $\mathcal{G}(t)$.

Specifically, for target node $i$ and anchor node $a$, $(i,a)\in\mathcal{E}(t)$ identifies the communication from $a$ to $i$. For any pair of target nodes $i$ and $j$ such that $(i,j)\in\mathcal{E}(t)$, then $(j,i)\in\mathcal{E}(t)$ and both nodes can communicate with each other. Moreover, two noisy range measurements $d_{ij}(t)$ and $d_{ji}(t)$ are taken by node $i$ and node $j$, respectively. Note that $d_{ij}(t)$ and $d_{ji}(t)$ may not be equal due to the use of different range sensors. While for a target node $i$ and an anchor node $a$ such that $(i,a)\in\mathcal{E}(t)$, only the range measurement $r_{ia}(t)$ is available to node $i$ and $(a,i)\notin\mathcal{E}(t)$. A target node $i$ is said to be connected to an anchor node $a$ in $\mathcal{G}(t)$ if there is a path of consecutive edges in $\mathcal{E}(t)$ that connects the two nodes. Given a target node $i$, let $\mathcal{T}_{i}(t)$ be the set of its neighboring target nodes, i.e., $\mathcal{T}_{i}(t)=\left\{j|j\in\mathcal{T},(i,j)\in\mathcal{E}(t)\right\}$ and $\mathcal{A}_{i}(t)$ is the set of neighboring anchor nodes, i.e., $\mathcal{A}_{i}(t)=\left\{a|a\in\mathcal{A},(i,a)\in\mathcal{E}(t)\right\}$. Thus, the set of range measurements available to the target node $i$ at time $t$ is given as

[TABLE]

II-B Coordinate alignment for cooperative localization

Let $p_{a}^{g}(t)\in\mathbb{R}^{3}$ be the global position of an anchor node $a$ and $p_{i}^{l}(t)\in\mathbb{R}^{3}$ be the local position of a target node $i$ at time $t$, whose local coordinate system is parameterized by a rotation matrix $R_{i}^{*}\in\text{SO}(3)$ which is defined as

[TABLE]

and a translation vector $T_{i}^{*}\in\mathbb{R}^{3}$. Clearly, the global position of the target node $i$ is expressed as $R_{i}^{*}p_{i}^{l}(t)+T_{i}^{*}.$ To localize the target node $i$, we aim to compute its coordinate parameters $(R_{i}^{*},T_{i}^{*})$ with noisy range measurements up to time $\bar{t}$, i.e.,

[TABLE]

where $e=(i,j)$ or $e=(i,a)$ indicates an edge and $\{\xi_{e}(t)\}_{t=1}^{\overline{t}}$ is a sequence of temporally uncorrelated with zero mean and the sequence $\{\xi_{e}(t)\}_{e\in\mathcal{E}(t)}$ is spatially uncorrelated at any time $t$, and $[1:\bar{t}]=\{1,\ldots,\bar{t}\}$. Given a pair of $R=[R_{1}^{\prime},\ldots,R_{n}^{\prime}]^{\prime}$ and $T=[T_{1}^{\prime},\ldots,T_{n}^{\prime}]^{\prime}$, the least squares estimate uses the quadratic loss function

[TABLE]

Our coordinate alignment problem for cooperative localization is formulated as a constrained optimization problem

[TABLE]

where each summand in the objective function is given by

[TABLE]

Under mild conditions, we show that the constrained optimization problem in (3) is solvable.

Proposition 1

If each target node is connected to an anchor node in the union graph $\bigcup_{t=1}^{\overline{t}}\mathcal{G}(t)$, then the constrained optimization problem in (3) contains at least an optimal solution.

Proof:

See Appendix -A.

In the sequel, we shall design explicit algorithms to solve the optimization problem in (3) by using projection technique.

III Localizing a Target Node with an Anchor Node

In this section, we consider the problem of localizing only a mobile target node with an anchor node. This is well motivated by localizing a GPS-denied AUV. Another AUV with known global position is deployed to serve as a communication and navigation aid (CNA) [19, 14]. They cooperatively work in the underwater and communicate with each other to obtain a series of range measurements, see Fig. 1. In this case, the minimum number of range measurements is $7$ [21].

To simplify notations of this section, let $p^{g}(t)\in\mathbb{R}^{3}$ be the global position of the anchor node, $p^{l}(t)\in\mathbb{R}^{3}$ be the local position of the GPS-denied target node and $r(t)\in\mathbb{R}$ be the range measurement between the two nodes at time $t$. Then, the information set for the target localization performed in the time interval $\bar{t}$ is collectively given by

[TABLE]

and the optimization problem in (3) is reduced as

[TABLE]

where the summand in the objective function is

[TABLE]

III-A Optimization problem reformulation using projection

To solve the optimization problem in (III), there are at least two challenges. The first is that $f_{t}(R,T)$ is non-convex, which usually is approximately solved by the convex relaxation [10, 9, 13, 22]. Here we solve it by expressing as the minimization of a convex function over a spherical surface. The second lies in the constraint set of a rotation group $\text{SO}(3)$, which fortunately can be explicitly solved as well.

One can show that $f_{t}(R,T)$ is the squared range between the point $Rp^{l}(t)+T$ and the spherical surface centered at $p^{g}(t)$ with a radius $r(t)$[9], see Fig. 2. That is,

[TABLE]

where $\mathcal{S}(t)$ is a spherical surface, i.e.,

[TABLE]

In view of (8), we obtain the following optimization problem

[TABLE]

Remark 1

Note that target localization is not instantaneous but performed in time interval $\bar{t}$. In [9], the so-called disk relaxation is adopted by relaxing the spherical surface $\mathcal{S}(t)$ into a closed ball $\mathcal{B}(t)=\{y|\lVert y-p^{g}(t)\rVert\leq r(t)\}$. This leads to an underestimated convex problem, and is useless here as $\text{SO}(3)$ is not convex.

Clearly, the two optimization problems in (III) and (10) are essentially equivalent in the sense that both achieve the same minimum value and the same optimal set of $(R,T)$. The good news is that the optimization problem (10) has favorable properties. First, its objective function is quadratically convex. Second, the newly introduced sets $\mathcal{S}(t),t\in[1:\overline{t}]$ are spherical surfaces which are not difficult to compute the associated Euclidean projection. In fact, given a vector $y\in\mathbb{R}^{3}$, its Euclidean projection onto a spherical surface $\mathcal{S}(t),t\in[1:\overline{t}]$ is explicitly expressed as

[TABLE]

To be specific, the projection of any matrix $\Omega\in\mathbb{R}^{3\times 3}$ onto $\text{SO}(3)$ is obtained by solving a constrained optimization problem, i.e.,

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm. In view of [23], $P_{\text{SO}(3)}(\Omega)$ is explicitly given as

[TABLE]

where $U$ and $V$ are obtained via the singular value decomposition of $\Omega$, i.e., $\Omega=U\Sigma V^{*}$, and

[TABLE]

Next, we shall design algorithms to effectively solve the optimization problem (10).

III-B Parallel projection algorithm

Once the target node has access the information set $\mathcal{I}(\bar{t})$ in (5), it solves the optimization problem (10) by a block coordinate descent algorithm [20] with parallel projections. We use master and worker to denote the order of updating per iteration. Specifically, one master is used to update $(R^{k},T^{k})$ and $\overline{t}$-parallel workers are responsible for simultaneously updating $y^{k}(t),t\in[1:\overline{t}]$. The superscript $k$ denotes the number of iterations for solving the optimization problem (10).

At the $k$-th iteration, each worker receives the latest update $(R^{k},T^{k})$ from the master, and then performs the following projection in a parallel way

[TABLE]

where $P_{\mathcal{S}(t)}(\cdot)$ is given in (11), and sends $y^{k}(t)$ to the master.

Once the master receives $(y^{k}(1),\ldots,y^{k}(\overline{t}))$, it solves the following constrained least squares optimization

[TABLE]

Proposition 2

The optimization problem in (14) is explicitly solved as

[TABLE]

where $P_{\text{SO}(3)}(\cdot)$ is given in (12), $\overline{p}^{l}=\frac{1}{\overline{t}}\sum\nolimits_{t=1}^{\overline{t}}{p^{l}(t)}$ and $\overline{y}^{k}=\frac{1}{\overline{t}}\sum\nolimits_{t=1}^{\overline{t}}{y^{k}(t)}$ are “mean” vectors of $\{p^{l}(t)\}$ and $\{y^{k}(t)\}$, and $P^{k}$ is their “correlation” matrix

[TABLE]

Proof:

See Appendix -B.

is interesting that (14) is closely related to the basic Procrustes problem [24] and can be found in its full version in [25]. For completeness, we also include a proof in Appendix. Finally, we summarize the above result in Algorithm 1.

Remark 2

*Instead of using a SDP initialization [13], we just randomly select a pair of $(R^{0},T^{0})$. Clearly, we can also adopt the same initialization to avoid getting into a bad local minimum.  *

Since the optimization problem in (10) is inherently non-convex, it cannot be guaranteed to converge to a global optimal solution. However, it at least sequentially reduces the objective function per iteration, and achieves a better solution. To exposit it, let $q:=(R,T,y_{1},\ldots,y_{\bar{t}})$ and $g(q):=\sum\nolimits_{t=1}^{\overline{t}}{\lVert Rp^{l}(t)+T-y(t)\rVert^{2}}$ be the decision variables and the objective function, respectively. We have the following result.

Proposition 3

Let $\{q^{k}\}$ be iteratively computed in Algorithm 1. Then, it holds that $g(q^{k})\leq g(q^{k-1}),\forall k$ and there exists a convergent subsequence of $\{q^{k}\}$.

Proof:

See Appendix -C.

Remark 3

In [13], a semidefinite programming (SDP) relaxation is firstly devised to find an initial estimate of $(R,T)$, which involves solving a SDP with $11$ equality constraints and the decision vector is a $17\times 17$ positive semi-definite matrix. Then, they solve the optimization problem in (III) by using the projection of the gradient $M:=\frac{\partial}{\partial R}\sum\nolimits_{t=1}^{\overline{t}}f_{t}(R,T)$ onto the tangent space of $\text{SO}(3)$, which is explicitly given as

[TABLE]

The discretized version is essentially gradient descent (GD) and given by $R^{k+1}=P_{\text{SO}(3)}\left(R^{k}-\alpha^{k}M_{T}(R^{k})\right),$ where $\alpha^{k}$ is a stepsize. Notably, they also explicitly state (without proof) that the SDP relaxation is important in providing a good initialization. Solving such a SDP and extracting a feasible $R^{0}\in\text{SO}(3)$ from the SDP’s solution inevitably increases the computation cost. Though Algorithm 1 is only randomly initialized, numerical results show that its localization accuracy is still comparable to that of [13].

*More importantly, the focus of the equivalent optimization problem in (10) allows us easily to devise a recursive algorithm to estimate $(R,T)$ in an online way (c.f. Section III-C) and generalize to the case of generic mobile sensor networks (c.f. Section IV). It is worthy mentioning that the approach in [13] currently only applies to a star topology.  *

III-C The recursive projection algorithm

While Algorithm 1 produces good results if $\overline{t}$ is moderately large, it does not exploit the sequential collection of the measurement, and the number of local intermediate variables $y(t)$ increases linearly with the number of range measurements. To resolve it, this subsection presents an approximate Recursive Projection Algorithm (RPA) which only performs one iteration whenever new measurement arrives.

At time $t$, suppose we have already obtained a prior estimate $(R(t-1),T(t-1))$ and collected a new measurement $\{p^{l}(t),p^{g}(t),r(t)\}$. Using this information, we shall recursively update the estimate of $(R,T)$ in an online way.

Similar to (13), we perform an online projection

[TABLE]

where $\mathcal{S}(t)$ is defined in (9). In comparison with (14), the projection operation for $y(t)$ is only performed once. Then, the new estimate of $(R,T)$ is set as follows

[TABLE]

which can be recursively computed.

Proposition 4

Let $\bar{y}(t),\bar{p}^{l}(t)$ and $P(t)$ be recursively computed by

[TABLE]

where $\bar{y}(0)=\bar{p}^{l}(0)=0\in\mathbb{R}^{3}$ and $P(0)=0\in\mathbb{R}^{3\times 3}$. Then, the optimization problem in (18) is solved by

[TABLE]

Proof:

See Appendix -D.

The recursive algorithm is summarized in Algorithm 2. In practice, we shall further adopt the idea of smoothing [26] to improve the algorithmic performance. Instead of solving (18), it is better to consider

[TABLE]

where $\acute{y}(k)=P_{\mathcal{S}(k)}(R(t-1)p^{l}(k)+T(t-1))$ if $k\in[t-b,t-1]$ and $\acute{y}(k)=y(k)$ if $k<t-b$. Here $b$ denotes the length of smoothing interval and indicates the tradeoff between computational cost and performance improvement. Clearly, Algorithm 2 corresponds to the special case $b=1$. Then, the optimization problem in (21) can be recursively solved. Let $\Delta(k)=\acute{y}(k)-y(k)$, which is zero if $k<t-b$, and compute $\bar{\Delta}(t)=\frac{1}{t}\sum_{k=1}^{t}\Delta(k)$. We solve it by replacing $\bar{y}(t)$ and $P(t)$ with $\bar{y}(t)+\bar{\Delta}(t)$ and

[TABLE]

in (20), respectively.

Since the localization problem is typically non-convex, we are unable to prove the asymptotic convergence of $(R(t),T(t))$. Jointly with (17) and (18), one may also use a discount factor $\alpha\in(0,1)$ to emphasize the importance of the latest range measurements, e.g.,

[TABLE]

and replace $1/t$ in (19) by $(1-\alpha)/(1-\alpha^{t})$.

IV Localizing Multiple Target Nodes in the Sensor Network

In this section, we are interested in the localization problem of multiple target nodes in the mobile sensor network with generic time-varying communication topology $\mathcal{G}(t)$. In [13], the SDP based approach can only deal with the network setting that the only one anchor is connected to all target nodes. Such a scenario gives a star topology, which is trivial to treat by using the results on the situation with one anchor and one target node. While for general mobile sensor networks, they leave it to future work. By using the approach in Section III, we are able to solve this problem, which is the focus of this section.

IV-A Optimization problem reformulation using projection

The loss function (4) introduces coupled summands, which makes the problem difficult. We shall use the projection idea in Section III to reformulate the optimization problem in (3). As in (9), define the spherical surfaces

[TABLE]

In view of (8), the loss functions in (2) are rewritten as

[TABLE]

With a slight abuse of notations, let

[TABLE]

Jointly with (22), the problem in (3) can be reformulated as

[TABLE]

where the summand is given by

[TABLE]

In the sequel, we shall design a block coordinate descent algorithm to solve the optimization problem in (24).

IV-B Parallel projection algorithms

Clearly, the objective function in (24) is quadratically convex. We only need to handle the non-convex constraints $\text{SO}(3)$ and spherical surfaces $\mathcal{S}(t),t\in[1:\overline{t}]$.

Now, we design a block coordinate descent algorithm [20] with parallel projections to solve (24). Specifically, given $(R^{k},T^{k})$, we update $y(t),t\in[1:\overline{t}]$ by

[TABLE]

which can be explicitly expressed as

[TABLE]

and the projection onto a spherical surface is given in (11).

Next, we shall update $(R,T)$ by fixing $y(t)=y^{k}(t)$, i.e.,

[TABLE]

To solve the above optimization problem, the major difficulty lies in the constraints of $\text{SO}(3)$. Two ideas are adopted.

IV-B1 Constrained least squares

The first idea is to solve an unconstrained least squares problem, i.e.,

[TABLE]

and then project $Z^{k}$ onto the constraints of $\text{SO}(3)^{n}$, i.e.,

[TABLE]

which is explicitly given in (12).

The remaining problem is how to effectively solve the least squares problem in (27). For this purpose, we represent $Z_{i}p_{i}^{l}(t)+T_{i}$ as a linear function of $\bm{x}_{i}$, where $\bm{x}_{i}\in\mathbb{R}^{12}$ is a column vector reshaping from $(Z_{i},T_{i})$. Specifically, denote

[TABLE]

where $\otimes$ denotes the Kronecker product, $\text{vec}(Z_{i})\in\mathbb{R}^{9}$ is a large vector by stacking all the columns of $Z_{i}\in\mathbb{R}^{3\times 3}$, and $I_{3}\in\mathbb{R}^{3\times 3}$ is an identity matrix. Then, it follows that

[TABLE]

and the objective function in (27) is rewritten as

[TABLE]

which clearly is quadratic in the decision vector ${\bm{x}}$.

For a graph $\mathcal{G}(t)$, we define a sparse block matrix $E(t)\in\mathbb{R}^{|\mathcal{E}(t)|\times n}\otimes\mathbb{R}^{3\times 12}$ over the graph for a compact form of $f({\bm{x}})$. Particularly, if $e=(i,j)\in\mathcal{E}(t)$ and $j\in\mathcal{T}_{i}(t)$, then the $(e,i)$-th block of $E(t)$ is $B_{i}(t)$ and the $(e,j)$-th block of $E(t)$ is $-B_{j}(t)$. If $e=(i,a)\in\mathcal{E}(t)$ and $a\in\mathcal{A}_{i}(t)$, then the $(e,i)$-th block of $E(t)$ is $B_{i}(t)$. All the unspecified blocks are set to be zero matrices with compatible dimensions. This implies that the objective function in (27) can be compactly expressed as

[TABLE]

Clearly, the minimizer of $f({\bm{x}})$ is simply given by

[TABLE]

To compute the above ${\bm{x}}^{ls}$, let $Q(t)=E(t)^{\prime}E(t)\in\mathbb{R}^{n\times n}\otimes\mathbb{R}^{12\times 12}$. Denote the $(i,j)$-th block of $Q(t)$ by $Q(t)_{(i,j)}\in\mathbb{R}^{12\times 12}$, it follows that

[TABLE]

and

[TABLE]

where $|\mathcal{T}_{i}(t)|$ and $|\mathcal{A}_{i}(t)|$ denote the cardinality of the sets $\mathcal{T}_{i}(t)$ and $\mathcal{A}_{i}(t)$ respectively, and

[TABLE]

Similarly, the $i$-th block of $E(t)^{\prime}{y}^{k}(t)$ is defined as $E(t)^{\prime}{y}^{k}(t)_{(i)}$ and given by

[TABLE]

Let $\widehat{y}_{i}^{k}(t)=\sum\nolimits_{j\in\mathcal{T}_{i}(t)}\left(y^{k}_{ij}(t)-y^{k}_{ji}(t)\right)+\sum\nolimits_{a\in\mathcal{A}_{i}(t)}{y^{k}_{ia}(t)}$, then

[TABLE]

Jointly with (29)-(31), the minimizer in (28) can be readily computed. If the graph $\mathcal{G}(t)$ is fixed, (28) can be cast as a sparse least squares problem, see e.g. [27] for details.

IV-B2 Jacobi iterative method

We can also solve the optimization problem in (26) by using the Jacobi iterative method [20]. Particularly, we compute $(R_{i},T_{i})^{k+1}$ by setting $(R_{-i},T_{-i})$ to be $(R_{-i},T_{-i})^{k}$, where $(R_{-i},T_{-i})=\{(R_{j},T_{j})\}_{j\in\mathcal{T},j\not=i}$, i.e.,

[TABLE]

where the objective collects all summands in the objective function of (26) containing the decision variables $(R_{i},T_{i})$, and $g_{i}(R,T,y^{k}(t))$ is given by

[TABLE]

Then, the optimization problem in (32) has a similar structure to that of (14), and can be solved as

[TABLE]

where $\overline{y}_{i}^{k}$ and $\overline{p}_{i}^{l}$ are two mean vectors and $P_{i}^{k}$ is a correlation matrix, i.e.,

[TABLE]

Different from (28), the Jacobi iterative method does not need to solve the least squares problem in (27), which may need to compute the inverse of a large matrix, i.e., $\sum_{t=1}^{\bar{t}}E(t)^{\prime}E(t)\in\mathbb{R}^{n\times n}\otimes\mathbb{R}^{12\times 12}$. Instead, we only need to use (33) to replace Step 5 in Algorithm 3.

IV-C Distributed implementation of Jacobi method for fixed graphs

Centralized algorithms are not scalable for the large network. If $\mathcal{G}(t)$ is fixed, the Jacobi iterative method can even be implemented in a distributed way, which is termed as DPPA and given in Algorithm 4.

Remark 4

By Proposition 2.3.1 [20], we can obtain the similar result as Proposition 3 for Algorithms 3-4. Take Algorithm 3 as an example. Let ${\bf y}=[y(1),\ldots,y(\bar{t})]$ and

[TABLE]

Then, there is a convergent subsequence of $\{(R^{k},T^{k},{\bf y}^{k})\}$. To elaborate it, we obtain from (27) that $g(Z^{k},T^{k+1},{\bf y}^{k})\leq g(R^{k},T^{k},{\bf y}^{k})$. Since the projection operator is non-expansive, it implies that $g(R^{k+1},T^{k+1},{\bf y}^{k})\leq g(Z^{k},T^{k+1},{\bf y}^{k})$. By (25), it holds that $g(R^{k+1},T^{k+1},{\bf y}^{k+1})\leq g(R^{k+1},T^{k+1},{\bf y}^{k})$. Combining the above, it finally yields that

[TABLE]

*The rest of proof follows exactly the same as that of Proposition 3.  *

V Numerical Experiments

In this section, we perform numerical experiments to validate the proposed algorithms in Python 2.7 environment on a MacBook Pro with 2.2 GHz Intel Core i7 CPU and 16GB DDR3. Open source packages such as Numpy 1.12.1 and cvxopt 1.1.9 are used for numerical computation. The experiments are implemented in both two dimensional space and three dimensional space. As there is no difference between the two cases, we only report results of the two dimensional case for visualization convenience.

V-A Experiment setup

For the two-node localization problem, the coordinate system of the target node is generated by a rotation matrix $R$ and a transformation vector $T$ as follow

[TABLE]

where the rotation angle $\theta\sim\mathcal{U}(0,2\pi)$, and $a,b\sim\mathcal{U}(-1,1)$ are randomly selected with uniform distributions. The target node and the anchor node are randomly moving in a square area $[1,9]\times[1,9]$. We also randomly select $p^{l}(t)$ and $p^{g}(t)$ such that $p^{g}(t)\in[1,9]\times[1,9]$ and $Rp^{l}(t)+T\in[1,9]\times[1,9]$.

Then their range measurements at time slot $t$ are generated by $r(t)=\|Rp^{l}(t)+T-p^{g}(t)\|+\xi$ where the random noise is $\xi\sim N(0,\sigma^{2})$. To quantify the noise level, define the signal-to-noise ratio (SNR) by

[TABLE]

where $d_{0}=4.1712$ is the average range of two nodes in the area $[1,9]\times[1,9]$. Clearly, a smaller SNR means a higher noise level. Our objective is to compute the coordinate system parameters under different signal-to-noise ratios by the proposed algorithms, which are denoted as $\widehat{R}$ and $\widehat{T}$. We are concerned with their relative errors

[TABLE]

For each target node $i\in\mathcal{T}$ in the multi-node localization problem, $\widehat{R}_{i}$ and $\widehat{T}_{i}$ are denoted as the same way as that in the two-node localization scenario. Similarly, all nodes are limited to the square area $[1,9]\times[1,9]$.

V-B Experimental results of the two-node localization problem

We compare the proposed PPA with the SDP based method [13]. Since the method in [13] is unable to deal with general multi-node situations, we only compare their algorithm for localizing one target node with one anchor node.

Numerical experiments are performed under two noise levels111Kindly note that the results of $\text{SNR}=10,20,30$ are consistent, we only report the case of $\text{SNR}=20$ for saving space. ($\text{SNR}=20$ and $\text{SNR}=80$) and three different $\overline{t}$. The results in Fig. 3 are obtained by averaging over $10^{4}$ independent simulations. Since rotation matrices are more difficult to estimate, we choose to report results mostly on rotations and only include the final results on translations for saving space. The green line corresponds to the use of the pure SDP, and the red line is the result of the PPA of Algorithm 1. The blue line is the result of the GD with the SDP initialization [13], i.e., SDP+GD, while the purple line is the result of the PPA with the SDP initialization, i.e., SDP+PPA. We also record the time used for running different algorithms. In Fig. 3(a), it takes 4.96e-01s to find the SDP based solution. To achieve the same relative rotation error, it only takes 2.19e-03s by using PPA. Moreover, it only takes 4.20e-03s for PPA to outperform the SDP+GD, whose running time is (4.96e-01+6.88e-03)s. We also observe that the SDP+PPA finally achieves the smallest relative rotation error. If the SNR is large, see Fig.3(b), the PPA cannot reduce the relative rotation error as small as that of the SDP due to the use of random initialization and the gap induced by the SDP relaxation decreases with SNR. However, the SDP+PPA performs much better than the SDP+GD, both in terms of running time and accuracy. In Table I, we include the final results, i.e., the number of iterations is set to $1000$, on the relative translation errors when $\bar{t}=20$. In summary, both Fig. 3 and Table I consistently validate the advantages of the PPA of Algorithm 1.

Next, the performance of the RPA of Algorithm 2 is shown in Fig. 4, which illustrates that the relative error of the rotation matrix essentially decreases with the number of range measurements. Due to the use of approximation in deriving the RPA of (17), it further induces performance degradation in comparison with the PPA. Note that the method in [13] is unable to write in a recursive form.

V-C Experimental results of the multi-node localization problem

In this subsection, we apply Algorithm 3 to a sensor network which contains $110$ target nodes and $4$ anchor nodes. All the target nodes are randomly deployed in a two dimensional area $[1,9]\times[1,9]$, and $4$ anchor nodes are located at $(2,2),(2,8),(8,2),(8,8)$ respectively. Each node moves randomly in a unit square centered at its initial position.

Two nodes can communicate only if their distance is within $1$, which clearly results in time-varying communication graphs. The SNR of each node is set to $\text{SNR}=100$. The localization results of target nodes at time slots $\overline{t}=5$ and $\overline{t}=25$ are presented in Fig. 5. We observe that the localization accuracy is improved when $\overline{t}$ increases and all target nodes are well localized.

At the time slot $\overline{t}=25$, we count the number of target-anchor range measurements for each target node, which is shown in Table II. One can observe that in our cooperative localization method, more than a half ($50.9\%$) of target nodes have never directly taken range measurements with respect to any anchor node. However, their positions can also be successfully localized by Algorithm 3 as shown in Fig. 5(b), which confirms the benefit of using cooperative methods.

Finally, we compare Algorithm 3 with the DPPA of Algorithm 4 in a fixed graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with $110$ target nodes and $4$ anchor nodes. Note that the DPPA is only applicable to a fixed graph. Define the average degree by

[TABLE]

which characterizes the edge density of a network. By varying the average degree and the SNR, we implement both algorithms using the range measurements in a period of time $\overline{t}$ ($\overline{t}$ is set from $5$ to $25$). The resulting coordinate alignment relative errors are presented in Fig.6, which illustrate that their performances are very close, and increasing the length of time interval $\overline{t}$ or the network density, both algorithms lead to better estimates. However, we recommend to use DPPA for a fixed graph as it involves simpler iterations and is a distributed version.

VI Conclusion

This work considers the cooperative localization as a coordinate alignment problem using range measurements. To align the coordinate of a target node with an anchor node, we present PPA and RPA respectively. Then, the algorithms are generalized to the case of multiple target nodes in a sensor network. The effectiveness of all algorithms have been validated by numerical experiments. The state-of-the-art works such as the SDP and the SDP+GD are also compared with our work, which confirms the advantages of the proposed algorithms.

-A Proof of Proposition 1

Proof:

We first prove that $f(R,T)$ is coercive [20] with respect to $T$, i.e.,

[TABLE]

where $\|T\|=(\sum_{i=1}^{n}\|T_{i}\|^{2})^{1/2}$. Suppose that $\|T\|\rightarrow\infty$, then there must exist some target node $i$ such that $\|T_{i}\|\rightarrow\infty$. We have two exclusive scenarios.

If $\bigcup_{t=1}^{\overline{t}}\mathcal{A}_{i}(t)$ is nonempty, e.g., there exists an anchor node $a$ such that $a\in\mathcal{A}_{i}(t)$ for some $t$, i.e., the component $f_{ia}^{\mathcal{A}}(t)$ exists in the objective function. Then, one can easily verify that $f_{ia}^{\mathcal{A}}(t)$ tends to infinity as $\|T_{i}\|\rightarrow\infty$. This implies that $\lim_{\|T_{i}\|\rightarrow\infty}f(R,T)=\infty$.

If $\bigcup_{t=1}^{\overline{t}}\mathcal{A}_{i}(t)$ is empty, the target node $i$ must connect to an anchor node via some target node $j$ with a nonempty $\bigcup_{t=1}^{\overline{t}}\mathcal{A}_{j}(t)$ in the union graph $\bigcup_{t=1}^{\overline{t}}\mathcal{G}(t)$ since otherwise, the target node $i$ is disconnected to anchor nodes. Particularly, let $(j_{0},j_{1}),\ldots,(j_{k-1},j_{k})\in\bigcup_{t=1}^{\overline{t}}\mathcal{E}(t)$ be the consecutive edges from node $i=j_{0}$ to node $j=j_{k}$. Suppose

[TABLE]

it follows from (2) that $\|T_{j_{0}}\|=\ldots=\|T_{j_{k}}\|=\infty$. Since $\bigcup_{t=1}^{\overline{t}}\mathcal{A}_{j}(t)$ is nonempty, it immediately implies that $\lim_{\|T_{j}\|\rightarrow\infty}f(R,T)=\infty$.

Overall, (35) is proved. Since the rotation group $\text{SO}(3)$ is compact and $f(R,T)$ is continuous, the rest of proof follows from the Weierstrass’ theorem [20, Proposition A.8].

-B Proof of Proposition 2

Proof:

For any fixed $R\in\text{SO}(3)$, it is obvious that $T=\overline{y}^{k}-R\overline{p}^{l}$ minimizes the objective function of (14) with respect to $T$. Let $T=\overline{y}^{k}-R\overline{p}^{l}$ in the objective function of (14). Then, it follows that

[TABLE]

where $c$ is independent of $R$ and is not explicitly given here.

Then, $R^{k+1}$ is obtained via the minimization problem

[TABLE]

where the second equality follows from the fact that $\|R-P^{k}\|_{F}^{2}=\text{trace}((R-P^{k})^{\prime}(R-P^{k}))=-2\text{trace}(R^{\prime}P^{k})+(P^{k})^{\prime}P^{k}+I$ for any $R\in\text{SO}(3)$.

-C Proof of Proposition 3

Proof:

By (13) and $y^{k-1}(t)\in\mathcal{S}(t),\forall t\in[1:\bar{t}]$, it follows $\|R^{k}p^{l}(t)+T^{k}-y^{k}(t)\|\leq\|R^{k}p^{l}(t)+T^{k}-y^{k-1}(t)\|,$ which implies that

[TABLE]

By (14), we obtain that $\sum\nolimits_{t=1}^{\overline{t}}{\lVert R^{k}p^{l}(t)+T^{k}-y^{k-1}(t)\rVert^{2}}\leq\sum\nolimits_{t=1}^{\overline{t}}{\lVert Rp^{l}(t)+T-y^{k-1}(t)\rVert^{2}}$ for all $R\in\text{SO}(3)$ and $T\in\mathbb{R}^{3}$. Since $R^{k-1}\in\text{SO}(3)$ and $T^{k-1}\in\mathbb{R}^{3}$, this implies that

[TABLE]

Thus, it holds that $g(q^{k})\leq g(q^{k-1})$. Since $\text{SO}(3)$ and $\mathcal{S}(t)$ are compact, it follows from (15) that $\{q^{k}\}$ is a bounded sequence. Thus, it contains a convergent subsequence.

-D Proof of Proposition 4

Proof:

Clearly, both $\bar{y}(t)$ and $\bar{p}^{l}(t)$ compute the time average of their associated vectors and can be expressed as

[TABLE]

Moreover, it holds that

[TABLE]

In fact, let $\tilde{y}(t)=y(t)-\bar{y}(t-1)$ and $\tilde{p}^{l}(t)=p^{l}(t)-\bar{p}^{l}(t-1)$. Then, it follows from (19) that

[TABLE]

The rest of proof follows directly from that of Proposition 2 and is omitted.