TWTOA Based Positioning in the Presence of Clock Imperfections
Abstract
This paper studies the positioning problem based on twoway timeofarrival (TWTOA) measurements in asynchronous wireless sensor networks. Since the optimal estimator for this problem involves difficult nonconvex optimization, we propose two suboptimal estimators based on squaredrange 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 concaveconvex procedure. Simulation results illustrate the high performance of the proposed techniques, especially for the DCP approach.
Index Terms– Positioning, twoway timeofarrival (TWTOA), clock imperfection, convex optimization, trust region subproblems, concaveconvex 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, selfposition 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 timeofarrival (TOA) technique provides a good estimate of the distance between two nodes for reasonable signaltonoise 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 leastsquares [8, 9, 10], squaredrange 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 nonideal 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 timedifferenceofarrival (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 socalled twoway timeofarrival (TWTOA) or timeofflight based technique, which is an elegant approach in removing the effect of the clock offset on range measurements [10]. Range estimates obtained via TWTOA are affected by the clock skew and a processing delay called the turnaround time [26]. A number of researchers have tackled the positioning problem or distance estimation based on TWTOA in fully or partially asynchronous networks [27, 28]. To improve the range estimation via TWTOA, 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 TWTOA 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 nonlineofsight errors.
In this study, we consider the positioning of a single target node based on TWTOA 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 turnaround time delay. As it is common in the literature, we assume that the reference node measures the turnaround time by a loop back test and transmits the estimate to the target node [30, 29]. The target node then computes the roundtrip delay based on an estimate of the turnaround time. Assuming an affine model for the clock of the oscillator, it is observed that the range estimation using the TWTOA 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 turnaround 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 turnaround times using a linear transformation and consequently obtain a nearoptimal estimator to find the location and clock skew of the target node. Both the optimal and nearoptimal 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 squaredrange leastsquares 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 concaveconvex 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 TWTOA 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 squaredrange measurements;

a suboptimal estimator formulated as the DCP that can be solved using a concaveconvex 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 network^{1}^{1}1The 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 TWTOA 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
(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 TWTOA 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 turnaround time in the th reference node, which is assumed to be fixed during the positioning process. The TWOTOA measurement is computed in the target local clock as
(2) 
where is the TWTOA measurement error, which we model it as a zeromean Gaussian with standard deviation , i.e., , and as the number of the TWTOA 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
(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
(4) 
We now combine (4) with (2) and obtain
(5) 
As mentioned, the approximation is good in the (reasonable) case when the reference nodes are equipped with accurate clocks.
Iii Maximum Likelihood Estimator
We define the measurement vector
(6) 
where
(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
(9) 
Then, the MLE is obtained as
(10) 
For the MLE formulated in (III) there are unknowns. Therefore, for low numbers of messages , the MLE problem can be illposed. 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
(11) 
Now, we collect in a vector . Next, we compute the pdf of as
(12) 
We then find an approximate MLE (AMLE)^{2}^{2}2We 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
(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).
Iv Proposed techniques
In this section, we propose two techniques based on squaredrange least squares and norm minimization of residuals. First, we divide both sides of (5) by (we safely assume that ) and express the model as
(14) 
where .
In the following, the model in (14) is employed in order to derive the proposed suboptimal estimators.
Iva SquaredRange measurement Least Squares (SQLS)
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
(15) 
where . Now, we apply a weighted least squares criterion to the model in (15) and obtain the following minimization problem:
(16) 
The problem in (16) can be expressed as a quadratic programming problem:
(17) 
where matrices , , and and vectors , , and are defined as
(18) 
with denoting a column vector of ones.
The problem in (IVA) 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 (IVA) is that there exists a such that [35]
(19) 
Under the conditions considered in (IVA), the solution to the problem of (IVA) is given by
(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]
(22) 
with representing the largest eigenvalue of [34]. In summary, the solution to (IVA) is obtained as follows:
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., SQLS, has better performance than the LLS approach, especially for low number of reference nodes.
IvB A concaveconvex 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:
(23) 
where with . Note that for high signaltonoise 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]
subject to  
(24) 
The nonconvex problem in (IVB) is reminiscent of a wellknown nonconvex problem called difference of convex functions programming (DCP) [40]. The general form of DCP is as follows:
(25) 
where and are both smooth convex functions for . A method to solve (IVB) 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 (IVB), linearize the concave function around and write the optimization problem in (IVB) as
(26) 
The convex problem in (IVB) can now be solved efficiently. Denote the solution of (IVB) as . Next we go for further improving the solution by convexifing (IVB) for the new point similar to the procedure employed for . This sequential programming procedure, called concaveconvex 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:
(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 (IVB), we get the following optimization problem:
subject to  
(28) 
where and . The optimization problem in (IVB), which is called second order cone programming (SOCP), can be efficiently solved. We call the corresponding CCCP as CCCPSOCP.
Note that if is not differentiable at , we can replace the term by a subgradient^{3}^{3}3Let 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, SQLS, and CCCPSOCP. To compute the complexity of the MLE, we assume that a good initial point is available, and an iterative algorithm such as the GaussNewton (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 SQLS, 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 SQLS twice. Thus the corresponding complexities are increased by a factor of two. Finally, the complexity of the CCCPSOCP can be computed as follows. Consider a general form of the SOCP problem as
(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 worstcase complexity of the problem in (V) can be computed as [43], where is an accuracy tolerance in solving the problem.
The complexity of the CCCPSOCP for every estimate can now be approximated as
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.
Method  Complexity 

MLE (GN initialized with good initial points)  
AMLE (GN initialized with good initial points)  
CCCPSOCP  
LLS  
SQLS 
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 turnaround times. The CCCPSOCP 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 CCCPSOCP 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 CCCPSOCP outperforms both the LLS and SQLS approaches in terms of the rootmeansquarederror.
Method  Time (ms) 

AMLE (GN initialized with true values)  32 
MLE (GN initialized with true values)  317 
LLS  0.74 
SQLS  6.2 
CCCPSOCP  976 
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 turnaround 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
(30) 
We compare the proposed techniques (CCCPSOCP and SQLS) 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, turnaround times, and clock skew), the LLS derived in Appendix A, and the CramérRao lower bound (CRLB) as derived in Appendix B. In the simulations, we assume that . We randomly initialize the CCCPSOCP 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 rootmeansquarederrors (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, CCCPSOCP, achieves good performance very close to the optimal estimator MLE and the CRLB, especially for low number of reference nodes and high signaltonoise ratios. From the figure, it is noted that the SQLS performance is worse than CCCPSOCP, but better than LLS, especially for a low number of reference nodes. As the number of reference nodes increases, the LLS and SQLS 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 CCCPSOCP 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:
(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 CCCPSOCP 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 CCCPSOCP approach. For large standard deviations of noise, which indicates the Gaussian measurement noise is dominant, the CCCPSOCP 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 CCCPSOCP 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 CCCPSOCP, we compute the residual , where is given by (23). It is observed that the CCCPSOCP approach converges very fast, approximately in three sequential updatings.
Vii Conclusions
In this paper, we have studied TWTOA 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 squaredrange 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 concaveconvex 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)):
(32) 
where , , , and are given in (IVA). 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]
(33) 
where is as in (IVA). The covariance matrix of can be computed as
(34) 
Note that for a large network, matrix can be illconditioned [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:
(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:
(36) 
Based on (A) and (A), we obtain a linear model as
(37) 
where
(38) 
The least squares solution to (37) is given by
(39) 
where the covariance matrix can be computed as
(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
(41) 
where denotes the signum function defined as
(42) 
Appendix B CramérRao Lower Bound (CRLB)
Considering the measurement vector in (6) with mean and covariance matrix where
(43) 
the elements of the Fisher information matrix can be computed as [33, Ch. 3]
(44) 
where
(45) 
From (B), can be obtained as follows:
(46) 
where
(47) 
After some calculations, the entries of the Fisher information matrix can be computed as follows: