Robust time-of-arrival localization via ADMM

Source localization (SL) refers to estimating the position of a signal-emitting target from measurements collected using multiple spatially separated sensors [1]. Owing to its great significance to a lot of location-based applications, such as emergency assistance [2], asset tracking [3], Internet of Things [4], and radar [5], the SL problem has received much attention in the literature over the past decades [6]. Depending on the measurements being used, methods for SL can be roughly divided into two groups: range- and direction-based. The category of range-based SL techniques, as the term range suggests, involves the use of a variety of distance-related location data, including time-of-arrival (TOA) [7, 8, 9, 10, 11, 12, 13, 14, 15, 16], time-difference-of-arrival [9, 12, 14, 17, 18, 19, 20], time-sum-of-arrival [12, 21, 22, 23, 24, 25, 26, 27], and received signal strength (RSS) [15, 16, 28], to determine the source position. Among these schemes, the time-based ones exploit the signal propagation time, while their RSS counterparts harness the signal power strength to acquire distance-type sensor observations. This dichotomy allows time-based SL solutions to excel in precision at the cost of more stringent synchronization requirements, while RSS-based approaches prioritize simplicity and inexpensiveness in implementation. In contrast, direction-based SL techniques utilize the source’s bearings relative to the sensors to fulfill the positioning task [17, 29]. They, like methods built upon the RSS measurement model, can dispense with clock synchronization. However, direction-based positioning mandates the installation of an antenna array at each sensor for angle-of-arrival estimation, resulting in a somewhat higher hardware cost compared to its range-based rivals.

TOA-based SL takes center stage in this contribution, given its stature as one of the most widely adopted and intensively examined distance-based positioning schemes. Such a choice is further justified by the widespread availability of smartphones as well as other mobile devices that possess communication and information processing capabilities, most of which can obtain TOAs without encountering significant technical challenges [6]. Under the Gaussian noise assumption, least squares (LS) has undoubtedly been the very thing one needs to deliver statistically meaningful estimation results. In this vein, a wealth of LS techniques have been proposed in the literature, including algebraic explicit and exact solutions [7, 8, 9], convex programming approaches [10, 11], and local optimization algorithms [12, 13, 14], to name but a few.

Statistical model mismatch will appear in real-world scenarios where the Gaussian noise assumption may be violated due to the presence of unreliable sensor data [30], greatly deteriorating the performance of the LS methodology [31]. To reduce the negative impact of non-Gaussian measurement errors on positioning accuracy, researchers focusing on TOA-based SL have explored the strategies of joint estimation of location coordinates and a balancing parameter [26, 27, 32, 33], worst-case LS criterion [33, 34, 35], and robust statistics [16, 28, 36, 37, 38, 39, 40, 41]. Here, we are particularly interested in the statistical robustification of the TOA-based LS position estimator. Countermeasures of this type have been attracting increasing attention in recent years, mainly because of their relatively low prior knowledge requirements and low complexity of implementation.

Consider a single-source localization system deployed in $H$-dimensional space. It consists of an emitting source that is to be located, with unknown position $\bm{x}\in\mathbb{R}^{H}$, and $L$ sensors with signal receiving capabilities and known positions $\{\bm{x}_{i}\in\mathbb{R}^{H}|i=1,...,L\}$. Synchronization is guaranteed not only among the sensors themselves but also between the source and each sensor, so that the TOAs of the emitted signal at all $L$ sensors can be estimated. Specifically, the TOA-based range measurements are modeled as [1]

where ${\|\cdot\|}_{2}$ denotes the $\ell_{2}$-norm and $e_{i}$ accounts for the observation uncertainty associated with the $i$th source-sensor path. Under the assumption that $\{e_{i}\}$ are uncorrelated zero-mean Gaussian processes with variances $\{\sigma_{i}^{2}\}$ and in a maximum likelihood (ML) sense [42], the TOA-based SL problem is mathematically stated as [12]

