Efficient Distributed Method for NLOS Cooperative Localization in WSNs

The accuracy of cooperative localization can be severely degraded in non-line-of-sight (NLOS) environments. Although most existing approaches modify models to alleviate NLOS impact, computational speed does not satisfy practical applications. In this paper, we propose a distributed cooperative localization method for wireless sensor networks (WSNs) in NLOS environments. The convex model in the proposed method is based on projection relaxation. This model was designed for situations where prior information on NLOS connections is unavailable. We developed an efficient decomposed formulation for the convex counterpart, and designed a parallel distributed algorithm based on the alternating direction method of multipliers (ADMM), which significantly improves computational speed. To accelerate the convergence rate of local updates, we approached the subproblems via the proximal algorithm and analyzed its computational complexity. Numerical simulation results demonstrate that our approach is superior in processing speed and accuracy to other methods in NLOS scenarios.


Introduction
Wireless-sensor-network (WSN) technology has rapidly developed because of its convenience and prospects. It can significantly improve living quality in many important fields, such as environment monitoring [1] and surveillance [2], vehicle tracking [3,4], exploration [5,6], and other sensing tasks [7]. According to Reference [1], WSNs extend the human ability "to monitor and control physical world". This means that WSNs provide a more intelligent way to link the physical world with humans. It is worth noting that the positions of sensor nodes are key information for WSNs to fulfil various tasks. Therefore, as a preliminary task, cooperative localization has aroused increasing interest, especially in indoor scenarios where satellite communications cannot be employed. In general, cooperative localization has two main categories: rangefree methods and range-based methods [8]. Rangefree methods are easier to implement, but their accuracy is lower than range-based methods [9]. The most common range-measurement techniques of range-based localization are based on received signal-strength indicator (RSSI) [10] and time of flight (TOF) [11]. The RSSI technique has a lower cost than the TOF technique, but the latter has higher accuracy. In this paper, we mainly study the localization problem based on range-based technology. Some related work in this field is listed as follows.
These distributed approaches are effective and bring a controlled computational burden. Convergence performance is generally guaranteed.
However, the aforementioned distributed algorithms work only in line-of-sight (LOS) scenarios. In some practical situations, such as in forests, cities, and indoor places, most connections between sensors are non-line-of-sight (NLOS) because of the obstacles in the direct paths of signal propagation. This phenomenon can severely degrade localization accuracy if the NLOS impact is not taken into consideration. In general, two approaches are applied to alleviate NLOS propagation in localization. The first approach distinguishes LOS and NLOS connections via prior information [30][31][32]. This approach needs other measurement techniques to provide NLOS information, such as Direction of Arrival (DOA). Hence, its estimation model is less suitable for cooperative localization methods based on range-based techniques. The second approach modifies the original model by weighting and adding constraints of range-measurement errors [33,34]. This kind of approach can be applied to more general scenarios. Some representative work is listed as follows. The work in Reference [35] is "the first one to address NLOS node localization in WSN". According to three different scenarios, the authors provided three modified models of the ML estimator by adding the upper and lower bounds of range measurements. These models were solved in standard SDP problems with heuristic estimation. The work in Reference [34] proposed an ESDP method and added constraints to improve robustness. In Reference [36], the authors introduced NLOS bias parameters to the estimation model for the source-node locations and turned the model into a SDP problem. In Reference [37], the authors proposed a three-block ADMM algorithm based on the model in Reference [36]. However, this method does not decrease the size of the original problem, and the authors also mention that "it is hard to distribute the calculation" because it involves matrix projection. The work in Reference [38] proposed a distributed algorithm based on the Huber estimation of Reference [39], but it did not provide theoretical proof for the convergence property.

Contributions
Up to now, most NLOS mitigation techniques applied to WSN cooperative localization still utilize centralized optimum algorithms. In this paper, we propose a parallel distributed algorithm based on a tight relaxation technique to both decrease computational complexity and improve accuracy in NLOS environments.
First, we propose a modified convex model based on projection relaxation, which relaxes the original nonconvex problem into its convex envelope. This model can be applied to tough situations where NLOS connections are unidentifiable in all range measurements. According to the bounds of range measurements, we formulated the problem in the form of projection distances. Then, we relaxed this formulation into the projection on convex sets by using Cauchy-Schwartz inequality. Compared with SDP, this approach improves the decomposition property of the convex model and the accuracy of the estimation results.
Second, we developed a parallel distributed algorithm to solve the convex model. We designed a consensus form to decompose the large-scale problem into numerous local subproblems and provide the relevant theoretical proof. The proposed consensus form is more suitable for an ADMM framework. It enables each node to solve each subproblem in parallel. We derived the concrete procedure of how to handle the problem in a parallel way. The distributed algorithm had much lower computational complexity than that in existing papers about NLOS cooperative localization.
Third, we further improved the convergence rate of local updates. The local subproblems are convex and nonquadratic differentiable, which makes them less appropriate for the Newton method and interior-point algorithm. Hence, with the guarantee of Lipschitz continuity, we propose an iterative algorithm to solve the untypical convex subproblems based on the proximal method.
The paper is organized as follows. Section 2 formulates the cooperative localization problem and the corresponding convex relaxation. Section 3 presents the consensus form and solves it in a distributed way. Section 4 derives the iterative method for local subproblems. In Section 5, the simulation results are reported. Section 6 concludes the paper.

Mathematic Model
The mathematic model of range-based cooperative localization is described as follows. As is shown in Figure 1, consider a sensor network consisting of N source sensors, of which locations x i ∈ R d , i = 1, 2, . . . , N are unknown, and M anchor sensors of which the location x k ∈ R d , k=N+1, . . . , N+M is known. d is the coordinate dimension. Source nodes are collected in the set S = {1, 2, . . . , N}, and anchors are collected in the A = {N + 1, N + 2, . . . , N + M} set. We denote the Euclidean distance between node i and j as d i,j , and the noisy range measurement as r i,j . In general, not every pair of nodes can communicate because the communication distance has an upper limit. We denoted this distance limit as r µ . The neighbor of node i is denoted as j ∈ N i if r i,j is available, i.e., N i = {j|r i,j ≤ r µ , ∀j ∈ S}, ∀i ∈ S ∪ A. We denoted the pairwise of source sensors as (i, j) ∈ Z S , and the pairwise between a source sensor and an anchor as (i, k) ∈ Z A . The distance between two nodes is defined as Considering the impact of NLOS propagation on distance estimation, we divided the range measurements into two sets. We used Z LOS (and, respectively, Z NLOS ) to denote the set of pairwise nodes in which connections between nodes are LOS (and, respectively, NLOS). Hence, the range measurements are defined as: where n i,j ∼ N(0, σ 2 i,j ) is the measurement noise following a zero-mean Gaussian distribution with variance σ 2 i,j , and ε i,j is the error of the NLOS measurement that is exponentially distributed with a mean parameter α i,j = α NLOS . The value of α NLOS depends on the NLOS propagation environment [40][41][42].
In most cases, prior information of NLOS connections, such as NLOS distribution, is unavailable. In this paper, our model was designed for such tough scenarios, i.e., NLOS connections being unidentifiable among all connections. First, since we could not distinguish which connections were NLOS, we assumed that all range measurements were NLOS. Model r i,j was modified as

Convex Relaxation
The method in Reference [35] provided two bounds of d i,j in NLOS environments. The upper bound is the lower bound is For the single-constraint case, the upper and lower bounds provide an annular feasible region. Heuristic localization estimation lies on the circle with a radius of (u i,j + l i,j )/2. Hence, the optimization problem of cooperative localization in an NLOS environment is modified as: This model takes NLOS impact into account and uses heuristic points (u i,j + l i,j )/2 to replace range measurements r i,j . In severe environments where most connections are NLOS, range measurements may contain much false information and, thus, cannot directly be used as distance metrics. In Model (6), the bounds can neutralize most NLOS errors in range measurements. However, Problem (6) is nonconvex. In this paper, we propose a tight relaxation method rather than SDP. We used projection relaxation to obtain the convex envelope of the original nonconvex problem. The theoretical derivation is given as follows.
First, we rewrote the fundamental part of Problem (6) in the form of a squared distance of the projection on a certain set by using Cauchy-Schwartz inequality.
The proof is given as follows: Then, function can be rewritten as where B = {b| b = d 0 } is a nonconvex set. Hence, the original problem can be written in a simpler form: where set B i,j is a spherical surface depending on the bounds, i.e., B i, A feasible approach is to find the convex envelope by relaxing constraint b = d 0 to b ≤ d 0 . Therefore, we obtained the convex form based on Function (8).
where C = {b| b ≤ d} is the convex hull of the non-convex set B. According to Equation (11), we could establish the convex model corresponding to Problem (10): where

