Nazanin Vafaei

Rita A Ribeiro

Luis M Camarinha-Matos

Normalization Techniques for Multi-Criteria Decision Making: Analytical Hierarchy Process Case Study

Keywords: Normalization, AHP, MCDM, Rank Reversal, Cyber-Physical Systems (CPSs)

Multi-Criteria Decision Making (MCDM) methods use normalization techniques to allow aggregation of criteria with numerical and comparable data. With the advent of Cyber Physical Systems, where big data is collected from heterogeneous sensors and other data sources, finding a suitable normalization technique is also a challenge to enable data fusion (integration). Therefore, data fusion and aggregation of criteria are similar processes of combining values either from criteria or from sensors to obtain a common score. In this study, our aim is to discuss metrics for assessing which are the most appropriate normalization techniques in decision problems, specifically for the Analytical Hierarchy Process (AHP) multi-criteria method. AHP uses a pairwise approach to evaluate the alternatives regarding a set of criteria and then fuses (aggregation) the evaluations to determine the final ratings (scores).

Introduction

Everybody makes decisions in their daily lives, as for example: "Should I take an umbrella today"? "Where should I go for lunch"? To make decisions we need access to information (or data) and to reach a decision we need to combine the data to obtain a final score for each candidate decision alternative (e.g. combining food prices and service of restaurants to recommend). The aim of Multi-Criteria Decision Making (MCDM) methods is to rate and prioritize a set of alternatives that best satisfy a given set of criteria [START_REF] Triantaphyllou | Multi-Criteria Decision Making Methods[END_REF]. Criteria are a set of requirements or independent attributes that have to be satisfied by several alternatives. Each criterion may be measured in different units, for example, degrees, kilograms or meters; but they all have to be normalized to obtain dimensionless classifications, i.e. a common numeric range/scale, to allow aggregation into a final score. Hence, data normalization is an essential part of any decision making process because it transforms the input data into numerical and comparable data, allowing using MCDM methods to rate and rank alternatives [START_REF] Jahan | A state-of-the-art survey on the influence of normalization techniques in ranking: Improving the materials selection process in engineering design[END_REF][START_REF] Nayak | Impact of Data Normalization on Stock Index Forecasting[END_REF].

In this work, the main research question that we address is: Which normalization technique is more suitable for usage with the AHP method?

The motivation for carrying out this work includes four interconnected issues: a) the importance of data normalization for decision problems where we need to fuse or aggregate data to obtain a final score per alternative; b) the reduced number of research studies available in this topic; c) continuation of previous work on suitability of normalization techniques for well-known MCDM methods (e.g. TOPSIS) [START_REF] Vafaei | Data Normalization Techniques in Decision Making: Case Study with TOPSIS Method[END_REF]; d) contributing to advances in Cyber Physical Systems [START_REF] Camarinha-Matos | Collaborative Systems for Smart Environments: Trends and Challenges[END_REF] research, where huge amounts of available data from heterogeneous sensors need to be fused (aggregated) to determine a combined view. Specifically, in this study we focus on the well-known AHP method because it is a well-known and widely used MCDM method [START_REF] Saaty | A scaling method for priorities in hierarchical structures[END_REF][START_REF] Saaty | The Analytic Hierarchy Process[END_REF][START_REF] Cheng | Evaluating attack helicopters by AHP based on linguistic variable weight[END_REF][START_REF] Gaudenzi | Managing risks in the supply chain using the AHP method[END_REF][START_REF] Zahedi | The analytic hierarchy processa survey of the method and its applications[END_REF][START_REF] Tuzkaya | An analytic network process approach for locating undesirable facilities: An example from Istanbul, Turkey[END_REF][START_REF] Saaty | Models, Methods, Concepts & Applications of the Analytic Hierarchy Process[END_REF][START_REF] Tzeng | Multiple Attribute Desicion Making: Methods and Applications[END_REF] but we plan to perform the same study for data fusion problems as well as other MCDM methods in the future, to determine which technique is more suitable for any decision problem that requires combining (fusing) data.

The Analytic Hierarchy Process (AHP) was introduced by Saaty [START_REF] Saaty | A scaling method for priorities in hierarchical structures[END_REF][START_REF] Saaty | The Analytic Hierarchy Process[END_REF] to solve unstructured problems in economics, social sciences, and management [START_REF] Cheng | Evaluating attack helicopters by AHP based on linguistic variable weight[END_REF]. AHP has been used in a vast range of problems from simple ones (e.g. selecting a school) to harder ones (e.g. in allocating budgets and energy domains) [START_REF] Cheng | Evaluating attack helicopters by AHP based on linguistic variable weight[END_REF]. When applying the AHP method, the decision maker is able to structure the decision problem and break it down into a hierarchical top-down process. Then, he/she performs a pairwise matrix comparison of criteria using a [START_REF] Triantaphyllou | Multi-Criteria Decision Making Methods[END_REF][START_REF] Jahan | A state-of-the-art survey on the influence of normalization techniques in ranking: Improving the materials selection process in engineering design[END_REF][START_REF] Nayak | Impact of Data Normalization on Stock Index Forecasting[END_REF][START_REF] Vafaei | Data Normalization Techniques in Decision Making: Case Study with TOPSIS Method[END_REF][START_REF] Camarinha-Matos | Collaborative Systems for Smart Environments: Trends and Challenges[END_REF][START_REF] Saaty | A scaling method for priorities in hierarchical structures[END_REF][START_REF] Saaty | The Analytic Hierarchy Process[END_REF][START_REF] Cheng | Evaluating attack helicopters by AHP based on linguistic variable weight[END_REF][START_REF] Gaudenzi | Managing risks in the supply chain using the AHP method[END_REF] scale (corresponding to semantic interpretations such has "A is much more important than B" regarding a criterion). After normalization, the priorities are determined using either Eigen vectors or a simplified version with weighted sum (SAW) [START_REF] Gaudenzi | Managing risks in the supply chain using the AHP method[END_REF][START_REF] Zahedi | The analytic hierarchy processa survey of the method and its applications[END_REF].

AHP involves five main steps [START_REF] Tzeng | Multiple Attribute Desicion Making: Methods and Applications[END_REF]:

Step 1: Decompose the problem into a hierarchical structure; Step 2: Employ pairwise comparisons. A pairwise comparison is the process of comparing the relative importance, preference, or likelihood of two elements (objectives) with respect to another element (the goal). Pairwise comparisons are carried out to establish priorities. Decision elements at each hierarchy level are compared pairwisely and then the reciprocal matrix is completed;

Step 3: Determine the logical consistency and if > 10% revise the pairwise classifications until the consistency index is below 10%. In the implementation of AHP, we may face with inconsistent judgment of input data that it may cause some bad effects on decision process. For example, A1 may be preferred to A2 and A2 to A3, but A3 may be preferred to A1. So, Saaty [START_REF] Saaty | The Analytic Hierarchy Process[END_REF] defined a measure of deviation from consistency that is called called a consistency index, as: C.I. = (λmax -N)/ (N-1), where N is the dimension of the matrix and λ is the largest eigenvalue of the matrix A. Then, Saaty calculated a consistency ratio (C.R.) as the ratio of the C.I. to a random index (R.I.) which is the average C.I. of sets of judgments (from a 1 to 9 scale) for randomly generated reciprocal matrices [START_REF] Saaty | The Analytic Hierarchy Process[END_REF]. Step 4: Estimate the relative weights by combining the individual subjective judgments. We can use the eigenvalue method to estimate the relative weights of the decision elements. In order to estimate the relative weight of the decision elements in a matrix, we can use A.Ⱳ=λmax .Ⱳ where W is the weight of criterion [START_REF] Tzeng | Multiple Attribute Desicion Making: Methods and Applications[END_REF]. Step 5: Determine the priority of alternatives by doing aggregation on relative weights which is obtained by combining the criterion priorities and priorities of each decision alternatives relative to each criterion. Since in our work we discuss the suitability of normalization techniques for the AHP method, we focus on step 4 and 5.

In this work we propose an assessment approach for evaluating five common normalization techniques (see Table 1), using an illustrative example solved with AHP method [START_REF] Triantaphyllou | Multi-Criteria Decision Making Methods[END_REF][START_REF] Jahan | A state-of-the-art survey on the influence of normalization techniques in ranking: Improving the materials selection process in engineering design[END_REF]. We choose AHP because it is a well-known and widely used MCDM method [START_REF] Saaty | A scaling method for priorities in hierarchical structures[END_REF][START_REF] Saaty | The Analytic Hierarchy Process[END_REF][START_REF] Cheng | Evaluating attack helicopters by AHP based on linguistic variable weight[END_REF][START_REF] Gaudenzi | Managing risks in the supply chain using the AHP method[END_REF][START_REF] Zahedi | The analytic hierarchy processa survey of the method and its applications[END_REF][START_REF] Tuzkaya | An analytic network process approach for locating undesirable facilities: An example from Istanbul, Turkey[END_REF][START_REF] Saaty | Models, Methods, Concepts & Applications of the Analytic Hierarchy Process[END_REF][START_REF] Tzeng | Multiple Attribute Desicion Making: Methods and Applications[END_REF] but we plan to perform the same study for other MCDM methods in the future. Our novel assessment approach calculating Pearson correlation for global weight of alternatives and Spearman correlation for rank of alternatives which are borrowed from [START_REF] Celen | Comparative Analysis of Normalization Procedures in TOPSIS Method: With an Application to Turkish Deposit Banking Market[END_REF] to determine mean values in order to ensure a more robust evaluation and selection of the best normalization technique in AHP. The novelty of this study is making adaptation between assessment process and AHP in order to find best normalization technique for AHP method. The next section presents the experimental study performed.

Relationship to Cyber-Physical Systems

Cyber-physical systems (CPS) involve merging computation and physical processes, often denoted as embedded systems [START_REF] Shi | A survey of Cyber-Physical Systems[END_REF]. In most CPS, the physical inputs and outputs are typically designed as a network of interacting elements. This conceptual model is tied to the notion of robotics and sensor networks and their usage has been increasing day by day [START_REF]Cyber-physical system[END_REF]. But CPS also inherits ideas from the areas of embedded and real-time systems. CPS have a broad scope of potential application in areas such as reactive interventions (e.g., collision avoidance); precision operations (e.g., robotic surgery and nano-level manufacturing); operation in dangerous or inaccessible environments (e.g., search and rescue, firefighting, and deep-sea exploration); complex systems coordination (e.g., air traffic control, war fighting); efficiency (e.g., zero-net energy buildings); and augmentation of human capabilities (e.g., healthcare monitoring and service delivery) [START_REF]Cyber-physical system[END_REF], to name a few.

There are some discussions on the relationship between Cyber-Physical Systems and Internet of Things [START_REF] Camarinha-Matos | Technological Innovation for the Internet of Things[END_REF][START_REF] Jeschke | Everything 4.0? in: Drivers and Challenges of Cyber Physical Systems[END_REF][START_REF]Normalization (statistics)[END_REF]. Camarinha and Afsarmanesh [START_REF] Camarinha-Matos | Collaborative Systems for Smart Environments: Trends and Challenges[END_REF] mention that "there is a growing convergence between the two areas since CPSs are becoming more Internet-based". For example, in smart car parking, data from the parking space is transferred to the car drivers with the help of CPS and IoT technologies. Data is collected from sensors, which are installed in the parking lot, and transferred to the data center to be processed with MCDM methods, to determine the ranking of alternatives (best parking spaces). The best parking spaces are provided to the car drivers to support them making more informed decisions. In the illustrative example section, we will compare several normalization techniques for usage with the AHP method to rank alternatives and support car drivers. The smart car parking example shows a robust relationship between cyber physical system (CPS), Internet of Thing (IoT) and multi-criteria decision making (MCDM) concepts.

