Evolutionary Optimization Strategy for Indoor Position Estimation Using Smartphones

: Due to their distinctive presence in everyday life and the variety of available built-in sensors, smartphones have become the focus of recent indoor localization research. Hence, this paper describes a novel smartphone-based sensor fusion algorithm. It combines the relative inertial measurement unit (IMU) based movements of the pedestrian dead reckoning with the absolute ﬁngerprinting-based position estimations of Wireless Local Area Network (WLAN), Bluetooth (Bluetooth Low Energy—BLE), and magnetic ﬁeld anomalies as well as a building model in real time. Thus, a step-based position estimation without knowledge of any start position was achieved. For this, a grid-based particle ﬁlter and a Bayesian ﬁlter approach were combined. Furthermore, various optimization methods were compared to weigh the different information sources within the sensor fusion algorithm, thus achieving high position accuracy. Although a particle ﬁlter was used, no particles move due to a novel grid-based particle interpretation. Here, the particles’ probability values change with every new information source and every stepwise iteration via a probability-map-based approach. By adjusting the weights of the individual measurement methods compared to a knowledge-based reference, the mean and the maximum position error were reduced by 31%, the RMSE by 34%, and the 95-percentile positioning errors by 52%.


Introduction
The automatic determination of a location in an earth-related reference system is essential in many areas of life and economy. This is the reason why global navigation satellite systems (GNSS), such as GPS, GLONASS, or Galileo, have become of high economic importance. The range of applications for automatic position determination extends from ascertaining the location of a person in the leisure sector (e.g., for hiking or sailing) up to locating people (e.g., for rescue and emergency forces) and tracking objects (e.g., goods and merchandise) to determine the position and orientation of air, water, and land vehicles for navigation and autonomous navigation purposes.
Whereas GNSS is typically used outdoors, its application is usually limited due to shading or attenuation as well as other propagation effects (fading, multipath) of the electromagnetic signals within buildings, underground, or in built-up areas ("indoor"). However, most of the aforementioned applications would also benefit from automatic indoor positioning [1,2]. As a result, numerous local positioning systems have been developed in the past to enable the automatic localization of people or objects in "GNSS-free" space [3,4]. In contrast to outdoor areas, however, there is no positioning technology analogous to GNSS. Instead, different technologies are used depending on the specific application.
For pedestrian indoor localization, smartphones have increasingly gained attention due to their distinctive presence in everyday life and the variety of available built-in sensors and communication technologies. and magnetometers for relative continuous localization [21,22]. Since IMU are embedded as Micro-Electro-Mechanical Systems (MEMS) in smartphones, they can be directly used for pedestrian dead reckoning. However, the sole use of MEMS-IMU for pose estimation is only feasible for a short period of time due to its drift behavior [23,24]. Alternatively, high-performance IMU sensors with high weight and power consumption could achieve high-quality measurements, but these are unsuitable for pedestrian localization.
Besides infrastructure-based and infrastructure-less systems, hybrid systems are becoming a standard for pedestrian localization with smartphones. In these systems, multisource information like IMU data and additional infrastructure-based systems (e.g., WLAN fingerprinting, BLE beacons, building model) are fused [25].
In [26], the hybrid multi-source information fusion of a small low-cost inertial navigation system (INS) with GPS, barometer, and a camera is described. In [27][28][29], the data of low-cost smartphone sensors (IMU, magnetometer, barometer) and approximate estimation methods (fingerprinting) based on communication technologies (WLAN, BLE) are combined with digital building models for pedestrian localization. In [30], fuzzy logic was used with an abstracted building model with 90 degree turns and BLE signals as landmarks. Furthermore, in ref. [31] the authors used fuzzy logic to consider the noisy and uncertain WLAN or BLE data. Moreover, robust multi-source information fusion positioning technology is shown in [32,33], while [34] identifies hybrid systems as most suited to mass-market applications and [35,36] explicitly focuses on IPS solutions applicable to smartphones. In the most recent work of [37], a survey of the different positioning techniques for collaborative localization is given.
While the localization in the outdoor areas is solved, there are numerous problems to solve in indoor areas. Therefore, to categorize the performance and compare different localization algorithms, the authors of [38] identified six key challenges of positioning. These challenges are accurate and continuous localization, scaling to multi-story buildings, signal bias adaptation at scale, upscaling to large numbers of measurements, identifying varying walking profiles, and observing signal delays. They also evaluated their localization algorithm based on these challenges for a shopping mall, including three multi-story buildings and a large open underground passageway.
To solve the problem of accurate and continuous localization, multi-source information fusion positioning from different sensors often uses Recursive Bayesian Estimation methods. Their basic idea is to estimate the current state of the system, taking into account all previous state measurements up to this point in time and represent the state estimation by a probability density function (PDF) [39,40]. Well-known forms of Recursive Bayesian Estimation algorithms are the Kalman filter (KF) [41] with its variants (Iterative Kalman Filter (IKF), Adaptive Kalman Filter (AKF)) and extensions (e.g., Extended Kalman Filter (EKF), Unscented Kalman Filter (UKF)) [42][43][44][45][46], or Sequential Monte-Carlo (SMC) methods [47].
The KF predicts the actual state based on the previous system state and linear functions, but the state transitions are often nonlinear. Therefore, the state transitions are approximated in the EKF with the tangent of a nonlinear function and Gaussian noise [39]. In contrast, the UKF uses a weighted, statistical, linear regression process for linearization. Thus, instead of the tangent, so-called sigma points are extracted from the input distribution and transformed using a nonlinear function.
The sigma points are based on the mean value and are symmetrically oriented using the covariance [46]. However, the KF and its variants are based on a Gaussian probability distribution. Alternatively, SMC filters, also known as PF [48], can be applied for state estimation since they can approximate any (non-Gaussian) distribution as well as linear and non-linear system models. The PF also often uses Gaussian functions to reflect the sensor uncertainties [49,50]. In particular, PF, in the context of (robot) localization also referred to as Monte Carlo Localization (MCL), has proven as a suitable method for multi-source information fusion for pedestrian localization, e.g., [8,12]. In [51][52][53][54], the digital building model was integrated into a PF to evaluate the particles against the building structure, while [55][56][57] used an abstracted discretized grid-based building representation.
Because of the digital building model's nonlinearities in combination with nonlinear pedestrian moving patterns and stochastically intricate representations of the measurement signals of the hybrid positioning systems, this work weights its information sources explicitly to improve the localization accuracy. Therefore, it aims to reproduce a system's PDF by weighted particles and reflect the sensor uncertainties with Gaussian functions, while a detailed grid-based building representation and a novel particle interpretation as localized grid cells are investigated. Thereby, the number of particles was drastically increased and the performance was enhanced via fast and efficient 2D matrix multiplications.
The input information sources (position data, sensor observations, building models) are weighted using optimization methods to determine the actual position with high accuracy. Therefore, classical optimization methods (e.g., hill climbing or gradient descend algorithms like the Nelder-Mead method) require good initial guesses for the global extreme estimation [58]. Otherwise, a convergence against local extrema may occur and the algorithm fails to estimate the global optimum. Classical optimization schemes like Gauss-Newton [59] or Levenberg-Marquart [60] show the problem that they need continuously differentiable functions for optimization, which are not available due to the complexity of the environment described in this work.
Thus, to avoid local minima in the search process, global optimization strategies can be used. One common strategy known for its robustness end flexibility is the GA [61][62][63][64][65]. GA can be assigned to heuristic optimization methods. These algorithms are used specifically if the functions that have to be optimized are nonlinear or the calculations are computationally intensive. This could be the case if the amount of data is either complex or extensive. Moreover, there are no restrictive requirements for the target function, which does not need to be differentiable or continuous. Furthermore, the GA's basic functionality is easy to understand, the algorithm can be applied even with little knowledge of the problem structure, and functions with numerous dependent and independent parameters can be optimized [66].