Distributed Framework
In order to approach the solution of convex optimization Problem (12) in a distributed way, we built the consensus form to formulate Problem (12) and design a parallel distributed algorithm based on ADMM. Problem (12) cannot be directly applied to parallel ADMM because it involves matrix calculation. We built the consensus form to decompose large-scale Problem (12) into N + M subproblems with independent variables. Then, we could design the ADMM framework to solve these subproblems in a parallel way. In this section, we prove the feasibility of decomposing Problem (12) and provide the concrete procedure of the parallel distributed algorithm.
The consensus form could independently represent the connections between each sensor and its neighbors in the whole network by duplicating a new vector of each node j at neighbor node i.
when we put the elements of neighbor set N i in corresponding vector N i in an ascending order, the relation between (z i ) j and x l is l = N i (j − 1). N i (j) means the jth element of N i . In other words, (z i ) j is the duplication of (j − 1)th neighbor node at node i. Therefore, Problem (12) has an equivalent form: wherex i is the linear function of x.
We denoted local objective function as Then, we could rewrite Problem (14) in another, simpler form: where R z is a linear space satisfying Matrix A i indicates the linear relation between local variable z i and global variable x , of which the form is where e l is the lth unit-row vector in R N+M , I d is an identity matrix of order d, and ⊗ is a Kronecker product. By introducing an indicator function G z (z) of R z , Problem (14) can be rewritten in a general consensus form as where z is the shorthand notation for z T 1 , z T 2 , . . . , z T N+M T , and A = [A 1 ; A 2 ; . . . ; A N+M ].
Before deriving the ADMM framework, we needed to build the augmented Lagrangian function of Problem (18) where λ is a vector consisting of Lagrangian multipliers, and c is the given penalty parameter. Indicator function G z (z) can be written as the sum of numerous subfunctions, namely, by decomposing linear space R z as the Cartesian product of (N + M) closed convex sets where T k is a matrix satisfying Thus, we can conclude that the objective function, constraints, and penalty term are separable for each sensor. Distributed optimization could be obtained by applying the ADMM method to Problem (18). The local updates are handled in an alternate way that can be separately optimized for variables z, x, and λ, as follows, for ∀i ∈ S ∪ A whereλ i = λ i /c is the scaled dual variable. z i −updates andλ i −updates can be independently carried out in parallel for each node i. The concrete implementation steps are shown as follows.
(1) Initialize variables z 0 i , x 0 ,λ 0 i for all nodes. (2) Each node updates its local variables z t+1 i according to the first step of Problem (21), corresponding to the action If we put indicate function G z i (z i ) in the constraints, Problem (22) can be written in an equivalent form: where γ i are the slack variables working as the quadratic penalty in the local cost function. Obviously, Problem (23) is untypically convex. The gradient of F i (z i ) is discontinuous. Hence, classical methods such as the interior point method and the Newton method are not appropriate to solve these problems. In Section 4, we further derive an iterative algorithm to approach the solution of Problem (23). This iterative algorithm has the same convergence rate with the interior point method.
(3) Broadcast local variables z t+1 i to all nodes of neighbor sets. (4) Each node computes variablesx t+1 i according to Problem (24) is actually easy to solve by averaging all components of z k+1 i of which the local indices correspond to global index l. The x-update step has a simpler solution: where k l is the number of local variable components corresponding to global components z l and A summary of the proposed parallel algorithm is illustrated as Algorithm 1 and Figure 2: Algorithm 1 ADMM_Projection (ADMM_P) algorithm for Problem (18).
Input: number of anchors M, location of anchors x k ∈ R 2 , k ∈ A, number of source sensors N, range measurements {r i,j , (i, j) ∈ Z NLOS ∪ Z LOS } ; Output: (z i ) 1 , i ∈ S 1.Initialization: t = 0; (1)Initialize local location of sensors x 0 and set dual variables λ 0 i to zero; (2) According to x 0 i and λ 0 i , compute the initial values:  (26).

