Abstract
An integral step in an ultra-wideband localization network installation is determining the positions of the fixed infrastructure nodes, the anchors. This process is time-consuming and usually requires specialized equipment. Additionally, it is difficult to achieve scalability, as any change or addition in the network requires a redetermination of the affected anchors. One can automate this process by utilizing the distance-measuring capabilities of the network infrastructure and employing a distributed position estimation algorithm, such as the consensus subgradient (CSG) algorithm. Yet, the CSG suffers from scalability issues due to high problem dimensionality and data-sharing bottlenecks in practical applications. Consequently, implementation in embedded devices is difficult. In this article, we propose a modification of this algorithm, the neighborhood CSG, which aims toward embedded implementation by local reduction of the problem dimensions without hindering the precision of the original CSG algorithm or its convergence rate.
1 INTRODUCTION
The goal in radio positioning networks is to determine the position of a user from their interaction with the network infrastructure. Naturally, it is essential that the positions of infrastructure nodes be known during the estimation. In the case of terrestrial systems, the infrastructure nodes are usually stationary during network operation (in which case, they are called anchors), and their positions are surveyed before network operation begins. The nodes are generally allowed to move, such as in global navigation satellite systems, where their movement is deterministic (Hegarty & Kaplan, 2005). In the following, only networks with stationary nodes (anchors) will be considered.
One of the typical methods of anchor position determination involves the use of laser total stations to survey the network area. This approach is tedious, time-consuming, and labor-intensive, especially in harsh areas; however, it achieves the most accurate results. Therefore, such an approach is suitable for networks with long-term placement. For temporary networks, the time needed for determination must be as short as possible, in order to keep the setup time substantially shorter than the expected network lifetime. Moreover, the usefulness of a network may diminish before the determination is completed; therefore, a faster, simpler method is needed. If anchors of such an ad hoc network are capable of distance measurements, e.g., by two-way ranging (TWR) (Neirynck et al., 2016) in ultra-wideband (UWB) networks, the distances between anchors can be used to estimate their positions by multilateration (Djaja-Josko, 2017; Koppanyi et al., 2018).
It is preferable and straightforward to collect the measurements from anchors into a single point in the network (e.g., a computer) and to then estimate the positions (centralized solution). The results would then be redistributed back to the anchors. The estimation algorithm could be any established algorithm used for tag positioning, with the modification that the anchor positions are estimated. Such approaches have been reported in relative positioning scenarios used in drone swarm positioning (Courtay et al., 2019; Kohlbacher et al., 2018). However, the distances used for estimation do not have to be measured exclusively between anchors, but can also be measured with respect to a “calibration” device, as shown by Djaja-Josko (2017), who placed a device at several locations inside the network and determined its position. The device then measured distances from the anchors within reach, and the anchor positions were estimated from the distances and known position of the calibration device during the measurements.
Centralized estimation methods require all anchors to be connected to a common network for ranging and position data-sharing (e.g., Wi-Fi, Ethernet). However, this requirement cannot be guaranteed, and anchors may directly communicate only with those anchors that are within range. Consequently, frequent data-sharing may introduce an excessive load for the network, along with high power consumption. Therefore, a distributed estimation is more suitable than a centralized distribution.
There are numerous methods for distributed network localization utilizing state probability, optimization, machine learning, or their combination. It is worth noting that methods primarily used for distributed user positioning can be modified for anchor network localization. Among probabilistic methods, a commonly adapted approach is maximum likelihood estimation (Copa et al., 2021; Erseghe, 2015; Meyer et al., 2012; Simonetto & Leus, 2014). Because the problem is non-convex, it is not guaranteed that any solver will find the global minimum instead of a local minimum. Simonetto and Leus (2014) applied convex relaxations in order to increase the likelihood of finding a solution closer to the optimum. The intermediate results are shared, and the alternating direction method of multipliers is employed to find the solution. The method proposed by Meyer et al. (2012) approximates the probabilities with their corresponding beliefs, which are subsequently shared and combined through the network and used as the basis for a particle filter.
Alternatively, the problem can be formulated in a weighted least-squares (WLSQ) sense. For special cases, the exact solution can be obtained via a normal equation (Chan & So, 2009; Lv et al., 2016), provided that a solution exists. In more generic cases, the WLSQ cost function may be minimized by applying the Gauss–Newton method (Soatti et al., 2014) or by machine-learning techniques such as support vector machine or k-means methods (Zhu & Kuh, 2007). This problem has also been approached as relative (formation) positioning with distributed WLSQ (Rajan et al., 2015).
Two similar distributed least-squares approaches were proposed by Koppanyi et al. (2018) and Calafiore et al. (2012). The former utilized the subgradient method with a constant step-size coefficient to update the estimate (all unknown anchor positions), which is then shared and combined until a network consensus is reached. The latter suggested that each node estimate only its own position with the gradient update equation. The value of the step-size coefficient is a result of the network consensus, where the initial step-size coefficients are computed in each node via the Barzilai–Borwein scheme (Calafiore et al., 2012). According to Shen et al. (2018), the anchors can also be localized by means of a variational Bayesian consensus extended filter, in which the state and its covariance are shared and combined among the nodes for a fixed number of iterations.
As we are focused on embedded devices with limited power and computational resources, the power and computational demands of the estimation method should be as low as possible. Consequently, we further focus on gradient-optimization-based methods, which require simple matrix operations only. Specifically, this paper focuses on extending the consensus subgradient (CSG) method (Koppanyi et al., 2018), with an emphasis on suitability for implementation in embedded hardware. We chose the CSG over the solution provided by Calafiore et al. (2012) because it provides individual anchors with more flexibility by allowing them to have different step-size coefficient values. However, the original CSG method does not scale well, as will be discussed in Section 3.3, and requires an impractically high data throughput between anchors for larger anchor networks. This paper aims to provide modifications to the CSG method that will decrease its memory and communication requirements, while maintaining its accuracy. The anchor devices used in this work utilize a UWB chip compliant with the IEEE802.15.4 standard (IEEE Computer Society, 2020; Sahinoglu et al., 2008), controlled by an on-board ARM Cortex-M4 microprocessor.
The implementation of this algorithm may extend the capabilities and utility of current positioning systems. The application of the proposed method is not limited to UWB systems. For example, systems with passive user positioning (Corbalán et al., 2019; Großwindhager et al., 2019; Ledergerber et al., 2015; Navratil et al., 2022) do not require an interconnected infrastructure or any central points for collecting measurements and performing the position estimation. Rather, the positioning is performed on the user side based on messages broadcasted by the infrastructure. Combining such systems with the proposed algorithm would allow their fully automatic deployment and operation. Furthermore, it would be possible to use systems for autonomous exploration, where infrastructure nodes are placed by autonomous vehicles. The vehicles would then benefit from the additional position information. Similar applications can be seen in ad hoc positioning networks deployed by emergency first responders. Alternatively, the algorithm could serve as a sanity check to detect and reconfigure anchors that have been misplaced or intentionally relocated.
The remainder of this paper is divided into several sections. Section 2 defines the network localization problem and formulates it in the form of an optimization problem. Next, Section 3 describes the CSG algorithm and proposes some modifications. Section 4 focuses on an evaluation of the proposed algorithms and their modifications. Section 5 discusses the results of the algorithms in simulated scenarios and one scenario from a real environment. Finally, Section 6 summarizes and concludes the article.
2 PROBLEM STATEMENT
Multilateration is a standard technique for position estimation of a node given several distances between the node of interest and nodes with known positions. To measure the distance d between UWB anchors, the TWR protocol is used. The TWR protocol is a method for measuring distances between devices with free running clocks (without synchronization) by exchanging several messages (Neirynck et al., 2016). By appropriately combining the arrival and departure timestamps, TWR is able to eliminate the mutual clock offset. In our devices, we have implemented the alternative double-sided TWR (ADS-TWR) variant, which also reduces errors induced by the mutual clock drift of the measuring devices. The specifics of the principle, distance calculation, and error analysis have been provided by Neirynck et al. (2016). Notably, Navrátil and Vejražka (2019) showed that the bias of nonlinear ADS-TWR is negligible in practical applications and that the variance of the estimator is similar to that of the linear symmetric double-sided TWR. For the remainder of this work, we assume that the distance measurements are affected by zero-mean additive white Gaussian noise.
In addition to the distance measurements, the initial conditions must be specified as well. It is possible to define the initial conditions in either a form of predetermined positions of a few anchors or a form of constraints on the positions (e.g., the x-coordinate can only be positive); only the former will be considered further.
At this point, it is convenient to represent the anchor network as an undirected graph G(V, E), where V is the set of anchors in the network (graph vertices) and E is the set of edges. Here, an edge el = (i, j) belongs to the set E if and only if there is a measured distance dij = dji between anchors i and j.
Next, let us define the estimation vector x. Vectors are assumed to be column vectors in the following calculations. The estimation vector x is a long vector composed of anchor positions ri:
1
where ne is the total number of anchors with unknown position. The position of each anchor is a three-dimensional Cartesian coordinate. The positions of several anchors are fixed by the initial conditions to ensure that the whole constellation is attached to a coordinate reference frame. We will refer to such anchors as fixed anchors, and nf will specify their total number. Anchors with their positions to be estimated will be referred to as estimated anchors, and ne indicates their total number in the network. It is obvious that the sum of the counts must be equal to the total number of anchors in the network, i.e., n = nf + ne.
For the sake of convenience, anchors are indexed such that the first ne anchors are the estimated anchors and anchors from ne +1 to n are the fixed anchors.
The cost function can be defined for each edge el ∈ E(G) of the measurement graph G as follows (Calafiore et al., 2012; Koppanyi et al., 2018):
2
where the function parameter 
is a set of fixed anchors, which specifies the function, and the operator ||·|| is a standard vector 2-norm defined by a scalar product as follows:
3
In Equation (2), the function also sums the cost for measurements between two fixed anchors. This does not pose a problem, as this cost is either zero or gives a constant offset to the function; thus, it does not change the shape of the function.
The cost function in Equation (2) differs from the function that is usually used for multilateration (Hegarty & Kaplan, 2005):
4
The cost function f does not require explicit computation of a vector length (involving the square root) and its inverse, which makes it computationally less expensive than fb. Consequently, the gradient of the cost function f is free of singularities that occur when any of the estimated distances are zero (i.e., when the two estimated positions are identical). Although such singularities are avoided in common multilateration scenarios, such conditions are more likely to occur for cases in which multiple anchors are estimated at once (e.g., when the position estimates are poorly initialized). Because the solver is intended to be implemented in embedded devices, only the former function f will be considered. For the sake of clarity, we will also assume that the set of fixed anchors Rf is implicitly assigned to the cost function f, which will not be explicitly stated in the equations.
The cost function f is a nonlinear sum of squares; therefore, it can be written as a product of two vectors or vector functions:
5
where 
holds the differences between the estimated distance and measurement, m is the number of measurements and therefore the number of edges is |E| = m. Each element gl(x) is then an individual difference:
6
The vector g is also called the vector of residuals or simply the residual vector.
We conclude this section with a statement of the optimization problem. We are optimizing the values of vector x (anchor positions ri), given the distance measurements and positions of fixed anchors, such that the value of the cost function f is minimized:
7
3 ALGORITHMS
The usual approach for solving Equation (7) would be to collect all of the measured distances at a single point and solve the problem at once (globally) using any of the available nonlinear least-squares (NLSQ) solvers, e.g., Gauss–Newton or Levenberg–Marquardt (LM) method. The resulting positions would then be distributed back to the anchors. Such an approach is further denoted as a global solver; in particular, the LM method will be described and further used as a reference.
As stated earlier, in the considered UWB network, the anchors do not have to be fully interconnected by data links (e.g., Ethernet or Wi-Fi); therefore, gathering the measurements to a single point may be difficult. Instead of utilizing a global solver, the problem shall be solved in a distributed manner by each anchor. Sections 3.2 to 3.4 cover the use of the distributed CSG algorithm and its modifications. Such solvers are referred to as distributed solvers.
3.1 Global Solvers
The LM algorithm is a well-known tool for finding the minima of NLSQ problems. This method combines the behaviors of Gauss–Newton and gradient-descent algorithms in order to speed up convergence and reduce the risk of divergence (Levenberg, 1944; Marquardt, 1963). The behavior is controlled by the value of the positive parameter λk. The algorithm is iterative, starting from the initial guess x0, with the update equation for the estimated vector xk at the k-th iteration as follows (Levenberg, 1944; Marquardt, 1963):
8
where the operator diag(·) sets all off-diagonal elements of a square matrix to zero and gk and Gk denote the residual vector and its Jacobian at point xk, respectively:
9
For the update (Equation (8)), we have used a modified version that accounts for the sensitivity of gk in the cardinal directions by augmenting the diagonal of the 
matrix (Fletcher, 1971).
Utilization of derivative-free (value) optimization methods, such as the Nelder–Mead method, has been considered as well. Nonetheless, the performance of such methods on problems with a high number of dimensions is generally poor (Gao & Han, 2010). A few initial experiments for Equation (7) confirmed that the Nelder–Mead method is not suitable, and therefore, this method was not further investigated in detail.
3.2 Distributed Solvers
In contrast to global solvers, distributed solvers search for a solution in multiple independent nodes in parallel. All of the solution candidates from each node are then shared though the network and combined. Every iteration consists of the update phase, where each node calculates a new solution candidate, and the consensus phase, where the candidates are shared and combined among the nodes. At the end of the estimation, the solution in each network node is identical (or at least sufficiently close), meaning that each estimate has converged to the common minimum and the network has reached consensus (Koppanyi et al., 2018).
The same concept is applied for anchor position estimation with the CSG algorithm, an application presented by Koppanyi et al. (2018). Each anchor i has its own estimation vector x[i], and each anchor runs its own instance of the CSG algorithm. During the update phase in the k-th iteration, the vector 
is used to compute the update vector u[i] with the subgradient step:
10
where 
is the step-size coefficient, which can be fixed or varying during the CSG iterations, 
is the partial residual vector at point 
, which has non-zero elements only for measurements performed by anchor i, and 
is the derivative of the partial residual vector evaluated at 
. The product 
is a subgradient of the cost function f.
Next, the consensus phase follows, where anchors share and mix their update vectors 
. Mixing is performed by gathering the received update vectors 
from the neighbors into a matrix 
(rows are formed from update vectors) and then multiplying by a weight matrix W:
11
where 
is obtained by taking the i-th row of the weight matrix W and taking only the elements that correspond to the neighbors of anchor i. The new update vector 
is then shared and mixed again. This process continues until the values of the update vectors converge, that is, 
. When this holds for all anchors, the network is said to have reached consensus. The value of the final update vector 
is then assigned to the estimation vector 
, and the algorithm continues with the estimation phase. The CSG iterations continue until the estimation vectors converge.
The values and properties of the weight matrix W have a direct influence on the convergence of the consensus phase and thus the algorithm. Therefore, the choice of its values is crucial. Koppanyi et al. (2018) used the Metropolis–Hastings metric (Boyd et al., 2004; Xiao & Boyd, 2004) to define W, whereas a convergence analysis was conducted by Johansson et al. (2008). We will use the same definition as Koppanyi et al. (2018), i.e., each element Wij of W is defined as follows:
12
where d(i) is the degree (number of measurements) of anchor i. Importantly, the row sums of W must all be equal to 1:
13
where 1 denotes a column vector of ones. This ensures that the consensus will converge to the average of the update vectors (Johansson et al., 2008). The definition in Equation (12) originates from random walks on Markov chains and rates the ability of a vertex to distribute information based on its degree. Intuitively, more information will flow through an anchor with a higher degree. When anchor i is connected to an anchor j with a higher degree, then the information of anchor i will not be distributed as well (in comparison to j), and some of the information of anchor i is retained (therefore, Wii > 0).
3.3 Dimension Reduction
The applicability and results of the CSG algorithm described above for anchor position determination were evaluated by Koppanyi et al. (2018) for simulated scenarios. While the algorithm gives promising results, several aspects may cause its implementation in embedded devices to be troublesome.
Firstly, the algorithm has a issue with scalability. The algorithm assumes that the total number of anchors is determined before anchor placement. This assumption arises from the principle of distributed algorithms, where each anchor estimates both its own position as well as that of all other anchors. This assumption also implies that the size of the estimation vector would be ne · ndim, where ndim is the number of dimensions of the positions (either two or three in practical scenarios).
Secondly, during the consensus phase, the estimation vector x[i] must be shared multiple times with other anchors. If the vector x[i] fits into a UWB message, convergence to the network mean could be slow. Moreover, only one message can be transmitted at a time through the UWB channel; therefore, some kind of arbitration is needed. This step could prove to be difficult to manage efficiently and would further increase the time demands of the algorithm.
To address the above issues and ensure that the algorithm is suitable for practical applications, the neighborhood CSG (N-CSG) modification has been devised. From the algorithm principle, it is apparent that an anchor cannot improve the position estimate of any non-neighboring anchor (without a distance measurement), because the corresponding rows of the Jacobian G[i] would be zero. The anchor only passes on such positions; moreover, the anchor has no use for these positions, and a knowledge of the neighboring anchor positions only is more than sufficient. In the neighborhood modification, the size of the estimation vector x[i] is reduced to hold only the position of anchor i and its neighbors, and the number of consensus iterations is limited. This modification substantially simplifies both the computation of the update vector 
and the consensus step. Furthermore, this modification improves the scalability of the algorithm, as its complexity no longer depends on the size of the whole network, but rather on the size of the anchor’s neighborhood.
Effectively, each anchor solves only a part of the problem, as it has only a relatively small portion of measurements available. Each individual anchor i then minimizes its partial cost function f[i]:
14
where E(i) is the set of edges of anchor i.
The mixing equation (Equation (11)) now takes the following form:
15
where 
is the weight matrix of anchor i and mi is the number of neighbors of anchor i.
The matrix of combined update vectors 
is constructed entirely from the standpoint of anchor i. Its columns correspond to the estimated positions of anchor i’s neighbors and the estimate of its own position, whereas the rows are update vectors received from the neighbors. However, the update matrix 
will surely contain zeros for neighbors that are not common for the two anchors. This direct consequence of the estimation vector reduction poses a problem, because the zeros pull all of the estimates toward the zero vector 0 in the consensus phase. The following example illustrates this effect.
Let us consider a portion of a network with four anchors, labeled from 1 to 4. Anchor 1 is a neighbor of anchors 2, 3, and 4, while anchor 3 is neighbors only with anchor 1. This subnet is shown in Figure 1. The update matrix 
for anchor 1 in the l-th iteration is as follows:
16
and the respective row of matrix W is as follows:
17
An example network to demonstrate the computation of 
As we can see in Equation (16), zero vectors do indeed appear in the matrix. These zero vectors remain constant during the consensus phase and are averaged with the remaining vectors. During mixing, the zero vectors will progressively bring the position estimates of anchors 2, 3, and 4 in the new update vector 
toward zero.
To avoid such behavior, the zeros should be replaced with the estimates of anchor 1, which will compensate for the estimation vector reduction to a certain degree. The choice of anchor 1’s estimates is based on convenience, as the estimates are guaranteed to be available (as discussed earlier). Nonetheless, we acknowledge that a heuristic approach (e.g., one that involves anchor degrees) might be more appropriate and might speed up the consensus. The resulting matrix 
from the example in Equation (16) filled with estimates instead of zeros yields the following:
18
3.4 Convergence Rate Improvement
By inspecting the CSG estimate update in Equation (10), one modification possibility can be noted, arising in the value of the step-size coefficient αk. In the original CSG algorithm, the step-size coefficient remains constant, αk = α0. Having a constant step-size coefficient will most likely result in oscillating (zigzag) behavior around a minimum as the estimate approaches the minimum. Dynamically adjusting the step-size coefficient αk can reduce these oscillations and increase the convergence rate. In the following, the step-size coefficient will be referred to simply as the step size.
A line-search method can be employed to find a more suitable value of the step size αk. Such methods focus on a smaller one-dimensional optimization problem, where the optimal αk is identified on the closed interval [0, α0] such that f(x + αk ∇f) is minimized, where ∇f is the derivative of f at x. In this work, two line-search methods are discussed: the backtracking line search (BLS) (Armijo, 1966) and two-way backtracking line search (TBLS) (Truong & Nguyen, 2020).
3.4.1 Backtracking Line Search
The BLS starts with an initial value of the step size, e.g., α = α0, and the descent direction d. Then, it verifies whether the Armijo condition (Armijo, 1966) is satisfied:
19
In other words, the reduction of the function value in the direction d must be greater than the reduction estimated by the product αγ ∇fTd, which is further scaled by the constant γ ∈ (0, 1).
If α does not satisfy the condition in Equation (19), it is decreased, i.e., multiplied by a constant factor β ∈ (0, 1):
20
This operation is repeated until an α value that satisfies the condition in Equation (19) is obtained. The final α value is then used for the estimation, that is, αk = α. The typical value for both β and γ is 1/2 (Armijo, 1966). For the CSG algorithm, the direction d is a subgradient of f at point 
, which is 
.
3.4.2 Two-Way Backtracking Line Search
The TBLS is an extension of the BLS (Truong & Nguyen, 2020). The BLS always starts from the largest step size allowed, α0, which is then sequentially reduced. The TBLS initiates the search with the previously found step-size value, as it assumes that the next optimal step size is similar to the previous one. The method finds the optimal value by both shrinking and expanding the step size in the interval [0, α0]. After initialization, the algorithm checks whether the current α satisfies the condition in Equation (19). If not, the step size is reduced by a factor of β in the same way as in the BLS. If the current step size satisfies Equation (19), then the TBLS attempts to increase the step size by a factor of β as follows:
21
The step size is increased until either the maximal step size (α0) is reached or the Armijo condition (Equation (19)) is violated. The resulting step size is the last value that satisfies both conditions.
3.4.3 Estimate-Update-Step Modification
An alternative option for improving the convergence rate is to use a different estimate update equation. For instance, it is possible to modify the LM method for distributed computations in the anchors to obtain a neighborhood consensus LM (N-CLM) approach. The update vector computation for the N-CLM method is as follows:
22
The consensus phase remains identical to that of the N-CSG, except that the gradient-update step in Equation (10) is replaced by Equation (22). This method should converge faster than the N-CSG, thanks to the inherent update step scaling, but at the cost of increased computational demands (computation of the matrix inverse). Special care should be taken when performing manipulations with values of 
, as the consensus phase will most likely change the estimate given by Equation (22) during the estimation phase.
4 EVALUATION
The performance obtained by the modifications of the CSG algorithm (N-CSG, N-CSG with BLS and TBLS) and N-CLM will be evaluated by applying the algorithms to three generated problem classes and one real problem class based on an anchor network spanning a floor of an office building. The results of the modified algorithms will be compared with the results of the original CSG algorithm and the LM global solver. The anchor positions will be solved in three dimensions. Firstly, the test cases and algorithm configuration will be described. Secondly, the method of evaluation will be detailed, and finally, the achieved results will be presented in Section 5.
Each problem consists of anchor placement, the selection of fixed anchors, distance measurements, and an initial position guess. Problems were generated by a numerical processing script, which could be configured to generate a problem with specific properties (number of anchors, topology, ranging error distribution). For the real, measured problem from an office building, the anchor placement and distance measurements were given, and only the initial guesses were generated.
Three generated problem classes have been considered, as depicted in Figure 2. In class “G” (Figure 2(a)), the anchors are placed in a square grid with even horizontal and vertical spacing (10 m). The neighborhood of an anchor was limited to its closest anchors in the horizontal, vertical, and diagonal directions. Thus, the maximum possible node degree (neighbor count) was eight. Such scenarios mimic the coverage of broad, relatively open areas, e.g., parking houses, warehouses.
Examples of the considered problem with varying selections of fixed anchors White circles indicate the true positions of estimated anchors, filled circles indicate fixed anchors, and lines represent measured distances. (a) GO — Grid with opposite corner selection (b) LC — A cut of ladder class with corner selection (c) RM — Random class with middle selection (d) MO — Measured class with opposite selection
For class “L,” the anchors were placed in a horizontal ladder-like shape, with two anchors in each column. The neighborhoods were created similarly as for class “G,” but the maximal degree was five. A ladder-like structure of anchors is typical for tunnels, roads, railway sidings, etc.
Class “R” cases are random with no predefined structure, and the neighborhoods are limited by a maximal distance (18 m).
The generated problem classes (G, L, and R) had 100 anchors in total; three of the anchors were selected as fixed (97 to be estimated). Fifty problem instances were generated each for classes “G” and “L.” Each instance had the same anchor placement but differed in the measurements and initial guesses. The initial guess 
for an anchor was obtained by perturbing its true position r[i] by a random error vector 
sampled from a uniform distribution in a range of -5 m to 5 m. The reasoning for this type of initial guess instead of the origin (zero) is to imitate placement guessing by using a floor plan, for example, placing the anchors in fairly imprecise guessed locations, rooms, or areas:
23
Similarly, the distance measurements were perturbed with random error. Ranging error was simulated as additive noise with a normal distribution with zero mean and 15-cm standard deviation, which is typical for TWR distance measurements (McElroy et al., 2014).
Furthermore, each problem instance (with the exception of class “R”) was solved for three different selections of fixed anchors, in order to capture the sensitivity of the algorithm to the initial boundary condition. The selection was either from one corner (labeled “C”), opposite corners (e.g., two in the bottom-right and one in the upper-left corner; labeled “O”), or from the middle of the network (labeled “M”). For instance, the “GC” problem class indicates a grid-like constellation of anchors, with the fixed anchors located in a corner. An example of different subclasses is shown in Figure 2. Varying the selection of fixed anchors may provide some insight into how algorithms behave on a given problem with different initial conditions. In practice, it can be expected that the position of only anchors that are close together can be easily determined (because of walls and other obstacles). Therefore, the “M” and “C” cases will be the most feasible in real situations. The “O” case is considered to test whether this initialization has a major impact on the accuracy of the estimation.
For the fully random problem class “R,” 100 instances were generated. In this case, the anchor placement was also generated for each instance, along with measurements and initial guesses. Likewise, each instance was solved with fixed anchors taken either from one corner or from the middle of the network (subclasses “C” and “M”). The subclass “O” was not implemented for the random instances because of the difficulty in maintaining a consistent selection of anchors from the opposite sides of randomly generated networks. An example of a randomly generated network (problem “RM”) is shown in Figure 2(c).
Additionally, the algorithms were tested on a problem class that involved a real-world anchor network and real distance measurements. We denote this test class as “Measured,” labeled by “M.” The constellation covers one floor of an office building and consists of 25 anchors installed under the ceiling. The anchors are connected via Ethernet to a central computer; the network connections also power the devices (power over Ethernet). Although the anchors feature an ARM Cortex-M4 processor, the actual computations were performed in post-processing, in order to maintain the ability to input identical data to all evaluated algorithms. Nonetheless, the distance measurements were collected by the anchor hardware, and the results were stored for later evaluation. The measurements were performed between anchors that both had a clear line of sight (LOS) and that had their LOS obscured (non-LOS [NLOS]) by an obstacle or a wall. Thus, the measurements were affected by both natural measurement errors and errors caused by propagation through objects. The NLOS errors are not compensated for, as these errors are difficult to determine because no prior knowledge of the environment is assumed. Therefore, the resulting position estimates are not expected to be perfect (absolute error below approximately 65 cm). The office anchor network is shown in Figure 2(d).
For the problem class “M,” 50 instances were evaluated; nonetheless, only the initial guesses were generated, as the measurements and true anchor positions were measured in the real environment. Likewise, each instance was solved for three selections of fixed anchors, either from one corner, from opposite corners, or from the middle. Because of the necessity of NLOS measurements, four anchors were chosen as fixed (as opposed to the other problems, in which only three anchors were fixed).
5 RESULTS AND DISCUSSION
In total, 650 instances of 11 problems were prepared to be solved by the six algorithms. Each solver was run for 500 iterations. The CSG and N-CSG (without line search) algorithms were configured to have a constant step-size coefficient α = 10−3. For the two N-CSG variants with a line search, the maximum step-size coefficient was set to α0 = 10−2. The consensus phase was limited to one iteration.
The performance of the algorithms for a given problem was evaluated by the root-mean-square (RMS) of the difference between the estimated and true anchor positions (estimation error). Such evaluations were performed for each of the tested algorithms applied to a given problem. The RMS of the estimation error (RMSE) εk at iteration k was calculated as follows:
24
where Xl denotes the l-th instance of problem X (e.g., GC for problem class G with fixed anchors in one corner), 
is the estimated position of the i-th anchor at the k-th iteration, ri is the true position of the i-th anchor, and n is the number of estimated anchors in the given problem (either 97 or 21).
The partial RMSEs 
were averaged to obtain the mean RMSE at iteration k as follows:
25
where I is the number of problem instances.
The RMS of the cost function ρ[X] value is computed in the same manner. Instead of the estimation error, the cost value f(xk) is used for the LM and CSG solvers, and 
is used for the neighborhood solvers (cost of only measured distances, for anchor i). The final values (after 500 iterations) of both the RMSE and RMS of cost for each combination of problem and solver are presented in Table 1. Naturally, the tabular results do not sufficiently express the speed of convergence; therefore, the average development of the RMSE and cost function RMS after each iteration is plotted in Figures 3 and 4, respectively. Separate plots for each problem are provided, in order to confirm the dependency of the convergence rate and final accuracy on the problem geometry itself.
Benchmark Results for Each Problem–Solver Combination
The RMS of the estimation error ε[X] is shown in meters in the left columns. The right columns show the RMS of the cost function value ρ[X] in m2.
Mean RMS of the estimation error across problem instances (a) Grid — Corner (GC) (b) Grid — Opposite (GO) (c) Grid — Middle (GM) (d) Ladder — Corner (LC) (e) Ladder — Opposite (LO) (f) Ladder — Middle (LM) (g) Measurement — Corner (MC) (h) Measurement — Opposite (MO) (i) Measurement — Middle (MM) (j) Random — Corner (RC) (k) Random — Middle (RM) (l) Legend
Mean RMS of the cost function value across problem instances (a) Grid — Corner (GC) (b) Grid — Opposite (GO) (c) Grid — Middle (GM) (d) Ladder — Corner (LC) (e) Ladder — Opposite (LO) (f) Ladder — Middle (LM) (g) Measurement — Corner (MC) (h) Measurement — Opposite (MO) (i) Measurement — Middle (MM) (j) Random — Corner (RC) (k) Random — Middle (RM) (l) Legend
As expected, the LM solver achieves the lowest value of the cost function among all of the solvers, as shown in Table 1. However, its accuracy is the worst according to the RMS of the estimation error. From Figure 3, it is apparent that the LM error magnitude stagnates above the 1-m level for all problems. This behavior is most likely due to the presence of measurement errors and the fact that the LM algorithm solves the problem as a whole. The LM solver seeks to find the minimum of the function for all variables at once, despite measurement errors. While other solvers consider only a portion of these errors, the LM solver considers them all.
The distributed solvers generally produced estimates with worse costs; yet interestingly, the estimates were much closer to the true anchor positions. The RMS of the estimation error of the distributed solvers is generally lower than that of the LM solver, generally staying below 60 cm. This property can be attributed to the usage of subgradients rather than gradients. Subgradients are affected only by a selection of measurements rather than by the complete set of measurements. The distributed algorithms presented within this article, except for the original CSG, reduce the RMSE at approximately the same rate throughout the iterations. The CSG is considerably slower, as shown in Figure 3, because the updated anchor positions have a relatively low weight in the network-wide consensus phase, in comparison to the majority of anchors that do not update the particular values of the estimation vector.
The ladder-like problem class “L” proved to be more challenging to solve, as all of the solvers experienced the highest estimation errors for this class. This difficulty arises from the nature of the problem class, in which every column consists of only two anchors that measure distances only with the adjacent column, as shown in Figure 2(b). Consequently, an ambiguity is introduced to the solutions. The problems are solved in three dimensions; thus, the set of possible positions given by two distances is a circle in a vertical plane, instead of a point. Such ambiguity worsens the results and can lead to oscillations or even divergence of the results, as shown in Figures 5(a) and 5(b), respectively. To resolve this ambiguity, distances between anchors from more distant rows or columns must be measured or additional fixed anchors must be included throughout the network. Another possibility is to introduce a height constraint, i.e., a constraint on the z-coordinates, if possible. Although the results are not satisfactory, the algorithms behaved in a stable manner and reduced the cost (Figures 4(d)–(f)) and estimation error (Figures 3(d)–(f)) in each iteration.
Selected results of distributed network determination
Rings indicate true positions, filled circles indicate fixed anchors, and crosses indicate estimated positions. (a) Problem LC, N-CSG TBLS solver (b) Problem LC, CSG solver (c) Problem MO, N-CSG TBLS solver (d) Problem GO, N-CSG TBLS solver (e) Problem RM, N-CSG TBLS solver
Concerning the selection of fixed anchors, the results suggest that for regular anchor placement (grid and ladder problems), it is negligibly better to select anchors from opposing corners of the network when possible (see Table 1). However, selecting anchors from one corner yields a smaller increase in estimation error than selecting fixed anchors in the middle. Interestingly, the choice of fixed anchors is more important for the global LM solver, whilst the effect on the distributed solvers is marginal. For random cases, it also appears that having fixed anchors in a corner is a better option, in regard to both estimation error and cost.
The “neighborhood” solvers performed worse for the real problem class “M” than for the other problem classes. The global LM solver and N-CLM solver had the highest estimation error whereas the CSG-based solvers performed much better (RMS of error below 1 m but above 0.5 m). As this problem class is the only class that has been actually measured, the measurements were affected by additional errors due to signal propagation through walls and other obstacles. Therefore, the errors are not guaranteed to have a zero mean. Moreover, the anchors were all installed at approximately the same height. Consequently, the solvers become more prone to range estimation errors, which manifest as uneven progress of the cost function values or even a slight increase throughout the iterations (see Figures 4(g)–(i)). Similarly, the RMS error values do not always monotonically decrease throughout the iterations, as shown in Figures 3(g)–(i).
The distributed solvers achieved the best results for problem class “M” when the fixed anchors were in the middle, rather than in the corners. It is likely that the measurement errors propagate and accumulate throughout the anchor network. By inspecting the network layout in Figure 2(d), one can observe that the paths to the boundary of the network are much shorter when the fixed anchors are selected from the middle, rather than from a corner. Therefore, the error accumulated along the way is smaller as well for the former.
The convergence of a solver is assessed by the development of its mean RMS of the cost function value ρ[X] and the first iteration from which its value becomes stable. Comparisons of the solver mean RMS errors ρ[X] for each problem are shown in Figure 4. As expected, the LM solver achieves the fastest convergence (in less than 50 iterations) for all of the problems, as indicated by a solid line without markers. Interestingly, apart from the ladder problem class, the N-CSG solvers (both with and without a line search) and the N-CLM solver quickly reduce the value of the cost function as well (in less than 100 iterations), whereas the CSG solver lags behind. For the ladder problem class, each of the distributed solvers converge at a similar rate to some of the local minima (this problem class is not fully determined).
For the random problem class, the behavior of the N-CSG solvers without a line search and with TBLS is compelling. The two algorithms exhibit oscillations in cost, unlike the others. Nonetheless, the cost oscillations do not result in RMSE oscillations; instead, the RMSE shows a gradual decrease. This behavior suggests that the solvers have encountered a descending “valley” in the cost function (which might only appear as a valley because of the use of subgradients). The solvers may oscillate between opposite slopes while proceeding in the descent direction.
Despite the oscillations, all N-CSG solvers and the N-CLM were able to rapidly reduce both the cost function and estimation error. However, they do not fully stop the decrease by the last iteration. Although there is still some improvement after 500 iterations, the RMSE reaches acceptable levels substantially earlier. Therefore, a suitable set of stopping conditions based on the cost function should be developed.
As shown by the estimation error in Figure 3 and cost in Figure 4, the N-CSG solvers with a line search performed better than the solver without a line search and achieved results comparable to those of the N-CLM. Figure 6 shows the development of the mean value of the step size α throughout iterations (mean among all anchors and all problem instances) for the two random problems. On average, the chosen value of the step-size coefficient is small in the beginning and then gradually increases to approximately 1.2 · 10−3 and 1.3 · 10−3 for N-CSG BLS and TBLS, respectively. By comparing this coefficient to the constant step size of N-CSG, α0 = 10−3 (displayed as a dotted line), one can see that having a constant step size indeed slows down the convergence. It is also worth mentioning that the step size converges at the same rate as the cost (Figures 4(j) and 4(k)), as expected. Both quantities change only slightly from the 100th iteration onward. This correlation could be used as a part of the stopping conditions.
Mean step sizes across all anchors and random problem instances
The dashed line represents the constant step-size value α0 = 10−3 used by the CSG and N-CSG solvers.
The above discussion and results suggest that the best distributed solvers for anchor network determination are N-CSG with a line search and N-CLM, as they achieved the lowest estimation error and cost value. The proposed solvers should be deployed on embedded hardware, for which the most suitable algorithm is the N-CSG with TBLS, as its convergence is comparable to that of N-CLM with lower computational costs and fewer numerical stability risks (no need to compute the inverse of a matrix). To conclude this section, we present several example results of the N-CSG TBLS solver for the problems “GO,” “LC,” “RM,” and “MM” in Figures 5(a) and 5(c)–(e).
6 CONCLUSION
In this article, the problem of anchor network localization has been described, along with distributed methods for its resolution. A neighborhood modification of CSG was proposed to mitigate scalability issues. The network nodes use only information concerning their neighborhood, i.e., regarding only themselves and the anchors to which they measure distance. As a result, each anchor estimates only the position of itself and its neighbors, rendering the calculations less computationally expensive. More importantly, the computational complexity depends on the size of the neighborhood (usually several anchors), rather than the whole anchor network (possibly tens or hundreds of anchors).
Both CSG and N-CSG solvers feature a constant step-size coefficient, which has a purpose similar to that in the gradient optimization method. The N-CSG algorithm was enhanced with line-search methods to speed up the convergence by finding a suitable value for the step-size coefficient in each iteration. The BLS and TBLS methods were implemented as line-search methods. As a last modification, the LM method was adjusted for distributed computations. The N-CLM solver was derived from the N-CSG via a modification to its estimate update equation.
Finally, the proposed solvers were tested on generated and measured problems, and their results were compared with results given by the original CSG solver and global LM solver. In the presence of measurement errors, the neighborhood algorithms outperformed both the LM and CSG solvers in terms of the estimation error, which was less than 50 cm for most of the problems considered. The estimation error was higher for the real problem, as its measurements were affected by both equipment measurement errors and errors from signal propagation through obstacles and walls.
Among the proposed solvers, N-CLM and N-CSG with TBLS achieved a cost function value convergence comparable to that of the global LM while achieving lower estimation errors. The N-CLM, however, is computationally more expensive and possibly less stable, as it requires the computation of a matrix inverse. Therefore, the N-CSG with TBLS was chosen as the most suitable solver for the anchor network determination problem among the considered algorithms, as it requires only basic matrix multiplication, whilst maintaining reasonable convergence properties.
HOW TO CITE THIS ARTICLE
Krška, J., & Navrátil, V. (2024). Distributed nonlinear least-squares solver for practical network determination. NAVIGATION, 71(3). https://doi.org/10.33012/navi.658
ACKNOWLEDGMENTS
The research presented in this paper was supported by the Grant Agency of the Czech Technical University in Prague, grant no. SGS22/169/OHK3/3T/13. The measurement hardware was developed in cooperation with RCD Radiokomunikace a.s. and supported by the Technology Agency of the Czech Republic, grant no. TE01020186.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.