Basic Localization Strategy
As described in Section 2, PF are often applied for pedestrian navigation in buildings. In a PF implementation, the system states are discretized by numerous independent weighted particles that approximate the PDF at each point in time. In the case of position estimation, the particles represent possible concrete positions in space.
After the initialization, the recursive update steps (2) to (4) are performed continuously for each filter epoch. As soon as the particles cluster around one point (convergence), they indicate the current position.
The present work has been built on the approach presented in [68,69], whose PFbased algorithm is briefly described in the following. In the initialization, numerous (500-2000) independent particles are randomly distributed in the walkable areas of a building and outside the walls. This step is followed by prediction sampling, which predicts the actual location based on PDR. Therefore, the observations of the built-in accelerometer, gyroscope, and magnetometer smartphone sensors are used. After that, the importance sampling weights each particle concerning the difference between the predicted and observed position in step three.
Finally, the weights are normalized in the resampling step and particles with low weights are replaced with new and randomly distributed ones. Thus, the PDF p(x k |Z k ) is Electronics 2021, 10, 618 5 of 23 estimated with the observations from step 0 to k, the current posteriori probability x i k , the individual weights w i k , and the Dirac delta function δ: Here, the weighting is defined as a two-dimensional standard normal distribution with (x i s , y i s ) as the predicted and (x k ,y k ) as the observed position as follows: (2) If a particle moves inside a wall-according to the building structure taken from a digital building model-it is deleted, or its weight is decreased drastically.
In addition, further existing infrastructure-based positioning systems such as WLAN fingerprinting, proximity estimates based on BLE beacons, and magnetic field anomaly maps can be integrated. For this purpose, probability density distributions (e.g., normal or uniform distribution) are assumed for each positioning information and fused with the PF estimate via a Mixture model [70]. Through this approach, the probability distribution of the sensor measurements of the PDR, BLE proximity and WLAN fingerprinting, and the magnetic field anomalies can be approximated and fused via the Bayesian Filter approach (described by Equation (3)) [39,71,72]. The likelihood p(z k |x k ) of the measurement model represents the system state x k and its correlation to the current observation z k at timestamp k. The a-priori probability p(x k |Z k−1 ) represents the information of the model and p(z k |Z k−1 ) is used for normalization. This filter fuses the respective position estimations of the involved localization techniques using multiplication and normalization to obtain a single position output: Subsequently, this recursive plausibility test removes some particles and the majority of the particles will group around the actual position. If the particles group in several clusters, strategies like the kernel density estimation with mixture models can help choose one group representing the actual position [73].