As mentioned earlier, the $\ell_{2}$-norm based estimator (2) is vulnerable to outlying (unreliable) data, which manifest themselves as range measurements immersed in non-Gaussian, potentially bias-like errors. Common causes of them in the localization context include non-line-of-sight and multipath propagation of signals, interference, sensor malfunction, and malicious attacks [5, 30, 31, 43]. To achieve resistance against the occurrence of outliers in such adverse situations, (2) was statistically robustified in [16, 28, 38, 39, 40, 41] by substituting the $\ell_{2}$ loss $(\cdot)^{2}$ with a certain fitting error measure less sensitive to biased $\{r_{i}\}$, i.e.,

where frequently utilized options for $f(\cdot)$ encompass the $\ell_{1}$ loss $|\cdot|$ [16, 39], the Huber loss:

with radius $R_{i}>0$ [16, 38], the Welsch loss $f_{\textup{W}}(\cdot)=\exp\left(-\tfrac{(\cdot)^{2}}{2\sigma_{\textup{W}}^{2}}\right)$ with parameter $\sigma_{\textup{W}}$ [40], and the $\ell_{p}$ loss ${|\cdot|}^{p}$ for $1\leq p<2$ [28, 41]. Below, let us provide a concise overview of several well-established optimization solutions available in the literature for addressing (3). They serve as competitors of our proposal, as we will discuss in Sections 3 and 4. By forming the associated epigraph problem [66]:

one may directly apply second-order cone programming (SOCP), a classic convex approximation technique, or difference-of-convex programming (DCP) [16] to address the $\ell_{1}$-minimization version of (3):

The epigraph form (1) plays a crucial role as the basis for deriving (3) and (3), which represent the SOCP-based solution and DCP-based solution for (6), respectively [16]. Smoothed approximations through the composition of natural logarithm and hyperbolic cosine functions [39], on the other hand, render the Lagrange-type neurodynamic optimization framework of Lagrange programming neural network (LPNN) [67] applicable to

where

with parameter $\nu>0$ is also known as (a.k.a.) the smoothed $\ell_{1}$ loss [68]. The idea of turning instead to the Huber convex underestimator for (1) based on the composite loss $f_{R_{i}}([-(\cdot)]^{+})$ was conceptualized in [38], where $[\cdot]^{+}=\max(\cdot,0)$ is the ramp function. Nonetheless, the emphasis of [38] is more on the extension of

to sensor network node localization, rather than single-source positioning. In [40], the challenging range-based Welsch loss minimization problem:

was approximated by

using the squared ranges, after which a half-quadratic (HQ) optimization algorithm was devised. The authors of [41] converted the $\ell_{p}$-minimization problem:

into an iteratively reweighted LS (IRLS) framework and applied sum-product message passing (MP) to solve the $\bm{x}$-update subproblem within each IRLS iteration in closed form.

While scaling poorly with the problem dimension is a well-known limitation of convex programming solutions whose implementation relies on interior-point methods, numerical realization of the neurodynamic approach in [39] typically takes several thousands of iterations to reach equilibrium [69]. Despite the low-complexity advantage of the HQ algorithm in [40], the squared-range formulation (11) it adopts inherently suffers from impaired statistical efficiency. The MP-assisted IRLS [41] is computationally attractive, yet the lack of a theoretical foundation to support convergence of it may affect the soundness of the technique. With these concerns in mind, we are motivated to explore new avenues for coping with (3), especially in a way less computationally demanding but more statistically efficient and theoretically complete compared to the existing solutions.

In this contribution, we employ the alternating direction method of multipliers (ADMM) as our algorithmic strategy. ADMM is an iterative optimization solver which works by breaking down the original problem into several more easily solvable subproblems and, technically, it can be seen as an effort to combine the advantages of dual ascent and augmented Lagrangian methods [44]. Having originated in the 1970s with roots dating back to the 1950s [44, 45], this simple but powerful algorithm has witnessed extensive use in the past few decades across various statistical and machine learning problems, including face recognition [46], subspace segmentation [47], portfolio optimization [48], image classification [49], and low-rank matrix recovery [25, 50, 51]. The key reasoning behind the adoption of ADMM in diverse applications of recent interest is that, while initially crafted for convex optimization, numerous extensions of it to the nonconvex setting still allow for lightweight yet reliable implementation and comprehensive theoretical analysis.

