The knowledge of environmental maps (i.e., map-awareness) can appreciably improve the accuracy of optimal methods for position estimation in indoor scenarios. This improvement, however, is achieved at the price of a significant complexity increase with respect to the case of map-unawareness, specially for large maps. This is mainly due to the fact that optimal map-aware estimation algorithms require integrating highly nonlinear functions or solving nonlinear and nonconvex constrained optimization problems. In this paper, various techniques for reducing the complexity of such estimators are developed. In particular, two novel strategies for restricting the search domain of map-aware position estimators are developed and the exploitation of state-of-the-art numerical integration and optimization methods is investigated; this leads to the development of a new family of suboptimal map-aware localization algorithms. Our numerical and experimental results evidence that the accuracy of these algorithms is very close to that offered by their optimal counterparts, despite their significantly lower computational complexity.

1. Introduction

Indoor localization systems have found widespread application in a number of different areas, both military and civilian [1]. In many cases, the maps of the environments where the users are supposed to lie (e.g., the floor plans of a given building or a spatial road network) are perfectly known and this form of a priori knowledge (dubbed map-awareness in the following) can be employed for improving localization accuracy. In particular, in the technical literature about localization systems map information is usually exploited to (a) improve the quality of position estimates generated by standard localization technologies; (b) set simple geometric constraints in the search space of fingerprinting systems (e.g., see [24]); (c) acquire an accurate knowledge of propagation of wireless signals on the basis of two-dimensional or three-dimensional ray-tracing techniques [5, 6]; (d) develop map-based statistical models for the measurements acquired by multiple anchors in a given indoor environment. It is important to point out the following:(i)Point (a) refers to the wide family of map-matching algorithms (see [7] and references therein), which includes both simple search techniques exploiting geometric and topological information provided by maps, and more refined techniques based on signal processing algorithms, such as Kalman filtering and particle filtering [8]. Distinct types of map-matching algorithms may offer a substantially different complexity-performance trade-off, but the improvement in localization accuracy is limited by the fact that map information is not exploited in the first stage of position estimation (i.e., in ranging). Note also that these algorithms have been mainly developed for transport applications (where they are employed to refine the estimated position provided by GPS or GPS integrated with dead-reckoning in outdoor environments). In principle, they can be also employed in indoor scenarios, provided that a proper map (e.g., a schematic representation of the possible user positions such as a Voronoi diagram) is available for the considered environment; for instance, map-matching methods based on particle filtering have been proposed in [913] for localization systems operating in indoor scenarios.(ii) Fingerprinting methods (see point (b)) make a limited use of the geometrical properties of the propagation scenario, since they mainly rely on a training database. In addition, the preliminary measurement campaign necessary to generate such a database can represent a formidable task in large buildings and/or time-varying scenarios.(iii) Ray-tracing techniques (see point (c)) can certainly outperform the above-mentioned techniques in terms of accuracy; however, real-time ray-tracing is often unfeasible because of its large computational complexity and the detailed knowledge it requires about various physical properties (e.g., electrical permittivity of walls and boundary conditions) of localization scenarios.(iv)Point (d) refers to the so-called map-aware statistical localization techniques, which have been developed for time of arrival (TOA), time difference of arrival (TDOA), and received signal strength- (RSS-) based indoor localization systems [14, Sec. 4.11.5], [15, 16]. Such techniques rely on the availability of map-aware statistical signal models and are characterized by the following relevant features: (a) they are able to compensate for the non-line-of-sight (NLOS) bias (i.e., for the extra delay or extra attenuation due to propagation through obstructions, mainly walls), which represents a major source of error in indoor localization; (b) their use requires measurement campaigns substantially less time consuming than those needed for fingerprinting techniques in the same scenarios. In addition, previous work illustrated in [15, 16] has evidenced that these models can be easily combined with optimal estimation techniques to derive novel localization algorithms.This paper represents a follow-up to [15], where new map-aware statistical models based on experimental measurements have been exploited to develop optimal map-aware minimum mean square error (MMSE) and maximum a posteriori (MAP) estimators. Our previous work has evidenced that these estimators may substantially outperform their map-unaware counterparts at the price, however, of a significant complexity increase, specially for large maps. In this paper, the problem of complexity reduction of map-aware MMSE and MAP estimators is tackled and the following novel contributions are developed:(1)Two novel algorithms for restricting the domain over which the agent position is searched for in map-aware estimation are devised. These algorithms, dubbed distance-reduced domain (DRD) and probability-reduced domain (PRD) in the following, can reduce the rate at which the complexity of map-aware estimators increases with map size without substantially affecting localization accuracy.(2)The application of specific mathematical tools (namely, cubature rules for numerical integration [17, 18] and direct-search methods [19, 20]), which can be exploited to reduce the implementation complexity of map-aware algorithms, is analysed.(3)Novel suboptimal map-aware algorithms resulting from the combination of the above-mentioned algorithms and mathematical tools are proposed and are compared, in terms of both root median square error (RMSE) and computational complexity (namely, overall number of floating point operations (FLOPs)), to maximum likelihood (ML) estimators (note that ML estimators are optimal for map-unaware localization and have been widely adopted for their limited computational complexity [21]). The following is worth pointing out:(i)Most of the performance results illustrated in this paper mainly rely on a set of data generated according to the statistical model developed in [15], which, in turn, is based on experimental measurements. This approach, which is commonly adopted in the technical literature (e.g., see [22, 23]), is motivated by the fact that experimental databases referring to localization systems have usually a limited size; this is mainly due to the time consuming tasks required by any measurement campaign in this field. However, some performance results evaluated on the basis of our experimental database are also shown.(ii)In our work only RSS-based localization algorithms are considered, even if the proposed approach can be easily extended to TOA and TDOA based systems, since all rely on similar statistical signal models [15].(iii)Our contribution focuses on the role played by map-awareness in the localization of still targets in static environments (consequently, the potential improvement deriving from the knowledge of mobility models is not taken into consideration).The remaining part of this paper is organized as follows. In Section 2, some models for indoor maps and measurements in map-aware RSS-based localization systems are illustrated. In Section 3, optimal and suboptimal localization algorithms based on the proposed models are developed, whereas in Section 4 various problems arising from their implementation are analysed and some solutions are proposed. In Section 5 various numerical results about the performance and complexity of the devised algorithms are illustrated. Finally, in Section 6, some conclusions are drawn.