Normalization

There are several definitions for data normalization, depending on the study domain. For example, in Databases, data normalization is viewed as a process where data attributes, within a data model, are organized in tables to increase the cohesion and efficiency of managing data. In statistics and its applications, the most common definition is the process of adjusting values measured on different scales to a common scale, often prior to aggregating or averaging them [START_REF]Normalization (statistics)[END_REF]. Many other definitions exist, depending on the context or study domain (see for example [START_REF]Normalization[END_REF]). Here we focus on normalization techniques for MCDM. In general, normalization in MCDM is a transformation process to obtain numerical and comparable input data by using a common scale [START_REF] Vafaei | Data Normalization Techniques in Decision Making: Case Study with TOPSIS Method[END_REF]. After collecting input data, we must do some pre-processing to ensure comparability of criteria, thus making it useful for decision modeling. Furthermore, in MCDM, normalization techniques usually map attributes (criteria) with different measurement units to a common scale in the interval [0-1] [START_REF] Pavlicic | Normalization affects the results of MADM methods[END_REF][START_REF] Etzkorn | Data Normalization and Standardization[END_REF]. Several studies on the effects of normalization techniques on the ranking of alternatives in MCDM problems have shown that certain techniques are more suitable for specific decision methods than others [START_REF] Celen | Comparative Analysis of Normalization Procedures in TOPSIS Method: With an Application to Turkish Deposit Banking Market[END_REF], [START_REF] Chakraborty | A Simulation Based Comparative Study of Normalization Procedures in Multiattribute Decision Making[END_REF][START_REF] Chakraborty | A Simulation Comparison of Normalization Procedures for TOPSIS[END_REF][START_REF] Chakraborty | Rank Similarity based MADM Method Selection[END_REF][START_REF] Milani | The effect of normalization norms in multiple attribute decision making models: a case study in gear material selection[END_REF][START_REF] Wang | Integration of correlations with standard deviations for determining attribute weights in multiple attribute decision making[END_REF][START_REF] Chatterjee | Investigating the Effect of Normalization Norms in Flexible Manufacturing Sytem Selection Using Multi-Criteria Decision-Making Methods[END_REF].

Chakraborty and Yeh [START_REF] Chakraborty | A Simulation Based Comparative Study of Normalization Procedures in Multiattribute Decision Making[END_REF] analyzed four normalization techniques (vector, linear max-min, linear max and linear sum) in the MCDM simple additive weight (SAW) method. They used a ranking consistency index (RCI) and calculated the average deviation for each normalization technique and concluded that the best normalization technique for SAW is the vector normalization. Further, the same authors analyzed the effects of those normalizations for order preference by similarity to ideal solution method (TOPSIS) by calculating ranking consistency and weight sensitivity of each normalization and proved that vector normalization technique is the best for implementing in TOPSIS method [START_REF] Chakraborty | A Simulation Comparison of Normalization Procedures for TOPSIS[END_REF]. The authors [START_REF] Chakraborty | A Simulation Comparison of Normalization Procedures for TOPSIS[END_REF] defined weight sensitivity as a method to analyze sensitivity level of different normalization procedures under different problem settings. They assumed same weights for attributes and then they increased their weights to find the sensitivity of the alternatives (normalization techniques) [START_REF] Chakraborty | A Simulation Comparison of Normalization Procedures for TOPSIS[END_REF].

Also, the result was further validated by Vafaei et al. [START_REF] Vafaei | Data Normalization Techniques in Decision Making: Case Study with TOPSIS Method[END_REF], who used Pearson and Spearman correlation coefficients to also conclude that the best normalization technique for TOPSIS method is the vector normalization.

In this work, we selected five (shown in Table 1) of the most promising normalization techniques [START_REF] Jahan | A state-of-the-art survey on the influence of normalization techniques in ranking: Improving the materials selection process in engineering design[END_REF][START_REF] Celen | Comparative Analysis of Normalization Procedures in TOPSIS Method: With an Application to Turkish Deposit Banking Market[END_REF] and analyzed their effect on the AHP method. In Table 1, each normalization method is divided in two formulas, one for benefit and another for cost criteria, to ensure that the final decision objective (rating) is logically correct, i.e. when it is a benefit criterion for high values it will correspond to high normalized values (maximization -benefit) and when it is a cost criterion high values will correspond to low normalized values (minimization -cost).

Summarizing, the aim of this study is to identify which normalization technique is best suited for the AHP method.

Table 1: Normalization techniques.

Normalization technique Condition of use Formula

Linear: Max (N1) [START_REF] Celen | Comparative Analysis of Normalization Procedures in TOPSIS Method: With an Application to Turkish Deposit Banking Market[END_REF] Benefit criteria 𝑛 𝑖𝑗 = 𝑟 𝑖𝑗 𝑟 𝑚𝑎𝑥

Cost criteria

𝑛 𝑖𝑗 = 1 -𝑟 𝑖𝑗 𝑟 𝑚𝑎𝑥 Linear: Max-Min (N2) [START_REF] Celen | Comparative Analysis of Normalization Procedures in TOPSIS Method: With an Application to Turkish Deposit Banking Market[END_REF] Benefit criteria 𝑛 𝑖𝑗 = 𝑟 𝑖𝑗 -𝑟 𝑚𝑖𝑛 𝑟 𝑚𝑎𝑥 -𝑟 𝑚𝑖𝑛

Cost criteria

𝑛 𝑖𝑗 = 𝑟 𝑚𝑎𝑥 -𝑟 𝑖𝑗 𝑟 𝑚𝑎𝑥 -𝑟 𝑚𝑖𝑛 Linear: sum (N3) [START_REF] Celen | Comparative Analysis of Normalization Procedures in TOPSIS Method: With an Application to Turkish Deposit Banking Market[END_REF] Benefit criteria

𝑛 𝑖𝑗 = 𝑟 𝑖𝑗 ∑ 𝑟 𝑖𝑗 𝑚 𝑖=1 Cost criteria 𝑛 𝑖𝑗 = 1 𝑟 𝑖𝑗 ⁄ ∑ 1 𝑟 𝑖𝑗 ⁄ 𝑚 𝑖=1
Vector normalization (N4) [START_REF] Jahan | A state-of-the-art survey on the influence of normalization techniques in ranking: Improving the materials selection process in engineering design[END_REF] Benefit criteria

𝑛 𝑖𝑗 = 𝑟 𝑖𝑗 √∑ 𝑟 𝑖𝑗 2 𝑚 𝑖=1 Cost criteria 𝑛 𝑖𝑗 = 1 - 𝑟 𝑖𝑗 √∑ 𝑟 𝑖𝑗 2 𝑚 𝑖=1
Logarithmic normalization (N5) [START_REF] Jahan | A state-of-the-art survey on the influence of normalization techniques in ranking: Improving the materials selection process in engineering design[END_REF] Benefit criteria 𝑛 𝑖𝑗 = ln (𝑟 𝑖𝑗 ) ln (∏ 𝑟 𝑖𝑗 ) 𝑚 𝑖=1

Cost criteria

𝑛 𝑖𝑗 =

-ln (𝑟 𝑖𝑗 ) ln (∏ 𝑟 𝑖𝑗 ) 𝑚 𝑖=1

𝑚 -1

Experimental Study with a Smart Car Parking Example

Here, we discuss the suitability of five normalization techniques for AHP with an illustrative example for smart car parking. This illustrative case consists of 3 criteria (C1, C2, C3), which correspond to time to park, distance, and size of the parking space, and 7 alternatives (A1, A2, …, A7), which correspond to candidate location sites for parking. Finding the best place for parking the car is the goal; C1 and C2 are cost criteria, where low values are better, and C3 is a benefit criterion, where high values are desirable. Following the AHP method we defined three pairwise comparison matrices for each criterion (example in Table 2) and then one pairwise comparison matrix between criteria. To these four matrices we applied the five normalization techniques, separately, to determine the ranking of alternatives and compare results. The pairwise decision matrix for criteria "time to park", after steps 1, 2, 3 of AHP, is shown in Table 2.

Table 2:

Pairwise Comparison matrix with respect to the time.

A1 A2 A3 A4 A5 A6 A7 A1 1 1/3 1/2 3 1/3 2 1 A2 3 1 1 4 1 3 1 A3 2 1 1 2 1/2 3 2 A4 1/3 1/4 1/2 1 1/4 1 1/3 A5 3 1 2 4 1 3 1 A6 1/2 1/3 1/3 1 1/3 1 3 A7 1 1 1/2 3 1 1/3 1
We started by testing the sum-based normalization (N3 in Table 1), the usual normalization technique for AHP [START_REF] Saaty | The Analytic Hierarchy Process[END_REF], because it ensures column sum per alternative is equal to one that is defined by Saaty [START_REF] Saaty | The Analytic Hierarchy Process[END_REF]. The other normalization techniques do not include this characteristic and the sum of the normalized values can be bigger than 1; hence, for comparison purposes we opted for re-normalizing the other four using N3. For illustrating the alternatives rating procedure we show the calculation for vector normalization of alternative A1 and the final results for all alternatives are shown in the Table 3 and4: The global weights of alternatives and ranking results for the four tested normalization techniques are shown in Table 5. We discarded the logarithmic normalization technique from our results because we obtained negative and infinite data (due to the characteristics of pairwise matrices), hence it is not usable (appropriate) for the AHP method. As it can be seen in Table 5, there is consensus on which normalization techniques is better for alternatives A2, A3, A4 and A5 (i.e. they all have the same ranking), but for the other alternatives there was no consensus. Since, it is not possible to distinguish which is the best normalization technique just by looking at the results, we used the evaluation approach proposed in [START_REF] Vafaei | Data Normalization Techniques in Decision Making: Case Study with TOPSIS Method[END_REF] to make the assessment. Hence, we calculated Pearson correlation and mean rs values [START_REF] Vafaei | Data Normalization Techniques in Decision Making: Case Study with TOPSIS Method[END_REF] with the global weights of alternatives and Spearman correlation with the ranks of alternatives to assess the suitability of the four tested normalization techniques for the AHP method. Table 6 displays that there exists complete consensus between Pearson and Spearman correlation's results and it is clear that the best normalization technique is N1 (linear: max) because it has the highest mean rs value (P=0.9606 & S=0.9524) and the worst one is N3 (linear sum) with the lowest mean rs value (P=0.9029 & S=0.8413).

𝑃 11= 𝑥 11 √∑ 𝑥 1𝑗 7 𝑗=1 = 1 √(1 2 )+(3 2 )+( 2 2 )+( 1 3 2 )+(3 2 )+( 1 2 2 )+(1) 2 =0

Table 5:

Global weight (G) and Ranking (R) of alternatives for the smart parking example. From the example we can conclude that linear max (N1) is the best normalization technique for the AHP method and linear sum (N3) is the worst one. It is interesting to note that the single normalization used in AHP (linear sum-N3) is the worst one from this comparison study. Although N1 is elected as the most suitable normalization technique it required a re-normalization with N3 because the sum of the normalized values has to be 1. Therefore, we may conclude that a combination of max-normalization (N1) with linear-sum (N3) seems the most appropriate for AHP.

N1 N2 N3 N4 G R G R G R G R A1 0.1972 2