Our main contribution in this work is to develop an ADMM approach for tackling the generalized robust TOA positioning problem (3). As the versatility of the cost function $f(\cdot)$ will only be embodied in one of the subproblems that amounts to computing the proximal operator of it, we are able to adapt $f(\cdot)$ to specific noise environments with ease. In an analytical manner, we prove that if the presented ADMM is convergent, the limit point to which it converges will satisfy the Karush-Kuhn-Tucker (KKT) conditions for the equivalent constrained reformulation of (3). We also show the possibility of insuring theoretical convergence of the ADMM, under several additional assumptions about $f(\cdot)$. Numerical examples are included to corroborate the positioning accuracy and computational simplicity superiorities of our solution over its competitors. It should be pointed out that compared with the ADMM in [12], our scheme holds independent significance in the following two aspects. First, the algorithm in [12] deals with the non-outlier-resistant formulation (2), whereas we aim at solving (3) with a robustness-conferring fitting error measure. Second, we provide analytical convergence results of the presented ADMM. The authors of [12], in contrast, incorrectly assumed the nonconvex constrained optimization reformulation of (2) to be convex, consequently omitting to probe into the convergence of their method under nonconvexity.

The rest of this contribution is organized as follows. The equivalent constrained reformulation of (3) and the ADMM solution to it are described in Section 2. The complexity and convergence properties of the proposed ADMM are discussed in detail in Section 3. Performance evaluations are conducted in Section 4. Ultimately, Section 5 concludes the article.

Let us begin by converting (3) into

where $\bm{d}=\left[d_{1},...,d_{L}\right]\in\mathbb{R}^{L}$ is a dummy vector of decision variables for the source-sensor distances, independent of $\bm{x}$.

We further rewrite (13) as

by introducing an auxiliary vector $\bm{\beta}=\big{[}\bm{\beta}_{1}^{T},...,\bm{\beta}_{L}^{T}\big{]}^{T}\in\mathbb{R}^{HL}$ into the characterization of range information [12]. Here, $\bm{\beta}_{i}\in\mathbb{R}^{H}$ is a unit vector that indicates the source’s direction with respect to (w.r.t.) the $i$th sensor, a.k.a. the direction vector of arrival [52].

The transformation from (13) to (14) bears some resemblance to the conversion from (13) to a formulation, where the Euclidean distances are re-expressed in a quadratic form. This conversion, in fact, is a widely adopted preprocessing step in the domain of range-based localization [10, 23], aimed at eliminating the cumbersome Euclidean norm terms in the initial formulation and simplifying the associated constrained optimization problem. In order to achieve an equivalent reformulation during this simplification process, it is crucial to consider the quadratic form of the Euclidean distance constraints along with the nonnegativity constraints on the distance-representing variables [10, 23]. In other words, it would be mathematically less rigorous to casually disregard these nonnegativity constraints without the backing of appropriate theoretical justifications, as demonstrated by the Appendix in [10] and Theorem 1 in [23]. We therefore follow this tradition in our transformation from (13) to (14).

For (14), we construct the following augmented Lagrangian with constraints:

where $\bm{\lambda}=\left[\bm{\lambda}_{1}^{T},...,\bm{\lambda}_{L}^{T}\right]^{T}\in\mathbb{R}^{HL}$ is a vector containing the Lagrange multipliers for (14a) and $\rho>0$ the augmented Lagrangian parameter. Splitting the primal variables into two parts, the ADMM for solving (14) consists of the following iterative steps:

where the iteration index is indicated by $(\cdot)^{(k)}$. To be specific, (16a) and (16b) sequentially minimize the augmented Lagrangian (2) w.r.t. the primal variables, after which (16c) updates the dual variables via a gradient ascent with step size $\rho$. We note that there are only two (namely, $\bm{x}$ and $(\bm{d},\bm{\beta})$) rather than three primal blocks in the ADMM governed by (16). By comparison, natural extension of the basic two-block ADMM to multi-block cases may result in divergence [53].

More detailed explanations for (16a) and (16b) are given as follows.

By ignoring the constant terms independent of $\bm{x}$, the subproblem (16a) can be simplified into an LS form

where

It has the following closed-form solution:

Similarly, another subproblem (16b) is re-expressed as

where

Exploiting the particular structure of (20), the optimal $\bm{\beta}$ can be obtained first as

Thus, (20) is reduced to

which is separable w.r.t. the partition of $\bm{d}$ into its $L$ elements and leads to the following $L$ subproblems:

With the help of the following proposition, we can simplify (24) to a certain degree to facilitate the calculations.

Proposition 1. Under mild assumptions about the robust loss and range observations that $f(\cdot)$ is an even function strictly increasing on the nonnegative semi-axis and $\{r_{i}\}$ are always nonnegative, coping with (24) will be equivalent to handling

Proof. See A.

It is worth noting that the assumptions made in Proposition 1 are actually rather common (see our justification in the following). In such cases (where Proposition 1 applies), solving (24) boils down to computing the proximal operator of $f(\cdot)$, namely,

complying with the definition of proximal mapping below.

Definition 1. The proximal mapping of a function $f:\mathbb{R}\rightarrow\mathbb{R}$ with parameter $\tau>0$ for any $b\in\mathbb{R}$ is [54]

Obviously, the computational simplicity of the proximal mapping procedure is crucial to the efficient update of $\bm{d}$ in each iteration of the proposed ADMM. In fact, the calculation of (27) can be done in a relatively unencumbered manner or in closed form for many choices of $f(\cdot)$ exhibiting outlier-resistance [54, 55, 56, 57]. As summarized in Table 1²²2For the purpose of completeness, we also include the non-robust case of $p=2$, in which the $\ell_{p}$-minimization estimator will coincide with (2), providing optimum estimation performance as long as the disturbances $\{e_{i}\}$ are independent and identically distributed (i.i.d.) Gaussian processes., this article restricts the scope of discussion to two typical options: (i) the $\ell_{p}$ loss ${|\cdot|}^{p}$ for $1\leq p<2$ and (ii) the Huber function $f_{R_{i}}(\cdot)$. The two corresponding instantiations of (3) are then (12) and

respectively. Our justification for them is given as follows. The $\ell_{p}$-minimization criterion with $1\leq p<2$ is known to show a considerable degree of robustness against outliers in a wide range of adverse environments [78, 75, 41, 55]. Closely resembling its $\ell_{1}$ (resp. $\ell_{2}$) counterpart for large (resp. small) fitting errors, the Huber loss minimization scheme offers another means of allowing somewhat controllability therebetween [58]. Most importantly, all these cost functions are compatible with our assumptions in Proposition 1.

We note that $a^{+}$ and $a^{-}$ in Table 1 are the unique positive root of $\tfrac{a-b}{\tau}+pa^{p-1}=0$ for $a\in[0,b]$ (when $b\geq 0$) and unique negative root of $\tfrac{a-b}{\tau}-p(-a)^{p-1}=0$ for $a\in[b,0]$ (when $b<0$), respectively, both of which can be conveniently found using simple root-finding algorithms (e.g., bisection or secant, among others). While these root-finding algorithms typically have only a constant time complexity of $\mathcal{O}(1)$ [55], the computational efficiency of the overall methodology may not always be optimal (e.g., they will lead to nested iterations in our $\ell_{p}$-norm-based ADMM use case). This discrepancy is evident in our CPU time comparisons, where we found that the $\ell_{p}$-norm-based ADMM method consistently required more runtime than the Huber-loss-based one. However, it is noteworthy that this increase in computational complexity does not significantly hinder the practicality of the corresponding ADMM algorithm. ADMM adopting the $\ell_{p}$ loss for $1<p<2$, even with the increased complexity, remains far more computationally efficient than several existing convex approximation approaches. Furthermore, in most cases, the Huber variant of our ADMM algorithm, which does not involve any iterative root-finding steps thanks to the closed-form proximal operator for Huber loss, already delivers decent performance. This means that there is usually limited need for invoking the $\ell_{p}$-norm version of the proposed ADMM technique (with iterative root-finding requirements) in typical TOA SL scenarios.

The steps of ADMM for dealing with the robust TOA SL problem are summarized in Algorithm 1, where $\delta$ is a user-defined tolerance to terminate the iterations.

The complexity and convergence properties of the proposed ADMM are discussed in detail in this section.