Notations. The probability density function (pdf) of a random vector evaluated at the point is denoted ; denotes the pdf of a Gaussian random vector having mean and covariance matrix , evaluated at the point ; denotes a square matrix having the arguments on its main diagonal and zeros elsewhere; denotes the cardinality of the set ; denotes the th element of the vector and denotes its size; denotes the floor of ; denotes the Euclidean projection of over a domain , that is, the point which minimises the Euclidean distance with ; denotes the centre of the rectangle ; () is the unit vector along the () axis.

2. Reference Scenario and Signal Model

In the following, we focus on a two-dimensional (2D) RSS-based localization system employing devices, called anchors, and whose positions are known, to estimate the position of an agent (see Figure 1). The mathematical models defined in [15] for the map, the wireless connectivity between the agent and the anchors, and the measurements acquired for position estimation are adopted; their essential features are summarised below.

2.1. Map Model

A rectangular uniform map, having support and area , is assumed in the following, since it provides a good approximation of the floor plans of typical buildings (see [15, eq. (1)]). This map consists of the union of nonoverlapping rectangles having their sides parallel to the axes of the adopted reference frame and representing the so-called rectangular covering of , so that (see Figure 1). Note that, generally speaking, the support of a rectangular map can be approximated by different coverings and there is no one-to-one mapping between the support of a (rectangular) map and its covering . However, various algorithms are already available in the technical literature to partition a generic polygon into a set of rectangles (e.g., see [24]) and, in particular, to generate a minimal nonoverlapping covering (MNC) [25, 26] of that is a partition consisting of a minimal number of rectangles for the considered domain. In the following, a MNC of , which is a partition consisting of a minimal number of rectangles for this domain, or a dense nonoverlapping covering (DNC) , which originates from splitting the MNC rectangles whose areas or side lengths exceed a given threshold, is considered, so thathere () is the number of rectangles in the selected MNC (DNC) of . As it will become clearer later, a MNC of can be very useful to reduce the computational load of some map-aware localization algorithms; some algorithms, however, perform better if a DNC is selected. In localization systems, is usually time invariant and such coverings can be computed offline resorting to a number of optimized (but complicated) algorithms (e.g., see [25, 26] and references therein); in our work, however, the simpler algorithm described in [27] is adopted, since it can handle maps containing “holes” (e.g., inaccessible areas of a building floor) and its input data can be easily extracted from the floor plans of typical buildings.

Note that our approach to map-aware localization requires not only the knowledge of a MNC (or, equivalently, a DNC) of the map support , but also that of the obstructions (e.g., walls) the map contains; this knowledge is exploited to evaluate the function , which represents the number of obstructions interposed between two arbitrary points and of . On the contrary, map-unaware localization systems, which are introduced later for comparison, rely on the knowledge of the bounding box of only.

2.2. Connectivity Model

