TW-TOA Based Positioning in the Presence of Clock Imperfections

# TW-TOA Based Positioning in the Presence of Clock Imperfections

Sinan Gezici, Senior Member, IEEE, and Erik G. Ström, Senior Member, IEEE
M. R. Gholami is with the ACCESS Linnaeus Center, Electrical Engineering, KTH–Royal Institute of Technology, SE-100 44 Stockholm, Sweden (e-mail: mohrg@kth.se)S. Gezici is with the Department of Electrical and Electronics Engineering, Bilkent University, Ankara 06800, Turkey (e-mail: gezici@ee.bilkent.edu.tr).E. G. Ström is with the Division of Communication Systems, Information Theory, and Antennas, Department of Signals and Systems, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden (e-mail: erik.strom@chalmers.se).This work was supported in part by the Swedish Research Council (contract no. 2007-6363) and in part by the European Commission in the framework of the FP7 Network of Excellence in Wireless COMmunications COMmunications # (contract no. 318306). Part of this work was presented in IEEE ICASSP 2013 [1].
###### Abstract

This paper studies the positioning problem based on two-way time-of-arrival (TW-TOA) measurements in asynchronous wireless sensor networks. Since the optimal estimator for this problem involves difficult nonconvex optimization, we propose two suboptimal estimators based on squared-range least squares and least absolute mean of residual errors. The former approach is formulated as a general trust region subproblem which can be solved exactly under mild conditions. The latter approach is formulated as a difference of convex functions programming (DCP), which can be solved using a concave-convex procedure. Simulation results illustrate the high performance of the proposed techniques, especially for the DCP approach.

Index Terms– Positioning, two-way time-of-arrival (TW-TOA), clock imperfection, convex optimization, trust region subproblems, concave-convex procedure.

## I Introduction

Location aware services are becoming vital requirements for many wireless systems. Due to some drawbacks of using GPS receivers at wireless nodes for some scenarios, self-position recovery has been proposed as an alternative approach and extensively investigated in the literature [2, 3, 4, 5, 6, 7]. Positioning based on range estimates between nodes is a popular technique in the literature. For synchronous networks, the time-of-arrival (TOA) technique provides a good estimate of the distance between two nodes for reasonable signal-to-noise ratios. A huge number of algorithms have been proposed in the literature to address the positioning problem based on range measurements, e.g., the maximum likelihood estimator [3], linear least-squares [8, 9, 10], squared-range least squares [11], projection onto convex sets [12, 13, 14], and convex relaxation techniques [15, 16, 17].

In asynchronous networks, the range estimate based on the TOA is highly sensitive to the clock imperfections. Therefore, the positioning accuracy can be considerably degraded in the presence of clock imperfection. In particular for an affine model describing the clock behavior, the accuracy of the positioning techniques based on the TOA measurements is affected by non-ideal clock offset and clock skew. The clock of a target node can be synchronized with a reference time (clock) using a synchronization technique using, e.g., the MAC layer time stamp exchanging, e.g., see,[18, 19, 20, 21, 22] and references therein. Motivated by pairwise synchronization techniques, the authors in [23] formulate a joint synchronization and positioning problem in the MAC layer. If the major delay is the fixed delay due to propagation through the radio channel, the joint position and timing estimation technique works well. The authors in [24] study the positioning problem in the presence of clock imperfection only due to clock offset. Considering the effect of imperfect clock on distance estimates in the physical layer, the authors in [25] investigate the positioning problem using time-difference-of-arrival (TDOA) in the presence of clock imperfections. The TDOA technique effectively removes the clock offset, but still suffers from the clock skew.

Another popular approach for estimating the distance between sensor nodes is to use a so-called two-way time-of-arrival (TW-TOA) or time-of-flight based technique, which is an elegant approach in removing the effect of the clock offset on range measurements [10]. Range estimates obtained via TW-TOA are affected by the clock skew and a processing delay called the turn-around time [26]. A number of researchers have tackled the positioning problem or distance estimation based on TW-TOA in fully or partially asynchronous networks [27, 28]. To improve the range estimation via TW-TOA, an effective technique based on a new clock counter scheme is proposed in [29]. The authors in [28] study the positioning problem in the presence of clock imperfections for a TW-TOA based technique and propose a linear least squares based approach to solve the problem. The proposed approaches work well in some scenarios, e.g., when there is a sufficient number of reference nodes at known positions. In general, the previously proposed approaches need modifications to be effectively applied to the positioning problem in which the clock skew and turn around times are also unknown. In addition, for practical applications the proposed algorithms may not be robust against outliers and non-line-of-sight errors.