From Algorithm 1, it is not difficult to find that the complexity of the ADMM is dominated by that of the $\bm{d}$-update steps, in which proximal mapping calculations are involved. Therefore, the total complexity of Algorithm 1 will be $\mathcal{O}(N_{\textup{ADMM}}L)$ if the loss function taken from Table 1 admits closed-form $\bm{d}$-updates and $\mathcal{O}(N_{\textup{ADMM}}KL)$ otherwise, where $N_{\textup{ADMM}}$ represents the iteration number of the ADMM and $K$ the number of searches involved in the root-finding process at each ADMM iteration. Empirically, a few steps of the bisection method applied in each $\bm{d}$-update and a few tens to hundreds of ADMM iterations will already be sufficient to yield an accurate source location estimate.

By comparison, based on interior-point methods, solving the SOCP problem derived from (1) [16]:

implementing $N_{\textup{CCCP}}$ concave-convex procedure (CCCP) iterations in tackling the DCP problem (1) [16]:

and realizing the ramp function based Huber convex underestimator (9) [38] result in $\mathcal{O}(L^{3.5})$, $\mathcal{O}(N_{\textup{CCCP}}L^{3.5})$, and $\mathcal{O}(L^{3.5})$ complexity, respectively [19]. The LPNN [39], the HQ technique [40], and the MP-assisted IRLS algorithm [41, 28] are three other representative TOA positioning schemes from the literature that adopt the strategy of statistical robustification as well. They lead to $\mathcal{O}(N_{\textup{LPNN}}L)$, $\mathcal{O}(N_{\textup{HQ}}KL)$, and $\mathcal{O}(N_{\textup{IRLS}}L)$ complexity, respectively, where $N_{\textup{LPNN}}$ is the number of iterations in the numerical implementation of LPNN, $N_{\textup{HQ}}$ the number of steps needed for the HQ algorithm to converge, and $N_{\textup{IRLS}}$ that of the IRLS iterations. $N_{\textup{HQ}}$ typically takes a value of several tens, as does $K$ and $N_{\textup{IRLS}}$, while $N_{\textup{LPNN}}$ is generally in the range of several thousands [40, 41, 69]. It turns out that the proposed ADMM with loss functions permitting closed-form proximal computations, the HQ algorithm in [40], and the IRLS in [41] are most computationally efficient.

With an appropriately selected value of $\rho$, Algorithm 1 is observed to be always convergent in our simulations. However, many existing theoretical analyses of ADMM convergence in nonconvex scenarios, such as [59], may not be directly applicable to our context of robust TOA positioning due to two primary reasons. First, the optimization paradigms (of which convergence analyses were carried out) in many of these established analytical works differ from the constrained reformulation of the robust TOA positioning problem we address in our research. For instance, in [59], any constraints on the optimization variables that are not linear constraints are treated as indicator functions and incorporated into the objective function of the optimization framework. In our problem (14), nevertheless, we explicitly impose unit-length constraints on the direction-indicating variables and nonnegativity constraints on the distance-representing variables, and even the equality constraint functions involved with (14a) are not affine as per their definitions. Such discrepancies, in turn, result in distinct variable-update strategies between the ADMM algorithms considered in the relevant existing studies and our research. We note that similar dilemmas have been faced by many practitioners, despite empirically sound ADMM use cases in their respective contributions [12, 37, 55, 60]. In what follows, we will present our own analytic proofs for the convergence of Algorithm 1 in lieu of pinning all hopes on the generally applicable results from the literature.

Establishing the equivalence between (14) and

the following proposition points out that (14b) can be disregarded and we may shift our focus to the simplification (31) instead in the convergence analysis:

Proposition 2. Under the same assumptions as in Proposition 1, the formulations (14) and (31) are equivalent to one another.

Proof. The proof of Proposition 2 is similar to that of Proposition 1. Please see B for the details.

Next, we derive the following theorem for the optimality of tuples produced by Algorithm 1:

Theorem 1. Let $\big{\{}\bm{x}^{(k)},\bm{d}^{(k)},\bm{\beta}^{(k)},\bm{\lambda}^{(k)}\big{\}}_{k=1,...}$ be the tuples of primal and dual variables generated by Algorithm 1. If