In our work, the coverage region of the th anchor (with ) is assumed to consist of a circle centred at whose radius depends on the transmitted power and signal propagation conditions; inside this region, if the agent is connected with the associated anchor, it acquires a single observation (in particular, a RSS measurement) for localization purposes by exchanging wireless signals with the anchor itself [15, Sec. II.B.] (if multiple measurements are acquired by the ith anchor, they are averaged or filtered in order to generate a single observation ). In practice, in the presence of harsh propagation conditions, the shape of the coverage region is likely to substantially differ from a circular shape, so that the number of observations available to the agent cannot be predicted theoretically; for this reason, the parameter is defined, where is the set of indices associated with the anchors truly connected with the agent. It is easy to show that if , then ; consequently, in a map-aware localization system the agent position has to be searched for inside the domain . In our work, the region is assumed to be well approximated by the rectangular region that consists of all the covering rectangles intersecting with the exact ; that is,whereand belongs to a MNC or a DNC set (in the following, unless explicitly stated, the symbol () denotes () or ()). Note that applying this connectivity model to the considered localization problem results in a preliminary domain reduction; in fact, instead of considering the or rectangles covering the whole map support, all the proposed localization algorithms restrict their search domain to rectangles approximating (2). In a map-unaware system, instead, the agent position is expected to belong toIn this case, it is assumed that for any the domain can be approximated by a square region whose centre is and whose side is ; this implies that (4) consists of a single rectangle.

Finally, we note that the sets and in (3) and the domains (2) and (4) are all available to the agent localization algorithm immediately after all RSS observations have been acquired.

2.3. Statistical Modelling of Observations

We consider a wireless localization system inferring the agent position from a set of RSS observations , with . In [15] it is shown that, given the trial agent position , the map-aware likelihood for the observation vector can be expressed as (see [15, Eq. (13)])Here is the Euclidean distance between the trial position and the th anchor anddenote the mean and the standard deviation, respectively, of the bias affecting the th observation and originating from signal propagation inside the obstructions, and represents the standard deviation of the bias-unrelated and position-dependent noise affecting the same observation. Note that the likelihood function (5) is influenced by the structure of the map support through the parameters and , since both these quantities depend on .

Map-unaware systems, instead, cannot rely on the knowledge of the function in bias modelling. Moreover, state-of-the-art map-unaware systems often employ a LOS/NLOS detection stage before localization, whose output () is the subset of containing the indices of the anchors detected to experience NLOS (LOS) conditions with the agent. Given and the trial agent position , the map-unaware likelihood (see [15, Eq. (15)])is adopted for the same observation vector considered in (5). Here and for (both are zero otherwise) and .

2.4. Calibration

The parameters appearing in both the map-aware and map-unaware likelihood functions (see (5) and (8), resp.) need to be extracted from a set of experimental measurements, as illustrated in [15]. In our work, the performance of the developed localization algorithms has been assessed on the floors of 3 different buildings of the University of Modena and Reggio Emilia (these buildings mainly host offices and laboratories). In each scenario, 15–20 measurement sites have been identified in order to acquire RSS data referring to almost 100 independent links. Then, the method illustrated in [15] has been exploited to extract the values of the unknown parameters appearing in (5) and (8); this has produced the following values:  m,  m,  m,  m,  m, , and  m for the map-aware model (5) and the parameters  m,  m,  m, and for the map-unaware model (8) (please see [15, Table I] for further details). It is important to note the following:(i)In our experimental campaign, we have also found that the same values of the considered parameters can accurately describe the statistical models to be adopted for RSS measurements in different scenarios; this result is motivated by the fact that actually similar propagation conditions are experienced in distinct buildings (which can be classified as light commercial in standard terminology), even if the areas of the considered propagation scenarios (floor plans) are significantly different (see Section 5).(ii)In principle, a preliminary measurement campaign is always required to adapt models (5) and (8) to different propagation scenarios. This need is common to various localization techniques, which may require a time consuming training phase.(iii)Map-aware solutions require a proper preprocessing of the map to compute its MNC or DNC and to evaluate the function (see Section 2.1); however, this step can be carried out offline and does not appreciably affect the overall computational complexity (see also Section 4.1).

3. Estimation Algorithms

In this Section, optimal map-aware (namely, MMSE and MAP) and map-unaware (namely, ML) localization algorithms are analysed first. Then, two novel heuristic (suboptimal) strategies are proposed to reduce the computational load required by optimal estimation algorithms.

3.1. Optimal Map-Unaware Estimator

In a map-unaware context (i.e., in the absence of prior information about the position of the agent), the ML approach can be adopted to generate the optimal estimate of [28], where is given by (8). In practice, this estimate can be evaluated aswhere is a single rectangle (see Section 2.2). The last expression deserves the following comments:(1)The cost function appearing in (9) exhibits a discontinuous behaviour, even if it does not contain map-aware functions (e.g., ), since and are assigned a constant value or zero depending on . Consequently, in principle, optimization methods involving gradients or Hessians cannot be used, unless the cost function is, at least approximately, continuous.(2)The evaluation of in (9) requires solving a constrained least squares problem with a quadratic regularization term (represented by the product appearing in the right-hand side of (9)) which penalizes trial positions affected by large noise and bias. Unluckily, standard optimization methods do not necessarily achieve the global minimum; when this occurs, large localization errors are found (note that some methods developed to mitigate the problem of local minima, like projection onto convex sets (POCS) and multidimensional scaling (MDS), cannot be straightforwardly applied in this case because of the significant complexity of our observation model).(3)The complexity of the cost function in (9) can be related to as .

3.2. Optimal Map-Aware Estimators

In a map-aware context, a MMSE or a MAP approach can be employed in position estimation. In particular, it is not difficult to show that the MMSE estimate of is given by [28]Substituting (5) in (10) yields, after some manipulation,whereThe last result deserves the following comments:(1)The MMSE estimation procedure, unlike its MAP and ML counterparts, does not involve an optimization step and, consequently, does not suffer from problem of local minima. However, it requires numerical multidimensional integration.(2)Our numerical results have evidenced that (a) the product appearing in (11) exhibits small variations inside the considered map if the anchor density is approximately constant and (b) the accuracy of MMSE estimation is negligibly influenced by the presence of this term for any such that (i.e., for any agent position far from the anchors). Then, discarding this term yields the approximate (and computationally simpler) variantof the MMSE estimator.(3)The computational load required by the evaluation of is (see Section 2.1), where is the number of segment intersection tests accomplished in processing a given observation; consequently, the complexity required by (12) is approximately , where is the average of evaluated over the observations collected in (note that is nonlinearly related to the shape of obstructions, the trial agent position, and the anchor positions).An alternative to the MMSE estimator is offered by the so-called MAP estimator, which can be expressed as [28] or, if (5) and Bayes’ rule are exploited, aswhereis the MAP cost function. This function deserves the following comments:(1)It is not differentiable, since it depends on . Consequently, its gradient/Hessian matrix cannot be evaluated analytically and, in principle, steepest descent methods cannot be used; in practice, such methods can be exploited if it is assumed that the cost function is approximately continuous.(2)For large , it is characterized by a better numerical stability than that of the function (12). In fact, for large , the argument of the function appearing in becomes small, so that the resulting numerical values of itself may quickly drop below machine precision. On the contrary, this problem is not experienced with (15), because of the presence of a natural logarithm in its expression.(3)Its computational complexity is similar to that of , which is .

3.3. Distance-Reduced Domain Map-Aware Estimator

The computational complexity of the MMSE and MAP estimators described above becomes unacceptable when the considered map is large. To mitigate this problem, a new technique for restricting the domain over which the agent position is searched for has been developed; it consists of the following three steps:(1) Map-Unaware Estimation. A raw estimate is generated by the map-unaware estimator (9), without performing any NLOS correction (this is equivalent to assuming and in (9)).(2) Domain Reduction. A portionof close to is extracted, where is a subset of the values of the index associated with the set of rectangles forming (2); in practice, is generated according to the heuristic criterionwhere is a threshold distance, is a fixed parameter, and is defined by (3). Note that the definition given above for entails that only the covering rectangles within standard deviations from are selected.(3) Map-Aware Estimation. The final estimate is generated using either the MMSE (13) or the MAP estimator (14) under the constraint that ; intuitively, since the subset is smaller than (2), the computational load required by this step is lower than that required by (13) or (14).Note that the proposed technique mitigates localization complexity by restricting the search domain on the basis of a distance criterion, in the sense that only a subset of rectangles, close to the raw estimate , is taken into consideration. For this reason, in the following, this procedure is dubbed distance-reduced domain MMSE (DRD-MMSE) or distance-reduced domain MAP (DRD-MAP), if a MMSE or a MAP estimator is employed in the last step, respectively. Note that DRD performance strongly depends on (a) the accuracy of the estimate and (b) the criterion adopted in generating the set (17) (if the proposed criterion fails to select the map portion where the agent truly lies, then a large error is made in step (3), where it is assumed that ). In principle, other heuristic criteria (e.g., based on selecting a fixed number of rectangles close to ) could be adopted; however, (17) turned out to provide the best accuracy/complexity trade-off among the “distance-based” criteria we tested in our computer simulations.

3.4. Probability-Reduced Domain Map-Aware Estimator

