Relational Positioning Method for 2D and 3D Ad Hoc Sensor Networks in Industry 4.0

Featured Application: Position location information plays a crucial role in the automated and adaptable operation of Industry 4.0, allowing sensor fault detection as well as operators localization. The development of 3D positioning techniques is crucial for massive sensor network operability imposed by the Internet of Things. Several network protocols depend on position location information for security and routing purposes. Alternative positioning techniques are necessary for new paradigms imposed by new network applications. The proliferation of 3D applications in Industry 4.0 demands alternative 3D location algorithms. Abstract: Industry 4.0, or smart factory, refers to the fourth industrial revolution, as a result of which it is expected that massive sensor deployment, high device connectivity, automation, and real-time data acquisition support improve industrial processes. Sensor fault detection and operator tracking allow adequate performance of the sensor network supporting Industry 4.0, and sensor localization has become crucial to enable sensor fault detection. Hence, the development of new localization methods is necessary for environments where GPS localization technology is unfeasible. Furthermore, position location information (PLI) is a crucial requirement for the deployment of multiple services and applications in wireless ad hoc and sensor networks supporting Industry 4.0. Three-dimensional (3D) scenarios are considered to extend the applicability of PLI-based services. However, PLI acquisition presents several challenges within any 3D ad hoc sensor network paradigm. For instance, conventional triangulation algorithms are often no longer applicable due to the lack of direct land-ﬁxed references and the inclusion of a third-dimension coordinate. In this paper, a PL relational searching algorithm suitable for 2D and 3D ad hoc environments was formulated in terms of a discretized searching space and a vertex weight assignment process according to distance relationships from anchor nodes to the node to be located on a feasibility space. The applicability of the proposed algorithm was examined through analytic formulation and simulation processes. Results show small distance errors (less than 1 m) are achievable for some scenarios. j,k for location estimation of node n 0 . This process is explained in detail in this section, and it is further developed when the 3D localization method is explained.


Introduction
On a daily basis, position location (PL) becomes a more important issue in wireless systems for security and law enforcement purposes, as well as for resource management, new location-based services (LBSs), and applications (LBAs). Thus far, most of the literature has been focused on planar scenarios. However, sensor networks and the Internet of Things (IoT) have become pervasive concepts that are expected to be used in billions of devices with multiple applications. As an example, the number of IoT-connected devices in 2018 was about 7 billion. Nevertheless, the value of IoT and sensor networks is often related to location information. Therefore, the importance of 3D localization becomes more relevant, as it enables a wide variety of new LBS and new LBA.
The concept of Industry 4.0 denotes the emergence of smart factories where automated and decentralized production is enabled by the Internet of Things (IoT). In Industry 4.0, all objects in the factory (i.e., sensors, tools, actuators, etc.) are meant to be connected and communicating with each other, in order to make possible automated and smart adaptive production. In smart factories, it is possible to dynamically adapt the factory operation depending on the operator's location. This implies the knowledge of the position of all the involved operators. Additionally, position location information is crucial for sensor fault detection, which ensures an adaptive factory operation. This creates the necessity for developing alternative position location methods for environments in which GPS systems or traditional indoor localization systems are not able to provide information, [1].
In modern network scenarios, position location is necessary to assist network tasks such as routing, resource management, and the delivery of location-aware services. However, in ad hoc and sensor networks, connectivity constantly changes due to variations in propagation conditions, and node activity or node motion, and these variations create changes in the network topology, [2,3]. In addition to global navigation satellite systems (GNSSs) that provide 3D geolocation for outdoor applications, most existing positioning techniques deal with planar location environments that are not sufficient for many applications, [4]. Although GNSS has been in place for many years with different technologies, currently, its availability and accuracy are often insufficient for indoor applications. It is recognized that 3D localization is necessary for several applications, such as in emergency call services in tall building environments, museum guides, warehouse stock management, crew monitoring, drone localization, sensor deployment, etc. [5]. Additionally, emerging applications, such as event monitoring for vehicular networks (VANET) [6] or special situations for the safety and security of people in smart cities, can be assessed by sensor networks, when accurate sensor location information is available.
In the IoT context, location information enables the implementation of new protocols in several functions of the network layer such as routing, sensor energy harvesting, and security. Similarly, it is responsible for enabling multiple applications such as asset tracking, connected vehicles, smart appliances, tracking of people by wearable devices, assisted living, etc. Therefore, innovative localization methods need to be developed in order to adapt to the new paradigms imposed by the massification of heterogeneous sensor nodes.
In this paper, we introduce a new relational searching localization method that is not based on any previous localization technique. Some of its advantages are that distances are estimated based on a hop count criterion, and a very low number of known location nodes or anchors are needed to provide sensor location estimation. These characteristics make it suitable for IoT environments, where low-equipped sensors with very low transmission ranges need to be considered. We present a location-searching method on a discretized space defined by available anchor nodes (ANs). The method used sole distance relationships between ANs and the node of interest, rather than mean square optimization, to develop likelihood indicators for vertices in the searching space according to those distance relationships. Since the process occurred in location systems, observations were prone to impairments due to errors in distance estimation. In addition, the absence of direct connectivity from nodes to ANs is observed in ad hoc routes that depart from a line-ofsight trajectory. The proposed methodology is adequate for 2D and 3D localization, and its performance was assessed throughout simulations assuming different node densities, transmission ranges, and the number of available ANs. Results under different scenarios show a remarkable performance of the method reaching a small distance error (less than 1 m) for some scenarios. Simulation results show that, compared with the 3D positioning algorithm introduced in [7], the proposed method achieves a significant improvement (approximately 50% improvement when combined with finer discretization).
The rest of the paper is as follows: In Section 2, a discussion on related research about localization in 2D and 3D situations is provided. In Section 3, we introduce the relational location-searching method proposed in the paper, starting with a 2D setup and extending it to the 3D scenario. In Section 4, results are presented with their corresponding discussion.
Lastly, Section 5 provides the conclusions.