Extended Localization Strategy
As an essential extension of the previously described basic approach, the particle representation was modified in this work from free-moving particles to a grid-based interpretation with fixed possible particle locations. For this, the building model given in vector format was transferred from its global Universal Transverse Mercator (UTM) coordinates to a local coordinate system and discretized to a raster grid. Thus, each particle is interpreted by a grid cell, which represents a quadratic element (or pixel) of a 2D matrix. This representation of the building model is shown in Figure 1a. Figure 1b illustrates a section of this model in case of an unknown starting position. Since the actual starting position is unknown, the presented PF does not initialize the possible position via random particles. Instead, each pixel of the initial matrix in Figure 1a (yellow cells) represents a start hypothesis of the current position (I 0 ), and only these pixel values can change. With this changed particle interpretation, the number of usable particles was drastically increased while keeping the computational load low via fast and efficient 2D matrix multiplications on the Graphical Processing Unit (GPU).
particles. Instead, each pixel of the initial matrix in Figure 1a (yellow cells) represents a start hypothesis of the current position ( ), and only these pixel values can change. With this changed particle interpretation, the number of usable particles was drastically increased while keeping the computational load low via fast and efficient 2D matrix multiplications on the Graphical Processing Unit (GPU). To reflect the building topology's valid positions, wall-free areas, e.g., rooms, floors, door areas, and stairs, are characterized as possible starting positions (yellow cells). In contrast, wall or outdoor areas are invalid starting positions (blue cells). In this exemplary cell array of 556 × 335 pixels with an edge length of 0.15 m, there exist 76.614 valid out of 186.260 possible initial positions. Only grid cells that are likely to be entered get a value of one. In contrast, inaccessible areas get zeros pixel values and waist-high obstacles like fixed desks or benches get a value in between to reflect the unlikelihood of passing through such objects.
The process of our position estimation algorithm is shown in Figure 2. In Figure 3, the grid-based building model and the related walkway structure (Figure 3a), as well as example maps of a BLE ( Figure 3b) and a WLAN fingerprinting-based radio map ( Figure  3c), and the magnetic field anomaly map (Figure 3d) are shown. All information, therefore, was converted into a grid cell representation of the corresponding PDF.
The walkway structure (Figure 3a) is an assumption of the paths in a building based on the skeletonization of the walkable indoor area. This map's highest values (yellow cells) represent areas with a high distance to the walls nearby inside the building.
The BLE (Figure 3b) and the WLAN (Figure 3c) fingerprinting radio map illustrate a position estimation based on a weighted k-nearest-neighbors (WKNN) fingerprinting (k=5) [74]. If a position can be estimated via BLE or WLAN-based fingerprinting, each of To reflect the building topology's valid positions, wall-free areas, e.g., rooms, floors, door areas, and stairs, are characterized as possible starting positions (yellow cells). In contrast, wall or outdoor areas are invalid starting positions (blue cells). In this exemplary cell array of 556 × 335 pixels with an edge length of 0.15 m, there exist 76.614 valid out of 186.260 possible initial positions. Only grid cells that are likely to be entered get a value of one. In contrast, inaccessible areas get zeros pixel values and waist-high obstacles like fixed desks or benches get a value in between to reflect the unlikelihood of passing through such objects.
The process of our position estimation algorithm is shown in Figure 2. In Figure 3, the grid-based building model and the related walkway structure (Figure 3a), as well as example maps of a BLE ( Figure 3b) and a WLAN fingerprinting-based radio map (Figure 3c), and the magnetic field anomaly map (Figure 3d) are shown. All information, therefore, was converted into a grid cell representation of the corresponding PDF.
The walkway structure ( Figure 3a) is an assumption of the paths in a building based on the skeletonization of the walkable indoor area. This map's highest values (yellow cells) represent areas with a high distance to the walls nearby inside the building.
The BLE (Figure 3b) and the WLAN (Figure 3c) fingerprinting radio map illustrate a position estimation based on a weighted k-nearest-neighbors (WKNN) fingerprinting (k = 5) [74]. If a position can be estimated via BLE or WLAN-based fingerprinting, each of these positions was convoluted with a symmetric 2D Gaussian function to reflect the position uncertainties based on previous measurements. The magnetic field anomaly map (Figure 3d) was estimated using the magnetic anomalies caused by ferromagnetic objects, in the present case, especially fire-protection doors. These anomalies influence the magnetometer if the smartphone is near these doors. Therefore, the areas of the fire-protection doors are convoluted with a symmetric 2D Gaussian function with a standard deviation (σ) of 1.0 m. The used values of σ were postulated in [5]. For the estimation of the positions, the measurement uncertainties of the PDR have to be considered. Therefore, the empirically determined probability distributions of the step size estimation (Figure 4a, σ = 0.15 m) and the heading estimation ( Figure 4b, σ = 40°) are fused to a resulting step probability distribution ( Figure 4c). The current position distribution matrix is then superimposed with the resulting probability distribution of step estimation at each position after a determined movement. This process corresponds to the discrete 2D convolution formulated in Equation (4) and predicts the actual movement. I is the image and I* represents the result after a discrete convolution with the convolution these positions was convoluted with a symmetric 2D Gaussian function to reflect the po-sition uncertainties based on previous measurements.
The magnetic field anomaly map (Figure 3d) was estimated using the magnetic anomalies caused by ferromagnetic objects, in the present case, especially fire-protection doors. These anomalies influence the magnetometer if the smartphone is near these doors. Therefore, the areas of the fire-protection doors are convoluted with a symmetric 2D Gaussian function with a standard deviation (σ) of 1.0 m. The used values of σ were postulated in [5]. For the estimation of the positions, the measurement uncertainties of the PDR have to be considered. Therefore, the empirically determined probability distributions of the step size estimation (Figure 4a, σ = 0.15 m) and the heading estimation ( Figure 4b, σ = 40°) are fused to a resulting step probability distribution ( Figure 4c). The current position distribution matrix is then superimposed with the resulting probability distribution of step estimation at each position after a determined movement. This process corresponds to the discrete 2D convolution formulated in Equation (4) and predicts the actual movement. I is the image and I* represents the result after a discrete convolution with the convolution The magnetic field anomaly map (Figure 3d) was estimated using the magnetic anomalies caused by ferromagnetic objects, in the present case, especially fire-protection doors. These anomalies influence the magnetometer if the smartphone is near these doors. Therefore, the areas of the fire-protection doors are convoluted with a symmetric 2D Gaussian function with a standard deviation (σ) of 1.0 m. The used values of σ were postulated in [5].
For the estimation of the positions, the measurement uncertainties of the PDR have to be considered. Therefore, the empirically determined probability distributions of the step size estimation (Figure 4a, σ = 0.15 m) and the heading estimation ( Figure 4b, σ = 40 • ) are fused to a resulting step probability distribution ( Figure 4c). The current position distribution matrix is then superimposed with the resulting probability distribution of step estimation at each position after a determined movement. This process corresponds to the discrete 2D convolution formulated in Equation (4) and predicts the actual movement. I is the image and I* represents the result after a discrete convolution with the convolution kernel k R n×n . The coordinates of a pixel are x and y, n is the kernel's size in the first dimension, while a is the center of the convolution kernel's size: Electronics 2021, 10, x FOR PEER REVIEW 8 of 23 kernel ℝ × . The coordinates of a pixel are x and y, n is the kernel's size in the first dimension, while a is the center of the convolution kernel's size: Due to the sensors' measurement uncertainties, the current position at building junctions could be determined wrong, and therefore, a false path would be taken. This process usually leads to a state where the determined position gets stuck in a wall. This scenario, known as the Kidnapped Robot Problem [39], is solved automatically by the algorithm without detecting a faulty state because the particles are not erased, and the grid cell values are not set to zero in a large environment around the current position. This is achieved by shifting the wall areas (blue areas in Figure 1a) for each step with the respective step length in pixel according to the movement's heading using Bresenham's line algorithm [75]. Using this algorithm, all position updates (step length and heading) are translated to pixel-wise shifts to reflect the raster-discretized translation.
All grid cell values, which were implausible to represent the current step's start, are attenuated by 90%. Therefore, the correct position can be estimated again if it was lost after only a few steps. This procedure reduces the values of wall areas and the areas which cannot be reached with the current step probability distribution. This weakening process is shown in Figure 5a after five steps.  Due to the sensors' measurement uncertainties, the current position at building junctions could be determined wrong, and therefore, a false path would be taken. This process usually leads to a state where the determined position gets stuck in a wall. This scenario, known as the Kidnapped Robot Problem [39], is solved automatically by the algorithm without detecting a faulty state because the particles are not erased, and the grid cell values are not set to zero in a large environment around the current position. This is achieved by shifting the wall areas (blue areas in Figure 1a) for each step with the respective step length in pixel according to the movement's heading using Bresenham's line algorithm [75]. Using this algorithm, all position updates (step length and heading) are translated to pixel-wise shifts to reflect the raster-discretized translation.
All grid cell values, which were implausible to represent the current step's start, are attenuated by 90%. Therefore, the correct position can be estimated again if it was lost after only a few steps. This procedure reduces the values of wall areas and the areas which cannot be reached with the current step probability distribution. This weakening process is shown in Figure 5a after five steps.  Due to the sensors' measurement uncertainties, the current position at building junctions could be determined wrong, and therefore, a false path would be taken. This process usually leads to a state where the determined position gets stuck in a wall. This scenario, known as the Kidnapped Robot Problem [39], is solved automatically by the algorithm without detecting a faulty state because the particles are not erased, and the grid cell values are not set to zero in a large environment around the current position. This is achieved by shifting the wall areas (blue areas in Figure 1a) for each step with the respective step length in pixel according to the movement's heading using Bresenham's line algorithm [75]. Using this algorithm, all position updates (step length and heading) are translated to pixel-wise shifts to reflect the raster-discretized translation.
All grid cell values, which were implausible to represent the current step's start, are attenuated by 90%. Therefore, the correct position can be estimated again if it was lost after only a few steps. This procedure reduces the values of wall areas and the areas which cannot be reached with the current step probability distribution. This weakening process is shown in Figure 5a after five steps.   The procedure of superimposing the probability distribution map (I Z−1 ) with the current step distribution (k CSD R n×n ) represents the prediction step. In contrast, the Hadamard products (element-wise multiplications ) of this result with the maps of walkway structure (I WS ), BLE (I BT ), WLAN (I W LAN ), and of the magnetic field anomalies (I MFA ), as well as of the weakening of implausible areas (I weak ) corresponds to the importance sampling or update step in the classical PF approach.
These calculations and the normalization are repeated recursively for each step. Therefore, the current position is determined from each previous position estimate and the current sensor data. This process and the respective probability distributions are shown in Figure 5b-d after 2, 5, and 20 steps. Furthermore, the video "S1 Extended Localization" of the Supplementary Materials illustrates the step-by-step positioning determination.
A maximum of three decimal places was allowed for each parameter to limit the parameter combinations' depth. For instance, the combination [0.153; 0.982; 0.349; 0.787; 0.622] is one of the possible individual parameter vectors.
In Equation (5), a pixel's coordinates are x and y, while a is the center of the convolution kernel's size (k CSD ).
To avoid zero-pixel values and ensure the further usability of cells, the distributions of the walkway structure, BLE, WLAN, and magnetic field anomalies are multiplied by a weighting factor w i (interval from 0 to 1) and an offset of +1 is added: To reflect the discrepancy between the possible body height, leg length, personal walking speed, and walking style between the examined pedestrian of [76] and the one in this work, an adjustment can be made. Therefore, the step length could be changed within the parameter optimization by multiplying a weighting factor 2·w 5 (interval from 0 to 1) to the estimated step length before the steps are discretized. Nevertheless, the minimal step length is equal to one cell's size and the maximal is equal to the doubled estimated step length. During the first steps or at building junctions, I Z will have multiple peaks, but only the highest cell value indicates the current position. The knowledge-based parameter combination is [1; 1; 1; 1; 0.5], which represents the equal trust in the fused position information of the walkway structure, BLE fingerprinting, WLAN fingerprinting, magnetic field anomalies, and the step length estimation.