In this study, we consider the positioning of a single target node based on TW-TOA measurements in the presence of clock imperfections. In this approach, a target node transmits a signal to a reference node located at a known position and the reference node responds to the received signal after an unknown turn-around time delay. As it is common in the literature, we assume that the reference node measures the turn-around time by a loop back test and transmits the estimate to the target node [30, 29]. The target node then computes the round-trip delay based on an estimate of the turn-around time. Assuming an affine model for the clock of the oscillator, it is observed that the range estimation using the TW-TOA measurement is affected by an unknown clock skew of the target node. Modeling the measurement errors as Gaussian random variables, we obtain the optimal estimator to find the clock skew, and the location of the target node, and the turn-around times for the reference nodes. The optimal estimator poses a high dimensional optimization problem and needs more than one distance estimate for every link to provide good estimates of the unknown parameters. We, then, omit the effect of the turn-around times using a linear transformation and consequently obtain a near-optimal estimator to find the location and clock skew of the target node. Both the optimal and near-optimal estimators for the positioning problem considered in this study are nonconvex and difficult to solve. Using some approximations, we obtain two suboptimal estimators. In the first approach, we apply a similar technique as considered in [11] for synchronized networks to the asynchronous scenario. Namely, we consider the squared-range least-squares approach and formulate the problem as a general trust region subproblem, which can be solved exactly under mild conditions. In the second approach, we minimize the residual errors based on the norm and arrive at a nonconvex problem in the form of the difference of convex functions programming (DCP). We then employ a concave-convex procedure to solve the problem. Note that the latter approach is robust against outliers. Simulation results indicate the high performance of the proposed techniques, especially the DCP.

In summary the main contributions of this study for TW-TOA based positioning in the presence of clock imperfections are:

• the MLE and an approximate MLE to find the location and clock skew of the target node;

• a suboptimal estimator based on general trust region subproblem for squared-range measurements;

• a suboptimal estimator formulated as the DCP that can be solved using a concave-convex procedure.

The remainder of the paper is organized as follows. Section II explains the signal model considered in this study. In Section III and IV the localization algorithms are studied. Complexity analyses of different approaches are discussed in Section V. Simulation results are presented in Section VI. Finally, Section VII makes come concluding remarks.

Notation: The following notations are used in this paper. Lowercase Latin/Greek letters, e.g., , denote scalar values and bold lowercase Latin/Greek letters represent vectors. Matrices are shown by bold uppercase Latin/Greek letters. is the by identity matrix. The operator is used to denote the expectation of a random variable (or vector). The norm of a vector is denoted by . The is a diagonal matrix with diagonal elements . For two matrices and , means is positive semidefinite. denotes the gradient of at . The set of all -vector with positive components are denoted by . We use to denote the Kronecker product.

## Ii System Model

Consider a two dimensional network111The generalization to a three dimensional network is straightforward, but is not explored in this study. with reference (anchor) nodes located at known positions , . Suppose that one target node is placed at unknown position . We assume that the target node estimates the distance to a reference node by performing a TW-TOA measurement. We further assume that the clock value of an imperfect clock follows an affine relation with the true (global) time  [31, 18, 19, 21, 22]. That is, the clock value of reference node is

 gi(t)≜wit+θi (1)

where is the skew and is the offset associated with the th node clock. Note that a perfectly synchronized clock has and . In practice, is a number close to 1. For convenience, we denote the target clock as , where .

A TW-TOA measurement between the target node and the th reference node is carried out as follows: (a) the target sends a message to the reference node at (global) time , (b) the message arrives at the reference node at time , (c) the reference node sends a return message at time , and (d) the return message arrives at the target node at time . Clearly, , where is the propagation speed and is the distance between the target and th reference node. Moreover, , where is the turn-around time in the th reference node, which is assumed to be fixed during the positioning process. The TWO-TOA measurement is computed in the target local clock as

 zki=12[g(tk4,i)−g(tk1,i)+nki]=wdic+wTi2+nki2, k=1,2,…,K (2)