then the limit $(\bm{x}^{\star},\bm{d}^{\star},\bm{\beta}^{\star},\bm{\lambda}^{\star})$ will satisfy the KKT conditions (a.k.a. the first-order necessary conditions [61]) for (31)³³3In cases where $f(\cdot)$ is non-differentiable at the origin, the condition (33b) should be substituted with the inclusion: $\bm{0}_{L}\in\hat{\partial}_{\bm{d}}(\bm{x}^{\star},\bm{d}^{\star},\bm{\beta}^{\star},\bm{\lambda}^{\star})$ [62], where $\hat{\partial}_{\bm{d}}(\cdot)$ denotes the generalized gradient of a function at $\bm{d}$ that takes into consideration the subdifferential calculus [63]. For simplicity’s sake, we confine our analyses here to the special case with a differentiable $f(\cdot)$. Nevertheless, it is worth pointing out that with slight modifications, the results can be readily extended to the scenarios in the presence of non-differentiability.:

where

is the associated Lagrangian of (31).

Proof. See C.

Remark 1. Theorem 1 does not guarantee the convergence of the proposed ADMM. Instead, it only reveals the optimality of solution when the algorithm converges. To further improve the analytical completeness, one may borrow the idea from [60] to relax the formulation (31) somewhat by introducing additional auxiliary variables and an extra quadratic penalty, whereby rigorously proving the convergence of the generated sequence is made possible. This will, however, give rise to degradation of the estimator’s performance [60] (i.e., the incompleteness of analytical convergence results can be seen as an acceptable trade-off for higher estimation accuracy). For avoiding such performance-impairing approximations to the source localization formulation, one would have to limit the scope of choice of $f(\cdot)$, as elaborated in the following statements.

Theorem 2. If $f(r_{i}-d_{i})$ is convex w.r.t. $d_{i},\forall i\in\{1,\ldots,L\}$, and for some constant $\bar{M}$, the gradient of $f(r_{i}-d_{i})$ satisfies that $\lvert\lvert\nabla_{d_{i}}f(r_{i}-d_{i})|_{d_{i}=d_{i,1}}\rvert+\lvert\nabla_{d_{i}}f(r_{i}-d_{i})|_{d_{i}=d_{i,2}}\rvert\rvert\leq\bar{M}{\lvert d_{i,1}-d_{i,2}\rvert}$, the sequence:

will be monotonically nonincreasing as the iterations (16a)–(16c) proceed, for a sufficiently large $\rho$.

Proof. See D.

Theorem 3. With the assumptions made about $f(\cdot)$ in Proposition 1 and Theorem 2, and an additional assumption that the gradient of $f(\cdot)$ at $d_{i}$ is Lipschitz continuous with a constant $M_{1}$, $\mathcal{L}_{\rho}\big{(}\bm{x}^{(k)},\bm{d}^{(k)},\bm{\beta}^{(k)},\bm{\lambda}^{(k)}\big{)}$ is bounded from below for sufficiently large $\rho$.

Proof. See E.

Corollary 1. $\big{\{}\mathcal{L}_{\rho}\big{(}\bm{x}^{(k)},\bm{d}^{(k)},\bm{\beta}^{(k)},\bm{\lambda}^{(k)}\big{)}\big{\}}_{k=1,...}$ is convergent under the above assumptions about $f(\cdot)$ and $\rho$.

Proof. This corollary is established via Theorems 2 and 3. $\square$

Equipped with these promising tools, we are finally enabled to present the following theorem, which asserts that Algorithm 1 can be theoretically guaranteed to converge with certain reasonable assumptions.

Theorem 4. The sequence $\big{\{}\bm{x}^{(k)},\bm{d}^{(k)},\bm{\beta}^{(k)},\bm{\lambda}^{(k)}\big{\}}_{k=1,...}$ is convergent, namely, (32) holds, under the same circumstances as in Corollary 1.

Proof. See F.

Corollary 2. Let us adopt the setting of Corollary 1. The solution obtained by Algorithm 1 corresponds to a KKT point of the nonconvex constrained optimization problem (31).

Proof. This corollary is established via Theorems 1 and 4. $\square$

In this section, we carry out computer simulations to evaluate the performance of Algorithm 1. Comparisons are made between the proposed ADMM and several existing TOA-based location estimators that are also built upon robust statistics. Table 2 gives a summary of these methods, and we note that the convex approximation techniques are all realized using the MATLAB CVX package [71]. To implement the LPNN [39], we invoke the MATLAB routine ode15s, a variable-step, variable-order solver relying on the formulas for numerical differentiation of orders 1 to 5 [72].