Optimization Strategies
For determining the weights w 1 to w 5 , used in the prediction step for fusion the position information of the different sources (Equation (5)), in general, optimization methods can be applied. For the implementation the procedure was to carry out test runs in buildings with ground truth points (GTPs) and then to estimate the weights with a root-mean-square error (RMSE)-based fitness metric shown in Equation (6). To obtain small values for the distances between the GTPs and the corresponding estimated positions, the product of the RMSE values' reciprocals for each involved test run is used. Furthermore, for better numerical handling, a scaling factor of 10 was chosen: Due to the high ambiguity and the unknown (mixed) stochastical distributions of the input variables, the optimization quickly converges to a local minimum when using classical search methods (e.g., gradient-based methods). Therefore, global optimization procedures are recommended to maximize the quality metric and avoid local extrema.
However, because of the building model's nonlinearities, the used 2D maps, and the cropped steps, jumps and plateaus of the metric values exist. This also makes the stochastic description of the weights difficult and no gradients can be easily derived analytically. Thus, gradient-based optimization methods are not applicable and only function evaluations appear expedient.
Hence, in this work a GA is investigated to avoid the convergence against local extrema and to find optimal weights despite the complex optimization problem. For evaluation and benchmarking, the GA is compared with four other local and global optimization strategies (Hill Climbing, Nelder-Mead method, Simulated Annealing, Particle Swarm) to illustrate the necessity for more complex and robust search techniques.