An alternative to the DRD approach is based on a sort of probabilistic criterion and consists of the following two steps:(1) Map-Aware Raw Estimation. Similarly to step (1) of the DRD algorithm, a portion (16) of is extracted. However, in this case, the integer set is generated by taking the values of the index associated with the largest elements of the setwhere is a fixed parameter, , is given by (5), is defined by (3), and is a proper integer set, whose target is reducing the cardinality of (and thus the overall computational load of the estimation algorithm). In practice, has been selected, where is a fixed real and positive parameter, so that only rectangles spaced by at least meters are considered (instead of all rectangles forming the search domain (2)). The effect of the choice expressed by (18) is to generate the reduced domain (16) by selecting the rectangles that, on the basis of their centers, exhibit the largest likelihood values.(2) Map-Aware Estimation. The final estimate is evaluated exploiting either the MMSE (13) or the MAP estimator (14) under the constraint that .The resulting estimation technique is dubbed probability-reduced domain MMSE (PRD-MMSE) or probability-reduced domain MAP (PRD-MAP), if a MMSE or a MAP estimator, respectively, is employed in the last step of the procedure illustrated above. It is important to point out that this approach offers some advantages with respect to its DRD counterpart, since (a) the criterion (18) is more closely related to the optimal MAP criterion than (17) and (b) all the PRD steps rely on map-awareness, that is, on the most accurate modelling we propose. However, it should be also taken into account the fact that the first step of the PRD approach may entail a significant complexity, since the evaluation of a map-aware likelihood is required.

4. Implementation of Estimation Algorithms

All the map-aware and map-unaware estimation algorithms illustrated in Section 3 require numerical integration and/or the application of optimization methods and involve nonlinear and nondifferentiable cost functions. For these reasons, their implementation raises various problems; some possible solutions are illustrated below.

4.1. Map-Aware Implementations

In this section, two different implementations for the evaluation of the MMSE estimate on the basis of (13) are proposed and analysed in detail. The first implementation, called , is based on the use of the so-called cubature formulas (e.g., see [17, 2931]), which approximate any multidimensional definite integral as a finite sum of terms depending on its integrand. In particular, in the following, any integration over the search domain is expressed as the sum of distinct integrals, each referring to a distinct rectangular domain (see (1)), and the cubature formulas illustrated in [30] are exploited for the evaluation of each of the integrals. Then, from (13) the approximate expressionis easily inferred, where is the th node referring to the th rectangle in (2), is the corresponding weight, , and , whereas and denote the th node and the th weight, respectively, of the cubature formula of order . It is worth mentioning that, in principle, the larger the size of is, the stronger the fluctuations of the integrand function to be expected over this domain are; consequently, high cubature orders (i.e., large values of ) should be selected for large rectangles, at the price, however, of an increase in the computational load. In our computer simulations, the heuristic formula has been adopted for selecting the number of nodes; here, denotes the function defined in [30] to map the degree of a polynomial into the number of integration nodes, is the polynomial degree heuristically associated with the (nonpolynomial) integrand functions appearing in (13), is the area of , whereas and are real parameters whose values are listed in Table 1. It is also important to point out that (a) all the sums appearing in (19) require the evaluation of the same quantities and (b) is noniterative and its approximate complexity is , where and denotes the average number of nodes selected for integrating over each of the rectangles. Further details about our implementation of can be found in [27, 32].

The second implementation, called , is based on the well-known trapezoidal integration rule; then, from (13), the approximate formulais easily inferred. Here, the nodes of the th integration formula correspond to the vertices of a regular grid extending over and characterized by a spacing (see Table 1), whereas the corresponding weights are evaluated aswith and . It is important to point out the following: (a) similarly to (19), the three sums in (20) require the evaluation of the same quantities ; (b) is noniterative and its approximate complexity is , where and () denotes the average number of integration points selected along the () direction.

In implementing the MAP estimator (14), the search over the space has been turned into searches over the distinct rectangles covering . This choice is motivated by the fact that (a) the constraint , unlike , cannot be directly formulated as a set of linear inequalities, as required by standard optimization techniques; (b) solving distinct optimization problems, each involving a search domain much smaller than , strongly reduces the probability of reaching local minima; (c) when iterative optimization algorithms are adopted over each of the rectangles , fast convergence can be usually achieved if their centres are selected as initial points. Note that in this case direct-search methods [19, 20] should be preferred to standard Newton-based or gradient-based methods, since the first class of methods usually admits simple implementations for constrained problems and does not require smooth cost functions. These considerations have led us to develop 3 different implementations, called , , and , of the MAP strategy (14). The first one (namely, ) is summarised in Algorithm 1 and is based on the use of the MATLAB fmincon routine to solve the optimization problem appearing in step (2). This implementation is characterized by various adjustable parameters; their values, which are listed in Table 1, have been selected to minimise the overall computational complexity.