where is the TW-TOA measurement error, which we model it as a zero-mean Gaussian with standard deviation , i.e., , and as the number of the TW-TOA measurements during the positioning process.

The unknown parameter either might be extremely small and can be neglected [20] or it needs to be estimated. One way to deal with the unknown parameter is to jointly estimate it along with the location of the target node [32]. It can also be estimated by reference node using a loop back test and is sent back to the target node [29]. In this study, we consider the latter approach. The estimate of is

 ~Tki=gi(tk3,i)−gi(tk2,i)+ϵki=wiTi+ϵki, k=1,2,…,K (3)

where we model the estimation error as .

In the sequel, we assume that the reference nodes are synchronized with a reference clock, e.g., via a GPS signal. Therefore and we can write , where

 ^Tki≜Ti+ϵki. (4)

We now combine (4) with (2) and obtain

 zki−w^Tki2=wdic+nki2−wϵki2. (5)

As mentioned, the approximation is good in the (reasonable) case when the reference nodes are equipped with accurate clocks.

In the following sections, we use the input data to obtain the optimal estimator based on models (2)-(4) or suboptimal estimators according to (5). The parameters and , , are considered as unknown nuisance parameters, while , , and are assumed to be known for .

## Iii Maximum Likelihood Estimator

We define the measurement vector

 m≜⎡⎢ ⎢ ⎢ ⎢ ⎢⎣m1m2⋮mK⎤⎥ ⎥ ⎥ ⎥ ⎥⎦, (6)

where

 mk≜[zk1zk2⋯zkN^Tk1^Tk2⋯^TkN]T. (7)

To obtain the MLE for joint estimation of the position and clock skew of the target node, the following optimization problem needs to be solved [33]:

 (8)

where is the probability density function (pdf) of vector indexed by the vector and . Since the TOA measurement errors are assumed to be independent and identically distributed random variables, the pdf of can be calculated from (2) and (4) as

 p(m;w,ta,x)=K∏k=1N∏i=1√2πσ2iexp(−2(zki−wTi/2−wd(x,ai)/c)2σ2i)√12πγ2iexp⎛⎝−(^Tki−Ti)22γ2i⎞⎠. (9)

Then, the MLE is obtained as

 [^xT ^w  ^ta]T =arg maxw∈R+;x∈R2;Ti∈R+p(m;w,ta,x) =arg minx∈R2;w∈R+;Ti∈R+K∑k=1N∑i=12σ2i(zki−wTi2−wdic)2+(^Tki−Ti)2γ2i. (10)

For the MLE formulated in (III) there are unknowns. Therefore, for low numbers of messages , the MLE problem can be ill-posed. To alleviate the difficulty for solving the optimal MLE in (III), we here investigate another approximate MLE based on the model obtained in (5).

Let us define a new set of measurements as

 yki≜zki−w^Tki2=wdic+nki2−wϵki2, i=1,2,…,N, k=1,2,…,K. (11)

Now, we collect in a vector . Next, we compute the pdf of as

 p(y;w,x) =K∏k=1N∏i=1√2π(σ2i+w2γ2i)exp⎛⎜⎝−2(zki−wdi/c−w^Tki/2)2(σ2i+w2γ2i)⎞⎟⎠. (12)

We then find an approximate MLE (AMLE)222We call the MLE in (III) as AMLE because it is based on the transformed set of measurements in (III) instead of the original measurements in (6). as

 [^xT ^w]T =arg maxw∈R+;x∈R2p(z;w,x) =argminx∈R2;w∈R+K∑k=1N∑i=12(σ2i+w2γ2i)(zki−w^Tki2−wdic)2+ln(σ2i+w2γ2i). (13)

It is observed that the search domain in the AMLE in (III) is limited to the location and the clock skew , thus a lower dimensional search compared to that of the MLE in (III).

It is also noted that both the MLE and AMLE formulations in (III) and (III) pose difficult global optimization problems. To avoid the drawbacks in solving these problems, we propose two suboptimal estimators in the next section.

## Iv Proposed techniques

In this section, we propose two techniques based on squared-range least squares and norm minimization of residuals. First, we divide both sides of (5) by (we safely assume that ) and express the model as

 zkiα−^Tki2=dic+nki2α−ϵki2, i=1,2,…,N, k=1,2,…,K (14)

where .