Genetic Algorithm
Since especially the GA is investigated in this work, it will be introduced first in this section.
The GA is a flexible and robust optimization strategy whose general structure is shown in Figure 6. The GA is based on the observation of nature in which biological evolution (following Darwin's thesis) has led to complex life forms optimally adapted to their environment. Therefore, it can be assumed that not only the results of biological evolution are optimal, but also its mechanisms [77]. Accordingly, a potential solution in the parameter space is abstracted as an individual in an artificial environment [78] and successful patterns are combined to achieve better results. Here, a specific weighting constellation is called an individual.
Due to the high ambiguity and the unknown (mixed) stochastical distributions of the input variables, the optimization quickly converges to a local minimum when using classical search methods (e.g., gradient-based methods). Therefore, global optimization procedures are recommended to maximize the quality metric and avoid local extrema.
However, because of the building model's nonlinearities, the used 2D maps, and the cropped steps, jumps and plateaus of the metric values exist. This also makes the stochastic description of the weights difficult and no gradients can be easily derived analytically. Thus, gradient-based optimization methods are not applicable and only function evaluations appear expedient.
Hence, in this work a GA is investigated to avoid the convergence against local extrema and to find optimal weights despite the complex optimization problem. For evaluation and benchmarking, the GA is compared with four other local and global optimization strategies (Hill Climbing, Nelder-Mead method, Simulated Annealing, Particle Swarm) to illustrate the necessity for more complex and robust search techniques.

Genetic Algorithm
Since especially the GA is investigated in this work, it will be introduced first in this section.
The GA is a flexible and robust optimization strategy whose general structure is shown in Figure 6. The GA is based on the observation of nature in which biological evolution (following Darwin's thesis) has led to complex life forms optimally adapted to their environment. Therefore, it can be assumed that not only the results of biological evolution are optimal, but also its mechanisms [77]. Accordingly, a potential solution in the parameter space is abstracted as an individual in an artificial environment [78] and successful patterns are combined to achieve better results. Here, a specific weighting constellation is called an individual. In the context of this work, the (μ, λ)-strategy (μ-parents and λ-children) is pursued [79], and as in biological evolution, a limited lifespan of individuals [78] is used. The maximum age of an individual is 5% of the number of generations and the population size is kept constant in every generation.
However, it is crucial that the optimum search is global and does not just converge to a local optimum. Therefore, inferior individuals also have a small probability of getting selected [80].
The genetic operator's selection, recombination, and mutation create the next generation from the random starting parameter space (the starting population). This process is repeated until most parameter combinations converge or a maximum number of genera- Figure 6. Structure of the Genetic Algorithm: The sequence of selection, recombination, mutation, and substitution is iteratively repeated until a stopping criterion is reached.
In the context of this work, the (µ, λ)-strategy (µ-parents and λ-children) is pursued [79], and as in biological evolution, a limited lifespan of individuals [78] is used. The maximum age of an individual is 5% of the number of generations and the population size is kept constant in every generation.
However, it is crucial that the optimum search is global and does not just converge to a local optimum. Therefore, inferior individuals also have a small probability of getting selected [80].
The genetic operator's selection, recombination, and mutation create the next generation from the random starting parameter space (the starting population). This process is repeated until most parameter combinations converge or a maximum number of generations has been reached. To achieve progress in the optimization, selecting the best individuals of a generation is required in each optimization step based on a quality metric [81].
Because the computation time for sorting is negligible in this work, and to prevent a convergence occurring directly at the beginning, a ranking selection is used in contrast to the frequently used recombination schemes of roulette wheel or competition selection. This selection applies the sorted rank instead of the absolute fitness values and determines whether an individual will be selected for recombination. A dominance of single individuals and convergence to a local extremum is thus prevented.
In this approach, individuals with rank 1 represent those with the lowest fitness value, whereas rank n describes the best [64]. To determine the probability of selection, [82] suggests the following procedure: the expected value E max is assigned (1 ≤ E max ≤ 2) to the highest-ranked individual. The value E min for the worst individual is estimated from E min = 2 − E max and the probability of an arbitrary individual E(i) to be chosen is calculated using rank(i). Thus E(i) is the linear index of the i-th element in the sorted array of length n: A real-valued genetic representation of each individual with five property vectors between 0 and 1 is used in this work. To take advantage of different recombination techniques, 80% of the next generation is created by intermediate and 20% by combining recombination. On the one hand, intermediate recombination is used in which each parameter of a new individual is represented as the average of the respective parent's parameters [79]. On the other hand, the combining recombination recomposes the details of different individuals and thus, the advantageous components of the parent individuals can be combined in the optimal case [83]. Here, individuals can also be combined with themselves. To avoid the resulting more homogeneous population and convergence without reaching the optimum, the individuals' parameters are mutated to prevent diversity loss. Thus, a local optimum can be left again and the search space is explored extensively.
Each child individual is created by recombination of two parents and subsequent mutation. This process concentrates the search space around promising areas of the model space. A mutation in the GA describes the arbitrary modification of the genetic material of an individual. For achieving this, the Gauss-mutation was used, which adds a Gaussian distributed random value u·mutation_rate (u ∼ N(0, σ)) to each component of the recombined child generation. Random changes with low probability secure the diversity within the population [84]. Due to this probability distribution, many mutations will only make small changes, but larger jumps are also possible [83]. Table 1 shows an overview of the features used for the GA in this work. All of the described procedures can replace the entire parent population with the new individuals, but this may result in the best individual's loss. The results of [85,86] had shown that elitism could significantly speed up the GA's performance, which can also help prevent the loss of good solutions once they are found. To achieve this and reach a better convergence, the worst offspring (20% of the population size) is replaced by the best parent individuals [64]. Furthermore, 20% of the parent individuals are chosen to be elite individuals [86] without any lifespan.

Hill Climbing
The Hill Climbing algorithm searches for the next maximum starting from the initial position (knowledge-based parameter combination) with decreasing step size (0.05; 0.02; 0.01; 0.005; 0.002; 0.001). The steps taken are identical to the Gauss-Newton algorithm (in contrast to the Levenberg-Marquardt algorithm, which uses a dynamic step size estimation) because minimizing the sum of least squares (therefore the RMSE) and maximizing the reciprocal generate the same step selection. Therefore, the step widths are always reduced if no neighboring position has a higher metric value than the current one. All neighbors with the respective step size are examined (initial values + [−step; 0; step]) for all five weighting parameters. Thus, there are 243 parameter combinations to be examined per step (3 5 possible permutations). However, the parameter values cannot exceed or fall below the limit values of 0 and 1. All combinations are stored in a matrix after calculation, which serves as a continuous comparison database to prevent the recalculation of already known parameter combinations.

Nelder-Mead Method
This algorithm (also known as the downhill simplex method) aims to minimize the metric value via function evaluations in the form of a simplex. Therefore, N + 1 function evaluations are calculated for the N parameters to be optimized. In each step, the point with the worst metric value is replaced by a new one [87][88][89]. The dynamic step size determination allows bypassing local minima at the beginning of the algorithm process before the algorithm's convergence occurs. Thus, local minima can be circumvented in some cases. To use this algorithm, the metric value is multiplied by "−1".

Simulated Annealing
The basic idea is to reproduce a cooling process such as that found in metallurgy. Here, the slow cooling ensures that stable crystals can form since the atoms have sufficient time to arrange themselves. Thus, a lower energy state is achieved, which corresponds to a local optimum. In the context of optimization methods, the temperature corresponds to the probability that an intermediate solution may deteriorate during the iterations. Thus, this optimization procedure can leave local optima and better values can be found [90,91].

Particle Swarm
This optimization procedure is similar to the GA but was derived from observing flocks of birds [92]. Individuals move in steps through the search space defined by the target function. Depending on the other individuals' metric values, the orientation of the particle movement, and thus, the values of the parameter values of the next individuals are determined. The particles will gather around the local maxima in the search area with the expectation to hit the global maximum. Problems arise when no particle is near the global maximum. Then the algorithm converges to the local maxima [93].

Experimental Evaluation and Data Collection
The experimental evaluation of the optimization strategy and the data collection were performed in a university building. The used smartphones were a Samsung Galaxy S7 (S7 SM-G930F, hereinafter referenced as S7) with an Exynos 8890 (eight cores) and an LG V30  Figure 7. All measurement runs are performed on the same floor with 107 GTPs, visualized with red dots, while the start is shown in black and the end of the track in blue. Eight test runs were performed with the S7 and five with the LG. In the case of the S7, five test runs were used to estimate the weights and the remaining three for evaluation, while in the case of the LG, the ratio was three to two runs. Each run consists of approximately 583 steps with an average walking speed of around 1.6 m/s. Because the first steps are needed to achieve a convergence of the particle values, the first 30 steps are excluded from the optimization process. Otherwise, the fastest possible convergence will dominate the overall RMSE estimation.
in blue. Eight test runs were performed with the S7 and five with the LG. In the case of the S7, five test runs were used to estimate the weights and the remaining three for evaluation, while in the case of the LG, the ratio was three to two runs. Each run consists of approximately 583 steps with an average walking speed of around 1.6 m/s. Because the first steps are needed to achieve a convergence of the particle values, the first 30 steps are excluded from the optimization process. Otherwise, the fastest possible convergence will dominate the overall RMSE estimation. The smartphones were held in hand in front of the body while walking through the measurement runs recording accelerometer, gyroscope, magnetometer, WLAN, and BLE using a data recording app. Every signal was logged with its timestamp, which was determined by the operation system. The app offers the opportunity to log the reference point when passing by as checkpoints via the current timestamp. The reference points sequence was the same during the measurement runs and each passing was logged with its timestamp. A run consists on average of 583 steps and a measurement duration of approximately 367 s. Each measurement run includes the passing of 157 checkpoints (some GTPs were used twice). The measurement runs were cut to size, so the first and last checkpoints' timestamps represent the start and end of all recorded signals of a specific run. Measurement data outside of this period was disregarded for further analysis.
It has to be noticed that the deviation analysis was estimated from the data of the runs and during movements. That is why [5] estimated up to 0.3 m deviations between an accurate centered foot on the GTP and the recorded steps. The smartphones were held in hand in front of the body while walking through the measurement runs recording accelerometer, gyroscope, magnetometer, WLAN, and BLE using a data recording app. Every signal was logged with its timestamp, which was determined by the operation system. The app offers the opportunity to log the reference point when passing by as checkpoints via the current timestamp. The reference points sequence was the same during the measurement runs and each passing was logged with its timestamp. A run consists on average of 583 steps and a measurement duration of approximately 367 s. Each measurement run includes the passing of 157 checkpoints (some GTPs were used twice). The measurement runs were cut to size, so the first and last checkpoints' timestamps represent the start and end of all recorded signals of a specific run. Measurement data outside of this period was disregarded for further analysis.
It has to be noticed that the deviation analysis was estimated from the data of the runs and during movements. That is why [5] estimated up to 0.3 m deviations between an accurate centered foot on the GTP and the recorded steps.
For better use of the individual signals, accelerometer, gyroscope, and magnetometer sensor data were converted with a spline interpolation to equidistant data points using a frequency of 100 Hz. Then the step detection was implemented according to [94]. For this, the magnitudes of the acceleration signals are determined. A step is only recognized if the acceleration values' variance in a sliding window with a width of 0.8 s exceeds a value of 0.8.
Since the smartphones were held in hand, a higher value for sig_thresh of 0.8 was empirically determined here compared to the original value of 0.6. Subsequently, the acceleration magnitude is smoothed by a moving-average filter with a width of 0.31 s and windowed peak detection with a width of 0.59 s is performed. The sequence of local maximum and minimum accelerations thus determined marks the steps. The step length is then determined using the mean absolute acceleration magnitude of a step, according to [76] (see Equation (8)).
K was determined to be 0.3179 based on an empirical investigation of one part of the data set of Wang [95]: The orientation was estimated using the magnetometer to determine the absolute magnetic north direction, while the accelerometer indicates the direction of the gravity vector to determine the down direction. The relative orientation changes were tracked via the gyroscope [96]. This work uses the step size (σ = 0.15 m) and heading (σ = 40 • ) uncertainties postulated in [5].
During the test runs, it was checked whether beacon signals were recognized for a step. If this is the case, the difference in signal strength between the measured beacon signal and all reference points which have received the recorded Beacon ID in the offline phase will be compared. The inverse squared signal strength differences between the recorded and the reference signals are used to estimate the weights for the position determination. If no match was found, the corresponding reference points for the current position determination are ignored. The position is then determined as the weighted average of all remaining reference positions. Here, in case of signal strength differences of 0, these values are changed to 0.1 to avoid impermissible weights. The median of the position deviations during the runs to the determined positions for registered BLE signals is 6.89 m for the S7 and 4.32 m for the LG.
The same procedure was used for positioning using WLAN fingerprinting. The WLAN fingerprinting provided a higher signal diversity with 109 different WLAN signal sources, whereas BLE uses three. That is why a reference point is only included in the positioning process if it has at least three signal overlaps with the currently registered WLAN signal. Finally, the position is determined from the current step's permissible reference points as a WKNN (k = 5). The median of the position deviations from the actual position to the determined positions for registered WLAN signals is 4.43 m for the S7 and 7.04 m for the LG.
A simple outlier determination was carried out to determine if magnetic field anomalies are present, which indicate characteristic building parts (especially fire doors). For this, a magnetic field anomaly was detected if the magnetic field magnitude has a greater or lower value than three median absolute deviations from the moving median. For this purpose, a window width of 4 s was used. Furthermore, a magnetic field disturbance is present if the magnetic field magnitude has a value that is more than twice the standard deviation from the mean value of the current run.

Parameter Optimization and Performance Analysis
After data collection, the weights for an optimized hybrid localization strategy based on the grid-based approach (Equation (5)) were estimated by means of the GA using the position estimates of the PF and the known position in each GTP as input data.
The results of the GA and for comparison of the different optimization strategies are shown in Table 2. The weights listed one after the other represent the walkway structure, BLE fingerprinting, WLAN fingerprinting, magnetic field anomalies, and the step length estimation. Although the particle swarm algorithm performs best on the LG data, it fails at the more complex data of the S7 and does not reach a good metric value. The Hill Climbing algorithm reaches the next local optimum but fails to find the global optimum in both datasets. The Nelder-Mead algorithm achieves slightly better results than the pure Hill Climbing but does also not find the possible global optimum. Simulated Annealing achieves good metric values, but overall, the GA outperforms the other algorithms. Although it needs the most computing time, the pure performance of estimating the best metric value is the optimization goal. This process of determining the optimal parameters must only be calculated once at the end of the offline phase.
In contrast to the offline phase, the estimated weights for position estimation are used in the online phase. Thus, a fast computation in the online phase is possible since only one convolution, the element-wise multiplications of the predefined sensor maps and the normalization must be calculated for each step. In Figure 8, the optimization process of the GA is shown. A substantial increase in fitness values can be seen at the beginning of the search during the first 20 generations.
Electronics 2021, 10, x FOR PEER REVIEW 15 of 23 After data collection, the weights for an optimized hybrid localization strategy based on the grid-based approach (Equation (5)) were estimated by means of the GA using the position estimates of the PF and the known position in each GTP as input data.
The results of the GA and for comparison of the different optimization strategies are shown in Table 2. The weights listed one after the other represent the walkway structure, BLE fingerprinting, WLAN fingerprinting, magnetic field anomalies, and the step length estimation. Although the particle swarm algorithm performs best on the LG data, it fails at the more complex data of the S7 and does not reach a good metric value.
The Hill Climbing algorithm reaches the next local optimum but fails to find the global optimum in both datasets. The Nelder-Mead algorithm achieves slightly better results than the pure Hill Climbing but does also not find the possible global optimum. Simulated Annealing achieves good metric values, but overall, the GA outperforms the other algorithms. Although it needs the most computing time, the pure performance of estimating the best metric value is the optimization goal. This process of determining the optimal parameters must only be calculated once at the end of the offline phase.
In contrast to the offline phase, the estimated weights for position estimation are used in the online phase. Thus, a fast computation in the online phase is possible since only one convolution, the element-wise multiplications of the predefined sensor maps and the normalization must be calculated for each step. In Figure 8, the optimization process of the GA is shown. A substantial increase in fitness values can be seen at the beginning of the search during the first 20 generations.   After the cutoff of the first 30 Steps, there are 147 checkpoints left. The evaluation runs' error distances are shown in Figure 9 for the S7 and Figure 10 for the LG. The specific benchmarks are shown in Tables 3 and 4, respectively. Considering all evaluation runs, the mean deviation values were reduced by about 35% for the S7 and 25% for the LG. Furthermore, the RMSE were reduced by about 39% and 26%. The RMSE is approximately 1.97 m for the S7 and 2.93 m for the LG on average. As shown in Figure 11, indicated with dashed lines, the third quartile positioning error was reduced by about 54% (3.57 m to 2.32 m) and 50% (5.35 m to 3.57 m). Moreover, the 95-percentile positioning error was reduced by 67% (6.68 m to 4.00 m) and 38% (7.87 m to 5.71 m) presented with dotted lines. Additionally, the maximum position error was reduced by 39% and 20%.
35% for the S7 and 25% for the LG. Furthermore, the RMSE were reduced by about 39% and 26%. The RMSE is approximately 1.97 m for the S7 and 2.93 m for the LG on average. As shown in Figure 11, indicated with dashed lines, the third quartile positioning error was reduced by about 54% (3.57 m to 2.32 m) and 50% (5.35 m to 3.57 m). Moreover, the 95-percentile positioning error was reduced by 67% (6.68 m to 4.00 m) and 38% (7.87 m to 5.71 m) presented with dotted lines. Additionally, the maximum position error was reduced by 39% and 20%.
Based on the determined optimal parameter values from the test runs, the RMSE values were calculated for all evaluation runs and these are shown in Figure 12. It illustrates the decrease of the deviations to the GTPs if the weighting parameters are optimized with the GA.         Based on the determined optimal parameter values from the test runs, the RMSE values were calculated for all evaluation runs and these are shown in Figure 12. It illustrates the decrease of the deviations to the GTPs if the weighting parameters are optimized with the GA.

Conclusions
In this work, an algorithm for smartphone-based pedestrian localization was shown, which combines the classical PDR with fingerprinting and a novel multiple hypothesis grid-based interpretation of the particles of the PF. The grid-based representations of the fingerprinting maps and the estimated step length are weighted with the optimization scheme of the GA to achieve higher position accuracy.

Conclusions
In this work, an algorithm for smartphone-based pedestrian localization was shown, which combines the classical PDR with fingerprinting and a novel multiple hypothesis grid-based interpretation of the particles of the PF. The grid-based representations of the fingerprinting maps and the estimated step length are weighted with the optimization scheme of the GA to achieve higher position accuracy.
The presented algorithm can deal with natural motion behavior, including dynamic changes of standing and moving with varying step length. This algorithm also does not require the artificial cut of curves to have 90-degree changes. It was also accelerated via GPU-based matrix multiplications to improve the vectorized particle computations.
With the presented probability-map-based approach, the sensor fusion algorithm can combine relative sensor-based movements, absolute fingerprinting-based position estimations, and digital building models in real-time to a step-based position estimation without knowledge of any start position. Moreover, it can deal with uncertain information, sensor noise, and building ambiguities. Other information sources such as brightness measurements, barometer signals, magnetic field inclination and declination anomaly maps, ultrasound, and GNSS can easily be included in the calculations as additional layers in a modular way. Thus, absolute position information is used via additional maps, and referential information is included with convolutions.
It is important to emphasize that in this novel particle interpretation the particles are not moving. In contrast, the particles are represented by fixed grid cells whose probabilities change with every new information source and stepwise iteration. Thus, smooth progress in the change of the position probability distribution indicates the current position, and several position theses can alternate without conflicting with each other.
During the reasoning stage, no particles are removed by walls, but the pixel-wise shifting of the wall maps in the respective direction of steps weakens the values of unlikely positions. The offline optimization of the weights with the GA has to be determined only once for a given environment. In contrast to classical, processing-intensive PF utilizing a large number (e.g., 5000) of freely moving particles, the grid-based implementation allows for seamless application on smartphones to enable real-time position estimation. Another advantage of the presented algorithm is the optionality of a starting position.
A limitation of the current implementation is the fingerprinting specific need of data premeasurement for BLE and WLAN signal maps. Hence, the weights have to be re-optimized when the environment changes (e.g., updated WLAN fingerprinting map, additional BLE beacons, etc.). However, this could be automated by triggering the optimization process after each change of the environment or the databases. Additionally, a denser fingerprinting grid and more BLE beacons would improve the results.
In future work, we will investigate the problems regarding the determination of translation for rotations without clearly recognizable steps or while walking backwards. Furthermore, the algorithm shows promising results and will thus be tested with other datasets, environments, different smartphones, and more subjects. It is therefore intended to use the datasets of the EvAAL Competitions of the last IPINs [97].
In summary, the implemented GA is beneficial for the position accuracy and the localization errors were reduced. In fact, the GA improves the overall mean position accuracy by 31%, as well as the maximum position error compared to the knowledge-based parameter combination using the same sensor data, digital building model, and particle representation.
The overall RMSE was reduced by 34%, whereas the third quartile positioning errors by 54% and 50%. Also, the 95-percentile positioning errors were reduced by 67% and 38%. This improvement was achieved by adjusting the weights of the individual measurement methods to an optimal constellation. Supplementary Materials: The following are available online at https://www.mdpi.com/2079-929 2/10/5/618/s1, Video S1 Extended Localization.
Author Contributions: J.G. performed the experiments, recorded and analyzed the data, conceived and designed the algorithms, prepared the original draft, and visualized the results; J.B. supervised this work, gave many ideas and valuable comments, and reviewed and edited the draft. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: All data included in this study are available upon request by contacting with the corresponding author.