Update the global variables
Update via (23) and Section 4 x Initialization Allocate corresponding value to each local variable

Solution to the Subproblem
Section 3 implies that the ADMM_P algorithm is an iterative method. Reference [27] guarantees that our algorithm converges to the same minimum as that in the centralized framework because its original counterpart (12) is convex. The algorithm summary indicates that the total computational complexity of ADMM_P almost relies on the local z−update (22). To satisfy the requirements on processing speed, we needed to design an effective algorithm that had both a fast convergence rate and low complexity. The Newton method and interior point method have these properties, but they were designed for smooth optimization problems [43]. Unfortunately, Problem (22) is a nonquadratic differential. Hence, we designed an iterative method based on the accelerated proximal algorithm, which is an analogous tool for nonsmooth optimization problems [21,44]. Given the structure of WSNs, we derived the proximal algorithm for source nodes and anchors, respectively.

Lipschitz Constant
The proximal algorithm needs the guarantee of Lipschitz continuity. First, we prove the premise in this subsection. To simplify notations, we defined the functions: where C i is Cartesian product of balls C i,l , l ∈ N i , and H i has the form of Function g(z i ) in Functions (27) is convex and differentiable, and its gradient is Obviously, this gradient is Lipschitz-continuous, and we could obtain its Lipschitz constant as follows where k g ≥ 1 is a constant.

Source Nodes
For source nodes, local Subproblem (22) has the equivalent simpler form min h(z i ) + g(z i ) the proximal mapping of function h(z i ) = 1 where µ = 1/L g is the step size, H † i is the pseudoinverse of H i , and P C i (H i z i ) is the orthogonal projection of point H i z i onto the convex set C i , i.e., Hence, we can obtain the accelerated proximal algorithm for local Subproblems (31), as follows, for k ≥ 1: This iterative procedure stops when the stopping criterion is satisfied. µ k is the step size. Here, we used the fixed step size µ k = µ = 1/L g . We used the estimation results of the last loop as the initial point of this subupdate, i.e., ξ i as the local estimation result.

Anchor Nodes
For anchor nodes (z i ) 1 = x i , i ∈ A, Problem (22) becomes an equality-constrained optimization problem for i ≥ N + 1, Our approach solving this problem was to eliminate the equality constraint via the equation where e d i , i = 1 or 2 is the ith unit vector in R d , I N i d is an identity matrix with order N i d, and 0 i , i = 1 or 2 is a zero matrix.
With Equation (36), we form eliminated optimization problem which is an unconstrained problem with variable β i . Then, we could derive the solution of equality-constrained Problem (35) according to Equations (34). Introducing the affine structure does not change the convexity and continuity of the functions. The gradient of theĝ(β i ) function is still Lipschitz-continuous. The gradient ofĝ(β i ) is Hence, the new Lipschitz constant could be obtained by the following inequality: By using the conclusion of Equation (32) to Equation (34),the accelerated proximal algorithm for the local subproblem of the anchor nodes is presented as follows:

Numerical Simulations
In this section, we present the simulation results to demonstrate the strengths and weaknesses of the ADMM_P algorithm. We considered the worst case, where NLOS connections are unidentifiable and the probability of these NLOS connections is up to 95%. We preliminarily considered a network that had 10 anchor nodes and 40 source nodes randomly distributed in a [−50 m, 50 m] × [−50 m, 50 m] square. Noisy range measurements were available when distances r i,j ≤ 30 m. The coordinate dimension is d = 2. Range measurements were generated according to Equations (1) and (2),where NLOS error ε i,j is exponentially distributed with mean parameter α NLOS . For a noisy environment, we assumed that the measurement noise was independent and identically distributed, i.e., σ i,j = σ n .
It is worth noting that ADMM_P is an iterative algorithm working in the parallel distributed framework, and that there are few authoritative publications about distributed algorithms for NLOS localization. Therefore, we only analyzed the convergence property for our algorithm. In the convergence analysis, we provide the estimation results of SDPH of Reference [35] as references.
All simulations were carried out in MATLAB. Local convex Problems (22) were handled via the MATLAB Optimization Toolbox with f mincon solver. Algorithm performance was evaluated via Cumulative Distribution Function (CDF), Root Mean Square Error (RMSE), and objective Function Value (FVAL): where p i is the true position, and (z i ) 1 is the estimated location of node i.

Performance Comparison
In this part, we verify the superiority of our method in the aspects of localization accuracy and computational complexity, respectively. We fixed the number of source nodes and anchors as N = 40, M = 10. The mean parameter of NLOS propagation is α NLOS = 3. Sensor nodes (including anchors) were randomly distributed in the square of [−50 m, 50 m] × [−50 m, 50 m]. To show the simulation results more precisely, we compared the ADMM_P algorithm with three other state-of-art methods in the same simulation environment. Figures 3 and 4 show the simulation results of the four available methods. As is shown, ADMM_P had higher accuracy than other three methods. In Figure 3, the CDF curve of ADMM_P is in the leftmost side of the three other methods. This means that the estimation results of ADMM_P were more concentrated around lower errors. Especially at point Error = 2.5 m, the CDF of our method improved by up to 50% compared with the ADMM_SF, which did not utilize any NLOS mitigation techniques. In Figure 4, we further analyzed the performance of the four methods in varying environments with different σ n and α NLOS , respectively. The simulation results show that ADMM_P has higher accuracy, and performance is more stable in different environments.

Complexity Analysis
We compared computational complexity in Table 1, in which we used N c to represent the connective numbers of pair nodes (i, j) ∈ Z S ∪ Z A . We could find that ADMM_P is more suitable for practical situations because it has much lower computational complexity and smaller problem size. The parallel distributed framework plays a key role in decreasing complexity.
The SDPH method in Reference [35] and the SDP method in Reference [36] were implemented in the centralized framework. The problem sizes of SDPH and SDP were, respectively, dN + 1/2 * N(N + 1) and dN + 1/2N(N + 1) + N 2 , both depending on the whole scale of the WSN N. The EDM method in Reference [37] applied a three-block ADMM, but problem size also depends on N. According to Table 1, we can see that the problem size of the ADMM_P algorithm was determined by the decomposed local minimization and could be reduced to d(N i + 1) +N i , whereN i represents the number of neighbors. Because of the limitation of r µ , the size of the neighbor set was much smaller than the whole scale of the WSN, i.e.,N i N. This conclusion can theoretically prove that our method can drastically reduce the precessing burden, no matter the problem size or computational complexity.
Our method decomposes the original large-scale problem into N + M small-scale subproblems via consensus form, and handles these subproblems in a parallel way. Hence, when compared with the centralized algorithm, our method can dramatically decrease computational complexity. In Sections 3 and 4, we could find that the total complexity of the proposed algorithm almost depends on solving the subproblems. Section 4 provides an iterative method to solve the subproblems, which are nonquadratic differentials. In this iterative method, with the guarantee of objective-function convexity, the optimal value is finite and attained at z * i . The derived method is known to converge at a rate of O(1/k 2 ) [44]: Subupdates require very simple operations, no matter if they are at anchor nodes or at source nodes. The inversion of matrix H i and H i Q i is easy to precompute because H i and Q i only consist of identity matrices and unit vectors. Projection on set C i can be decomposed into projections on sets C i,j , j ∈ N i , which involves very limited dimension d ≤ 3 and is easy to compute. Therefore, the local update given by Equations (34) and (41) has much lower complexity than the three other algorithms. However, because every node shares the processing burden and plays the same role in the parallel framework, the communication cost of each node is higher than that of centralized frameworks.