In the following, the model in (14) is employed in order to derive the proposed suboptimal estimators.

### Iv-a Squared-Range measurement Least Squares (SQ-LS)

In this section, we assume that the measurement noise is small compared to . Then, taking the square of both sides of (14) and dropping the small terms yield

 (zkiα)2+(^Tki)24−zki^Tkiα ≃1c2(xTx−2aTix+∥ai∥22)+νki, (15)

where . Now, we apply a weighted least squares criterion to the model in (15) and obtain the following minimization problem:

 minimizex∈R2; α∈R+ K∑k=1N∑i=11d2i(α2σ2i+γ2i)(1c2xTx−2c2aTix−(zki)2α2+zki^Tkiα+1c2∥ai∥22−(^Tki)24)2. (16)

The problem in (16) can be expressed as a quadratic programming problem:

 minimizey  ∥W1/2(Ay−b)∥22 subject to yTDy+2fTy=0 (17)

where matrices , , and and vectors , , and are defined as

 W=IK⊗diag(1d21(α2σ21+γ21),…,1d2N(α2σ2N+γ2N)), A≜⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣1c2−2c2aT1−(z11)2z11^T11⋮⋮⋮⋮1c2−2c2aTN−(z1N)2z1N^T1N⋮⋮⋮⋮1c2−2c2aT1−(zK1)2zK1^TK1⋮⋮⋮⋮1c2−2c2aTN−(zKN)2zKN^TKN⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, b≜⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣−1c2∥a1∥22+(^T11)24⋮−1c2∥aN∥22+(^T1N)24⋮−1c2∥a1∥22+(^TK1)24⋮−1c2∥aN∥22+(^TKN)24⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, D≜IK⊗diag(0,1,1,0,1), f≜1K⊗[−12  0  0 −12  0]T, y≜[∥x∥22 xT α2 α]T. (18)

with denoting a column vector of ones.

The problem in (IV-A) minimizes a quadratic function over a quadratic constraint. This type of problems are called a generalized trust region subproblem [34] and can be solved exactly. It has also been known that the general trust region subproblem has zero duality gap and the optimal solution can be extracted from the dual solution [35, 36, 34]. A necessary and sufficient condition for to be optimal in (IV-A) is that there exists a such that [35]

 (ATWA+μD)y∗=(ATWb−μf), (y∗)TDy∗+2fTy∗=0, (ATWA+μD)⪰0. (19)

Under the conditions considered in (IV-A), the solution to the problem of (IV-A) is given by

 y(μ)=(ATWA+μD)−1(ATWb−μf). (20)

In such a situation to find , we simply replace (20) into constraint , i.e.,

 (21)

where the interval consists of all such that . The interval is given by [11]

 I=(−1/μ1,∞), (22)

with representing the largest eigenvalue of  [34]. In summary, the solution to (IV-A) is obtained as follows:

• Use a bisection search to find a root of , say . Note that is a strictly decreasing function with respect to [34].

• Replace into (20) to obtain .

• Estimate the unknown parameters as and , with denoting the th to the th elements of vector .

Note that since the weighting matrix depends on the unknown distance and , we first replace with the identity matrix and find an estimate of the location and as described above. Then, we reconstruct the distance considering the estimate as and form a new approximate weighting matrix. This approach can be continued for a number of iterations, however, as we have observed through simulations, after two updatings, the estimation accuracy improves only slightly via additional iterations.

Another estimator based on a linear least squares (LLS) approach obtained in Appendix A can be alternatively applied to the model in (15). Note that the algorithm derived in Appendix A is similar to the one proposed in [28], except the correction technique introduced in this study. As will be observed in the simulations section the proposed approach in this section, i.e., SQ-LS, has better performance than the LLS approach, especially for low number of reference nodes.

### Iv-B A concave-convex procedure (CCCP)

In this section, we take the norm minimization of residual errors into account and propose a technique to solve the positioning problem. Namely, based on (14), we consider the following norm minimization problem:

 minimizex∈R2; α∈R+ ∥r∥1 (23)

where with . Note that for high signal-to-noise ratios (low standard deviations of noise), the and minimization approaches have similar performance [37]. Moreover, the norm minimization in (23) is robust against outliers [37]. The optimization problem in (23) can be written (in the epigraph form) as [37, 38, 39]

 minimizex∈R2;α∈R+;t∈RN K∑k=1N∑i=1tki subject to zkiα−^Tki/2−di/c≤tki zkiα−^Tki/2−di/c≥−tki. (24)