An SL system with $H=2$ and $L=8$ is considered. Unless otherwise mentioned, the locations of the source and sensors are randomly generated inside an origin-centered 20 m $\times$ 20 m square area, in each of the Monte Carlo (MC) runs. The positioning accuracy metric is the root-mean-square error (RMSE), defined as

where $N_{\textup{MC}}$ denotes the total number of MC runs and is fixed at $3000$ here, and $\tilde{\bm{x}}^{\{j\}}$ represents the estimate of the source position, ${\bm{x}}^{\{j\}}$, in the $j$th MC run. User-specified parameters of the ADMM are set as $\delta={10}^{-5}$ and $\rho=5$⁴⁴4It is essential to emphasize that the optimal selection of $\rho$ depends on both the problem scale and the specific setup of the loss function, i.e., careful consideration is required to ensure its optimality. For instance, when all variables in (14) are increased by a factor of 10, the last term of the augmented Lagrangian (2) will increase by a factor of 100, while the first two terms will increase by a factor of only 10 when using, e.g., the $\ell_{1}$ loss function. To maintain a balance between the two parts of the augmented Lagrangian, it is necessary to decrease $\rho$ to some extent, optimizing the numerical performance of ADMM. Nonetheless, for the sake of simplicity, our simulations exclusively employ a fixed $\rho$-setting with $\rho=5$, though this value may not be optimal in every scenario. We anticipate a thorough exploration of the impacts of $\rho$ to be a topic for future research, as explained in Section 5.. The tunable parameter $\sigma_{\textup{W}}$ of the Welsch loss is adaptively chosen according to the Silverman’s heuristic [73], in the same way as shown in [40]. The smoothing parameter $\nu$ is set to 0.1 to ensure feasibility of the ode15s solver.

Taking into account the presence of outliers, we follow [41, 75] to use the class of $\alpha$-stable distributions for modeling $\{e_{i}\}$ in (1). Stable processes are well known for their suitability to characterize skewness and heavy-tailedness. Except for a few special instances, members of the stable distribution family do not have an explicit expression of probability density function (PDF). Instead, their PDF $p(z)$ is implicitly described through the inverse Fourier transform of the characteristic function $\Phi(t;\alpha,\zeta,\gamma,\mu)$:

where the detailed analytical parameterization of $\Phi(t)$ can be found in [76]. As one may see, there are four parameters defining the family. The stability parameter, $0<\alpha\leq 2$, controls the tails of the distribution. Generally speaking, the smaller the value of $\alpha$, the heavier the tails and the more impulsive the random variable being modeled. $\mu$ determines the location. The skewness parameter, $-1\leq\zeta\leq 1$, is a measure of asymmetry. In the simplest case where $\zeta=0$, the distribution becomes symmetric about its mean (which is, $\mu$, when $\alpha>1$) and degenerates into the so-called symmetric $\alpha$-stable ($S\alpha S$) distribution. By contrast, the distribution is said to be right-skewed (resp. left-skewed) for $\zeta>0$ (resp. $\zeta<0$). $\gamma>0$ is the scale parameter, measuring to what extent the distribution spreads out (similar to the variance of the normal distribution). For illustration purposes, Fig. 1 plots the PDF functions for several representative choices of $\alpha$ and $\zeta$. As in our case, the stable-distributed range measurement noise $\{e_{i}\}$ is assumed to be i.i.d. with stability, skewness, scale, and location parameters $\alpha$, $\zeta=0$, $\gamma$, and $\mu=0$, respectively.

Because the variance of the stable distribution is undefined for $\alpha<2$, we introduce the concept of generalized signal-to-noise ratio (GSNR) from [41] to quantify the relative noise level:

Furthermore, the square root of the trace of the MC-approximated Cramér-Rao lower bound matrix (termed RCRLB) [77] is included in the comparison results to offer a benchmark for the accuracy of different robust position estimators in i.i.d. non-Gaussian noise. As it is in general difficult to work out a consistent schedule for adjusting the Huber radius under stably modeled non-Gaussianity, we simply follow [74] to assign a fixed value of 1 to $R_{i}$.