ADMM_P Properties in Different Noisy Environments
In this scenario, we fixed the number of anchors and the distribution parameter of NLOS errors to show localization accuracy and convergence property in different noisy environments. The mean parameter of NLOS error was fixed as α NLOS = 3. Ten anchors were randomly deployed in the region of 100 × 100 m. As is shown in Figure 5a,b, we present the simulation results with σ 2 n = 0.01 × 100 and σ 2 n = 0.1 × 100. We found that ADMM_P mitigates the effect of NLOS range measurements and provides better estimation accuracy. SDPH is solved with the interior-point method under a centralized framework, so its simulation results are only provided as accuracy references in Figure 5. From the figure, we can see that ADMM_P has good convergence performance in different noisy environments and higher accuracy than SDPH. In Figure 5b, we further analyzed the convergence performance of ADMM_P. Given that the local objective functions are based on heuristic solutions,the function value cannot converge to an extremely low level. Hence, the curve tendencies in Figure 5b still demonstrate a stable convergence of ADMM_P.

ADMM_P Properties in Different NLOS Environments
In this scenario, we fix the number of anchors and the noise level to study ADMM_P performance with different NLOS error distributions. The measurement noise variance was fixed as σ 2 n = 0.01 × 100. Other simulation settings were the same as those in Scenario B. In Figure 6a,b, we compared the performance of ADMM_P with α NLOS = 3, 4, 5. As we can see, although RMSE slightly increased as mean parameter α NLOS grew, the accuracy of ADMM_P is still better than SDPH. This also indicates that the proposed algorithm effectively mitigates NLOS error. Simulation results in Figure 6b verify the stable convergence performance of ADMM_P as varying NLOS error distribution.

Anchor Placement
In this scenario, we deployed the anchors at the perimeter of the region. We set σ 2 n = 0.01 × 100 and α NLOS = 3. We compared localization accuracy when the anchor placement was random or fixed. In Reference [13], the author proposed that an appropriate anchor-placement deign could further improve localization accuracy. Hence, we placed the anchors around the perimeter of the region to avoid getting crowded points, so that the higher dimensional projection of anchors could contain more sensor nodes. The corresponding simulation results are shown in Figure 7b, from which we can see that the fixed anchor placement improves the localization performance of the ADMM_P algorithm in the same NLOS environment. Furthermore, in the scenario of fixed anchors, the accuracy of the ADMM_P algorithm is higher than that of the SDPH algorithm, even if the former employs the less anchors.

Varying Numbers of Anchors
Simulations in the last subsection show that a predesigned placement of anchors has a positive influence on localization performance. In this scenario, we evaluated the performance of the ADMM_P algorithm by varying the number of anchors. Taking the simulation results of Figure 7 into consideration, we fixed 10 anchors on the region boundary and randomly deployed the remaining anchors inside the region. Other simulation settings were the same as those in the last subsection. Node distribution is shown in Figure 8a, and the corresponding simulations are shown in Figure 8b, which shows that proper anchor placement and an increase in anchor number both lead to an improvement of localization performance.

Conclusions
In this paper, we proposed an efficient distributed algorithm for NLOS cooperative localization. We employed tight convex relaxation for the heuristic model, and approached it in a parallel ADMM framework to offer powerful computation ability. The distributed framework can drastically reduce the computational burden since it decomposes a large-scale problem into numerous small-scale problems, which similarly retains computational time when the number of sensors grows. Furthermore, we derived the accelerated proximal gradient method to improve the convergence rate of local subproblems. Simulation results demonstrate that our method efficiently alleviates the influence of NLOS propagation and has fairly good convergence performance.
In future works, we plan to study the level of influence that step size and penalty parameters could have on the performance of local subproblems and an ADMM framework. How to establish the mathematic model to find the optimal parameters needs further study as well. Moreover, the lower and upper bounds of range measurements provide a good feasible region, which is nevertheless nonconvex. We are interested in finding better convex relaxation rather than utilizing a heuristic solution as the optimal solution.