The nonconvex problem in (IV-B) is reminiscent of a well-known nonconvex problem called difference of convex functions programming (DCP) [40]. The general form of DCP is as follows:

 minimizex f0(x)−g0(x) subject to wi(x)−gi(x)≤0, i=1,…,M (25)

where and are both smooth convex functions for . A method to solve (IV-B) is to sequentially solve the problem. That is, we first approximate the concave function with a convex one by an affine approximation. Let us consider a point in the domain of the problem in (IV-B), linearize the concave function around and write the optimization problem in (IV-B) as

 minimizex f0(x)−g0(xj)−▽g0(xj)T(x−xj) subject to  wi(x)−gi(xj)−▽gi(xj)T(x−xj)≤0. (26)

The convex problem in (IV-B) can now be solved efficiently. Denote the solution of (IV-B) as . Next we go for further improving the solution by convexifing (IV-B) for the new point similar to the procedure employed for . This sequential programming procedure, called concave-convex programming (CCCP), continues for a number of iterations. The convergence of the CCCP to a stationary point has been shown in the literature, e.g., [40, 41] and references therein. Fig. 1 shows an example of the CCCP approach for a DCP, where and are as follows:

 f(x)=∥x−a1∥2+∥x−a2∥2+∥x−a3∥2, g(x)=∥2x−a4∥2, (27)

with , and . In this figure, the original nonconvex problem is transferred to a convex problem using a linear approximation of the nonconvex problem around .

Applying the CCCP technique to the problem in (IV-B), we get the following optimization problem:

 minimizex∈R2;α∈R+;t∈RN K∑k=1N∑i=1tki subject to zkiα−hTi,jx−bji,k−tki≤0 1c∥x−ai∥2−zkiα+^Tki2−tki≤0 (28)

where and . The optimization problem in (IV-B), which is called second order cone programming (SOCP), can be efficiently solved. We call the corresponding CCCP as CCCP-SOCP.

Note that if is not differentiable at , we can replace the term by a subgradient333Let be a nonempty set in . A vector is a subgradient of a function at if for all  [42]. of at .

## V Complexity analysis

In this section, we study the complexity of the proposed techniques in terms of floating point operations (flops) and also running time in Matlab. We compare the complexity of the MLE, LLS, SQ-LS, and CCCP-SOCP. To compute the complexity of the MLE, we assume that a good initial point is available, and an iterative algorithm such as the Gauss-Newton (GN) method converges to the global minimum after a number of iterations. Of course, finding a good initial point for the MLE is a challenging problem and this study also aims to tackle it. For the problem at hand the complexity of the MLE for every Newton step can be computed as . For the AMLE, the complexity for every Newton step can be computed as . The corresponding LLS needs an order of to implement. For the SQ-LS, we need to use a bisection search to solve (21), which is the most complex part of the algorithm. Suppose the bisection search takes steps, then the total cost of the the proposed approach can be approximated as . In the simulations, we have observed that the bisection search algorithm usually takes 20 to 30 iterations to find the optimal value of . Note that we need to run the LLS and SQ-LS twice. Thus the corresponding complexities are increased by a factor of two. Finally, the complexity of the CCCP-SOCP can be computed as follows. Consider a general form of the SOCP problem as

 minimizex∈Rn cTx subject to ∥Aix+bi∥2≤cTix+di, i=1,…,m, ∥x∥2≤R (29)

where , and . Note that the constraint on the norm of ensures the strong convexity of the centering problem in the barrier approach [37]. The worst-case complexity of the problem in (V) can be computed as  [43], where is an accuracy tolerance in solving the problem.

The complexity of the CCCP-SOCP for every estimate can now be approximated as

 O((KN)3.5log1/ϵ).

As mentioned before we need to solve the problem in steps, hence the total cost is . As we observe, a small number of updatings, usually three, , is enough to obtain the solution. Table  I summarizes the complexity of the different approaches.

We have also measured the average running time of different algorithms for a network consisting of 6 reference nodes as considered in Section VI. In the simulations, we set and . The algorithms have been implemented in Matlab on a MacBook Pro (Processor 2.3 GHz Intel Core i7, Memory 8 GB 1600 MHz DDR3). The MLEs are implemented by Matlab function fminsearch initialized with the true values of the target position, the clock skew, and turn-around times. The CCCP-SOCP is implemented by the CVX toolbox [44]. We use three updatings to get an estimate.