Require: covering of (2); (15).
(1) for     do
(4) end for
(5) ,

The second implementation (namely, ) is summarised in Algorithm 2 and is based on discretisation of the search space and on the use of iterative direct-search methods (see [19, Sec. 3]). In practice, in its th iteration, for each rectangle of the selected covering, the cost function (15) is evaluated at the vertices of a rectangular grid characterized by spacing (see lines (4)-(5)); then, on the basis of the resulting values, only the most promising rectangles are saved for the next iteration (see line (8)), in which a halved step size is adopted (see line (9)). Note that both the computational complexity and the accuracy of depend on the initial and the final values of the step size (denoted by and , resp.), so that a proper trade-off has to be achieved in the selection of these parameters (the values we adopted are listed in Table 1).

Require: covering of (2); (15); , and (Table 1).
(1) ,
(2) while     do
(3)  for     do
(4)   ,
(6)   ,
(7)  end for
(8)   indices of the smallest elements of
(10) end while
(11) ,

The third implementation (namely, ) is summarised in Algorithm 3 and is based on a different direct-search method, known as the compass algorithm [20, Sec. 8.1]. Similarly to , exploits a discretisation of the search space and an iterative reduction of the step size; however, the grid it employs is not regular. In fact, in each rectangle, the domain is explored moving from a “current trial position” on the basis of given displacement vectors or (see line (6)) in order to identify the direction along which the cost function decreases (see line (8)). Even in this case, the overall computational complexity and the accuracy of this implementation depend on the initial and final values of the step size (see Table 1), which need to be carefully selected.

Require: covering of (2); (15); and (Table 1).
(1) for     do
(3)  ,
(4)  while     do
(6)   for d     do
(8)    if   and   then
(9)     ,
(10)     bImproved True, break
(11)    end if
(12)   end for
(13)    only if bImproved = False
(14)  end while
(15) end for
(16) ,

It is not difficult to show that the overall computational complexity of each of the MAP implementations proposed above is , where denotes the overall number of times the cost function (15) is evaluated. Unluckily, cannot be easily related to the other parameters of the proposed MAP implementations because of its nonlinear dependence on such parameters; consequently, a more accurate estimate of the computational burden required by the proposed implementations cannot be provided.

Let us focus now on the problem of combining the implementation of the MMSE and MAP estimators illustrated above with the DRD and PRD techniques described in Section 3.3. As far as the DRD technique is concerned, we note that all the proposed MAP/MMSE implementations can be employed in its last step. In the following, however, the implementations of DRD-MMSE and DRD-MAP techniques are always based on and , respectively; for this reason, they are denoted by   and , respectively. It is important to point out the following:(1)The algorithms and are not intrinsically iterative. However, they can be easily modified to include a mechanism for iteratively updating the value of the parameter in order to improve their robustness against inaccuracies in the initial setup (see Algorithm 4).(2)The complexity of and is approximately , where denotes the overall number of evaluations of the MMSE integrand (for ) or that of the MAP cost function (for ), and represents the overall number of evaluations of the ML cost function.

Require: covering of (2); () estimator; constants and (Table 1).
(1) , computed assuming
(3) for (see (17))
(4) while     do
(5)  ,
(6) end while
(7) () evaluated in the reduced map

In our work, and have been also selected for the last step of the PRD technique; the resulting implementations are denoted by and , respectively, in the following. It is worth mentioning that(1)both and are not iterative;(2)their complexity is , where and denote the overall number of evaluations of the MMSE integrand (for ) and that of the MAP cost function (for ) in the first step and in the second step, respectively (see Section 3.4);(3)the values of parameters and need to be carefully selected.Our numerical results (see Section 5) have evidenced that whatever option for map-aware estimation is selected, the parameter (i.e., the average number of intersection tests accomplished in processing each observation) appearing in the complexity estimate of all the implementations developed in this paragraph plays a significant role in penalizing them with respect to their map-unaware counterparts. To reduce , the following two techniques have been devised:(1) Segment Partitioning (SP). This consists of an offline phase, where (a) is partitioned in multiple cells using a nonuniform rectangular grid (in a way that each anchor lies on a vertex of the grid) and (b) the obstruction segments falling inside each cell are identified. Then, when is evaluated, only the segments of the cells intersecting with the segment are taken into account, thus reducing .(2) Precomputing (PC). This is based on an offline evaluation of the function over points (uniformly) distributed inside each rectangle of the available MNC or DNC; this requires storing or integers in a table, where denotes the average number of connected anchors. This technique is “memory hungry” for large maps; however, the precomputed table replaces that memorising information about the position of obstruction segments and can be stored efficiently using a limited number of bits for the integer representation of the precomputed values of . Moreover, this memory requirement is often more than counterbalanced by the reduction of computational complexity (since in this case).In our computer simulations, the first technique has been adopted for all the map-aware estimators. The second technique, instead, has been used in the first step of PRD-MAP estimation only; the resulting implementation is dubbed . This choice is motivated by the fact that this estimator can be easily adapted to the use of precomputed quantities, since it evaluates the map-aware likelihood (8) at a set of fixed points (see Section 3.4).