Related Research
Traditional localization methods using trilateration are commonly used in 2D for the line of sight (LOS) wireless scenarios, resulting in a good performance in the presence of small distance estimation errors [8]. However, those methods do not deal with the 3D multi-hop wireless sensor scenario, in which one generally finds environments with no LOS between nodes and anchors or access points (AP). Thus, range estimations from anchor nodes (with known coordinates) to the node to be located can be based on the aggregation of hop lengths in the connecting paths, which leads to an additive distance error, [9], impairing the node-to-anchor distance estimation process.
In ref. [10], the authors propose a beacon-based geometric localization method for 2D sensor networks, where a node initially chooses two reference points, then they create what they call the node's residence area by the range intersection of the two reference points. Subsequently, the node differentiates its estimated position within the residence area using a third reference point, by approximating the radius and half-length of the chord with respect to one of the reference points. Then, the radius and half-length of the chord are used to estimate the sagitta of an arc. Finally, the node estimates its position using the radius, half-length of the chord, and sagitta of an arc. There are several reported works for 3D localization, some of them proposed for indoor environments, where the position of a node is estimated by relating the actual observation with a set of collected observation data previously taken by a set of sensors deployed strategically in the indoor scenario. In ref. [11], the authors provide a survey of indoor position location systems updating the new transmission technologies.
In ref. [12], the authors present a survey for indoor localization techniques in the context of IoT. In this paper, the authors discuss different indoor localization systems focusing on user device tracking. They develop an evaluation framework to compare existing indoor localization techniques focusing on the IoT demands. They divide the existing indoor localization techniques into two classifications: monitor-based localization (MBL) and device-based localization, (DBL). In the MBL, a set of anchor nodes passively obtains the position of the user or entity connected to the reference node. In the DBL, the user device uses some reference nodes (RNs) to obtain its relative location. Their evaluation framework focuses on several IoT relevant parameters such as cost, service availability, energy efficiency, accuracy, scalability, etc.
Outdoor 3D localization has been commonly conducted through GNSS systems, [4]; however, there is a need to develop alternative 3D PL techniques suitable for places where satellite coverage is not available or limited. In [13], the authors present a localization system for 3D ad hoc networks using aerial sensor nodes equipped with UWB devices to provide 3D positioning capabilities. The proposed technique in [13] firstly uses a local coordinate system that is defined using sensor nodes. Secondly, a calculation of the sensor's coordinates is performed based on a regular trilateration technique. Finally, the sensor nodes progressively become anchor nodes as their locations are estimated. The method is evaluated in a single hop scenario. The method relies on nodes upgrading with UWB devices and antennas to provide 3D positioning capability. The system performance depends on the goodness of distance estimation. In [14], the authors propose a received signal strength indicator (RSSI) range-based 3D localization method for wireless sensor networks, where they formulate a degree of co-planarity and combine it with the traditional DV-Distance algorithm. Then, the set of four anchor nodes with the best degree of coplanarity (the flattest) is selected to calculate the node-of-interest location using simple quadrilateration. After a node calculates its position, it is promoted as anchor node assistant to improve the positioning coverage. Finally, a quasi-Newton method is used to optimize the positioning result.
In ref. [15], the authors propose a sequence localization algorithm based on the use of a 3D Voronoi diagram. The Voronoi diagram is used to divide the working space depending on the real beacon nodes (known location nodes) and construct a rank sequence table of virtual beacon nodes (determined after partition). This algorithm uses RSSI estimations between nodes as a reference to correct the measured distances and to fix the location sequence of unknown nodes. Then, the algorithm commands a weighted location estimate with N valid virtual beacon nodes by a normalization process of rank correlation coefficients.
In a different approach, the authors of [7] propose a method that is based on the convenient deployment of access points (APs) with known coordinates and on the Manhattanization of the network space. A set of algebraic relationships is defined for the node's coordinates based on distances from APs to the node of interest. In addition, a maximum a posteriori probability (MAP) criterion is used to improve estimated distances and to improve the performance of the algebraic method. This method is evaluated in multi-hop scenarios, and the distance between nodes is estimated based on the hop count criterion [16]. The main contribution is a new method for location estimation suitable for 3D network spaces. Table 1 summarizes some characteristics of the 3D localization methods explained above, as well as some evaluation performance and scenarios setup. Table 1. Summary of outdoor 3D PL methods.