Starting off, Figs. 2 and 3 show the convergence behavior of $\ell_{p}$-ADMM ($p$ taking different values from 1 to 2) and Huber-ADMM in a single MC run for $\alpha=1.5$ and GSNR $=20$ dB. We should also note that for reproducibility, in this test $L=8$ sensors are evenly placed on the perimeter of the afore-defined square region and the source is deterministically deployed at $\bm{x}=[2,3]^{T}$ m. Both $\ell_{p}$-ADMM and Huber-ADMM are observed to rapidly decrease the objective function and converge to a point close to the true source position, within the first few tens of iterations. In our following simulation studies, the value of $p$ required by $\ell_{p}$-ADMM and $\ell_{p}$-IRLS will be set to the optimal one hinging on $\alpha$ [79], in a way analogous to [41, 75].

Fig. 4 depicts the RMSE versus GSNR $\in[17,25]$ dB while fixing $\alpha$ at $1.5$. It is seen that the two proposed ADMM approaches deliver the lowest level of positioning errors for all GSNRs. Especially, the RMSE performance of $\ell_{p}$-ADMM is the closest to the RCRLB benchmark (with a small gap of about half a meter). Among the six competitors, $\ell_{p}$-IRLS produces the best results. This reconfirms the findings reported in [41], that $\ell_{p}$-IRLS can be fairly statistically efficient in impulsive noise if $p$ is optly chosen. The direct nonlinear optimization based $\ell_{1}$-minimization techniques, $\ell_{1}$-DCP and $\ell_{1}$-LPNN, are slightly inferior to $\ell_{p}$-IRLS due to model mismatch but still capable of outperforming the remaining schemes, whose estimation performance looks relatively poorer. $\ell_{1}$-LPNN will not be comparable to $\ell_{1}$-DCP as the GSNR decreases from 25 dB, attributed to the smoothed approximations introduced in (7). Nonetheless, the former does not have tightness issues from which the traditional convex relaxation methods $\ell_{1}$-SOCP and Huber-CUE suffer, nor is it quite as statistically inefficient as Welsch-HQ.

Setting the GSNR to 20 dB and the impulsiveness-controlling parameter as $\alpha=1.5$, we plot the RMSE versus the number of sensors used for localization, $L\in[5,9]$, in Fig. 5. Again, the performance of the two ADMM solutions is optimal. Although none of the estimators can attain RCRLB, our algorithms are still the closest ones to it, particularly in the scenarios with a large number of sensors. Fig. 6 further plots the RMSE as a function of $\alpha\in[1.1,1.9]$ for $L=8$ and GSNR $=20$ dB. In this scenario, once more, $\ell_{p}$-ADMM stands out as superior compared to other approaches.

In Table 3, we list the per-sample average CPU time of the eight location estimators in the above simulations. It is seen that the running time can vary by several orders of magnitude. Huber-ADMM that allows for closed-form proximal mapping of $f(\cdot)$ and Welsch-HQ are the second and third fastest, respectively, taking only slightly more time than the state-of-the-art method $\ell_{p}$-IRLS. Because of the extra bisection procedure incorporated into the ADMM algorithm, $\ell_{p}$-ADMM may not be as computationally efficient as Huber-ADMM in actual CPU time comparisons. Nonetheless, its problem-solving process still takes less time than those of the traditional convex approximation techniques. In general, the simulation results agree with our complexity analysis in Section 3.

This article focused on the problem of outlier-resistant TOA SL. Based on the ADMM, we presented an iterative algorithm to handle the statistically robustified positioning formulation. Each iteration of our ADMM consists of two primal variable minimization steps and a dual variable update, all of which can be effortlessly implemented provided that the loss function relied on allows for convenient proximal computations. What is more, we proved that the ADMM will converge to a KKT point of the nonconvex constrained optimization problem under certain conditions, and verified that the LICQ holds at the point for nontrivial SL configurations. The superiority of the devised scheme over a number of robust statistical TOA SL methods in terms of positioning accuracy in impulsive noise and computational simplicity was demonstrated through simulations.

Future research endeavors could involve (i) deriving analytical convergence results under a milder setting of $\rho$, (ii) embarking on a more detailed exploration of the influences of $\rho$, and (iii) extending ADMM to the novel use case of robust direction-based SL with outliers.