Conclusion

Normalization is the first step of any decision making process to transform data in different units into a common scale and comparable units. In this study we tested five common normalization techniques to assess their suitability for the AHP MCDM method. The tests showed that the logarithmic normalization technique (N5) is not usable in the AHP method because it can result in zero or infinite values in the normalized data, which is not acceptable to use in the method. Further, since AHP requires the columns of the pairwise matrices to sum up 1, the techniques: linear max, linear max-min and vector normalization techniques had to be re-normalized with linear sum (N3) before being compared. To assess the suitability of the normalization techniques for AHP we used Pearson and Spearman correlation and mean rs values; the results showed that the best normalization technique is N1 (linear: max) combined with N3 (linear-sum) to ensure the sum is 1, while the worst one is N3 alone.

In a previous work we did the same assessment study for TOPSIS and in the future we plan to extend it to other well-known MCDM methods, with the aim to support decision makers by recommending the most suitable normalization techniques for usage with each MCDM method.

Table 3 :

 3 Normalization results for vector normalization technique for cost criteria.

		P1	P2	P3	P4	P5	P6	P7	Average
	P1	0.7974 0.8390 0.8091 0.5991 0.8227	0.6524	0.7583	0.7540
	P2	0.3922 0.5169 0.6182 0.4655 0.4681	0.4786	0.7583	0.5283
	P3	0.5948 0.5169 0.6182 0.7327 0.7341	0.4786	0.5165	0.5988
	P4	0.9325 0.8792 0.8091 0.8664 0.8670	0.8262	0.9194	0.8714
	P5	0.3922 0.5169 0.2365 0.4655 0.4681	0.4786	0.7583	0.4737
	P6	0.8987 0.8390 0.8727 0.8664 0.8227	0.8262	0.2748	0.7715

Table 4 :

 4 Re-normalization results for vector normalization technique for cost criteria.

		A1	A2	A3	A4	A5	A6	A7	Average
	A1 0.1659 0.1814 0.1695 0.1304 0.1769 0.1393 0.1598	0.1605
	A2 0.0816 0.1118 0.1295 0.1013 0.1007 0.1022 0.1598	0.1124
	A3 0.1238 0.1118 0.1295 0.1595 0.1578 0.1022 0.1089	0.1276
	A4 0.1941 0.1901 0.1695 0.1886 0.1864 0.1764 0.1938	0.1856
	A5 0.0816 0.1118 0.0495 0.1013 0.1007 0.1022 0.1598	0.1010
	A6 0.1870 0.1814 0.1828 0.1886 0.1769 0.1764 0.0579	0.1644
	A7 0.1659 0.1118 0.1695 0.1304 0.1007 0.2012 0.1598	0.1485
	sum	1	1	1	1	1	1	1	1

Table 6 :

 6 Pearson correlation between global weights & Spearman correlation between ranks of alternatives for each normalization technique.

				0.1925	2	0.1505	4	0.1693	2
	A2	0.0681	6	0.0634	6	0.0762	6	0.1165	6
	A3	0.1143	5	0.1161	5	0.0993	5	0.1297	5
	A4	0.2469	1	0.2658	1	0.2876	1	0.1755	1
	A5	0.0460	7	0.0291	7	0.0749	7	0.1101	7
	A6	0.1765	3	0.1869	3	0.1598	2	0.1450	4
	A7	0.1509	4	0.1462	4	0.1517	3	0.1538	3

Acknowledgements. This work was partially funded by FCT Strategic Program UID/EEA/00066/203 of Computational Intelligence Group of CTS/UNINOVA.


Slavisa Tomic
email: s.tomic@campus.fc.unl.pt

Marko Beko
email: mbeko@uninova.pt

Rui Dinis
email: rdinis@fct.unl.pt

Milan Tuba
email: tuba@ieee.org

A WLS Estimator for Target Localization in a Cooperative Wireless Sensor Network

Keywords: Angle-of-arrival (AoA), generalized trust region sub-problem (GTRS), received signal strength (RSS), wireless sensor network (WSN)

This paper addresses target localization problem in a cooperative 3-D wireless sensor network (WSN). We employ non-traditional methodology which merges distance and angle measurements, respectively withdrawn from the received signal strength (RSS) and angle-of-arrival (AoA) information. Based on RSS measurement model and effortless geometry, a novel non-convex estimator according to the weighted least squares (WLS) criterion is obtained, which closely approximates the maximum likelihood (ML) estimator for small noise. It is shown that the devised estimator is appropriate for distributed implementation. Following the squared range (SR) approach, we propose a suboptimal SR-WLS estimator according to the generalized trust region sub-problem (GTRS) framework, to estimate the locations of all targets in the WSN. According to our simulations, the new estimator has excellent performance in a great variety of considered settings, in which the effectiveness of fusing two radio measurements is confirmed.

Introduction

Accurate information about sensor's location is a key component in many practical applications. Wireless localization algorithms usually rely upon range measurements extracted from angle-of-arrival (AoA), received signal strength (RSS), time-of-arrival (ToA) information, or a combination of them [START_REF] Patwari | Location Estimation in Sensor Networks[END_REF].

In [START_REF] Yu | 3-D Localization Error Analysis in Wireless Networks[END_REF], a hybrid methodology which fuses distance and angle measures has been considered. Both estimators proposed in [START_REF] Yu | 3-D Localization Error Analysis in Wireless Networks[END_REF] deal with the non-cooperative target localization problem in a 3-D space: linear least squares (LS) and optimization based. The former one is a fairly simple estimator, and the later one uses Davidson-Fletcher-Powell algorithm [START_REF] Fletcher | Practical Methods of Optimization[END_REF]. The authors in [START_REF] Wang | Hybrid RSS/AOA Emitter Location Estimation Based on Least Squares and Maximum Likelihood Criteria[END_REF] derived an LS and a maximum likelihood (ML) estimator for a hybrid scheme that fuses RSS difference (RSSD) and AoA measurements. The authors in [START_REF] Wang | Hybrid RSS/AOA Emitter Location Estimation Based on Least Squares and Maximum Likelihood Criteria[END_REF] employed non-linear constrained optimization to estimate the unknown location from multiple RSS and AoA measurements. In [START_REF] Gazzah | Selective Hybrid RSS/AOA Weighting Algorithm for NLOS Intra Cell Localization[END_REF], a selective weighted LS (WLS) estimator for mixed RSS/AoA localization problematic was presented. The target position was established by employing loaded distances from the two closest anchor measures, together with the serving base station AoA measurement. Another WLS estimator for a non-cooperative localization problem for the case where the transmitted power is not known was proposed in [START_REF] Chan | Hybrid Localization of an Emitter by Combining Angle-of-Arrival and Received Signal Strength Measurements[END_REF]. Nonetheless, similar as the method proposed in [START_REF] Gazzah | Selective Hybrid RSS/AOA Weighting Algorithm for NLOS Intra Cell Localization[END_REF], the WLS method has been derived without requiring information about the statistical properties of RSS and AoA measurements, which might lead to significant performance degradation in practice. Also, the authors in [START_REF] Chan | Hybrid Localization of an Emitter by Combining Angle-of-Arrival and Received Signal Strength Measurements[END_REF] investigated a small-scale wireless sensor network (WSN) only, with extremely low noise power.