3D PL Method Principle Accuracy
A 3D ad hoc localization system using aerial sensor nodes [13].
Traditional triangulation using UWB equipped aerial nodes.
In terms of standard deviation σ. They achieve 0.24 m σ. An improved 3D localization algorithm for the wireless sensor network [14].
Traditional quadrilateration combined with an optimization method.
Localization error is normalized to r 0 . Achieving error of 0.35 r 0 for 17 ANs.
3D localization algorithm based on Voronoi diagram and rank sequence in wireless sensor network [15].
Sequence localization using a Voronoi table to divide space. Then, construct a rank sequence for virtual beacon nodes.
Localization error is normalized to r 0 . Achieving 0.12 r 0 for 8 ANs.
Algebraic Manhattan restrictions.
In terms of mean Euclidean distance error, normalized to the Tx range r 0 . They achieve 0.8 r 0 for 4 anchors.
In ref. [17], the authors present a survey of multidimensional scaling (MDS) based localization techniques. These are range-free localization techniques that use proximity information instead of actual distance measurements. MDS keeps the important information while reduces multidimensional data into a lower-dimensional space. In this paper, the authors present an overview of the different variants of outdoor and indoor MDS-based localization systems. They also present some 3D MDS-based localization methods for WSN and discuss how some proposed methods for 2D are not effective in 3D environments, as complexity increases due to several factors, such as node deployment, increment of anchors needed for location calculation, connectivity, and physical obstructions.
As explained before, in this paper, we present a different positioning approach based on distance relationships and on a likelihood function. The relational positioning method formulates a searching equation based on a likelihood function that depends on a weights proportion where the points closest to the node to be located will have greater weight according to an assignment based on the ratio of distance relationships between ANs and the node of interest.
For the sake of simplicity, we describe the method both for planar and 3D scenarios in the following sections. Note that in this paper, we present the conceptualization of a new 3D localization method and its evaluation in different network scenarios. However, concerns about the implementation of the method in terms of the network layer perspective (e.g., control messages, message passing, etc.) are not necessary to assess the performance of the proposed algorithm.

Relational Positioning Method
In the current positioning system, we propose to determine the location of a node n 0 in a given area that is assumed to be a subset of a region R, defined by the convex hull contained by a number of ANs. Although our main concern is with 3D localization, for explanatory purposes, we initially consider a planar 2D scenario in which we discretize the region according to an immaterial regular square lattice with a basic separation unit L. We also assume that the locations of ANs and the node to be located coincide with some vertices v in the lattice. In the absence of further information, the uncertainty of the location of the node equates to log(N c ), where N c denotes the number of cross points in the lattice, and the node of interest is considered to be equally likely to be located at any of those cross points of the lattice. The proposed method allows us to reduce the location uncertainty through likelihood proximity estimates.