We run the algorithms for 500 realizations of the network and compute the average running time in ms. The results are shown in Table II. Considering the complexity analysis and average running time in Tables I and II, respectively, we can conclude that the proposed approach has reasonable complexity and running time. Although CCCP-SOCP takes a longer amount of time than MLE, it does not need a good initial point. While for the MLE with an arbitrary initial point, the algorithm may converge to a local minimum resulting in a large position error. As we will see in the next section, the CCCP-SOCP outperforms both the LLS and SQ-LS approaches in terms of the root-mean-squared-error.

## Vi Numerical results

In this section, we evaluate the performance of the proposed approaches through computer simulations. We consider a m by m area and a number of reference nodes that are located at fixed positions , and . In the simulations, we pick the first reference nodes, i.e., . One target node is randomly distributed inside the area. The turn-around time is set to ms. The clock skew is assumed to be unknown and is set to PPM, i.e., . Such a value for clock skew is common for a practical oscillator. For example for a center carrier frequency MHZ and a frequency offset KHZ, we can write

 T=1fc∓Δf≈1fc(1±Δffc)=1fc(1±0.0001). (30)

We compare the proposed techniques (CCCP-SOCP and SQ-LS) with the MLE and AMLE in (III) and (III), respectively, (which are implemented by Matlab function fminsearch [45] initialized with the true values of the target location, turn-around times, and clock skew), the LLS derived in Appendix A, and the Cramér-Rao lower bound (CRLB) as derived in Appendix B. In the simulations, we assume that . We randomly initialize the CCCP-SOCP inside the network and we also set . To simulate the range measurements and estimates of turn around times, we use models (4) and (2), respectively. To implement the bisection search, we consider an interval defined by and and investigate if the zero crossing of in (21) occurs in the interval. To check if the solution lies in the interval, we simply check the sign of at and . No change in sign means that the solution lies outside of the current interval. For initialization, we set and . If the solution of is not found in the interval, we change the interval as and . If the solution lies in an interval, we bisect the interval and investigate which subinterval contains the solution.

Fig. 2 shows the root-mean-squared-errors (RMSEs) of location estimates for different approaches for various numbers of reference nodes. In the simulations, we set . It is observed that the proposed approach, CCCP-SOCP, achieves good performance very close to the optimal estimator MLE and the CRLB, especially for low number of reference nodes and high signal-to-noise ratios. From the figure, it is noted that the SQ-LS performance is worse than CCCP-SOCP, but better than LLS, especially for a low number of reference nodes. As the number of reference nodes increases, the LLS and SQ-LS show similar performance. It is also observed that for large numbers of reference nodes, the least squares based approach shows better performance compared to the CCCP-SOCP approach for the low standard deviation of noise.

Next, we study the effects of NLOS measurements on the performance of estimators. We assume that a range measurement can be affected by NLOS errors with probability 0.2. For every NLOS measurement, we add a uniform noise to the measurements as follows:

 zki=w(dic+Ti2)+uki+nki2 (31)

where we assume that . The uniform distribution is commonly used to model NLOS error, e.g., [46, 47, 13].

Fig. 3 depicts the performance of different approaches in NLOS conditions for . It is observed that the CCCP-SOCP achieves high performance compared to the other approaches, especially for low standard deviations of noise, and it is robust against outliers as expected. For small , the dominant perturbation is outlier disturbance and consequently the MLE derived in this study is not optimal, explaining why the MLE is worse than the CCCP-SOCP approach. For large standard deviations of noise, which indicates the Gaussian measurement noise is dominant, the CCCP-SOCP seems to outperform the MLE. This can be explained by the fact that the MLE is only guaranteed to be asymptotically optimal, i.e., for low noise standard deviation or large number of measurements. Note that we have employed the MLE computed in (III) and (III) to study their robustness against NLOS conditions. It may be possible to derive an MLE to deal with NLOS measurements if the distribution of outliers is known.

Finally, we study the convergence of CCCP-SOCP through simulations. Fig. 4 depicts the convergence speed of the proposed approach for 50 random initializations. In the simulations, we set . For every estimate given by CCCP-SOCP, we compute the residual , where is given by (23). It is observed that the CCCP-SOCP approach converges very fast, approximately in three sequential updatings.