All of the above approaches examine non-collaborative localization problem only, where the location of a single target, which communicates with anchors exclusively, is established at a time. Contrarily to the above mentioned approaches, here, the target localization problem in an extensive WSN is considered, where the number of anchors is insufficient and the communication range of all sensors is limited (e.g., to prolong sensor's battery life). Hence, only some targets can directly communicate with anchors and sensor cooperation is required in order to acquire adequate amount of information to carry out the localization. Through employing the RSS spreading model and straightforward geometry, we first develop a new local non-convex estimator based on the WLS criterion that closely approximates the local ML one in the case where noise is low. Next, by following the squared range (SR) approach, we propose a suboptimal SR-WLS estimator based on the generalized trust region sub-problem (GTRS) framework, which is possible to solve exactly with a bisection method [START_REF] Beck | Exact and Approximate Solutions of Source Localization Problems[END_REF]. To the best of the authors' knowledge, distributed localization algorithms for hybrid RSS/AoA systems in cooperative WSNs are yet to be published.

Relationship to Cyber-Physical Systems

Cyber-physical systems (CPSs) have recently attracted much attention from both the educational and engineering public. These schemes are displaying an enormous potential in interacting with the physical world and creating its management more efficient. Integrating the computation and communication competencies into the modules of the physical surrounding [START_REF] Koubaa | Taxonomy of Fundamental Concepts of Localization in Cyber-Physical and Sensor Networks[END_REF] enables this. CPSs represent nowadays the new peer group of networks and embedded structures.

The placement of CPS raises numerous difficulties, where the most significant topic that has activated massive quantity of study is localization [START_REF] Koubaa | Taxonomy of Fundamental Concepts of Localization in Cyber-Physical and Sensor Networks[END_REF]. Basically, localization goals at estimating the location of CPS modules. In CPSs where data are closely related with the surroundings and the location at which they are produced, localization is a fundamental task.

Collaboration of a large number of scattered sensors in a WSN can be seen as a CPS. In such systems, localization is of a crucial importance, since a system may be configured to react locally to the variations within sensor records; hence, accurate determination of the location where the deviations arise is the key. Furthermore, services based on location-awareness represent a key component of countless wireless structures nowadays. By exploring the synergies between computational and physical components we can form smart environments which offer improved safety and efficiency in everyday life, e.g. smart parking, assistance for elderly or people with disabilities, monitoring of storage conditions and goods, etc.

Problem Statement

We examine a large-scale WSN with 𝑁 anchors and 𝑀 targets, that may also be viewed In favor of ease of expression, we describe a matrix 𝑿 = [𝒙 1 , ⋯ , 𝒙 𝑀 ] (𝑿 ∈ ℝ 3×𝑀 ) as the matrix of all unknown target locations. These unknown locations are estimated with the means of a hybrid methodology which merges distance and angle measures.

Throughout this work, we imply that the distance measurements are extracted from the RSS measurements, because ranging based on RSS does not oblige additional hardware [START_REF] Patwari | Location Estimation in Sensor Networks[END_REF]. The path loss between two sensors is given by:

𝐿 𝑖𝑗 = 𝐿 0 + 10𝛾 log 10 ‖𝒙 𝑖 -𝒂 ̂𝑗‖ 𝑑 0 + 𝑛 𝑖𝑗 , ∀(𝑖, 𝑗) ∈ ℇ 𝒜 ∪ ℇ 𝒯 (1)
(see [START_REF] Rappaport | Wireless Communications: Principles and Practice[END_REF], [START_REF] Sichitiu | Localization of Wireless Sensor Networks with a Mobile Beacon[END_REF]) where 𝐿 0 represents the path loss value at a reference distance 𝑑 0 (‖𝒙 𝑖 -𝒂 ̂𝑗‖ ≥ 𝑑 0 ), 𝛾 is the path loss exponent (PLE), 𝒂 ̂𝑗 is the 𝑗-th neighbor of the 𝑖-th target (𝒂 𝑗 if 𝑗 ∈ 𝒜, 𝒙 𝑗 if 𝑗 ∈ 𝒯), and 𝑛 𝑖𝑗 is the log-normal shadowing considered as a Gaussian random variable, i.e., 𝑛 𝑖𝑗 ~𝒩(0, 𝜎 𝑛 𝑖𝑗 2 ). The target/target path loss measures are assumed to be symmetric 1 . Fig. 1 gives a representation of a target and anchor locations within 3-D. Looking at Fig. 1, 𝒙 𝑖 = [𝑥 𝑖1 , 𝑥 𝑖2 , 𝑥 𝑖3 ] 𝑇 and 𝒂 𝑗 = [𝑎 𝑗1 , 𝑎 𝑗2 , 𝑎 𝑗3 ]

𝑇 represent the coordinates of the 𝑖th target (not known) and the coordinates of the 𝑗-th anchor (known) respectively, whereas 𝑑 𝑖𝑗 , 𝜙 𝑖𝑗 and 𝛼 𝑖𝑗 respectively stand for the range, azimuth and elevation angle amongst the 𝑖-th target and the 𝑗-th anchor. According to the measurement model [START_REF] Patwari | Location Estimation in Sensor Networks[END_REF] the ML distance between two sensors is obtained as follows [START_REF] Patwari | Location Estimation in Sensor Networks[END_REF]:

𝑑 ̂𝑖𝑗 = 𝑑 0 10 𝐿 𝑖𝑗 -𝐿 0 10𝛾 , ∀(𝑖, 𝑗) ∈ ℇ 𝒜 ∪ ℇ 𝒯 . (2)
By applying simple geometry, azimuth and elevation angle measurements can be modeled respectively as [START_REF] Yu | 3-D Localization Error Analysis in Wireless Networks[END_REF]:

𝜙 𝑖𝑗 = tan -1 ( 𝑥 𝑖2 -𝑎 𝑗2 𝑥 𝑖1 -𝑎 𝑗1 ) + 𝑚 𝑖𝑗 , ∀(𝑖, 𝑗) ∈ ℇ 𝒜 , (3)
and

𝛼 𝑖𝑗 = cos -1 ( 𝑥 𝑖3 -𝑎 𝑗3 ‖𝑥 𝑖 -𝑎 𝑗 ‖ ) + 𝑣 𝑖𝑗 , ∀(𝑖, 𝑗) ∈ ℇ 𝒜 , (4)
According to (1), ( 3) and ( 4), we can obtain the ML estimator as [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]:

𝑿 ̂= arg min 𝑿 ∑ 1 𝜎 𝑖 2 [𝜃 𝑖 -𝑓 𝑖 (𝑿)] 2 3|ℇ 𝒜 |+|ℇ 𝒯 | 𝑖=1 , (5)
𝒇(𝑿) = [ ⋮ 𝐿 0 + 10𝛾 log 10 ‖𝒙 𝑖 -𝒂 ̂𝑗‖ 𝑑 0 ⋮ tan -1 ( 𝑥 𝑖2 -𝑎 𝑗2 𝑥 𝑖1 -𝑎 𝑗1 ) ⋮ cos -1 ( 𝑥 𝑖3 -𝑎 𝑗3 ‖𝒙 𝑖 -𝒂 𝑗 ‖ ) ⋮ ] , 𝝈 = [ ⋮ 𝜎 𝑛 𝑖𝑗 ⋮ 𝜎 𝑚 𝑖𝑗 ⋮ 𝜎 𝑣 𝑖𝑗 ⋮ ]
.

The non-convex LS estimator in [START_REF] Gazzah | Selective Hybrid RSS/AOA Weighting Algorithm for NLOS Intra Cell Localization[END_REF] does not have a closed-form solution. Nonetheless, we will show that it can be transformed into GTRS framework and solved efficiently.

The Proposed SR-WLS Estimator

Assuming that the initial location estimations of the targets, 𝑿 ̂(0) , are given, the ML in (5) can be solved locally per target, resorting merely to the data collected from its neighbors and employing an iterative methodology. Consequently, target 𝑖 updates its location estimation in each iteration, 𝑡, by solving the below local ML estimator:

𝒙 ̂𝑖 (𝑡+1) = arg min 𝒙 𝒊 ∑ 1 𝜎 𝑗 2 [𝜃 𝑗 -𝑓 𝑗 (𝒙 𝑖 )] 2 , ∀𝑖 ∈ 𝒯 3|ℇ 𝒜 𝑖 |+|ℇ 𝒯 𝑖 | 𝑗=1 , (6)
where ℇ 𝐴 𝑖 = {𝑗: (𝑖, 𝑗) ∈ ℇ 𝐴 } and ℇ 𝑇 𝑖 = {𝑘: (𝑖, 𝑘) ∈ ℇ 𝑇 } represent the set of all anchor and all target neighbors of the target 𝑖 respectively. Given 𝑿 ̂(0) , for the case where noise power is low, out of (1) we can write:

𝜆 𝑖𝑗 ‖𝒙 𝑖 -𝒂 ̂𝑗‖ 2 ≈ 𝑑 0 2 , ∀𝑖 ∈ 𝒯, ∀𝑗 ∈ ℇ 𝒜 𝑖 ∪ ℇ 𝒯 𝑖 , (7)
where 𝜆 𝑖𝑗 = 10 𝐿 0 -𝐿 𝑖𝑗 5𝛾 . Similarly, from ( 3) and ( 4) we respectively get:

𝒄 𝑖𝑗 𝑇 (𝒙 𝑖 -𝒂 𝑗 ) ≈ 0, ∀𝑖 ∈ 𝒯, ∀𝑗 ∈ ℇ 𝒜 𝑖 , (8)
and

𝒌 𝑇 (𝒙 𝑖 -𝒂 𝑗 ) ≈ ‖𝒙 𝑖 -𝒂 𝑗 ‖ cos(𝛼 𝑖𝑗 ) , ∀𝑖 ∈ 𝒯, ∀𝑗 ∈ ℇ 𝒜 𝑖 . (9)
where 𝒄 𝑖𝑗 = [-sin(𝜙 𝑖𝑗 ), cos(𝜙 𝑖𝑗 ), 0] 𝑇 and 𝒌 = [0,0,1] 𝑇 . To grant more significance to nearby connections, bring into play weights, 𝒘 = [ √ 𝑤 𝑖𝑗 ], where 𝑤 𝑖𝑗 = 1 -𝑑 ̂𝑖𝑗 ∑ 𝑑 ̂𝑖𝑗 ∀𝑖∈𝒯,∀𝑗∈ℇ 𝒜 𝑖 ∪ℇ 𝒯 𝑖 ⁄ , and substitute ‖𝒙 𝑖 -𝒂 ̂𝑗‖ in (9) with 𝑑 ̂𝑖𝑗 described in [START_REF] Yu | 3-D Localization Error Analysis in Wireless Networks[END_REF].

According to the WLS criterion and ( 7), ( 8) and ( 9) every target updates its location through minimization of the below estimator:

𝒙 ̂𝑖 (𝑡+1) =
) 10
The above WLS problem is not convex and does not have a closed-form solution. Still, (10) may be expressed as a quadratic programming problem of which the global solution is found effortlessly [START_REF] Beck | Exact and Approximate Solutions of Source Localization Problems[END_REF]. Applying the replacement 𝒚 𝑖 = [𝒙 𝒊 𝑻 , ‖𝒙 𝑖 ‖ 2 ] 𝑇 , ∀𝑖 ∈ 𝒯, equation ( 10) is modified as:

𝒚 𝑖 (𝑡+1) = arg min 𝒚 𝑖 ‖𝑾(𝑨𝒚 𝑖 -𝒃)‖ 2 subject to 𝒚 𝑖 𝑇 𝐃𝐲 i + 2𝑙 𝑇 𝒚 𝑖 = 0, (11)
where 𝑾 = 𝑰 𝟑 ⨂diag(𝒘), and ⨂ denotes the Kronecker product, 𝑰 𝑞 represents the 𝑞 dimensional identity matrix, and

𝑨 = [ ⋮ ⋮ -2𝜆 𝑖𝑗 𝒂 ̂𝑗 𝑇 𝜆 𝑖𝑗 ⋮ 𝒄 𝑖𝑗 𝑇 ⋮ 𝒌 𝑇 ⋮ ⋮ 0 ⋮ 0 ⋮ ] , 𝒃 = [ ⋮ 𝑑 0 2 -𝜆 𝑖𝑗 ‖𝒂 ̂𝑗‖ 2 ⋮ 𝒄 𝑖𝑗 𝑇 𝒂 𝑗 ⋮ 𝒌 𝑇 𝒂 𝑗 + 𝑑 ̂𝑖𝑗 cos(𝛼 𝑖𝑗 ) ⋮ ] , 𝑫 = [ 𝑰 3 𝟎 3×1 𝟎 1×3 0 ] , 𝒍 = [ 𝟎 3×1 1 2 ],
i.e., 𝑨 ∈ ℝ 3|ℇ 𝒜 𝑖 |+|ℇ 𝒯 𝑖 |×4 , 𝒃 ∈ ℝ 3|ℇ 𝒜 𝑖 |+|ℇ 𝒯 𝑖 |×1 , 𝑾 ∈ ℝ 3|ℇ 𝒜 𝑖 |+|ℇ 𝒯 𝑖 |×3|ℇ 𝒜 𝑖 |+|ℇ 𝒯 𝑖 | , and 𝟎 𝑝×𝑞 is the 𝑝 × 𝑞-dimensional matrix of all zeros.

Both the objective and the restraint functions in [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF] are quadratic. This type of problem is identified as GTRS [START_REF] Beck | Exact and Approximate Solutions of Source Localization Problems[END_REF]; it may be solved exactly by a bisection method [START_REF] Beck | Exact and Approximate Solutions of Source Localization Problems[END_REF].

Assuming that 𝒞 is the set of colors of the nodes (in order to coordinate the network [START_REF] Tomic | Distributed RSS-Based Localization in Wireless Sensor Networks Based on Second-Order Cone Programming[END_REF]), Algorithm 1 encapsulates the suggested distributed SR-WLS algorithm. Lines 5-7 are carried out concurrently by all targets 𝑖 ∈ 𝒞 𝑐 , which might reduce the realization time of the approach. Information interchange occurs exclusively at Line 7, when targets send their location updates 𝒙 ̂𝒊 (𝑡+1) to their neighbors. Because 𝒙 ̂𝒊 (𝑡+1) ∈ ℝ 3 , one concludes that the new approach obliges a transmission of a maximum 3 × 𝑇 max × 𝑀 real values. The worst case computational complexity of the proposed method is linear, i.e., 𝑇 max × 𝑀 × 𝒪 (𝑁 max × max 𝑖 {3|ℇ 𝒜 𝑖 | + |ℇ 𝒯 𝑖 |}), where 𝑁 max denotes the maximum allowed number of iterations of the bisection method. In the further text, we will refer to Algorithm 1 as "SR-WLS". for all 𝑖 ∈ 𝒞 𝑐 (simultaneously) do 5.

Gather 𝒂 ̂𝑗, ∀𝑗 ∈ ℇ 𝒜 𝑖 ∪ ℇ 𝒯 𝑖

6.