4.2. Map-Unaware Implementations

Before illustrating our implementations of the ML estimator in (9), it is important to mention that if the constraint was neglected, standard algorithms (e.g., the Nelder-Mead simplex method [33]) could be employed to solve the optimization problem it involves. Unluckily, such a constraint plays an important role, specially in large maps, since (a) minima different from the global one are substantially less likely to be contained in the domain than in the whole space and (b) the map-unaware likelihood (8) is not defined for , as already explained in [15, Sec. II.C.]. For these reasons, in the 3 different implementations of the ML estimator (9) proposed below, the constraint is always accounted for. The first implementation, called , relies on the MATLAB fmincon routine, which is employed to solve the optimisation problem of (9). The second one, dubbed , is based on partitioning into subrectangles (whose sides have a length not exceeding a threshold ) and on running the MATLAB fmincon routine over each of the subrectangles; in this case, the final estimate is selected among candidates according to a minimum-cost criterion (note that a similar approach has been adopted in Algorithm 1). Finally, the third implementation, called , employs the compass algorithm (see Algorithm 3) under the assumption that and . Note that the complexity of each of the proposed implementations is , where is the number of times the ML cost function is evaluated. The values of the parameters selected for each of them are listed in Table 1.

Finally, a diagram summarising all the implementations developed in this section is provided in Figure 2.

5. Simulation Results

In this section, the implementations described in Section 4 are compared in terms of accuracy and complexity. In our analysis, the three different maps, shown in Figure 3, are considered. Such maps are characterized by different sizes and obstruction densities (see Table 2), by a specific deployment of anchors (whose positions are identified by green squares in Figure 3), and by a similar number of anchors per unit area (1 anchor each  m2), so that a fair comparison to the data referring to different maps is ensured. In addition, in all cases, the anchor coverage radius (see Section 2.2) is  m (with ).

The remaining part of this section is organized as follows. In Section 5.1, the method employed to estimate the computational complexity of the proposed implementations is described. In Section 5.2, the structure of the software simulator developed to assess both the accuracy and the complexity of our implementations is illustrated. Finally, in Sections 5.3 and 5.4, various numerical results based on statistical models and experimental measurements, respectively, are analysed.

5.1. Evaluation of the Computational Complexity of Different Implementations

Map-aware implementations are expected to outperform their map-unaware counterparts at the price of a larger computational complexity. However, in this case, an accurate assessment of computational complexity is a challenging task, mainly because of the adaptive and data-dependent nature of the proposed algorithms; for instance, the term (appearing in all the formulas summarised in Table 3) cannot be easily related to the size of the search domain or other implementation-specific parameters. For this reason, the computational complexity of each implementation has been assessed by estimating the average amount of FLOPs required in our computer simulations; this task has been accomplished under the following assumptions:(1)Basic algorithmic steps (e.g., basic checks and assignments in Algorithms 14 for map-aware estimation and LOS/NLOS detection in map-unaware estimation) play a negligible role.(2)The computational complexity of each implementation mainly depends on the average number of evaluations of integrands (in MMSE estimators) or of cost functions (in MAP and ML estimators); as already pointed out in Section 4, the complexity of such evaluations depends, in turn, on and . For this reason, in our simulations, counters for , , and have been used in order to estimate , , and , respectively (see Table 3).(3)Since distinct types of operations are characterized by different computational weights, the expression has been adopted, where is the number of FLOPs of type performed in each evaluation of integrands or cost function and is the associated computational weight. In particular, in our simulations, the weights have been assigned to square roots, exponentials, logarithms, divisions, multiplications, additions, and comparisons, respectively.Our approach is expected to provide statistically meaningful results about only if a large number (in practice, at least ) of statistically independent measurements are available. Acquiring such a large number of experimental data would require a burdensome experimental campaign involving at least 150 measurement sites (so that independent links would be available). To circumvent this problem, the computational complexity of the devised algorithms has been assessed mainly through computer simulations, following the approach always adopted in related work (e.g., see [22, 23]). However, some numerical results based on the available measurements are also shown (see Section 5.4) to support our conclusions.