Two-Dimensional (2D) Localization
For explanatory purposes, we first consider a planar scenario, and later the method is extended to a 3D scenario. A vertex in the lattice at cross point (j,k), with coordinates (x v j,k , y v j,k ), is denoted as v j,k and it is meant to be δ i,(j,k) Manhattan units apart from AN i with coordinates (x i , y i ). The Manhattan distance between AN i and vertex v j,k is defined where L is the Manhattan unit. Thus, the node of interest n 0 is assumed to be at the unknown vertex v j,k in the discretized region, and it will have an associated set {δ i,(j,k) }, i = 1, 2, 3, of Manhattan distances [18] to the corresponding AN i , (see Figure 1). This assumption is not necessarily true in real environments; however, the searching method perceives this as distance estimation noise and makes a search of the most likely vertex v j,k for location estimation of node n 0 . This process is explained in detail in this section, and it is further developed when the 3D localization method is explained. In order to construct a location likelihood matrix ρ = [ j,k W = [wjk] , where we assign an initial weight wj,k = 1 to each convex hull containing the ANi, and otherwise, wj,k = 0, me Actual Manhattan distances from ANs to n 0 at vertex v j,k , may not be available. Although in an ad hoc network, they can be estimated by multiplying the mean length of a hop and the hop count in the shortest path linking the concerned node and the AN [16]. Estimated distances from AN i to n 0 , ∆ i,(j,k) , may depart from the corresponding Manhattan distances δ i, (j,k) ; this is ∆ i,(j,k) = δ i,(j,k) + e i, (j,k) , where the term e i,(j,k) accounts for an error associated with the estimation method, and it is often related to the node density and routing algorithm [16]. Note that the actual Manhattan distance between any pair of anchors AN i (i = 1,2,3) is determined from the known AN i (i = 1, 2, 3) coordinates, as stated above.
In order to construct a location likelihood matrix ρ = [ρ j,k ], we define an initial matrix W = [w jk ], where we assign an initial weight w j,k = 1 to each vertex v j,k lying inside the convex hull containing the AN i , and otherwise, w j,k = 0, meaning that each vertex can equally be a candidate for the node n 0 location. Weights w j,k will iteratively be changed, allowing to correspondingly define a location likelihood matrix The weights w j,k can further be modified, bearing in mind the relative distances to different AN nodes. For instance, let ∆ i,(j,k) and ∆ m,(j,k) be estimated distances from n 0 to AN i and AN m , respectively. As many interference elements such as physical barriers, multipath phenomena, moving objects and measurement errors may impair distance estimates, our approach is based on the hop count criterion; that is, ∆ i , (j,k) is proportionally approximated to the number of hops in the route connecting AN i to node n 0 assumed to be at vertex v j,k . Although we introduce some quantizing errors, making distance estimation dependent upon the hop count criterion leads to the proposed method becoming more robust to propagation impairments than the previously mentioned methods.
We define an ANs pair (i,m) proximity ratio as Then, for those vertices in the region of interest (that is, for w i,k = 0), their weights w i,k are incremented according to w i,k = w i,k + 1, for those locations (j, k) where ζ i,m ≤ 1. In this form, the weights are heavier for those nodes in closer proximity to AN i , and likelihood indexes ρ j,k are accordingly modified, denoting the more likely location for n 0 . With these likelihoods ρ j,k , we can use entropy to measure the uncertainty of a vertex, and thus, we can focus on those vertices with less uncertainty. Entropy for the area/space discretized according to L can be expressed in this case as follows: Note that the likelihood ρ j,k assignment process is conducted a number of times given by the combinations N C 2 for all the pairs of ANs (where N is the number of available ANs). As explained before, the proposed position location method diminishes the uncertainty or entropy by modifying the location likelihood of the nodes in the feasibility space using distance relationships.
In Up to this stage, we worked with three ANs in order to keep the infrastructure requirements as low as possible. However, note that the feasibility vertex set can be reduced if more ANs are available. In Figure 3, we increase the number of ANs in the scenario illustrated in Figure 2, and it is possible to see that by applying an additional relation as ζ1,4 ≤ 1, the feasible vertex set is reduced. Recall that j,k, is the probability (likelihood) of vertex vj,k, to be the position of n0, where vertex probabilities ρj,k, change after every relationship index ζi,m is considered for weight assignment; this implies that the entropy function HL(v) also changes after every stage. In Figure 4, we can see how entropy behaves after every ζi,m is applied for the representative example explained in Figures 2 and 3. Up to this stage, we worked with three ANs in order to keep the infrastructure requirements as low as possible. However, note that the feasibility vertex set can be reduced if more ANs are available. In Figure 3, we increase the number of ANs in the scenario illustrated in Figure 2, and it is possible to see that by applying an additional relation as ζ 1,4 ≤ 1, the feasible vertex set is reduced. Up to this stage, we worked with three ANs in order to keep the infrast requirements as low as possible. However, note that the feasibility vertex set reduced if more ANs are available. In Figure 3, we increase the number of ANs scenario illustrated in Figure 2, and it is possible to see that by applying an add relation as ζ1,4 ≤ 1, the feasible vertex set is reduced. Recall that j,k, is the probability (likelihood) of vertex vj,k, to be the positio where vertex probabilities ρj,k, change after every relationship index ζi,m is conside weight assignment; this implies that the entropy function HL(v) also changes afte stage. In Figure 4, we can see how entropy behaves after every ζi,m is applied representative example explained in Figures 2 and 3. Recall that ρ j,k, is the probability (likelihood) of vertex v j,k, to be the position of n 0 , where vertex probabilities ρ j,k, change after every relationship index ζ i,m is considered for weight assignment; this implies that the entropy function H L (v) also changes after every stage. In Figure 4, we can see how entropy behaves after every ζ i,m is applied for the representative example explained in Figures 2 and 3.

Three-Dimensional (3D) Localization
In the case of 3D scenarios, the node to be located is meant to be within the 3D convex hull volume V defined by the ANs. For instance, in the case of four equidistant Ans, the convex hull coincides with a regular tetrahedron, as is illustrated in Figure 5a. As explained before, a discretization of the searching space is conducted with vertices vj,k,l in a three-dimensional regular grid with step L. The weight wj,k,l assignment process is extended to 3D, where (i,k,l) stands for the coordinates in the searching space, though it remains based on distance relationships ζi,m. As an example, consider the case in which n0 is contained within a regular tetrahedron with 4C2 = 6 (where 4C2 is the combinatorial coefficient) relative distance comparisons.
It is anticipated that performance improvements can be achieved when the number of ANs is increased in the same volume or area as more distance relationships, e.g., 5C2, from ANs to n0 can be used. Therefore, we obtain a different weight assignment from that obtained for four ANs in the tetrahedron. As stated before, a discretization process is performed first. Secondly, the weight assignment process is carried out according to distance relationships ζi,m. Finally, the searching equation as that shown in (2) can be reformulated based on the distance sets Δi, (j,k,l), and δi,(j,k,l), where sub-indices 1 ≤ i ≤ m refer

Three-Dimensional (3D) Localization
In the case of 3D scenarios, the node to be located is meant to be within the 3D convex hull volume V defined by the ANs. For instance, in the case of four equidistant Ans, the convex hull coincides with a regular tetrahedron, as is illustrated in Figure 5a.

Three-Dimensional (3D) Localization
In the case of 3D scenarios, the node to be located is meant to be within the 3D convex hull volume V defined by the ANs. For instance, in the case of four equidistant Ans, the convex hull coincides with a regular tetrahedron, as is illustrated in Figure 5a. As explained before, a discretization of the searching space is conducted with vertices vj,k,l in a three-dimensional regular grid with step L. The weight wj,k,l assignment process is extended to 3D, where (i,k,l) stands for the coordinates in the searching space, though it remains based on distance relationships ζi,m. As an example, consider the case in which n0 is contained within a regular tetrahedron with 4C2 = 6 (where 4C2 is the combinatorial coefficient) relative distance comparisons.
It is anticipated that performance improvements can be achieved when the number of ANs is increased in the same volume or area as more distance relationships, e.g., 5C2, from ANs to n0 can be used. Therefore, we obtain a different weight assignment from that obtained for four ANs in the tetrahedron. As stated before, a discretization process is performed first. Secondly, the weight assignment process is carried out according to distance relationships ζi,m. Finally, the searching equation as that shown in (2) can be reformulated based on the distance sets Δi, (j,k,l), and δi,(j,k,l), where sub-indices 1 ≤ i ≤ m refer As explained before, a discretization of the searching space is conducted with vertices v j,k,l in a three-dimensional regular grid with step L. The weight w j,k,l assignment process is extended to 3D, where (i,k,l) stands for the coordinates in the searching space, though it remains based on distance relationships ζ i,m . As an example, consider the case in which n 0 is contained within a regular tetrahedron with 4 C 2 = 6 (where 4 C 2 is the combinatorial coefficient) relative distance comparisons.
It is anticipated that performance improvements can be achieved when the number of ANs is increased in the same volume or area as more distance relationships, e.g., 5 C 2 , from ANs to n 0 can be used. Therefore, we obtain a different weight assignment from that obtained for four ANs in the tetrahedron. As stated before, a discretization process is performed first. Secondly, the weight assignment process is carried out according to distance relationships ζ i,m . Finally, the searching equation as that shown in (2) can be reformulated based on the distance sets ∆ i, (j,k,l) , and δ i,(j,k,l) , where sub-indices 1 ≤ i ≤ m refer to the corresponding AN i , and (j,k,l) refers to the coordinate (x, y, z) in a tridimensional space.
As stated before, uncertainty reduction helps to decrease the searching space, allowing lower processing time, required processing capacity, and power consumption. The final location will be determined by an optimal searching process over the feasible region R* = {v j,k,l , w j,k,l =0} where a minimum-error searching equation is formulated as follows: where ρ j,k,l = w j,k,l / ∑ j,k,l w j,k,l and 1 ≥ i ≥ 4.
A vertex v j,k,l will be the solution to the localization problem, if such vertex gives the minimum difference between its associated distance set δ i,(j,k,l) and the estimated distances from ANs to node n 0 , ∆ i,(j,k,l) , multiplied by the inverse of ρ j,k,l = w j,k,l / ∑ j,k,l w j,k,l .
Note that the region of interest R can have different shapes and sizes depending on the available ANs and their locations. Additionally, it is worth noting that the inclusion of additional ANs can improve localization accuracy. In the following section, we deduce that uncertainty is reduced significantly when the proposed searching method is combined with a finer discretization which helps to achieve significant improvements. We also evaluate how the weight assignment process can reduce uncertainty. Additionally, the performance of the proposed method is assessed under different scenarios.
We summarize the 3D relational positioning method with the following set of steps: 1.
Discretize the working operational area or space, obtaining v j,k,l ; 2.
Calculate the Manhattan distance from AN i to v j,k,l, δ i,(j,k,l) ; 3.
Obtain the route from ANs to the node of interest n 0 ; 4.
Calculate distance from AN i to n 0 based on the number of hops in the route ∆ i,(j,k,l) ; 5.
Modify vertex weight w j,k,l for every pair of ANs, as follows: If ζ i,m < 1, then the node is closer to AN i than to AN m , hence every vertex weight w j,k,l closer to AN i is increased by 1. Otherwise, increase by 1 every vertex weight w j,k,l closer to AN m ;

Simulations and Results
In this section, we evaluate the performance of the proposed method based on the average quadratic error criterion. Note that influences from protocols of the network layer or other layers are not considered in the evaluation, because results may be affected and such protocols change for different technologies, and therefore, we would have no clear trends of the performance of the method. First, we evaluate the PL relational searching method under two-dimensional network scenarios with three and four ANs.

The 2D Evaluation
To evaluate the performance of the proposed PL searching method, and without loss of generality, we set ANs' locations in an equilateral triangle defined by ANs coordinates as follows: AN 1 (0,0), AN 2 (d,0), AN 3 (d cos π 3 , d sin π 3 ), when three ANs are available (as in Figure 1), and as in a regular square AN 1 (0,0), AN 2 (d,0), AN 3 (d,d), AN 4 (0,d), when four ANs are available. Nodes are deployed randomly according to a uniform pdf in the working area or space; that is, nodes' coordinates are selected as x~U(0,d) and y~U(0,d). The node to be located is also randomly selected.
In the following analysis, we evaluate the average quadratic error of the Euclidean distance between the estimated position and the true position ε, and we normalize it to the AN separation distance d, which is related to ANs' coordinates and defines the size of the network scenario, i.e., ε/d. Error normalization is conducted in order to relate results under different scenarios (2D and 3D) in terms of the same scenario constructive variable d (associated with ANs coordinates). This is carried out in order to keep a fair comparison of performance between different scenarios.
In the following simulations, we analyze the positioning method under different network conditions. For the 2D scenarios, we perform a basic analysis considering different values of the coverage radius r 0 and number of nodes n , in order to prove the feasibility of the method for 2D networks, and we also show that location estimation is improved when the number of ANs is increased.
For instance, in Figure 6, we can observe the performance of the relational searching method for two different 2D scenarios with the same working area, where three and four ANs are considered. In this figure, we can see that a consistent improvement from about 7% to 18% is achieved when we increment the number of ANs from three to four. Additionally, we can see a slight improvement on the algorithm performance as we increase the number of nodes from n = 30 to n = 40, which can be explained by an expected improvement on routes obtained from ANs to the node to locate. As the hop count method is used to estimate distances, better routes may imply better distance estimations. Note that there is an increment of the normalized positioning error ε/d when the maximal transmission range r 0 is increased from r 0 = 2.5 m to r 0 = 5 m. This can also be explained because worse distance estimations from ANs to the node to be located are obtained (as the hop count distance estimation is selected) when r 0 is increased. The performance obtained for r 0 = 2 m, which is worse than the performance for r 0 = 2.5 m can be explained due to the low node density.
In the following analysis, we evaluate the average quadratic error of the Euclidean distance between the estimated position and the true position ε, and we normalize it to the AN separation distance d, which is related to ANs' coordinates and defines the size of the network scenario, i.e., ε/d. Error normalization is conducted in order to relate results under different scenarios (2D and 3D) in terms of the same scenario constructive variable d (associated with ANs coordinates). This is carried out in order to keep a fair comparison of performance between different scenarios.
In the following simulations, we analyze the positioning method under different network conditions. For the 2D scenarios, we perform a basic analysis considering different values of the coverage radius r0 and number of nodes n, in order to prove the feasibility of the method for 2D networks, and we also show that location estimation is improved when the number of ANs is increased.
For instance, in Figure 6, we can observe the performance of the relational searching method for two different 2D scenarios with the same working area, where three and four ANs are considered. In this figure, we can see that a consistent improvement from about 7% to 18% is achieved when we increment the number of ANs from three to four. Additionally, we can see a slight improvement on the algorithm performance as we increase the number of nodes from n = 30 to n = 40, which can be explained by an expected improvement on routes obtained from ANs to the node to locate. As the hop count method is used to estimate distances, better routes may imply better distance estimations. Note that there is an increment of the normalized positioning error ε/d when the maximal transmission range r0 is increased from r0 = 2.5 m to r0 = 5 m. This can also be explained because worse distance estimations from ANs to the node to be located are obtained (as the hop count distance estimation is selected) when r0 is increased. The performance obtained for r0 = 2 m, which is worse than the performance for r0 = 2.5 m can be explained due to the low node density. In order to consider more realistic scenarios than those evaluated in Figure 6, the AN3 coordinates are modified. Note that when changing x3, irregular triangle scenarios are generated. In Table 2 we can observe these results in terms of the normalized localization error ε/d, for different maximal transmission range r0, n = 40, and different x3 values. In Table 2, we can infer that the location estimation given by ε/d stays stable as the analyzed  In order to consider more realistic scenarios than those evaluated in Figure 6, the AN 3 coordinates are modified. Note that when changing x 3 , irregular triangle scenarios are generated. In Table 2 we can observe these results in terms of the normalized localization error ε/d, for different maximal transmission range r 0, n = 40, and different x 3 values. In Table 2, we can infer that the location estimation given by ε/d stays stable as the analyzed shape scenario is changed. Even more, we can find some cases in which the normalized error decreases, as better relationships might be found.

The 3D Evaluation
Let us consider a 3D scenario defined by a regular tetrahedron, as illustrated in Figure 5a, having four ANs with coordinates AN 1 (0,0,0), AN 2 (d,0,0), AN 3 (d cos π 3 , d sin π 3 ,0), and AN 4 ( d 2 , d 2 tan π 6 ,d). Routes from ANs to node n 0 are obtained according to Dijkstra's routing algorithm [19]. As stated before, the estimated distance between AN i and node n 0 , ∆ i, (j,k,l) is approximated by the number of hops in the route connecting them. In addition, the maximum coverage radius is initially set to r 0 = 2 m, and the number of nodes in the network is set as follows: for the 2D scenario, we randomly deploy 70 to 90 nodes in a dxd square meter area, and as the working area is a triangle defined by ANs coordinates, the proportional number of nodes within the triangle is n = 30 to 40 nodes. For the 3D scenario, 300 to 400 nodes are also randomly deployed in a dxdxd cube meter volume, and as the working space is a tetrahedron defined by ANs coordinates, the proportional number of nodes in such tetrahedron is n = 43 to 57.
For the 3D assessment, we extended the variety of scenarios analyzed, varying more network conditions (ANs locations, coverage radius r 0 , and the number of nodes) and under different values of the space discretization unit L. Note that L defines the vertices v j,k,l on the discretized space and thus the searching points related to distances δ i,(j,k,l) , in Equation (2). Figure 7 shows the performance of the relational searching method for a 3D scenario, as shown in Figure 5a, with AN 1 (0,0,0), AN 2 (d,0,0), AN 3 (d cos π 3 , d sin π 3 ,0) and AN 4 ( d 2 , d 2 tan π 6 ,d), where d = 10 m. In this figure, we analyzed the performance of the method under different coverage ranges, r 0 , and the number of nodes n in the working space. In this analysis, we can observe that as r 0 increases, the normalized error ε/d increases; however, the Euclidean quadratic error ε shows stability, as it remains close to 2 m. In this figure, we can also observe that by increasing the number of nodes in the space, we can achieve a significant improvement for r 0 = 2 m. This is because better routes can be found if more nodes are in the network, and this is more significant when r 0 is small. Even that the transmission range r 0 is selected in the range r 0 = [2 m 5 m] for the following evaluations, it can be larger than 5 m; however, the normalized localization error is very close to that obtained for r 0 = 5 m. Simulations were performed for r 0 = 6 m, 7 m, 8 m and n = 43 nodes, obtaining ε/d = 0.2236, 0.2235 and 0.2526, showing error saturation, which is a desirable behavior.
In Figure 8, results from a different perspective are shown; on the x-axis, the average Euclidean quadratic distance error ε (m) is presented, and on the y-axis, the values of the cumulative distribution function describing the percentage of cases when ε is under a given distance are presented. In Figure 8, results from a different perspective are shown; on the x-axis, the average Euclidean quadratic distance error ε (m) is presented, and on the y-axis, the values of the cumulative distribution function describing the percentage of cases when ε is under a given distance are presented.
Simulations in Figure 8 were performed for a 3D scenario with four ANs, (n = 43, L = 1). Note that results in Figure 8 are not normalized to d. In this analysis, we can observe that for r0 = 2.5 m, approximately 100% of the calculations are under 3 m of error. The worst analyzed case is when r0 = 2 m, which can be explained because distance estimations from anchors to the n0 are based on the hop count, and for a very small r0, we achieve longer routes, resulting in bad distance estimations. However, this can be alleviated by increasing the number of nodes in the space, as previously observed in Figure 7. From these results, we can infer that the best hop distance estimation is achieved for r0 = 2.5 m. We can also infer that for r0 < 2.5 m, the distance might be overestimated, and for r0 > 2.5 m, the distance might be underestimated.  In Figure 8, results from a different perspective are shown; on the x-axis, the average Euclidean quadratic distance error ε (m) is presented, and on the y-axis, the values of the cumulative distribution function describing the percentage of cases when ε is under a given distance are presented.
Simulations in Figure 8 were performed for a 3D scenario with four ANs, (n = 43, L = 1). Note that results in Figure 8 are not normalized to d. In this analysis, we can observe that for r0 = 2.5 m, approximately 100% of the calculations are under 3 m of error. The worst analyzed case is when r0 = 2 m, which can be explained because distance estimations from anchors to the n0 are based on the hop count, and for a very small r0, we achieve longer routes, resulting in bad distance estimations. However, this can be alleviated by increasing the number of nodes in the space, as previously observed in Figure 7. From these results, we can infer that the best hop distance estimation is achieved for r0 = 2.5 m. We can also infer that for r0 < 2.5 m, the distance might be overestimated, and for r0 > 2.5 m, the distance might be underestimated.  Figure 8 were performed for a 3D scenario with four ANs, (n = 43, L = 1). Note that results in Figure 8 are not normalized to d. In this analysis, we can observe that for r 0 = 2.5 m, approximately 100% of the calculations are under 3 m of error. The worst analyzed case is when r 0 = 2 m, which can be explained because distance estimations from anchors to the n 0 are based on the hop count, and for a very small r 0 , we achieve longer routes, resulting in bad distance estimations. However, this can be alleviated by increasing the number of nodes in the space, as previously observed in Figure 7. From these results, we can infer that the best hop distance estimation is achieved for r 0 = 2.5 m. We can also infer that for r 0 < 2.5 m, the distance might be overestimated, and for r 0 > 2.5 m, the distance might be underestimated.

Simulations in
In practice, it is difficult to guarantee uniform deployment of ANs, as in Figure 5a. Therefore, evaluations were carried out for regular and non-regular scenarios. Simulations were conducted for different shapes of the tetrahedron base, and the results obtained were similar to those reported in Figures 7 and 8. These results are presented in Table 3. Note that moving only AN 3 coordinates gives enough variability to the scenario setup, as the tetrahedron's base is changed. In these simulations, coordinate x 3 was varied, and coordinate y 3 was fixed in order to keep the same tetrahedron volume to keep the comparison fair. Simulation parameters were set as n = 43, d = 10 m, AN 1 (0,0,0), AN 2 (d,0,0), AN 3 (x 3 , d sin π 3 , 0), AN 4 ( d 2 , d 2 tan π 6 , h), and h = 10 m. However, when the tetrahedron height h is increased, we observe a slight degradation of the proposed PL method, which is due to the increased size of the searching space. In Table 4, these results are presented for different transmission ranges, n = 43, d = 10 m and AN 1 (0,0,0), AN 2 (d,0,0), AN 3 (d cos π 3 , d sin π 3 , 0), AN 4 ( d 2 , d 2 tan π 6 , h). Note that the proposed method shows robustness under volume increments; that is, for a volume increment of 30% (from h = 10 m to h = 13 m), the normalized error ε/d increment ranges from 3.5% to 17%. In Figure 9, we analyzed the localization method for different values of the discretization unit L. As we anticipated, a finer discretization grid helps to reduce the normalized error ε/d significantly. In this figure, we can observe that when the discretization unit is decreased to L = 1 2 m, improvements ranging from 39% to 45% can be achieved. As we continue to decrease L, ε/d diminishes in a lower proportion. With this analysis, we demonstrate that a finer discretization leads to a lower distance error (close to 1 m) for L = 1 2 m, and to 0.8 m for L = 1 3 m. However, as we decreased the discretization step L, we increased the number of vertices in the searching space, increasing the uncertainty. In Table 5, we analyzed the entropy function for the 3D scenario illustrated in Figure 5a, for L = 1 m and L = 1 2 m. The entropy values for L = 1 m and L = 1/2 m without the proposed method on the discretized space are H L=1 = 5.04 and H L=1/2 = 6.9735, respectively. Note that the method helps to significantly decrease the uncertainty.
In Figure 10, we present the performance comparison for a 3D scenario with four to seven ANs. This analysis was carried out for n = 43 randomly deployed nodes in the tetrahedron illustrated in Figure 5b, with AN 1 to AN 4 deployed as stated before and AN 5 , AN 6 , AN 7 deployed at the centroid of planes conformed by the corresponding sets (AN 1 , AN 2 , AN 3 ), (AN 2 , AN 3 , AN 4 ), and (AN 1 , AN 2 , AN 4 ), with d = 10. As anticipated, when we increase from 4 ANs to 5 ANs, we can see improvements from approximately 5% to 12% for all the cases reported. Additionally, when we increase from 5 ANs to 6 ANs, improvements from 2% to 6.6% are achieved. Finally, when we increase from 6 ANs to 7 ANs, improvements from 1.6% to 2.84% are observed. Simulations were performed by moving AN 5 on the plane conformed by the set (AN 1 , AN 2 , AN 3 ), and it was observed that the achieved improvement decreases as the AN 5 approaches AN 1 , AN 2 , or AN 3 . A similar experiment was conducted regarding AN 6 and AN 7 , and similar behavior was observed. For these scenarios, the entropy was also calculated, and results were very similar to those reported in Table 4 for L = 1. Note that when the number of ANs increases, there is an improvement in the location. However, for a given reachability ratio r 0 , the further increase in the number of ANs provides marginal improvement. However, as we decreased the discretization step L, we increased the number of vertices in the searching space, increasing the uncertainty. In Table 5, we analyzed the entropy function for the 3D scenario illustrated in Figure 5a, for L = 1 m and L = m. The entropy values for L = 1 m and L = 1/2 m without the proposed method on the discretized space are HL=1 = 5.04 and HL=1/2 = 6.9735, respectively. Note that the method helps to significantly decrease the uncertainty. In Figure 10, we present the performance comparison for a 3D scenario with four to seven ANs. This analysis was carried out for n = 43 randomly deployed nodes in the tetrahedron illustrated in Figure 5b, with AN1 to AN4 deployed as stated before and AN5, AN6, AN7 deployed at the centroid of planes conformed by the corresponding sets (AN1, AN2, AN3), (AN2, AN3, AN4), and (AN1, AN2, AN4), with d = 10. As anticipated, when we increase from 4 ANs to 5 ANs, we can see improvements from approximately 5% to 12% for all the cases reported. Additionally, when we increase from 5 ANs to 6 ANs, improvements from 2% to 6.6% are achieved. Finally, when we increase from 6 ANs to 7 ANs, improvements from 1.6% to 2.84% are observed. Simulations were performed by moving AN5 on the plane conformed by the set (AN1, AN2, AN3), and it was observed that the achieved improvement decreases as the AN5 approaches AN1, AN2, or AN3. A similar experiment was conducted regarding AN6 and AN7, and similar behavior was observed. For these scenarios, the entropy was also calculated, and results were very similar to those reported in Table 4 for L = 1. Note that when the number of ANs increases, there is an improvement in the location. However, for a given reachability ratio r0, the further increase in the number of ANs provides marginal improvement.   Thus far, we have focused on evaluations for indoor scenarios, which are frequently used on Industry 4.0 applications. However, there are some applications in smart industries and sensor networks in which outdoor environments have to be considered. Since the proposed relational positioning method does not rely on range estimation, there is no need to consider signal impairments associated with outdoor scenarios. However, the scenario size d, as well as the transmission range r0, have to be increased in order to assess outdoor networks. As shown in Figure 11, we evaluated the performance of the relational positioning method for d = 20 m, 30 m, and for transmission ranges r0 from 6 m to 12 m for the basic set up, n = 43, L = 1, and regular deployment of ANs (AN1(0,0,0), Thus far, we have focused on evaluations for indoor scenarios, which are frequently used on Industry 4.0 applications. However, there are some applications in smart industries and sensor networks in which outdoor environments have to be considered. Since the proposed relational positioning method does not rely on range estimation, there is no need to consider signal impairments associated with outdoor scenarios. However, the scenario size d, as well as the transmission range r0, have to be increased in order to assess outdoor networks. As shown in Figure 11, we evaluated the performance of the relational positioning method for d = 20 m, 30 m, and for transmission ranges r 0 from 6 m to 12 m for the basic set up, n = 43, L = 1, and regular deployment of ANs (AN 1 (0,0,0), AN 2 (d,0,0), AN 3 (d cos π 3 , d sin π 3 ,0) and AN 4 ( d 2 , d 2 tan π 6 ,d)). Results in Figure 11 show that the estimation error does not grow in the same proportion as the size of the network d, as maximal estimation errors for d = 20 m and d = 30 m are in the order of 4.1 m and 4.9 m, respectively, for r 0 = 12 m. Note that as stated before, finer discretization, as well as additional ANs, can help to improve the location estimation. Thus far, we have focused on evaluations for indoor scenarios, which are frequently used on Industry 4.0 applications. However, there are some applications in smart industries and sensor networks in which outdoor environments have to be considered. Since the proposed relational positioning method does not rely on range estimation, there is no need to consider signal impairments associated with outdoor scenarios. However, the scenario size d, as well as the transmission range r0, have to be increased in order to assess outdoor networks. As shown in Figure 11, we evaluated the performance of the relational positioning method for d = 20 m, 30 m, and for transmission ranges r0 from 6 m to 12 m for the basic set up, n = 43, L = 1, and regular deployment of ANs (AN1(0,0,0), AN2(d,0,0), AN3(d cos , d sin ,0) and AN4( , tan ,d)). Results in Figure 11 show that the estimation error does not grow in the same proportion as the size of the network d, as maximal estimation errors for d = 20 m and d = 30 m are in the order of 4.1 m and 4.9 m, respectively, for r0 = 12 m. Note that as stated before, finer discretization, as well as additional ANs, can help to improve the location estimation.  It is important to observe that as the method relies on the routes connecting ANs to n 0 , better route algorithms might help to obtain better positioning estimations. Note that proactive routing protocols might render better positioning estimations, as routes are updated periodically. However, reactive (on-demand) routing algorithms might lead to large positioning errors if the location estimation is calculated based on old routing data.

Conclusions
A position location-searching method suitable for both 2D and 3D sensor networks and based on relative distance relationships was proposed. This method can be used for indoor and outdoor scenarios, and it was demonstrated that its accuracy can be improved without the need to increase the infrastructure requirements. Since in Industry 4.0, traditional indoor and outdoor localization methods are not suitable due to the lack of reliability of distance estimation, the proposed method provides a reliable alternative to solve the need for position location information. The proposed method is robust, as it does not rely on the goodness of estimated distances. Its performance was assessed for multi-hop environments using a low-cost distance estimation method (hop count) considering several ANs for 3D networks. The 3D PL searching method is based on the space discretization once the searching space has been defined and enclosed by a set of ANs and on a vertex weight assignment due to distance relations. It was shown through simulations that the PL searching method improves as we increase the number of ANs in the operational area/space. However, using only four ANs combined with a small discretization step (L = 1 3 m) is sufficient to achieve the best performance reported. Additionally, marginal improvements can be achieved when the number of nodes in the network increases, as better routes and therefore better distance estimations become available.

Conflicts of Interest:
All the authors of the present paper declare no conflict of interest.