## Vii Conclusions

In this paper, we have studied TW-TOA based positioning in an asynchronous network. Since the optimal ML estimator is highly nonconvex and difficult to solve, we have obtained two efficient suboptimal estimators for the problem under some approximations and conditions. The first method is based on squared-range least squares which can be solved exactly under mild conditions. The second approach is derived by replacing the norm minimization of residuals by an norm minimization, which in turn can be formulated as difference of convex programming (DCP). We have used a concave-convex procedure to solve the resulting DCP. Simulation results show the high performance of the proposed techniques, especially the DCP approach. It has also been observed through simulations that the DCP approach is robust against NLOS errors.

## Appendix A Linear least squares (LLS)

In this section we obtain an LLS estimator similar to [9, 27]. We consider the following linear model (originated from (15)):

 b=Ay+ν, (32)

where , , , and are given in (IV-A). We assume that has full column rank. A necessary condition for this is that .

The unconstrained weighted least squares solution to (32) is given by [33]

 ^y=(ATWA)−1ATWb. (33)

where is as in (IV-A). The covariance matrix of can be computed as

 C^y=(ATWA)−1. (34)

Note that for a large network, matrix can be ill-conditioned [27]. Then, we can use a regularization technique to resolve the drawback in the least squares solution [37, 27].

We can further improve the location estimate by applying a correction technique similar to [9, 27]. We consider the following relations:

 [y]1=∥x∥22+ξ1, [y]4=α2+ξ4, [y]2=x1+ξ2, [y]3=x2+ξ3, [y]5=α+ξ5, (35)

where is the estimation error. Assuming small estimation errors, we take the squares of both sides of last three equations in (A) and obtain the following expressions:

 [y]22≃x21+2x1ξ2, [y]23≃x22+2x2ξ3, [y]25≃α2+2αξ5. (36)

Based on (A) and (A), we obtain a linear model as

 h=Bθ+Pξ, (37)

where

 B=⎡⎢ ⎢ ⎢⎣111100010001⎤⎥ ⎥ ⎥⎦,P=⎡⎢ ⎢ ⎢⎣1000102x1000002x2000002α0⎤⎥ ⎥ ⎥⎦ h=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣[y]1+[y]4[y]22[y]23[y]24⎤⎥ ⎥ ⎥ ⎥ ⎥⎦,θ=[x21 x22 α2]T. (38)

The least squares solution to (37) is given by

 ^θ=(BTC−1^θB)−1BTC−1^θh, (39)

where the covariance matrix can be computed as

 C^θ=PC^yPT. (40)

To compute the matrix , we use the estimate of obtained in (33) instead of unknown vector .

Finally the location estimate can be obtained as

 ~xi=sgn([y]i+1)√|[^θ]i|,i=1,2, (41)

where denotes the signum function defined as

 sgn(x)={1if x≥0;−1if x<0. (42)

The covariance matrix of the estimate in (41) can be obtained similar to [27].

## Appendix B Cramér-Rao Lower Bound (CRLB)

Considering the measurement vector in (6) with mean and covariance matrix where

 μ=[f(d1c+T12)…f(dNc+TN2) T1…TN]T, C=diag(σ214,…,σ2N4,γ21,…,γ2N), (43)

the elements of the Fisher information matrix can be computed as [33, Ch. 3]

 Jnm=[J]nm=[∂μK∂ψn]TC−1K[∂μK∂ψm],n,m=1,2,…,N+3, (44)

where

 ψn=⎧⎨⎩xn,if n=1,2w,if n=3Tn,if n>3. (45)

From (B), can be obtained as follows:

 [∂μK∂ψn]=1K⊗[∂μ1∂ψn … ∂μN∂ψn]T,n=1,2,…,N+3, (46)

where

 ∂μi∂ψn=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩wxn−a1,ncd(ai,x),if n=1,2, i≤Ndic+Ti2,if n=3, i≤N0,if n=1,2,or 3, i>Nw2,if n>3, i≤N1,if n>3, i>N. (47)

After some calculations, the entries of the Fisher information matrix can be computed as follows:

 J11=4Kw2N∑i=1(x1−ai,1σicd(ai,x)