5.2. Description of the Developed Simulator

A software simulator has been developed (the MATLAB code used to generate the numerical results is publicly available (together with other results not included in this paper for space limitations) at http://frm.users.sourceforge.net/publications.html, in accordance with the philosophy of reproducible research standard [34]) to test the proposed implementations in the harsh propagation conditions characterizing typical indoor environments [15]. In each simulation run, the following steps are accomplished sequentially:(1) Generation of the Agent Position. The true agent position is drawn from the selected map according to the rule , where for and elsewhere (in other words, the uniform statistical map model of [15] has been adopted).(2) Selection of the Connected Anchors. The th anchor is deemed to be “connected” with the agent only if the following conditions are jointly satisfied: (a) , (b)  m2, and (c)  m. In addition, it is assumed that each agent-anchor link is characterized by a “failure probability” (due, e.g., to synchronization or connection problems between the agent and the anchor itself).(3) Generation of the Noisy Observations. The noisy observation vector is drawn from the map-aware likelihood (5) according to rule ; the associated sets and are randomly generated assuming the presence of errors with probability (to account for the presence of a nonideal LOS/NLOS detector).(4) Position Estimation. All the estimator implementations discussed in Section 4 are run sequentially; they are all fed by and (2) (map-aware estimators ) or (4) (for map-unaware ones). Their output is either a position estimate or a failure indication (occurring, e.g., when optimization routines do not converge).In our simulations, steps (1)–(4) have been repeated for at least times in order to generate accurate estimates of the RMSE (which is more robust to outliers than root mean square error) and of the mean number of FLOPs for each implementation.

5.3. Accuracy and Complexity Comparison

Figure 4 shows the accuracy (in terms of RMSE ) of the implementations discussed in Section 4 for the 3 considered maps. These results show the following:(1)The accuracy achieved by the DRD-MMSE (DRD-MAP) and PRD-MMSE (DRD-MAP) estimator is similar to that of MMSE (MAP) estimators (e.g., for map #3).(2)The implementations and provide a 10%–15% performance improvement with respect to , , and (e.g., and for map #3).(3)Map-awareness significantly improves estimation accuracy (up to 250%), as already evidenced in [15] (e.g., and for map #3).It is important to point out that the estimation accuracy of the proposed estimators is mainly related to (a) the average spatial density of obstruction segments (i.e., to the ratio ) and (b) the average spatial density of anchors (i.e., to the ratio ). Since the propagation scenario is similar in the three maps (light commercial buildings) and a similar anchor density has been adopted for the sake of a fair comparison (see Table 2), the RMSEs provided by the considered estimators are similar over the three maps and do not necessarily increase with (e.g., map #3 is larger than map #2, but a lower RMSE is achieved for map #3).

Note also that Figure 4 does not unveil some important aspects of the considered estimation algorithms, which emerge, instead, from the related box plot of Figure 5 describing the distribution of the RMSE outcomes for map #3. Here, the box represents the 75th percentile, the line within the box denotes the median, the black dashed whiskers are upper adjacents, and the cross markers represent the outliers. From Figure 5, the following can be easily inferred:(1)The behaviour of all the estimators is characterized by the presence of multiple RMSE outliers associated with  m as well as of upper adjacents and values of 75th percentiles extending well beyond their median values. This is due to the statistical model of Section 2, which reproduces harsh propagation conditions degrading RSS measurements in practical systems. Note that nonlinear filtering (e.g., particle filtering) could be employed to reject, at least in part, outliers when tracking moving agents.(2)The implementations and are characterized by the smallest dispersion of estimation errors (in fact, upper adjacents are at 7.3 m for both).Some numerical results about the complexity of the estimators implementations are shown in Figure 6 (note that a logarithmic scale is adopted for since complexities range roughly from to FLOPs). These results lead to the following conclusions:(1)The MAP and MMSE estimators based on the DRD and PRD techniques are characterized by a computational complexity 4–11 times lower than that of their counterparts not employing such domain reduction techniques, especially for large maps (e.g., and for map #3). This result can be related to the mean ratio between the areas (16) and (2), that is, to the amount of “domain reduction.” In fact, our simulations have evidenced that, in the case of map #3, on average, the DRD and PRD-based estimators restrict their search domain to 13% and 4% only, respectively, of the whole domain . This makes the complexities of DRD and PRD estimators much closer to that of map-unaware estimators (at the price of a negligible sacrifice in accuracy).(2)Even if all map-aware implementations exhibit similar accuracy, they may require significantly different complexities. In fact, the parameter ranges from (for the most demanding map-aware algorithm) to