𝒙 ̂𝑖 (𝑡+1) ← solve (

Broadcast 𝒙 ̂𝑖 (𝑡+1) to 𝒂 ̂𝑗, ∀𝑗 ∈ ℇ 𝒜 𝑖 ∪ ℇ 𝒯 𝑖 8.

end for 9. end for 10. 𝑡 ← 𝑡 + 1 11. until 𝑡 < 𝑇 max

Performance Evaluation

A set of simulation results is presented here with the intention of assessment of the performance of the new approach in the view of the estimation precision and convergence. All of the considered approaches were solved via MATLAB. In order to demonstrate the benefit of fusing two radio measurements versus traditional localization systems, we include also the performance results of the proposed method when only RSS measurements are employed, called here "SR-WLSRSS".

In order to produce the radio measures, models (1), ( 3) and ( 4) are employed. A random arrangement of sensors within a box with the edge span 𝐵 = 20 m in each Monte Carlo (𝑀 𝑐 ) trial was considered. If not declared otherwise, the PLE was set to 𝛾 = 3, the reference distance to 𝑑 0 = 1 m, the reference path loss to 𝐿 0 = 40 dB, and the maximum number of steps in the bisection procedure was 𝑁 𝑚𝑎𝑥 = 30. In practice however, it is almost impossible to perfectly estimate the value of the PLE. Thus, in order to account for a realistic measurement model mismatch and test the robustness of the considered approaches to flawed knowledge of the PLE, the true PLE for each connection was chosen from a uniform distribution on an interval 𝛾 𝑖𝑗 ∈ [2.7, 3.3], ∀(𝑖, 𝑗) ∈ ℇ 𝒜 ∪ ℇ 𝒯 . Finally, the zero estimate of the targets' locations, 𝑿 ̂(0) , was assumed to be at the intersection of the big diagonals of the cube area. The main performance criterion used is the normalized root mean square error (NRMSE), defined as

NRMSE = √∑ ∑ ‖𝒙 𝑖𝑗 -𝒙 ̂𝑖𝑗 ‖ 𝑀 𝑀 𝑐 𝑀 𝑗=1 𝑀 𝑐 𝑖=1
, where 𝒙 ̂𝑖𝑗 represents the estimation of the real location of the j-th target, 𝒙 𝑖𝑗 , in the i-th 𝑀 𝑐 trial. From these figures, we can notice that the output from considered approaches improves as 𝑡 grows, as anticipated. Also, it can be seen from Fig. 2 that the performance of all algorithms improves considerably as more anchors are included into the network. This behavior is anticipated, because when 𝑁 is increased more reliable information is introduced in the network. Furthermore, Fig. 3 reveals that both methods need a somewhat elevated number of repetitions to converge for augmented 𝑀. Nonetheless, the estimation precision of the considered algorithms improves when more targets are added in the network. Moreover, from Figs. 2 and3 one can perceive that all major changes in the performance for the considered algorithms occur in the first few iterations (𝑡 ≤ 15), and that the performance gain is insignificant afterwards. This result is very important because it shows that our approach necessitates a low number of signal broadcast, which might augment the exploitation productivity of the spectrum, a valuable resource for wireless communications. It also shows that our algorithm is energy efficient; the communication stage is considerably more demanding (in terms of energy) than the data processing one [START_REF] Patwari | Location Estimation in Sensor Networks[END_REF], [START_REF] Tomic | Distributed RSS-Based Localization in Wireless Sensor Networks Based on Second-Order Cone Programming[END_REF].

Conclusion

Here, the hybrid RSS/AoA target localization problem in a collaborative 3-D WSN was addressed. By fusing information from RSS and AoA, we devised a new problem formulation based on a nonconvex WLS. We showed that the derived novel estimator is suitable for distributed implementation, and we presented a novel non-centralized procedure founded on the GTRS framework. This new algorithm is light in terms of computational complexity, and it can provide accurate localization in a variety of scenarios, having exhibited excellent performance both in view of the estimation precision and convergence. Moreover, the simulation results confirmed the effectiveness of combining RSS and AoA radio measurements in comparison with conventional RSS localization system, displaying a notable enhancement of the estimation precision.

  as a connected graph, 𝒢(𝒱, ℇ), with |𝒱| = 𝑀 + 𝑁 vertices and |ℇ| edges (connections), where | • | represents the cardinality of a set. The set of targets and the set of anchors are respectively denoted as 𝒯 (|𝒯| = 𝑀) and 𝒜 (|𝒜| = 𝑁), and their locations are denoted by 𝒙 1 , 𝒙 2 , ⋯ , 𝒙 𝑀 and 𝒂 1 , 𝒂 2 , ⋯ , 𝒂 𝑁 , (𝒙 𝑖 , 𝒂 𝑗 ∈ ℝ 3 , ∀𝑖 ∈ 𝒯 and ∀𝑗 ∈ 𝒜), respectively. Due to limited communication range, 𝑅, two sensors, 𝑖 and 𝑗, can interchange data if and only if they are inside the communication range of each other. The sets of all target/anchor and target/target edges are defined as ℇ 𝒜 = {(𝑖, 𝑗): ‖𝒙 𝑖 -𝒂 𝑗 ‖ ≤ 𝑅, ∀𝑖 ∈ 𝒯, ∀𝑗 ∈ 𝒜} and ℇ 𝒯 = {(𝑖, 𝑘): ‖𝒙 𝑖 -𝒙 𝑘 ‖ ≤ 𝑅, ∀𝑖, 𝑘 ∈ 𝒯, 𝑖 ≠ 𝑘}, respectively.

Fig. 1 .

 1 Fig. 1. Representation of a target and anchor locations within 3-D.

  where 𝜽 = [𝑳 𝑇 , 𝝓 𝑇 , 𝜶 𝑇 ] 𝑇 (𝜽 ∈ ℝ 3|ℇ 𝒜 |+|ℇ 𝒯 | ), with 𝑳 = [𝐿 𝑖𝑗 ] 𝑇 , ∀(𝑖, 𝑗) ∈ ℇ 𝒜 ∪ ℇ 𝒯 ,𝝓 = [𝜙 𝑖𝑗 ] 𝑇 , ∀(𝑖, 𝑗) ∈ ℇ 𝒜 , 𝜶 = [𝛼 𝑖𝑗 ] 𝑇 , ∀(𝑖, 𝑗) ∈ ℇ 𝒜 , and

Algorithm 1 -

 1 The proposed distributed SR-WLS algorithm Require: 𝑿 ̂(0) , 𝑇 max , 𝒞, 𝒂 𝑗 , ∀𝑗 ∈ 𝒜 1. Set: 𝑡 ← 0 2. replicate 3. for 𝑐 = 1, … , 𝒞 do 4.

Fig. 2 .

 2 Fig. 2. NRMSE set against 𝑡 assessment, when 𝑀 = 𝟓𝟎, 𝑹 = 𝟔. 𝟓 𝐦, 𝝈 𝒏 𝒊𝒋 = 𝟑 𝐝𝐁, 𝝈 𝒎 𝒊𝒋 = 𝟔 𝐝𝐞𝐠, 𝝈 𝒗 𝒊𝒋 = 𝟔 𝐝𝐞𝐠, 𝜸 𝒊𝒋 ∈ [𝟐. 𝟕, 𝟑. 𝟑], 𝜸 = 𝟑, 𝐁 = 𝟐𝟎 𝐦, 𝑳 𝟎 = 𝟒𝟎 𝐝𝐁, 𝒅 𝟎 = 𝟏 𝐦, 𝑴 𝒄 = 𝟓𝟎𝟎.

Figs. 2

 2 Figs. 2 and 3 illustrate the NRMSE versus 𝑡 performance when 𝑅 = 6.5 m, 𝜎 𝑛 𝑖𝑗 = 3 dB, 𝜎 𝑚 𝑖𝑗 = 6 deg, and 𝜎 𝑣 𝑖𝑗 = 6 deg, for 𝑀 = 50, 𝑁 = 20 and 𝑁 = 30, and 𝑁 = 20, 𝑀 = 50 and 𝑀 = 60, respectively.From these figures, we can notice that the output from considered approaches improves as 𝑡 grows, as anticipated. Also, it can be seen from Fig.2that the performance of all algorithms improves considerably as more

Fig. 3 .

 3 Fig. 3. NRMSE versus 𝑡 comparison, when 𝑁 = 20, 𝑅 = 6.5 m, 𝜎 𝑛 𝑖𝑗 = 3 dB, 𝜎 𝑚 𝑖𝑗 = 6 deg, 𝜎 𝑣 𝑖𝑗 = 6 deg, 𝛾 𝑖𝑗 ∈ [2.7, 3.3], 𝛾 = 3, B = 20 m, 𝐿 0 = 40 dB, 𝑑 0 = 1 m, 𝑀 𝑐 = 500.

Fig. 4

 4 Fig.4illustrates the NRMSE versus 𝜎 𝑛 𝑖𝑗 (dB) assessment, when 𝑁 = 20, 𝑀 = 50, 𝑅 = 6.5 m, and 𝑇 max = 30. In this figure, one can perceive the performance degradation of the proposed algorithm as the quality of the RSS measurement reduces, as projected. Nonetheless, we can see from the figure that the deterioration in the performance is lower than 15 % for the proposed algorithm, which is relatively low for the considered noise range. Finally, Fig.4confirms the effectiveness of measurement fusion, showing a gain of almost 1 m when the hybrid system is employed in comparison with the traditional RSS system.

Fig. 4 .

 4 Fig. 4. NRMSE versus 𝜎 𝑛 𝑖𝑗 (dB) assessment, for 𝑁 = 20, 𝑀 = 50, 𝑅 = 6.5 m, 𝜎 𝑚 𝑖𝑗 = 1 deg, 𝜎 𝑣 𝑖𝑗 = 1 deg, 𝛾 𝑖𝑗 ∈ [2.7, 3.3], 𝛾 = 3, 𝑇 max = 30, 𝐵 = 20 m, 𝐿 0 = 40 dB, 𝑑 0 = 1 m, 𝑀 𝑐 = 500.

Acknowledgments. This work was partially supported by Fundação para a Ciência e a Tecnologia under Projects PEst-OE/EEI/UI0066/2014, and PEst-OE/EEI/LA0008/ 2013 (IT pluriannual founding and HETNET), PEst-OE/EEI/UI0066/2011 (UNINOVA pluriannual founding), ADIN PTDC/EEI-TEL/2990/2012, COPWIN PTDC/EEI-TEL/1417/2012 and PTDC/EEITEL/6308/2014-HAMLeT, as well as the grants SFRH/BPD/108232/2015, SFRH/BD/91126/2012 and Ciência 2008 Post-Doctoral Research grant.


Ondrej Kachman
email: ondrej.kachman@savba.sk

Marcel Balaz
email: marcel.balaz@savba.sk

Effective Over-the-Air Reprogramming for Low-Power Devices in Cyber-Physical Systems

Keywords: reprogramming, low-power devices, over-the-air programming, firmware similarity, differencing algorithm, update agent, power consumption

Cyber-physical systems often include sensor devices in their structure. These devices may require firmware updates once deployed and these updates must be energy efficient for battery powered, physically inaccessible sensors. The problem of energy saving reprogramming can be split to four tasks -making old and new firmware versions more similar, generating small delta files using differencing algorithms, propagating delta files and applying updates at the end devices. This paper describes existing approaches dealing with this problem, analyzes their power consumption and introduces new optimizations for differencing algorithms. A new approach is presented that requires no external flash memory, device reboot or complex update agent at the sensor device.

Introduction

As the self-programming flash memories and wireless communication entered the world of embedded systems, over-the-air reprogramming became possible. It was especially examined for wireless sensor networks (WSNs), where hundreds of small sensor devices may be deployed in a large area. Firmware preloaded on devices in WSN is usually developed under test conditions and may not work properly once the devices are deployed. Over-the-air reprogramming made it easier to update firmware on the devices, but transferring the whole firmware image to many devices is energy inefficient as sending 1 bit wirelessly can consume the same energy as 1000 instructions [START_REF] Levis | Maté: A Tiny Virtual Machine for Sensor Networks[END_REF]. This is critical for wireless low-power devices powered by batteries that require firmware update.

The problem of effective over-the-air reprogramming motivates the research question of this paper. How do over-the-air updates influence energy consumption of the updated devices and what can be done to improve the energy efficiency of these updates. The task is also to identify partial problems in this area and their solutions. We present our contribution to differencing algorithms and show how it can help to perform firmware updates effectively.

Relationship to Cyber-Physical Systems

Cyber-physical systems (CPS) are a promising direction in automated interaction between the physical and the virtual world. Based on information collected from physical world, CPS can properly respond with programmed actions. Interaction between the two worlds is provided by intelligence of CPS, which may be created by 3D displays, sensors, actuators, cameras and other devices. This links CPS with previously mentioned WSNs and low-power devices. Sensors and actuators often create an important part of CPS and face the same challenges [START_REF] Wu | From wireless sensor networks towards cyberphysical systems[END_REF], for example limited battery life, network traffic processing with limited computing capabilities or the adaptation to environment. Over-the-air reprogramming for small devices within CPS should be effective in terms of energy consumption, speed and security.

Energy can be wasted by excessive amount of update data disseminated through the network to end devices. Many protocols have been developed for WSNs [START_REF] Stathopoulos | A Remote Code Update Mechanism for Wireless Sensor Networks[END_REF] [4] [START_REF] Aschenbruck | Selective and Secure Over-The-Air Programming for Wireless Sensor Networks[END_REF]. However, WSNs usually consist of the same or similar devices and serve to collect data. CPS may consist of many devices with different resources and tasks that may be networked by different technologies, for example WLAN, Bluetooth or GSM [START_REF] Shi | A Survey of Cyber-Physical Systems[END_REF]. When embedding over-the-air programming capabilities into the low-power device's firmware, developers must consider the right networking technology and protocols in order to reduce the amount of energy consumed by the firmware updates. Another way of energy wasting can be caused by the firmware code. Modern microcontrollers usually have firmware code stored in a program flash memory. Firmware updates are applied from delta files that encode the differences between older and newer versions. Some approaches proposed alteration to compilers [START_REF] Zhang | Towards Energy-Efficient Code Dissemination in Wireless Sensor Networks[END_REF] [8] and linkers [START_REF] Koshy | Remote Incremental Linking for Energy-Efficient Roprogramming for Sensor Networks[END_REF] in order to make the versions more similar. Compiler alterations led to more instructions, thus worse execution times and some linker alterations resulted in fragmented program memory. Calls to functions in different parts of flash memory can lead to worse energy consumption as activation of different flash regions requires additional energy [START_REF] Pallister | A high-level model of embedded flash energy consumption[END_REF]. When considering many different devices and platforms in CPS, compiler and linker alterations would be required for every different platform. Developers can instead work with the product of compilers -object files, which usually have the standard executable and linkable (ELF) format, and work with sections and relocation entries [START_REF] Dong | R2: Incremental Reprogramming Using Relocatable Code in Networked Embedded Systems[END_REF].

Update speed is important for CPS too. Ineffective update dissemination can slow down network traffic and long downtime of important devices may cause CPS to malfunction. The problem of network traffic can be addressed by effective protocols, Quality of Service mechanisms (QoS) and reduced amount of update data using delta files generated by the differencing algorithms that compare the old and new firmware versions. The delta files are later processed by the end devices and their decoding followed by application of the update may greatly influence the update speed, so their structure and encoding must be taken into consideration by the firmware developers.

The security problem of over-the-air updates concerns CPS as well as any other systems of collaborating computational devices. A device with rogue firmware could infect the network and compromise the system. Network protocols disseminating the update data must implement some security mechanisms and the delta files should include an integrity check to prevent them from tampering.

State of the Art and Related Work

The area of effective over-the-air reprogramming can be split into different problems. Based on our analysis of the existing research in this area, we split it to four partial problems: (1) Firmware version similarity, (2) Delta generating differencing algorithms, (3) Update dissemination, and (4) Update application.

Overview of over-the-air update procedures is shown in Fig. 1. The following subsections analyze state of the art for each of the listed problems. One more subsection on flash memory energy consumption modeling is included as a base for our energy analysis of over-the-air updates.

Fig. 1 Over-the-air updates overview

Firmware Version Similarity

Researchers agree, that firmware updates are more energy effective when data that must be sent through the network are minimal, so the devices save energy on the network communication. First step towards the smaller update data is through improvement in the similarity of two consecutive firmware versions. Considering multiplatform CPS, ELF files are the best way to improve firmware similarity. Authors of [START_REF] Shafi | No-reboot and Zero-Flash Over-the-air Programming for Wireless Sensor Networks[END_REF] identified 4 ELF file properties that caused larger delta file generation -branches, global variables, indirect addressing and relative jumps. Preserving addresses allocated for functions, variables and constants between the firmware versions results in a smaller difference between the linked files. Some of these problems can be addressed by changes to object file relocation entries [START_REF] Dong | R2: Incremental Reprogramming Using Relocatable Code in Networked Embedded Systems[END_REF] [START_REF] Dong | R3: Optimizing relocatable code for efficient reprogramming in networked embedded systems[END_REF]. Tools altering the ELF files to improve the similarity of firmware images can help the differencing algorithms to perform significantly better.

Delta Generating Differencing Algorithms

The differencing algorithms that generate delta files (also called deltas or patches) for low-power devices must generate delta as small as possible and encode it in such a way, that the end device can decode it and apply the update easily. Basic approaches used block level comparison [START_REF] Jeong | Incremental Network Programming for Wireless Sensors[END_REF], which was fast but resulted in unnecessarily large deltas. These deltas transferred whole blocks of fixed size that have changed. Soon, byte level differencing algorithms and delta files were introduced [START_REF] Hu | Reprogramming with Minimal Transferred Data on Wireless Sensor Network[END_REF] [16] [START_REF] Mo | An efficient differencing algorithm based on suffix array for reprogramming wireless sensor networks[END_REF]. These algorithms detected common sequences in the firmware images and new data. The difference information were encoded into delta files using usually two basic operations -COPY and ADD, each with different attributes. Basic format encoded into file is following:

 <COPY><old address><new address><number of bytes>  <ADD><new address><number of bytes><first byte> … <last byte> Every attribute may cost different amount of bytes. Cost for every attribute of an operation adds up to operation's total cost within a delta file, usually in bytes or words.

Update Dissemination

The delta dissemination stage greatly depends on the network protocol. One of the most notable protocols for WSNs updates is Deluge [START_REF] Hui | The Dynamic Behavior of a Data Dissemination Protocol for Network Programming at Scale[END_REF], included in TinyOS, running on many wireless sensor nodes. As mentioned in section 2, CPS consist of different platforms and the network topology may be dynamic. Standard protocols can be used on the higher layers of the CPS. Supported by the QoS mechanisms, update data should reach end devices quickly. If these low power devices cannot support standard the communication protocols (like the most common TCP/IP protocol) due to their computational restraints, they should be provided with border routers that would communicate with them using some lightweight protocol [START_REF] Wu | From wireless sensor networks towards cyberphysical systems[END_REF].

Update Application

Updates are applied at the end devices using update agents [START_REF] Jurković | Remote Firmware Update for Constrained Embedded Systems[END_REF]. The update agents decode the received delta file and execute decoded operations. Depending on implementation, the new firmware image can be reconstructed in an external memory and loaded into the program memory after reboot [START_REF] Dong | R2: Incremental Reprogramming Using Relocatable Code in Networked Embedded Systems[END_REF], or it can be applied on the fly, when the update agent edits the program memory directly and runs new firmware right after all the operations from the delta file are executed [START_REF] Shafi | No-reboot and Zero-Flash Over-the-air Programming for Wireless Sensor Networks[END_REF]. This saves the requirement of reboot and firmware loading, but may result in device malfunction if the update is not applied properly or interrupted. The Rollback option from external memory should be therefore present on these devices if possible.

Flash Memories Energy Consumption Modeling

Creating an accurate model for energy consumption of flash memories is not an easy task. There are many parameters that must be taken into account. The energy model proposed and validated in [START_REF] Pallister | A high-level model of embedded flash energy consumption[END_REF] is focused on embedded flash memories. The flash memory controller activates different regions of flash. These regions are created by 2 k bytes and activation of each region consumes Ek energy. Formal representation:

𝑖 → 𝑗 = ∑ 𝐸 𝑘 𝑁(𝑖,𝑗) 𝑘=0 (1)
i and j represent memory address. Term N(i,j) represents the largest changed region -2 N(i,j) bytes, 𝑁(𝑖, 𝑗) = ⌊log 2 (𝑖 ⊕ 𝑗)⌋ . The authors of [START_REF] Bazzaz | An Accurate Instruction-Level Energy Estimation Model and Tool for Embedded Systems[END_REF] created and validated an instruction-level energy model for embedded systems. They define a memory access energy consumption of a process as Ememory:

𝐸 𝑚𝑒𝑚𝑜𝑟𝑦 = 𝐸 𝑐𝑛𝑡𝑟 + 𝐸 𝐹𝑙𝑎𝑠ℎ + 𝐸 𝑆𝑅𝐴𝑀 (2)
Effective Over-the-Air Reprogramming for Low-Power Devices 289

Ecntr is energy consumed by the memory controller. Eflash is energy consumed by the flash memory write, read and erase operations. ESRAM is energy consumed by write and read operations in RAM memory. This basic model serves us for theoretical evaluation of our differencing algorithm optimizations.

Research Contribution and Innovation

In our research of update energy consumption, apart from delta file size, we also take into account encoding of a delta file and update application. We propose optimizations of the differencing algorithms and delta file encoding for updates without external flash. We also create an energy consumption model that can evaluate our solution.

Differencing Algorithms for Updates without Use of an External Memory

Low-power devices can be updated without use of an external memory. Selfprogrammability allows the device to rewrite its program memory on the fly. Some low-power devices used in CPS for control may not have an external memory, so they have to be updated without it if needed. On the fly update does not require copying of identical firmware parts like reconstruction in the external memory. Update application in device's program memory therefore requires different delta file encoding.

Delta file encoding.

Updates executed on the fly require different delta file encoding. The traditional way that the reconstructs firmware mixes ADD and COPY operations.

Fig. 2 Delta file format for updates without external memory usage

Now, COPY operations must be carried out first before new data are added. We propose the delta file to be split into three parts -header, operation data and integrity check data (Fig. 2). The header consists of two numbers, the number of COPY and the number of ADD operations. The operation data encode addresses and bytes for each operation. Integrity check data encode CRC-16 code that is checked prior to the firmware update in order to prevent corrupted delta files from updating firmware.

Optimizing differencing algorithms. All analyzed differencing algorithms work with full firmware images. For two images of length n and m, the space complexity in the worse cases is O(n x m) [START_REF] Hu | Reprogramming with Minimal Transferred Data on Wireless Sensor Network[END_REF], in the better cases O(n + m) [START_REF] Dong | R3: Optimizing relocatable code for efficient reprogramming in networked embedded systems[END_REF] [START_REF] Mo | An efficient differencing algorithm based on suffix array for reprogramming wireless sensor networks[END_REF]. If we apply xor operation on these images and only compare non-matching sections of a new image to an old image, it is possible to improve algorithm's space complexity, which also results in better execution time. Let k be the number of non-matching bytes between n and m bytes of compared images. Space complexity is reduced to O(n + k) as k is less than m.

Non-matching segments are built using COPY and ADD operations. The generated operations may be sometimes so small that they can be merged into preceding or following operation. For ADD operation with cost x bytes (without data bytes) and COPY operation with cost y bytes, we propose following optimizations:

 If COPY copies at most y bytes, it can be merged into preceding ADD  If two consecutive COPY operations copy less than 2x bytes, they can be merged into ADD operation  If result of COPY operations merging are two consecutive ADD operations, they can be merged into one ADD operation

Modeling Energy Consumption of the Update

We define energy consumed by flash write operation as λflash(write). This is energy consumed by the erase and write operations for one page. Let s be the flash page size, addr the address we start writing data to and n the number of bytes we write. NAND memories can only be written one page at time, so the number P of pages we must write in order to store n bytes is:

𝑃 = ⌊(𝑎𝑑𝑑𝑟 + 𝑛) ÷ 𝑠⌋ -⌊𝑎𝑑𝑑𝑟 ÷ 𝑠⌋ + 1 (3)
Using this formula, we can also determine the number of pages we read from, but flash memory can read bytes directly. If energy needed to perform each read operation is λflash(read), then reading n bytes and writing them during COPY and ADD operations will result in flash memory energy consumption:

𝐸 𝑓𝑙𝑎𝑠ℎ = 𝑛 • 𝜆 𝑓𝑙𝑎𝑠ℎ(𝑟𝑒𝑎𝑑) + 𝑃 • 𝜆 𝑓𝑙𝑎𝑠ℎ(𝑤𝑟𝑖𝑡𝑒) (4)
Ecntr can be calculated based on equation [START_REF] Levis | Maté: A Tiny Virtual Machine for Sensor Networks[END_REF]. We omit small region activations, for example reading a file byte by byte, and take into account only jumps from location we read to page we write and back (2P jumps total). This usually causes activation of more regions. If jump between two addresses i and j requires activation of R regions, jumping between these addresses during writing of P pages will consume energy:

𝐸 𝑐𝑛𝑡𝑟 = 2𝑃 ∑ 𝐸 𝑘 𝑅 𝑘=0 (5)

Experimental Results and Discussion

We implemented the differencing algorithm called Delta Generator (DG) with proposed optimizations that generates delta files in the proposed format. The evaluation was done for 7 firmware change cases on ATmega32U4 microcontroller with 32KB of self-programming flash memory. The page size of this memory is 128 bytes. We compare our solution to R3diff differencing algorithm [START_REF] Dong | R3: Optimizing relocatable code for efficient reprogramming in networked embedded systems[END_REF], one of the best differencing algorithms that generate deltas for firmware updates in an external memory. For this purpose we created tool that calculated the amount of required read and write operations encoded in deltas based on equation [START_REF] Stathopoulos | A Remote Code Update Mechanism for Wireless Sensor Networks[END_REF]. We also calculated the approximate number of region activations (equation ( 5)) that update process required. Results of our analysis are shown in Table 1.

The experimental results show upgrade over existing solution -R3diff. The amount of read operations is roughly the same, but the proposed optimizations reduced the amount of page writes. The calculation of activated flash regions is not precise, byte by byte readings were excluded from the calculations. We were not yet able to determine exact power consumption of read and write operations for the flash memory used. The exact numbers could better show how not only size, but also encoding influences the energy efficiency of over-the-air updates. We also have not yet evaluated similarity improving techniques. We provided basic model that is a base for future improvements.

Conclusion and Further Work

We identified partial problems in the area of over-the air updates. In addition to this, we created the basic energy consumption model of reprogrammed memories and investigated the effect of delta file encoding on update energy performance. The current model is basic. It should be upgraded to a more detailed model and experimentally evaluated with real data. The executed experiments have shown good results. The vision is to create a precise model that takes into account similarity improvements, differencing algorithms, energy consumption of wireless interfaces, energy consumption of update agents and updated firmware. These are the main problems that influence effectivity of over-the-air updates. Perfection of solutions to these problems in combination with good energy consumption model can lead to fast, secure and energy effective updates for low-power devices in cyber-physical systems.

Table 1

 1 Comparison of R3diff [13] differencing algorithm to proposed Delta Generator (DG)

Acknowledgment. This work has been supported by Slovak national project VEGA 2/0192/15.


Luis Irio
email: l.irio@campus.fct.unl.pt

Rodolfo Oliveira

Real-Time Estimation of the Interference in Random Waypoint Mobile Networks

Keywords: Interference Estimation, ad hoc Networks, Mobility

It is well known that the stochastic nature of the interference deeply impacts on the performance of emerging and future wireless communication systems. In this work we consider an ad hoc network where the mobile nodes adopt the Random Waypoint mobility model. Assuming a time-varying wireless channel due to slow and fast fading and, considering the dynamic path loss caused by the node's mobility, we start by characterizing the interference caused to a receiver by the moving nodes positioned in a ring. Based on the interference distribution, we evaluate two different methodologies to estimate the interference in real-time. The accuracy of the results achieved with the proposed methodologies in several simulations show that they may be used as an effective tool of interference estimation in future wireless communication systems, being the main contribution of this work.

Introduction

Interference is an important metric in the future generation of wireless communication systems because the traditional single transmitter and receiver model is being progressively replaced by a different approach, where multiple nodes may transmit simultaneously for a single or even multiple receivers.

The interference in wireless mobile networks, and particularly its characterization, is important for many applications. In most wireless mobile scenarios the characterization of the interference is a non-trivial task. While several authors model the interference in non-mobile networks [START_REF] Win | A Mathematical Theory of Network Interference and Its Applications[END_REF], the assumption nodes' mobility introduces a novel degree related with the time-varying nature of nodes' positions. The works already published approaching a formal description of the interference in mobile networks are mainly focused on modeling. The use of statistics describing the level of mobility of the interferers in the modeling process was considered in [START_REF] Yarkan | Impact of Mobility on the Behavior of Interference in Cellular Wireless Networks[END_REF][START_REF] Zhang | Interference dynamics in MANETs with a random direction node mobility model[END_REF][START_REF] Gong | Interference and Outage in Mobile Random Networks: Expectation, Distribution, and Correlation[END_REF]. [START_REF] Yarkan | Impact of Mobility on the Behavior of Interference in Cellular Wireless Networks[END_REF] models the aggregate interference caused by static interferers, being considered that the nodes' mobility only causes a time-varying displacement with respect to the different non-mobile cells. [START_REF] Zhang | Interference dynamics in MANETs with a random direction node mobility model[END_REF] admits a mobile scenario where the nodes adopt the Random Direction mobility model. [START_REF] Gong | Interference and Outage in Mobile Random Networks: Expectation, Distribution, and Correlation[END_REF] considers that the mobile nodes adopt the Random Waypoint mobility model (RWP), but only the interference power from the nearest interferer to the receiver is considered. [START_REF] Irio | Aggregate Interference in Random Waypoint Mobile Networks[END_REF] proposes an interference model for ad hoc mobile networks where the nodes move in accordance with the RWP and all the contributions of the nodes located within a defined region are considered.

This work starts by characterizing the distribution of the interference caused to a receiver by multiple moving nodes located in a ring, considering path loss and slow and fast fading. Based on the interference distribution, we evaluate two different methodologies to estimate the interference in real-time. The major contribution of this work is the identification of a method to estimate the aggregate interference in random waypoint mobility networks, leading to accurate results when used in realtime.

The next section describes the main contributions of this work. Section 3 presents the general assumptions. In Section 4 the distribution of the interference values obtained through simulation is approximated by known distributions in order to identify possible approximations. Section 5 describes two estimation methodologies as well as the real-time estimates obtained through simulation and finally the conclusions are present in Section 6.

Relationship to Cyber-Physical Systems

Recently, Cyber-physical Systems have attracted much attention from the academic community. These systems are mainly focused on the link between computation and physical processes in terms of their reciprocal interaction. Instead of considering standalone physical devices, Cyber-Physical Systems adopt an integrated network of multiple physical devices to enrich the interactions and cooperation between the devices and the virtual worlds available through computation.

Recent advances in wireless communications systems and distributed wireless networks have supported a plethora of innovation in Cyber-Physical Systems. Significant progresses have been observed in mobile ad hoc networks and wireless sensor networks. Examples of Cyber-Physical Systems include mobile robotics and mobile sensors or actuators.

Our work contributes to the development of mobile Cyber-Physical Systems, by studying interference phenomena in mobile wireless networks formed without a central coordinator. By characterizing the interference caused by multiple mobile nodes, the wireless communication process can be improved. Consequently, mobile Cyber-Physical Systems may benefit when higher throughput or reliability is needed.

Basically, we show the impact of the mobility, in terms of average velocity of the cyber-physical devices, in the interference caused to a central receiver. We characterize the interference power at the receiver taking into account the specifics associated with the propagation and mobility scenario. In this way, we contribute to the advance of Cyber-Physical Systems, by proposing an effective solution to estimate the interference, which may be used for different purposes ranging from wireless energy harvesting to the improvement of the wireless communication system.

System Description

Mobility Assumptions

This work considers that the nodes move in accordance with the RWP mobility model [START_REF] Johnson | Dynamic Source Routing in Ad Hoc Wireless Networks[END_REF]. In a RWP model all nodes are firstly placed in a random position ,ݔ( .)ݕ ,ݔ( )ݕ is sampled from an uniform distribution denoted by ݔ ∈ [0, ܺ ௫ ] and ݕ ∈ [0, ܻ ௫ ]. ,ݔ( )ݕ denotes the starting point, and the following procedure is the definition of the ending point ( ݔ ᇱ , ݕ ᇱ ), which is uniformly selected as the starting point (i.e. ݔ ᇱ ∈ [0, ܺ ௫ ] and ݕ ᇱ ∈ [0, ܻ ௫ ] ). Afterwards a node samples the velocity ݒ ∈ [ܸ , ܸ ௫ ] from an uniform distribution, which is adopted to travel from the starting point to the ending point.

After arriving at the ending point ݔ( ᇱ , ݕ ᇱ ), a node selects the duration of a pause (ܶ ) during which it remains stopped at the ending point. After the time ܶ , a node selects another value for the velocity to travel to a different ending point. After arriving at the ending point, a node repeats the same procedure as many times as parameterized in the mobility simulations.

Considering that E[]ܮ represents the expected distance between two random points and E[ܸ ௪ ] represents the expected velocity of the nodes without considering pause, the expected velocity of the nodes considering pause is given by

E[ܸ] = ൭ E[]ܮ Eൣܸ ௪ ൧ ିଵ E[]ܮ + E[ܶ ] ൱ , (1)
where E[ܸ ௪ ] and E[]ܮ are defined in [START_REF] Bettstetter | The node distribution of the random waypoint mobility model for wireless ad hoc networks[END_REF] as

Eൣܸ ௪ ൧ = ቀ ೌೣ ି ୪୬( ೌೣ / ) ቁ , E[]ܮ ≈
521.405 m, and E[ܶ ] represents the expected value of the pause duration.

Network Scenario

The network scenario considered in this work assumes a RWP scenario where ݊ nodes travel in a rectangular region with area ܺ ௫ × ܻ ௫ . The network model considered in this work assumes that a static node ܰ is placed in the center of the considered scenario (located at (ܺ ௫ 2 ⁄ , ܻ ௫ 2 ⁄ )). ܰ is a receiver of the mobile transmitters. The objective of this paper is the study (in terms of statistics and estimation) of the interference caused to ܰ by the transmitting nodes ሼܰ ଵ , ܰ ଶ , … , ܰ ሽ positioned within the interference region, i.e. the mobile transmitters located within the ring bounded by the smaller circle of radius ܴ and the larger circle of radius ܴ , represented in Figure 1. Knowing that the interference depends on the distance between the transmitter and receiver, a circular model was chosen. The parameters describing the network and the mobility conditions are described in Table 1. 20 m/s

Radio Propagation Assumptions

This subsection describes the radio propagation scenario considered in this work.

The aggregate interference power received by the node ܰ located in the centre is expressed by

ܫ = ܫ ಲೃ ೃ ୀଵ , (2)
where ܫ is the interference caused by the node ݅ , and ݊ ோ ோ is the number of transmitters positioned in the interference area ߨ(ܴ ଶ െ ܴ ଶ ). The interference power ܫ is given by

ܫ = ܲ ்௫ ߰ ݎ ିఈ , (3)
where ܲ ்௫ is the transmitted power of the node ݅ (ܲ ்௫ = 10 ଷ mW). ߰ denotes the fading occurring in the channel between the node ܰ and a transmitter ݅ and ݎ is the distance between the transmitter ݅ and the receiver. Finally, ߙ denotes the path-loss coefficient. In this works we consider that the transmitters do not adjust its transmitting power (i.e. no power control is considered). The fading ߰ includes the small-scale fading and shadowing effects. The smallscale fading effect is assumed to be distributed in accordance with a Rayleigh distribution, which is represented by

݂ )ݔ( = ݔ ߪ ଶ ݁ ି௫ మ ଶఙ അ మ , (4)
where ݔ is the envelope amplitude of the received signal, and 2ߪ ଶ is the mean power of the multipath received signal. ߪ = 1 is adopted in this work.

Regarding the fading effect, we have assumed that it follows a Lognormal distribution

݂ క )ݔ( = 1 √2ߨ ߪ క ݔ ݁ ି(୪୬(௫)ିఓ) మ ଶఙ మ , (5)
where ߪ క is the shadow standard deviation when ߤ = 0. The standard deviation is usually expressed in decibels and is given by ߪ కௗ = 10ߪ క ln(10) ⁄ . For ߪ క → 0, no shadowing results. Although [START_REF] Irio | Aggregate Interference in Random Waypoint Mobile Networks[END_REF] appears to be a simple expression, it is often inconvenient when further analyses are required. Consequently, [START_REF] Abdi | On the utility of gamma PDF in modeling shadow fading (slow fading)[END_REF] has shown that the log-normal distribution can be accurately approximated by a gamma distribution, defined by

݂ క )ݔ( = 1 Γ(ߴ) ൬ ߴ ߱ ௦ ൰ ణ ݔ ణିଵ ݁ ି௫ ణ ఠ ೞ , (6)
where

ߴ is equal to 1 ቀ݁ ఙ మ െ 1ቁ ൗ
and ߱ ௦ is equal to ݁ ఓ ඥ(ߴ + 1) ߴ ⁄ . Γ(. ) represents the Gamma function. The probability distribution function of the fading ߰ is thus represented by

݂ ట )ݔ( = 2 Γ(ߴ) ൬ ߴ ߱ ௦ ൰ ణାଵ ଶ ݔ ణିଵ ଶ ܭ ణିଵ ቌඨ 4 ߴ ݔ ߱ ௦ ቍ , (7)
which is the Generalized-K distribution, where ܭ ణିଵ (. ) is the modified Bessel function of the second kind.

Characterization of the Interference Distribution

Following the assumptions considered in the previous section, several simulations were performed considering two different mobility scenarios:

• Regarding the propagation conditions, we have considered the following scenario:

• Radio scenario -ߙ = 2 and ߪ కௗ = 3 dB.

During the simulations, the interference power sensed by the node ܰ was sampled every second in order to compute its Cumulative Distribution Function (CDF). The samples acquired in 1000 different simulations, totaling a sample set of ݈ = 3 × 10 samples, were also used to determine the parameters of a set of different probability density functions (PDFs) using a maximum-likelihood (ML) fitting methodology. For each one of the considered PDF ݂, an average logarithm likelihood was defined as follows

݃ ො = 1 ݈ ln ݔ(݂ |Θ) ୀଵ , (8)
where Θ represents the parameters of the PDF and ݔ represents each individual sample. ML was used to maximize the likelihood in order to determine Θ, which is described as follows

Θ ொ = argmax ݃ ො(Θ; ݔ ଵ , … , ݔ ) . (9)
Figure 2 represents the CDFs computed with the parameters obtained in (9) for the Generalized Extreme Value (GEV) and Gamma distributions. As illustrated, the fitting obtained with the GEV distribution presents a better approximation for the two mobility scenarios. Because of this observation, the estimation methods proposed in the next section assume that the interference distribution follows a GEV distribution.

Interference Estimation

This section assumes that the aggregate interference can be approximated by a GEV (Generalized Extreme Value) distribution, being its PDF represented by

;ݔ(݂ ߪ, ߛ, ߤ) = 1 ߪ )ݔ(ݐ ఊାଵ ݁ ି௧(௫) , (10)
where

)ݔ(ݐ = ቐ ቀ1 + ߛ ݔ െ ߤ ߪ ቁ ିଵ ఊ ൗ , ߛ ് 0 ݁ ି(௫ିఓ) ఙ ⁄ , ߛ = 0 . (11)
A Maximum Log-likelihood estimator (MLE) and a Probability Weighted Moments (PWM) estimator are introduced in the next subsections, in order to be used in realtime to estimate the aggregate interference. Hereafter, we denote the elements of an interference sample set by ߯ = ܺ ଵ , ܺ ଶ , … , ܺ . We also consider the ordered sample set, which is denoted by ܺ ଵ, ⋯ ܺ , .

Log-Likelihood Estimator

The log-likelihood function for a sample set ߯ = ሼܺ ଵ , … , ܺ ሽ of i.i.d GEV random variables is given by

log ܮ (ߪ, ߛ, ߤ) = െ݉ log ߪ െ ቀ 1 ߛ + 1ቁ ∑ log ቀ1 + ߛ ܺ െ ߤ ߪ ቁ ୀଵ െ ∑ log ቀ1 + ߛ ܺ െ ߤ ߪ ቁ ିଵ ఊ ൗ ୀଵ , (12)
under the condition 1 + ߛ ିఓ ఙ > 0 . The MLE estimator (ߪ, ෝ ߛ ො, ߤ) for (ߪ, ߛ, ߤ) is obtained by maximizing (12).

PWM Estimator

As described in [START_REF] Greenwood | Probability weighted moments: definition and relation to parameters of several distributions expressable in inverse form[END_REF], the PWM of a random variable ܺ with distribution function )ܺ(ܨ = ܲ(ܺ )ݔ are the quantities

ܯ ,,௦ = Eൣܺ ൫(ܨX)൯ ൫1 െ (ܨX)൯ ௦ ൧ , (13)
for real , ݎ and ݏ values. For the GEV distribution, [START_REF] Hosking | Estimation of the generalized extreme-value distribution by the method of probability-weighted moments[END_REF] shows that Eൣܺ ൫(ܨX)൯ ൧ can be written as

ܯ ଵ,, = 1 ݎ + 1 ൜ߤ െ ߪ ߛ [1 െ ݎ( + 1) ఊ Γ(1 െ ߛ)]ൠ, (14)
with ߛ < 1 and ߛ ് 0 .The PMW estimators (ߪ, ෝ ߛ ො, ߤ) of the Generalized Extreme Value distribution parameters (ߪ, ߛ, ߤ) may be computed through the following system of equations

ە ۖ ۖ ۔ ۖ ۖ ۓ ܯ ଵ,, = ߤ െ ߪ ߛ (1 െ Γ(1 െ ߛ)) ܯ2 ଵ,ଵ, െ ܯ ଵ,, = ߪ ߛ Γ(1 െ ߛ) (2 ఊ െ 1) ܯ3 ଵ,ଶ, െ ܯ ଵ,, ܯ2 ଵ,ଵ, െ ܯ ଵ,, = 3 ఊ െ 1 2 ఊ െ 1 , (15)
in which ܯ ଵ,, is replaced by the unbiased estimator proposed in [START_REF] Landwehr | Probability weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles[END_REF] ܯ

ଵ,, = 1 ݉ ൭ෑ ݆ െ ݈ ݉ െ ݈ ୀଵ ൱ ܺ , ୀଵ . (16)

Simulation Results

Figure 3 presents the simulation results obtained for the same scenario adopted in Figure 2. The "Simulation" curve represents the CDF obtained with the entire set of samples (3 × 10 samples). To apply the MLE and the PWM estimators in real-time we have considered a sample set ߯ of length ݉ = 100 samples. The estimators were computed for 20 different sets of samples, thus 20 different CDFs were computed (one per set). The CDF presented in Figure 3 is the average of the 20 CDFs computed for each sample set. All of the algorithms to find the estimators were solved by using the MATLAB.

Regarding the accuracy of the proposed estimators, both MLE and PWM present high accuracy. As a final remark, the results presented in Figure 3 validate the proposed estimation methodologies, being the PWM estimator more adequate for the real-time estimation due to its higher accuracy. Finally, we highlight that approximate results were observed for smaller sample set sizes using the PWM estimator, and similar results may be achieved using only ݉ = 10 samples per sample set, which is a remarkable low number of samples. The results are not show in the paper due to lack of space.

Conclusions

In this work we consider an ad hoc mobile network where the nodes move in accordance with the Random Waypoint mobility model. Assuming a time-varying wireless channel due to slow and fast fading and, considering the dynamic path loss due to the mobility of the nodes, we start by characterizing the interference distribution caused to a receiver by the mobile interferers located in a ring. The simulation results confirmed that the distribution of the aggregated interference may be accurately approximated by a Generalized Extreme Value distribution. Based on the interference distribution, two different methodologies (MLE and PWM) were assessed to estimate the interference in real-time. The accuracy of the results achieved with the proposed methodologies show that they may used as an effective tool of interference estimation in future wireless communication systems. Moreover, the low number of required samples constitutes one of the advantages of the proposed PWM estimator, even when the samples are highly correlated.

Fig. 1 .

 1 Fig. 1. Interference sensed by ܰ due to the mobile interferers located in the annulus area ߨ൫ܴ ଶ െ ܴ ଶ ൯.

Mobility scenario 1 -

 1 ܸ = 5 m/s, ܸ ௫ = 20 m/s, and ܶ = 0 s, representing an average node's velocity E[ܸ] = 10.82 m/s; • Mobility scenario 2 -ܸ = 5 m/s, ܸ ௫ = 20 m/s, and ܶ = 300 s, representing an average node's velocity E[ܸ] = 1.50 m/s.

Fig. 2 .

 2 Fig. 2. Interference CDF.

Fig. 3 .

 3 Fig. 3. Simulation results obtained with the MLE and PWM estimators.

Table 1 .

 1 Parameters adopted in the simulations.

	ܺ ௫ ܻ ௫ Simulation time ܸ ܸ ௫	1000 m 1000 m 3000 s 5 m/s	100 0 s ; 300 s 20 m ܴ 120 m ݊ ܶ ܴ

Acknowledgments. The authors gratefully acknowledge financial support from the Portuguese Science and Technology Foundation (FCT/MEC) through the project ADIN-PTDC/EEITEL/2990/2012 and the grant SFRH/BD/108525/2015.