Optimization Techniques in the Localization Problem: A Survey on Recent Advances

Abstract

Optimization is a mathematical discipline or tool suitable for minimizing or maximizing a function. It has been largely used in every scientific field to solve problems where it is necessary to find a local or global optimum. In the engineering field of localization, optimization has been adopted too, and in the literature, there are several proposals and applications that have been presented. In the first part of this article, the optimization problem is presented by considering the subject from a purely theoretical point of view and both single objective (SO) optimization and multi-objective (MO) optimization problems are defined. Additionally, it is reported how local and global optimization problems can be tackled differently, and the main characteristics of the related algorithms are outlined. In the second part of the article, extensive research about local and global localization algorithms is reported and some optimization methods for local and global optimum algorithms, such as the Gauss–Newton method, Genetic Algorithm (GA), Particle Swarm Optimization (PSO), Differential Evolution (DE), and so on, are presented; for each of them, the main concept on which the algorithm is based, the mathematical model, and an example of the application proposed in the literature for localization purposes are reported. Among all investigated methods, the metaheuristic algorithms, which do not exploit gradient information, are the most suitable to solve localization problems due to their flexibility and capability in solving non-convex and non-linear optimization functions.

1. Introduction

Optimization is a mathematical discipline which provides tools to select the best element from a set of possible alternatives based on some given criteria. In general, an optimization problem is one of maximizing or minimizing a real function by methodically selecting input values from a permitted set and figuring out the output values of the function itself. Optimization problems can be classified as continuous or discrete, and to solve them, several optimization techniques and approaches have been developed. These techniques can be divided into exact methods, approximate algorithms, metaheuristics, and greedy algorithms [1]. Each of them has pros and cons, and the user must choose the proper and most appropriate technique to solve the specific problem.

Nowadays, optimization, sometimes also referred to as mathematical programming, is used in every scientific field where it is necessary to deal with mathematical problems characterized by high complexity, constraints, interdependencies among variables, and a large space of solution [2]. Finance [3], chemical engineering [4], mechanical engineering [5], electrical engineering [6], and system and database design [7] are only a few examples where optimization has been applied in the last decades. In the field of localization, optimization methods have been widely used to solve localization problems pertaining to robots [8], smartphones [9], permanent magnets [10], vehicles [11], underground structures [12], urban utilities [13], and so on. By referring to this field, this paper aims to investigate and collect the main used optimization methods in the localization literature, such as the Gradient Descent method, Levenberg–Marquardt method, Genetic Algorithm (GA), Particle Swarm Optimization (PSO), etc., by reporting the concept of the used technique and an example of its application in the considered field. Furthermore, the paper also aims to provide in the initial part a general introduction to the classic optimization problem, providing the mathematical formulation, and, in the second part, a clear mathematical formulation of each investigated technique.

The rest of the paper is organized as follows. In Section 2, firstly, the single objective (SO) optimization problem is mathematically defined; then, the local vs global cases are presented, and at the end of the section, the multi-objective optimization (MO) problem with the related approaches is reported. In Section 3, some of the most interesting and important applications of the optimization methods in the field of localization are summarized and classified; for each of them, a brief explanation of the method and the inspired idea are reported. In Section 4, the conclusions are reported.

2. A Brief Overview of the Optimization Problem

2.1. Definition

Optimization methods are used to find a set of inputs to an objective function that results in the evaluation of the minimum or maximum of the function. Mathematically, the optimization problem is defined as follows [14]:

$$\begin{array}{lll}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}:& {\mathrm{f}}_0\left(\mathrm{x}\right)& \\ \mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h} \mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}:& {\mathrm{g}}_{\mathrm{j}}\left(\mathrm{x}\right)\le 0& \mathrm{j}=1,\dots ,\mathrm{m}\\ & {\mathrm{h}}_{\mathrm{k}}\left(\mathrm{x}\right)=0& \mathrm{k}=1,\dots ,\mathrm{p}\end{array}$$

(1)

where ${\mathrm{f}}_0\left(\mathrm{x}\right)$ is a function, called an objective function, and defined as ${\mathrm{f}}_0:{\mathbb{R}}^{\mathrm{n}}\to \mathbb{R}$, $\mathrm{x}={\left({\mathrm{x}}_1;{\mathrm{x}}_2;\dots ;{\mathrm{x}}_{\mathrm{n}}\right)}^{\mathrm{T}}$ is the optimization vector variable (or the design variables) of the problem, and ${\mathrm{g}}_{\mathrm{j}}\left(\mathrm{x}\right)$ and ${\mathrm{h}}_{\mathrm{k}}\left(\mathrm{x}\right)$ are the inequality and equality functions, respectively, that are constraint functions.

If for any $\mathrm{z}$ with ${\mathrm{g}}_1\left(\mathrm{z}\right)\le 0,\dots ,{\mathrm{g}}_{\mathrm{m}}\left(\mathrm{z}\right)\le 0$ and ${\mathrm{h}}_1\left(\mathrm{z}\right)=0,\dots ,{\mathrm{h}}_{\mathrm{k}}\left(\mathrm{z}\right)=0$, the relation ${\mathrm{f}}_0\left({\mathrm{x}}^{\mathrm{*}}\right)\le {\mathrm{f}}_0(\mathrm{z})$ is verified, then the vector ${\mathrm{x}}^{\mathrm{*}}$ is the optimal solution of the problem and it has the smallest objective value among all vectors that satisfy the constraint functions. The maximum problem of a function $\mathrm{f}(\mathrm{x})$ can be solved considering the minimum problem of −$\mathrm{f}(\mathrm{x})$.

The objective function and the constraint functions can be linear or non-linear and an optimization problem can be called a linear program if the objective $({\mathrm{f}}_0,\dots ,{\mathrm{f}}_{\mathrm{m}})$ and constraint functions satisfy the linear condition that can be expressed, for example, as follows:

$${\mathrm{f}}_{\mathrm{s}}\left(\mathsf{\alpha }\mathrm{x}+\mathsf{\beta }\mathrm{y}\right)={\mathsf{\alpha }\mathrm{f}}_{\mathrm{s}}\left(\mathrm{x})+\mathsf{\beta }{\mathrm{f}}_{\mathrm{s}}(\mathrm{y}\right)$$

(2)

where the variable s $x,y$ $\in {\mathbb{R}}^n$, and the coefficients $\alpha ,\beta \in \mathbb{R}$. If the equality in (2) is not verified, the problem is a non-linear program [14]. The problems can also be classified as convex optimization problems if the following condition is satisfied:

$${\mathrm{f}}_{\mathrm{s}}\left(\mathsf{\alpha }\mathrm{x}+\mathsf{\beta }\mathrm{y}\right)\le {\mathsf{\alpha }\mathrm{f}}_{\mathrm{s}}\left(\mathrm{x})+\mathsf{\beta }{\mathrm{f}}_{\mathrm{s}}(\mathrm{y}\right)$$

(3)

where $\mathsf{\alpha }+\mathsf{\beta }=1$, $\mathsf{\alpha }\ge 0$, and $\mathsf{\beta }\ge 0$. The convex optimization can be considered a generalization of a linear programming because any linear program is also a convex optimization problem [14].

The solution to the problem is found inside the searchable design space that is subjected to the side constraints defined as ${\mathrm{x}}_{\mathrm{i}\mathrm{L}}\le {\mathrm{x}}_{\mathrm{i}}\le {\mathrm{x}}_{\mathrm{i}\mathrm{U}}$ with $\mathrm{i}=1,\dots ,\mathrm{n}$, where variables are limited between the upper and lower bounds [15].

Optimization problems can have some or all of the design variables restricted to integer or discrete values, and they are referred to as integer or discrete optimization problems, respectively.

To find a solution for a class of optimization problems in relation to a given accuracy, a specific algorithm is needed; most of them are able to deal with the side constraints and the constraints function separately because the side constraints are managed directly by the algorithms themselves; unconstrained problems exist, but they can still have side constraints [15].

One or more equality and/or inequality constraints, with or without side restrictions, characterize a constrained optimization problem. For a generic problem, an equality constraint can only either be violated or satisfied, while an inequality constraint can be violated, active, or satisfied. An active inequality constraint, defined as ${\mathrm{g}}_{\mathrm{j}}\left(\mathrm{x}\right)\le 0$, is said to be active (or tight) when at a design point ${\mathrm{x}}^{\mathrm{*}}$, it is satisfied as an equality, i.e., ${\mathrm{g}}_{\mathrm{j}}\left(\mathrm{x}\right)=0$ [15].

Optimization techniques, or algorithms, use a specific procedure that is able to find the solution among design variable values that results in the best objective function value, and at the same time, it satisfies all the constraint functions (equality, inequality) and the side constraints; more than one optimum can be found, and they are referred to as local or relative optima.

2.2. A Comparison between Local and Global Extrema

The extreme minimum or extreme maximum of the objective function over the whole input search space is known as a global optimum. In general, finding a local optimum may be rather simple, whereas finding the global optimum may be more challenging, in particular for non-linear functions [16]. As an example, Figure 1 shows a generic function $\mathrm{f}\left(\mathrm{x}\right)$ which has three minima: two of them are of the local type, and one is the global minimum of the problem.

It follows that optimization problems and algorithms can be classified into two types: local optimization and global localization algorithms. The former type of algorithm looks for local optimum solutions in a particular area of the search space, whereas the latter looks for the best solution for problems with multiple local optimums. Furthermore, there could be multiple local optimal points for an objective function, and if there is only one, that local optimum point is referred also to as the global optimum point. In case a global optimum does not exist for an objective function, the problem cannot be classified as an optimization problem.

In general, optimization problems are typically classified as multimodal or unimodal in terms of their modality. When a problem has more local optima solutions in addition to the global optimum, it is said to be multimodal; when it has only the global optimum as its local optima solution, it is said to be unimodal. Multimodal problems can sometimes be trickier and therefore more complex to solve [17].

The authors in reference [18] classify optimization problems as reported in Figure 2. It is highlighted that non-convex functions do not need to be multimodal because a non-convex function can have only one local optimum, and it follows that this is also the global optimum.

Local and global optimization algorithms deal with different problems; a local optimization technique is utilized when the objective function includes a single optimum or when the region of the global optimum is known; when the objective function response surface’s structure is unknown or the function has local optimum points, a global optimization procedure is employed [19].

In order to discover the extrema of the function in a certain region of the search space or to come as close to them as possible, a local optimization method is developed to cross the considered region looking for the minima (or maxima) [14].

Global optimization techniques are suitable for situations where there are few variables and no necessity for fast computation time. They work with a single candidate solution or a population of them, from which new candidate solutions are generated iteratively and assessed using a predetermined function to see whether the newly generated variables are getting better or worse [14].

2.3. Local Optimization Algorithms

Most local optimization algorithms exploit the gradient information of the function to find the optimum solution. They work well for solving problems with lots of design variables, and fine-tuning their settings is rather simple. Furthermore, they are considered efficient compared to the number of function evaluations (or iterations) needed to find the solution. Among the drawbacks are that they can only locate a local optimum, they have issues solving discrete optimization problems, they require algorithms that are difficult to implement efficiently, and they suffer from numerical noise [15].

There are many types of gradient-based algorithms, and each of them is characterized by a different logic or idea to determine the search direction. In general, gradient-based algorithms typically use a two-step iterative process to reach the optimum: the first one is to exploit the gradient information to find a search direction $\mathrm{S}$ in which to move, and the second step is moving in the direction $\mathrm{d}$; the last one is called one-dimensional or a line search and also provides the optimum step size, ${\mathsf{\alpha }}_{\mathrm{k}}$. The two steps are repeated until no more progress is obtained and mathematically, they can be expressed as follows:

$${\mathrm{x}}_{\mathrm{k}}={\mathrm{x}}_{\mathrm{k}-1}+{\mathsf{\alpha }}_{\mathrm{k}}\mathrm{d}$$

(4)

where $\mathrm{k}$ identifies the number of the $\mathrm{k}$-th-step [20]. The gradient information is not always easily accessible in complex problems, and techniques such as Finite Difference Gradient Calculations, Automatic Differentiation, or Analytic/Semi-Analytic Gradient Calculations can be used to compute the gradient. Moreover, the first one provides gradient information that is accurate at the working precision; the second one provides only an approximation of the gradient, giving an accuracy that is related to the selected step size; the last one can be applied to linear finite element codes, and it is inexpensive [15].

2.4. Global Optimization Algorithms

Many problems can have multiple local optima (like the problem represented in Figure 1), and local globalization algorithms are able to find one of them without guaranteeing or saying if it is a local or global minimum point because all minima satisfy the Karush–Kuhn–Tucker (KKT) condition, which is a mathematical tool to verify if a solution is a minimum or not, but without specifying if it is a global one.

The found minima depend on the starting point, giving the first met minima as a result. To find other minima in the design space, a multi-start approach can be adopted, which starts from different points [21].

Global optimization algorithms have a far higher possibility of locating the global or near-global optimum. It is preferable to speak of these algorithms as having global properties because, in general, no algorithm can ensure convergence on a global optimum.

They fall into one of two categories: metaheuristic algorithms and deterministic algorithms [22].

2.4.1. Metaheuristic Algorithms

Metaheuristic (or stochastic) algorithms can be defined as computational intelligence paradigms which are especially used for solving complex optimization problems [23]. When precise optimization methods are unable to provide outcomes, they can provide satisfactory results. Metaheuristic algorithms can be classified into evolution-based, swarm intelligence-based, physics-based, and human behavior-related algorithms [24].

This class of optimization techniques is suitable for complex non-linear and discontinuous problems where classical optimization techniques might fail [25].

These algorithms, in contrast to the local methods, use a certain number of design points, which are often referred to as a population of individuals; they allow for determining the optimum without adopting any gradient information. These techniques frequently draw inspiration from various natural events, and some of their advantages include a higher likelihood of locating a global or nearly global optimum, robustness, and simplicity in implementation; moreover, they are also suitable for solving discrete optimization problems. On the other hand, some disadvantages include high computing costs, the limited size of the problem, weak constraint-handling capabilities, and difficulties in tuning problem-specific parameters [15].

These days, PSO and the GA are the two most often used metaheuristic algorithms. Additionally, evolutionary algorithms include genetic programming, Ant Colony Optimization (ACO), Simulated Annealing (SA), Differential Evolution (DE), Evolutionary Programming, Harmony Search, and others [22].

2.4.2. Deterministic Algorithms

Deterministic algorithms aim to find the global solution of an optimization problem by guaranteeing that the found solution is the real global minimum (or maximum) of the problem with a predefined tolerance. These algorithms usually deal with specific classes of problems, such as linear programming (LP), mixed-integer linear programming (MILP), non-linear programming (NLP), and mixed-integer non-linear programming (MINLP). An example of a deterministic technique is the DIRECT algorithm [26], which can find the global minimum of a multivariate function with simple bounds; it is a modification of the conventional Lipschitzian method where the need to specify a Lipschitz constant is eliminated.

In general, deterministic approaches can obtain advantages by exploiting the analytical properties of a problem to find the global solution; however, over the years, heuristic methods have shown to be more efficient and flexible than deterministic ones [18].

2.5. From Single to Multi-Objective Optimization

The case described by Equation (1) is referred to as a single objective optimization problem because the function to be optimized is only one. However, in many applications, the optimization of more than one objective function at the same time is requested, and this case, it is referred to as a multi-objective (MO) optimization problem. It is mathematically described as follows:

$$\begin{array}{lll}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}:& {\mathrm{f}}_{\mathrm{m}}\left(\mathrm{x}\right)& \mathrm{m}=1,\dots ,\mathrm{M}\\ \mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h} \mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}:& {\mathrm{g}}_{\mathrm{j}}\left(\mathrm{x}\right)\le 0& \mathrm{j}=1,\dots ,\mathrm{J}\\ & {\mathrm{h}}_{\mathrm{k}}\left(\mathrm{x}\right)=0& \mathrm{k}=1,\dots ,\mathrm{K}\end{array}$$

(5)

where $\mathrm{m}$ is the $\mathrm{m}$-th function among the $\mathrm{M}$ functions to be optimized.

In MO optimization problems, two spaces are defined: the first is the decision variable space of the solution vector that contains all possible solutions $\mathrm{x}$ of the problem; the second is the multi-dimensional space of the objective function vector that contains the evaluations of each objective function; each $\mathrm{x}$ solution belonging to the decision variable space corresponds to a point in the objective function space. The spaces for a problem with a vector solution $\mathrm{x}$ with three components and two membership functions are represented in Figure 3 [27].

If every objective function and solution region in a MO problem are convex, the problem is said to be convex too, and the convexity is crucial to solving the optimization problem.

The methods to find the solution to MO problems can be classified into two different types: the Pareto method and the scalarization method [28]. The former is used when performance indicators and desired solutions are distinct, and in such a way, the Pareto approach is applied to generate a compromise solution, or trade-off, that can be represented as a Pareto optimal front (POF). Unlike in the latter, a performance indicator component method is exploited to create a scalar function that is integrated into the fitness function [29].

The process of finding the optimal solution can be assisted by using metaheuristic algorithms such as GA, PSO, ACO, and so on.

2.5.1. Pareto Method

From a mathematical point of view, a MO problem with $\mathrm{n}$ objective functions solved by using the Pareto method is expressed as follows:

$$\begin{array}{l}{\mathrm{f}}_{1,\mathrm{o}\mathrm{p}\mathrm{t}}=\min {\mathrm{f}}_1(\mathrm{x})\\ \begin{array}{l}{\mathrm{f}}_{2,\mathrm{o}\mathrm{p}\mathrm{t}}=\min {\mathrm{f}}_2(\mathrm{x})\\ \dots \\ \begin{array}{l}\dots \\ {\mathrm{f}}_{\mathrm{n},\mathrm{o}\mathrm{p}\mathrm{t}}=\min {\mathrm{f}}_{\mathrm{n}}(\mathrm{x}).\end{array}\end{array}\end{array}$$

(6)

The method keeps the elements of the solution vectors independent and exploits the concept of dominance with which a solution can be considered a Pareto optimal solution (POS) or a non-dominated solution if around the considered solution there is no way of improving any objective without degrading at least one other objective; on the contrary, a solution is said to be a Pareto-dominated solution (PDS) if it can be improved without worsening at least one other objective. As an example, Figure 4 shows the solutions in the space of the objective functions where POSs are represented by red dots, whilst the PDSs are represented by red squares. Furthermore, in agreement with the representation in Figure 4, it is possible to give the following definitions [29]:

• Anchor point: the best value provided by an objective function.
• Utopia point: a point which is obtained by intersecting the best values provided by the objective functions.
• Pareto optimal front (POF): the line where the POSs lay.

Once the POF has been obtained, the MO problem is solved by computing which of the solutions on the POF has the shortest Euclidean distance. The following equation is used:

$${\mathrm{d}}_{\mathrm{E}}=\mathrm{m}\mathrm{i}\mathrm{n}\sqrt{{\left({\displaystyle \frac{{\mathrm{Q}}_1-{\mathrm{Q}}_1^{\mathrm{*}}}{{\mathrm{Q}}_{1\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}}}\right)}^2+{\left({\displaystyle \frac{{\mathrm{Q}}_2-{\mathrm{Q}}_2^{\mathrm{*}}}{{\mathrm{Q}}_{2\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}}}\right)}^2}$$

(7)

where $({\mathrm{Q}}_1^{\mathrm{*}},{\mathrm{Q}}_2^{\mathrm{*}})$ are the coordinates for the utopia point related to ${\mathrm{f}}_1$ and ${\mathrm{f}}_2$, respectively, $\left({\mathrm{Q}}_1,{\mathrm{Q}}_2\right)$ are the coordinates of the solution point on the POF, and ${\mathrm{Q}}_{1\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}$ and ${\mathrm{Q}}_{2\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}$ are normalization factors based on the minimum values of each function.

In general, by considering a MO problem with two objective functions, four types of problems can be obtained and are given as follows:

$$\mathrm{Problem} 1: \left\{\begin{array}{l}{\mathrm{f}}_{1,\mathrm{o}\mathrm{p}\mathrm{t}}=\min {\mathrm{f}}_1(\mathrm{x})\\ {\mathrm{f}}_{2,\mathrm{o}\mathrm{p}\mathrm{t}}=\min {\mathrm{f}}_2(\mathrm{x})\end{array}\right.,$$

(8)

$$\mathrm{Problem} 2: \left\{\begin{array}{l}{\mathrm{f}}_{1,\mathrm{o}\mathrm{p}\mathrm{t}}=\min {\mathrm{f}}_1\left(\mathrm{x}\right)\\ {\mathrm{f}}_{2,\mathrm{o}\mathrm{p}\mathrm{t}}=\max {\mathrm{f}}_2\left(\mathrm{x}\right)\end{array}\right.,$$

(9)

$$\mathrm{Problem} 3: \left\{\begin{array}{l}{\mathrm{f}}_{1,\mathrm{o}\mathrm{p}\mathrm{t}}=\max {\mathrm{f}}_1\left(\mathrm{x}\right)\\ {\mathrm{f}}_{2,\mathrm{o}\mathrm{p}\mathrm{t}}=\min {\mathrm{f}}_2\left(\mathrm{x}\right)\end{array}\right.,$$

(10)

$$\mathrm{Problem} 4: \left\{\begin{array}{l}{\mathrm{f}}_{1,\mathrm{o}\mathrm{p}\mathrm{t}}=\max {\mathrm{f}}_1(\mathrm{x})\\ {\mathrm{f}}_{2,\mathrm{o}\mathrm{p}\mathrm{t}}=\max {\mathrm{f}}_2(\mathrm{x})\end{array}\right..$$

(11)

Each of the four problems gives rise to a different shape of the POF and they are shown in Figure 5.

The objective function space for a problem with three functions can be represented in a three-dimensional space, but for a problem with a number of objective functions greater than three, the space cannot be represented.

To find non-dominated solutions, specific algorithms are used such as the continuously updated method reported in [29].

2.5.2. Scalarization Method

The scalarization method includes each function of the MO problem into a single scalar fitness function, practically transforming the MO problem into a single objective problem. Different types of scalarization can be defined and some basic approaches are the Weighted Global Criterion Method (WGCM), Weighted Sum Method (WSM), Exponential Weighted Criterion (EWC), or the Weighted Product Method (WPM) [30].

From a conceptual point of view, the WSM is the easiest way to implement the scalarization method. Each function of the MO problem is summed by adopting weights; it gives rise to a new total objective function $\mathrm{F}\left(\mathrm{x}\right)$ which is given as follows:

$$\mathrm{F}\left(\mathrm{x}\right)={\mathrm{w}}_1{\mathrm{f}}_1\left(\mathrm{x}\right)+{\mathrm{w}}_2{\mathrm{f}}_2\left(\mathrm{x}\right)+\dots +{\mathrm{w}}_{\mathrm{i}}{\mathrm{f}}_{\mathrm{i}}(\mathrm{x})+\dots +{\mathrm{w}}_{\mathrm{n}}{\mathrm{f}}_{\mathrm{n}}(\mathrm{x})$$

(12)

where $\mathrm{i}$ is the $\mathrm{i}$-th function of the original MO problem, ${\mathrm{w}}_{\mathrm{i}}$ represents the weight of the $\mathrm{i}$-th function, and ${\mathrm{f}}_{\mathrm{i}}$ is the $\mathrm{i}$-th function of the original MO problem composed of $\mathrm{n}$ objective functions [29]. The values of the weights are chosen in agreement with the performance priority: a large weight is given to functions with higher priority and vice versa for functions with lower priority. The following are three examples of how weights can be set:

• Equal weights: ${\mathrm{w}}_{\mathrm{i}}=1/\mathrm{n}$.
• Rank Order Centroid (ROC) weights: ${\mathrm{w}}_{\mathrm{i}}={\displaystyle \frac{1}{\mathrm{n}}}\sum _{\mathrm{k}=\mathrm{i}}^{\mathrm{n}}{\displaystyle \frac{1}{\mathrm{k}}}$.
• Rank Sum (RS) weights: ${\mathrm{w}}_{\mathrm{i}}={\displaystyle \frac{2(\mathrm{n}+1-\mathrm{i})}{\mathrm{n}(\mathrm{n}+1)}}$.

To provide a balance among functions, normalization is needed, and it can be obtained by dividing each objective function by its own RMS value ${\mathrm{f}}_{\mathrm{i}}^{\mathrm{R}\mathrm{M}\mathrm{S}}$ [31]. The final expression is given as follows:

$$\mathrm{F}\left(\mathrm{x}\right)=+{\displaystyle \frac{{\mathrm{w}}_1{\mathrm{f}}_1\left(\mathrm{x}\right)}{{\mathrm{f}}_1^{\mathrm{R}\mathrm{M}\mathrm{S}}}}+{\displaystyle \frac{{\mathrm{w}}_2{\mathrm{f}}_2\left(\mathrm{x}\right)}{{\mathrm{f}}_2^{\mathrm{R}\mathrm{M}\mathrm{S}}}}+\dots +{\displaystyle \frac{{\mathrm{w}}_{\mathrm{n}}{\mathrm{f}}_{\mathrm{n}}\left(\mathrm{x}\right)}{{\mathrm{f}}_{\mathrm{n}}^{\mathrm{R}\mathrm{M}\mathrm{S}}}}$$

(13)

In Equation (13), the sign of each objective function must agree with the request of minimization (minus) or maximization (plus).

The WSM can be considered as a particular case of the WGCM approaches, in which the scalarization can be obtained by adopting one of the following definitions [30]:

$$\mathrm{F}\left(\mathrm{x}\right)=\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\omega }}_{\mathrm{i}}{\left[{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)\right]}^{\mathrm{p}}$$

(14)

$$\mathrm{F}\left(\mathrm{x}\right)=\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left[{\mathsf{\omega }}_{\mathrm{i}}{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)\right]}^{\mathrm{p}}$$

(15)

where $\mathrm{p}$ is a parameter which can enhance the minimization of the function, and weights must satisfy $\sum {\mathsf{\omega }}_{\mathrm{i}}=1$.

The EWC method is defined as follows [32]:

$$\mathrm{F}\left(\mathrm{x}\right)=\sum _{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{e}}^{\mathrm{p}{\mathsf{\omega }}_{\mathrm{i}}}-1\right){\mathrm{e}}^{\mathrm{p}{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)},$$

(16)

whilst the WPM is given as follows [33]:

$$\mathrm{F}\left(\mathrm{x}\right)=\prod _{\mathrm{i}}^{\mathrm{n}}{\left[{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)\right]}^{{\mathsf{\omega }}_{\mathrm{i}}}.$$

(17)

2.5.3. Summary of the Optimization Problems

In summary, from the point of view of the number of objectives, optimization problems can be divided into two types: SO and MO problems. The main difference in the formulation of the problem is that in the former, only one objective is considered, whereas in the latter, there are at least two objectives. The MO case can be considered an extension of the SO case, and allows the optimization of more functions at the same time that belong to the same problem and have the same design variables.

Finally, the MO problem can be solved by using either the Pareto method or one of the scalarization approaches. The first one deals with the optimization problem by solving each fitness function separately, whereas the second ones provide a function that is a single figure of merit that englobes all the objectives of the problem in one function. The main advantage of the Pareto approach is that each fitness function is solved independently without any regard for the other objectives, and the final solution is extracted by considering the solution that has the shortest distance from the utopia point. On the other hand, scalarization transforms the MO problem into the SO type, but has the drawback of setting the weights of each objective properly in order to correctly set their priority.

3. Optimization Algorithms Used in the Field of Localization

In this section, all optimization methods in the field of localization are reported. Generally, a description of the algorithm, its concept, and the mathematical formulation are reported; then, a real and meaningful application for localization purposes is described.

3.1. Newton’s Method, Gradient Descent, and Gauss–Newton Method

Newton’s Method and the Gradient Descent (GD) algorithm are the first methods which were studied in the field of optimization to deal with local problems, and they have a similar structure but with the difference that the GD algorithm is a first-order method and Newton’s method is a second-order type. The latter has the disadvantage of a higher computational cost because it requires the calculation of the Hessian matrix. At the end of the paragraph, an extension of Newton’s method called the Gauss–Newton method is presented.

3.1.1. Newton’s Method

Newton’s algorithm is a classical method and is an unconstrained algorithm that is derived from a second-order Taylor series expansion of the objective function starting from an initial design point ${\mathrm{x}}_0$. The function can be expressed as follows [15]:

$$\mathrm{f}\left(\mathrm{x}\right)\approx \mathrm{f}\left({\mathrm{x}}_0\right)+\nabla \mathrm{f}{\left({\mathrm{x}}_0\right)}^{\mathrm{T}}\left(\mathrm{x}-{\mathrm{x}}_0\right)+{\displaystyle \frac{1}{2}}{\left(\mathrm{x}-{\mathrm{x}}_0\right)}^{\mathrm{T}}\mathrm{H}({\mathrm{x}}_0)(\mathrm{x}-{\mathrm{x}}_0)$$

(18)

where $\mathrm{x}$ is the variable of the problem, and $\mathrm{H}\left({\mathrm{x}}_0\right)$ is the Hessian matrix that contains the second-order gradient information of the objective function. By differentiating Equation (18) with respect to $\mathrm{x}$ and setting the result equal to zero according to the KKT conditions, it gives rise to the following iterative formula for the current design point [34]:

$${\mathrm{x}}_{\mathrm{n}+1}={\mathrm{x}}_{\mathrm{n}}-{\displaystyle \frac{{\mathrm{f}}^{\prime }({\mathrm{x}}_{\mathrm{n}})}{{\mathrm{f}}^{\prime \prime }({\mathrm{x}}_{\mathrm{n}})}}$$

(19)

where ${\mathrm{f}}^{\prime }$ and ${\mathrm{f}}^{\prime \prime }$ are the first and the second derivative, respectively, which can be expressed as ${\mathrm{f}}^{\prime }\left({\mathrm{x}}_{\mathrm{n}}\right)=\nabla \mathrm{f}\left({\mathrm{x}}_{\mathrm{n}}\right)$ and ${\mathrm{f}}^{\prime \prime }\left({\mathrm{x}}_{\mathrm{n}}\right)=\mathrm{H}{\left({\mathrm{x}}_{\mathrm{n}}\right)}^{-1}$. The step size cannot be changed, and no one-dimensional search is required because it is provided by $-\mathrm{H}{\left({\mathrm{x}}_{\mathrm{n}}\right)}^{-1}\nabla \mathrm{f}\left({\mathrm{x}}_{\mathrm{n}}\right)$.

The method has a quadratic rate of convergence, requiring only a single step to obtain the optimum; however, the estimation of the second-order derivative requires a high computational cost and often the method can result in being unpractical. To overcome this drawback, some variation of gradient-based methods can be adopted, which makes use only of the first-order gradient information.

3.1.2. Gradient Descent Algorithm

The GD algorithm is a first-order optimization algorithm that finds a local minimum (or maximum) by using an iterative process. The function to be optimized must be differentiable and convex. However, it is also possible to apply the algorithm to quasi-convex functions, but with the risk of being stuck at saddle points.

The GD algorithm is based on this iterative equation [35]:

$${\mathrm{x}}_{\mathrm{k}+1}={\mathrm{x}}_{\mathrm{k}}-\mathsf{\alpha }\nabla \mathrm{f}({\mathrm{x}}_{\mathrm{k}})$$

(20)

where $\mathrm{f}$ is the function to be optimized, ${\mathrm{x}}_{\mathrm{k}}$ is the input variable at the $\mathrm{k}$-th step, and $\mathsf{\alpha }$ is a parameter which controls the step size. It is highlighted that the value of $\mathsf{\alpha }$ controls the behavior of the method because, for example, the smaller the learning rate, the longer the GD convergence time, and a too-small $\mathsf{\alpha }$ could make the algorithm reach the maximum number of iterations too early without finding the right solution; on the other hand, a too-large $\mathsf{\alpha }$ could make the process diverge or just jump around to the minimum (or maximum) of the problem. The value $\mathsf{\alpha }\nabla \mathrm{f}({\mathrm{x}}_{\mathrm{k}})$ in (20) is subtracted when the minimum is looked for, whereas it is summed when the maximum is the desired output.

The GD algorithm deals with unconstrained problems, and it can be extended to also handle constraints by adopting a specific approach called the penalty function method, which is able to transform an unconstrained problem into a constrained one [20]. The method adds a penalty function to the original objective function that is composed of a penalty parameter multiplied by a measure of violation of the constraints. The latter is zero when there is no violation and greater than zero when violations are obtained. Considering a generic problem $\mathrm{m}\mathrm{i}\mathrm{n}=\mathrm{f}(\mathrm{x})$ subjected to ${\mathrm{c}}_{\mathrm{i}}(\mathrm{x})\le 0$, it is transformed into a series of unconstrained minimization problems, given as follows [36]:

$$\min {\mathsf{\Phi }}_{\mathrm{k}}\left(\mathrm{x}\right)=\mathrm{f}\left(\mathrm{x}\right)+{\mathsf{\sigma }}_{\mathrm{k}}\sum _{\mathrm{i}}\mathrm{g}({\mathrm{c}}_{\mathrm{i}}\left(\mathrm{x}\right))$$

(21)

where $\mathrm{g}\left({\mathrm{c}}_{\mathrm{i}}\left(\mathrm{x}\right)\right)={\max \left(0,{\mathrm{c}}_{\mathrm{i}}\left(\mathrm{x}\right)\right)}^2$ is the exterior penalty function, ${\mathsf{\sigma }}_{\mathrm{k}}$ is the penalty coefficient, and $\mathrm{k}$ represents the iteration. In each step of the algorithm, the penalty coefficient is increased, the unconstrained problem is solved, and the provided result is used in the next iteration as the initial guess. In this manner, the solution of the original constrained problem is obtained by the convergence of the unconstrained one.

In general, this method has been used in many different scientific fields, and, in particular, one of them where it has been widely adopted is Machine Learning (ML), in which the GD algorithm has three variations: the Batch GD, Stochastic GD, and Mini-Batch GD algorithms. These three different approaches manage the training process of a neural network in different ways: in the first one, the algorithm calculates the gradient to correct the weights by considering the complete training dataset in each iteration, in the second one, the gradient is calculated by considering only one sample of the dataset in each iteration, and in the last case, which is an intermediate approach between the two aforementioned approaches, the dataset is split into small subsets (batches) and the gradients are computed for each sub-batch in each iteration [37].

3.1.3. Gauss–Newton Method

The Gauss–Newton method is adopted to solve convex non-linear least-squares problems, i.e., a class of problems where it requested to minimize a sum of squared function values. Given an objective function which is twice differentiable, it has a form given as follows:

$$\mathrm{f}\left(\mathrm{x}\right)={\displaystyle \frac{1}{2}}\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{f}}_{\mathrm{i}}{\left(\mathrm{x}\right)}^2$$

(22)

and its gradient is defined as follows:

$$\nabla \mathrm{f}\left(\mathrm{x}\right)=\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)\nabla {\mathrm{f}}_{\mathrm{i}}(\mathrm{x}).$$

(23)

The search direction of the Gauss–Newton method is as follows [14]:

$$\Delta {\mathrm{x}}_{\mathrm{g}\mathrm{n}}=-{\left(\sum _{\mathrm{i}=1}^{\mathrm{m}}\nabla {\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)\nabla {\mathrm{f}}_{\mathrm{i}}{\left(\mathrm{x}\right)}^{\mathrm{T}}\right)}^{-1}\left(\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)\nabla {\mathrm{f}}_{\mathrm{i}}(\mathrm{x})\right).$$

(24)

where it is assumed that the inverse exists. The difference with the classic Newton’s method is that in this case, the direction is obtained by neglecting the second term from the Hessian. This result can be obtained by considering the first order of Taylor’s series ${\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}+\mathrm{v}\right)={\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)+\nabla {\mathrm{f}}_{\mathrm{i}}{\left(\mathrm{x}\right)}^{\mathrm{T}}\mathrm{v}$, which gives rise to the following:

$$\mathrm{f}\left(\mathrm{x}+\mathrm{v}\right)={\displaystyle \frac{1}{2}}\sum _{\mathrm{i}=1}^{\mathrm{m}}{\left({\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)+\nabla {\mathrm{f}}_{\mathrm{i}}{\left(\mathrm{x}\right)}^{\mathrm{T}}\mathrm{v}\right)}^2$$

(25)

where $\mathrm{v}=\Delta {\mathrm{x}}_{\mathrm{g}\mathrm{n}}$ is the needed exact value to minimize this approximation; furthermore, it is highlighted that the Gauss–Newton step direction can be estimated by solving a linear least-squares problem.

3.1.4. The Integer Least-Squares Problem

A particular case among least-squares problems is when the design variables are integers and these types of problems are named integer-least-squares (ILS) problems; they arise, for example, to solve the integer frequency ambiguity in Global Navigation Satellite Systems (GNSSs) for real-time applications, or in other fields such as communications, cryptography, and lattice design [38].

The pseudorange and carrier phase are extracted from GNSS measurements, and only the latter allows pose estimation with high precision in the order of the millimeter and centimeter level. However, the unknown integer part of the carrier phase, i.e., the frequency integer ambiguity, needs to be solved, and it can be performed by adopting the geometry-based model, which relies on the resolution of an ILS problem [39]. Such a problem can be defined as follows [38]:

$$\underset{\mathrm{a}\in {\mathbb{Z}}^{\mathrm{n}}}{\min } \left({\left(\mathrm{a}-\widehat{\mathrm{a}}\right)}^{\mathrm{T}}{\mathrm{Q}}_{\widehat{\mathrm{a}}}^{-1}\left(\mathrm{a}-\widehat{\mathrm{a}}\right)\right)$$

(26)

where $\mathrm{a}$ is the double-differenced integer ambiguity vector, $\widehat{\mathrm{a}}$ is the integer estimation of $\mathrm{a}$, and ${\mathrm{Q}}_{\widehat{\mathrm{a}}}$ is the symmetric positive $\mathrm{n}\times \mathrm{n}$ variance–covariance matrix.

To solve the problem represented by Equation (22), specific algorithms, such as the least-squares AMBiguity Decorrelation Adjustment (LAMBDA) algorithm or its variations, have been developed to reduce the computational resolution complexity of the problem [38]. They are composed of two parts: the reduction process, and the discrete search process. In this review, they are not investigated further because these algorithms are considered numerical techniques adopted to simplify the resolution of an optimization problem.

3.1.5. Integer Frequency Ambiguities Resolution

The authors in reference [39] propose a single-epoch triple-frequency real-time kinematic positioning method that exploits the information provided by the GNSS to obtain high-precision positioning in urban environments.

In their article, they propose the Random Sample Consensus (RANSAC) that is adopted after the ILS is solved to detect and exclude incorrectly fixed ambiguities of the Extra-Wide-Lane (EWL), Wide-Lane (WL), and original frequencies. The method deals with a triple-frequency integer ambiguity, and the double-differenced pseudorange ${\mathrm{P}}_{\mathrm{i}}$ and carrier phase ${\mathrm{L}}_{\mathrm{i}}$ observations for the $\mathrm{i}$-th frequency are given as follows:

$$\begin{array}{l}{\mathrm{P}}_{\mathrm{i}}=\mathsf{\rho }+{\mathrm{M}}_{\mathrm{i}}+{\mathsf{\epsilon }}_{\mathrm{i}}\\ {\mathrm{L}}_{\mathrm{i}}=\mathsf{\rho }+{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{N}}_{\mathrm{i}}+{\mathrm{m}}_{\mathrm{i}}+{\mathsf{ϵ}}_{\mathrm{i}} \end{array}$$

(27)

where $\mathsf{\rho }$ is the double-differenced geometric range; ${\mathrm{M}}_{\mathrm{i}}$ and ${\mathrm{m}}_{\mathrm{i}}$ are the multipath errors in the pseudorange and carrier phase, respectively, ${\mathsf{\lambda }}_1$ is the wavelength of the carrier; ${\mathrm{N}}_{\mathrm{i}}$ is the integer ambiguity, and ${\mathsf{\epsilon }}_{\mathrm{i}}$ and ${\mathsf{ϵ}}_{\mathrm{i}}$ are random noises of the pseudorange and carrier phase, respectively. By considering three frequencies, i.e., $\mathrm{i}=\left\{1,2,3\right\}$, as presented in reference [39], it is possible to define two robust double-differenced geometric ranges, ${\tilde{\mathsf{\rho }}}_1$ and ${\tilde{\mathsf{\rho }}}_2$, which allows for the computation of the EWL and WL ambiguities that are given, respectively, as follows:

$${\mathrm{N}}_{(0,-\mathrm{1,1})}=\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\left(\left({\displaystyle \frac{{\mathrm{L}}_3}{{\mathsf{\lambda }}_3}}-{\displaystyle \frac{{\mathsf{\lambda }}_2}{{\mathsf{\lambda }}_2}}\right)-\left({\displaystyle \frac{{\tilde{\mathsf{\rho }}}_1}{{\mathsf{\lambda }}_3}}-{\displaystyle \frac{{\tilde{\mathsf{\rho }}}_1}{{\mathsf{\lambda }}_2}}\right)\right)$$

(28)

$${\mathrm{N}}_{(1,-\mathrm{1,0})}=\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\left(\left({\displaystyle \frac{{\mathrm{L}}_1}{{\mathsf{\lambda }}_1}}-{\displaystyle \frac{{\mathsf{\lambda }}_2}{{\mathsf{\lambda }}_2}}\right)-\left({\displaystyle \frac{{\tilde{\mathsf{\rho }}}_2}{{\mathsf{\lambda }}_1}}-{\displaystyle \frac{{\tilde{\mathsf{\rho }}}_2}{{\mathsf{\lambda }}_2}}\right)\right).$$

(29)

Equations (28) and (29) allow for the definition of the following system of equations:

$$\left\{\begin{array}{l}{\mathrm{N}}_{(1,-\mathrm{1,0})}={\mathrm{N}}_1-{\mathrm{N}}_2\\ {\mathrm{N}}_{(0,-\mathrm{1,1})}={\mathrm{N}}_3-{\mathrm{N}}_2\\ {\mathrm{L}}_1=\mathsf{\rho }+{\mathsf{\lambda }}_1{\mathrm{N}}_1+{\mathrm{m}}_1+{\mathsf{ϵ}}_1\\ {\mathrm{L}}_2=\mathsf{\rho }+{\mathsf{\lambda }}_2{\mathrm{N}}_2+{\mathrm{m}}_2+{\mathsf{ϵ}}_2\\ {\mathrm{L}}_3=\mathsf{\rho }+{\mathsf{\lambda }}_3{\mathrm{N}}_3+{\mathrm{m}}_3+{\mathsf{ϵ}}_3\end{array}\right.$$

(30)

and by adopting an ILS method, it is possible to estimate $\mathsf{\rho }$, ${\mathrm{N}}_1$, ${\mathrm{N}}_2$, and ${\mathrm{N}}_3$ of the problems.

3.2. Levenberg–Marquardt

The Levenberg–Marquardt algorithm (LMA), or the Damped least-squares (DLS) method, is used to solve non-linear optimization problems. The LMA exploits a parameter $\mathsf{\lambda }$ such that it can switch between the GD algorithm and the Gauss–Newton algorithm; practically, the LMA interpolates between the Gauss–Newton algorithm and the method of GD according to the current prediction accuracy, but at the expense of a certain convergence rate [40]. The algorithm needs the initial guess and, if it is too far, then the error might be large, and the algorithm might not be able to give a correct global solution due to the presence of many local minima, and it might provide a local minimum which is not necessarily the global minimum. To obtain a better correct initial guess for the LMA, it is possible to use EAs in the initial phase because they do not suffer from the initial value problem.

Considering a non-linear least-squares problem, it can be expressed mathematically as follows:

$$\mathrm{f}\left(\mathrm{x}\right)={\displaystyle \frac{1}{2}}\sum _{\mathrm{j}=1}^{\mathrm{m}}{\mathrm{r}}_{\mathrm{j}}^2\left(\mathrm{x}\right)={\displaystyle \frac{1}{2}}{‖\mathrm{r}\left(\mathrm{x}\right)‖}^2$$

(31)

where $\mathrm{x}$ is the vector input, $\mathrm{f}$ is the objective function, ${\mathrm{r}}_{\mathrm{i}}$ is the residual and it is defined as a function such that ${\mathbb{R}}^{\mathrm{n}}\to \mathbb{R}$ with the assumption of $\mathrm{m}\ge \mathrm{n}$, and $\mathrm{r}\left(\mathrm{x}\right)={(\mathrm{r}}_1\left(\mathrm{x}\right),{\mathrm{r}}_2\left(\mathrm{x}\right),\dots ,{\mathrm{r}}_{\mathrm{m}}{\left(\mathrm{x})\right)}^{\mathrm{T}}$ is the residual vector and is defined such that ${\mathbb{R}}^{\mathrm{m}}\to {\mathbb{R}}^{\mathrm{n}}$.

The Levenberg–Marquardt algorithm is a modification of the Gauss–Newton method, and it is based on a local linearization of the residuals as follows [40,41,42]:

$${\mathrm{r}}_{\mathrm{m}}\left(\mathrm{x}+\mathsf{\delta }\mathrm{x}\right)\approx {\mathrm{r}}_{\mathrm{m}}\left(\mathrm{x}\right)+{\mathrm{J}}_{\mathrm{m}\mathsf{\mu }}\mathsf{\delta }\mathrm{x}$$

(32)

where ${\mathrm{J}}_{\mathrm{m}\mathsf{\mu }}=\mathrm{J}\left(\mathrm{x}\right)={\displaystyle \frac{\partial {\mathrm{r}}_{\mathrm{m}}}{\partial {\mathrm{x}}_{\mathsf{\mu }}}}$ is the Jacobian Matrix, and $\mathsf{\delta }\mathrm{x}$ is an infinitesimal of $\mathrm{x}$. Next, the Gauss–Newton technique recursively updates in accordance with the following:

$$\mathsf{\delta }\mathrm{x}=-{\left({\mathrm{J}}^{\mathrm{T}}\mathrm{J}\right)}^{-1}\nabla \mathrm{f}=-{\left({\mathrm{J}}^{\mathrm{T}}\mathrm{J}\right)}^{-1}{\mathrm{J}}^{\mathrm{T}}\mathrm{r}.$$

(33)

If the initial guess of the Gauss–Newton method is near the optimum of $\mathrm{f}$, it will usually converge quickly; however, this is not always the case, and the method can take very large and uncontrolled steps, failing to converge. To deal with this issue and control the steps in a useful way, Levenberg [41] and Marquardt [42] each proposed the damping of the ${\mathrm{J}}^{\mathrm{T}}\mathrm{J}$ matrix by a diagonal cut-off, modifying (33) as follows:

$$\mathsf{\delta }\mathrm{x}=-{\left({\mathrm{J}}^{\mathrm{T}}\mathrm{J}+\mathsf{\lambda }{\mathrm{D}}^{\mathrm{T}}\mathrm{D}\right)}^{-1}\nabla \mathrm{f}$$

(34)

where ${\mathrm{D}}^{\mathrm{T}}\mathrm{D}$ is a positive–definite, diagonal matrix that is the relative scale of the parameters and $\mathsf{\lambda }$ is a damping parameter adjusted by the algorithm. If $\mathsf{\lambda }$ is large, the steps are small along the gradient direction, and vice versa when $\mathsf{\lambda }$ is small. As the method is reaching a solution, $\mathsf{\lambda }$ is chosen to be small and the method converges quickly via the Gauss–Newton method.

The implementation of the Levenberg–Marquardt algorithm follows these iterative steps [40]:

• Based on the current parameter values, update the function and Jacobian values (if necessary);
• Update the scaling matrix ${\mathrm{D}}^{\mathrm{T}}\mathrm{D}$ and damping parameter $\mathsf{\lambda }$;
• Propose a parameter step $\mathsf{\delta }\mathrm{x}$ and evaluate the function at the new parameter values $\mathrm{x}+\mathsf{\delta }\mathrm{x}$;
• Accept or reject the parameter step depending on whether the cost has decreased at the new parameters;
• Stop if any of the desired convergence criteria of the algorithm are met or if the limit of function evaluations or iterations has been exceeded.

The methods for the choice strategy of the damping parameter $\mathsf{\lambda }$ and the scaling matrix ${\mathrm{D}}^{\mathrm{T}}\mathrm{D}$ require specific methods (as described in [40]) which are not further investigated in this article.

Magnetic Localization Technique Based on the Levenberg–Marquardt Algorithm

The authors in reference [10] propose a method for the localization of a permanent magnet based on the LMA. The permanent magnet is modelled as a dipole, and the magnetic field is measured by an array of $4\times 4$ magnetometers (16 in total). The general optimization problem can be expressed as follows:

$$\mathrm{f}\left(\mathrm{a},\mathrm{b},\mathrm{c},\mathsf{\theta },\mathsf{\phi }\right)=\sum _{\mathrm{i}=1}^{\mathrm{N}}{\left|{\widehat{\mathrm{B}}}_{\mathrm{i}}-{\mathrm{B}}_{\mathrm{i}}\right|}^2$$

(35)

where $\mathrm{f}$ is the objective function, $\mathrm{a}$, $\mathrm{b},$ and $\mathrm{c}$ are the coordinates of the permanent magnet position, $\mathsf{\theta }$ and $\mathsf{\phi }$ identify the orientation of the permanent magnet, $\mathrm{N}$ is the number of the $\mathrm{i}$-th magnetometer, ${\widehat{\mathrm{B}}}_{\mathrm{i}}$ is the measured magnetic field, and ${\mathrm{B}}_{\mathrm{i}}$ is a function representing the model of the system based on the magnetic dipole theory, and it includes the coordinates and the orientation angles.

The problem can be solved by using the LMA, but the issue is the need for a good initial guess to obtain a solution because when a too-far initial point is provided, divergence can occur. To deal with this problem, a solution could be to combine the PSO and LMA methods: the former is used to calculate an initial guess, and the second is used to calculate the local minima because it results in a faster method. The drawback of this approach is that it is not suitable for real-time applications, and the authors of paper [10] deal with it by proposing to use only the LMA and a properly modified model of Equation (35). In the first step, the new model neglects the orientation to find only the coordinates of the permanent magnet; then, in the second step, the orientation is also obtained. The modified model is expressed as follows:

$$\mathrm{f}\left(\mathrm{a},\mathrm{b},\mathrm{c}\right)=\sum _{\mathrm{i}=1}^{\mathrm{N}}{\left|{\widehat{\mathrm{B}}}_{\mathrm{i}}-{\displaystyle \frac{{\mathrm{B}}_{\mathrm{T}}}{{\left|{\mathrm{R}}_{\mathrm{i}}\right|}^3}}{\mathrm{A}}_{\mathrm{i}}{\widehat{\mathrm{m}}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\right|}^2$$

(36)

where ${\mathrm{B}}_{\mathrm{T}}$ is a scalar quantity that represents the magnetic strength, ${\mathrm{R}}_{\mathrm{i}}$ is the vector from the center of the permanent magnet to the $\mathrm{i}$-th sensor, ${\mathrm{A}}_{\mathrm{i}}$ is a matrix that is a function of the magnet position, and ${\widehat{\mathrm{m}}}_{\mathrm{a}\mathrm{v}\mathrm{g}}$ is the averaged orientation vector, expressed as follows:

$${\widehat{\mathrm{m}}}_{\mathrm{a}\mathrm{v}\mathrm{g}}={\displaystyle \frac{\sum _{\mathrm{i}=1}^{\mathrm{N}}\frac{{\tilde{\mathrm{m}}}_{\mathrm{i}}}{\left|{\tilde{\mathrm{m}}}_{\mathrm{i}}\right|}}{\left|\sum _{\mathrm{i}=1}^{\mathrm{N}}\frac{{\tilde{\mathrm{m}}}_{\mathrm{i}}}{\left|{\tilde{\mathrm{m}}}_{\mathrm{i}}\right|}\right|}}$$

(37)

where ${\tilde{\mathrm{m}}}_{\mathrm{i}}={\displaystyle \frac{{\left|{\mathrm{R}}_{\mathrm{i}}\right|}^3}{2{\mathrm{B}}_{\mathrm{T}}}}({\mathrm{A}}_{\mathrm{i}}-\mathrm{I}){\widehat{\mathrm{B}}}_{\mathrm{i}}$, and $\mathrm{I}$ is the one-diagonal matrix. The orientation vector is obtained by applying the following equation:

$$\widehat{\mathrm{m}}={\displaystyle \frac{{\left|{\mathrm{R}}_{\mathrm{i}}\right|}^3}{2{\mathrm{B}}_{\mathrm{T}}}}\left({\mathrm{A}}_{\mathrm{i}}-\mathrm{I}\right){\mathrm{B}}_{\mathrm{i}}$$

(38)

The authors claim that in such a way, a more suitable optimization problem for the LMA is obtained, and for this application, the method always provides an acceptable solution independently from the distance of the initial guess; moreover, the speed of computing is improved.

3.3. Lagrange Multipliers Method and Lagrange Duality

3.3.1. Lagrange Multipliers Method

The Lagrange multipliers method is an analytic method for finding the local maxima and minima of a function subject to equality constraints. The concept is to reformulate the original problem, obtaining a new function called the Lagrangian function, which allows for solving a new problem without constraints. Given a problem with an objective function $\mathrm{f}\left(\mathrm{x}\right)$ subjected to the equality constraint $\mathrm{h}\left(\mathrm{x}\right)=0$, the Lagrangian function is expressed as follows [43]:

$$\mathrm{L}\left(\mathrm{x},\mathsf{\lambda }\right)\equiv \mathrm{f}\left(\mathrm{x}\right)+\mathsf{\lambda }\mathrm{h}(\mathrm{x})$$

(39)

where $\mathsf{\lambda }$ is the Lagrangian multiplier. The solution is obtained by finding the stationary point of $\mathrm{L}$ by setting ${\displaystyle \frac{\partial \mathrm{L}\left(\mathrm{x},\mathsf{\lambda }\right)}{\partial \mathrm{x}}}=0$ and ${\displaystyle \frac{\partial \mathrm{L}\left(\mathrm{x},\mathsf{\lambda }\right)}{\partial \mathsf{\lambda }}}=0$, which gives rise to the following system:

$$\left\{\begin{array}{l}{\displaystyle \frac{\partial \mathrm{f}\left(\mathrm{x}\right)}{\partial \mathrm{x}}}+\mathsf{\lambda }{\displaystyle \frac{\partial \mathrm{h}\left(\mathrm{x}\right)}{\partial \mathrm{x}}}=0\\ \mathrm{h}\left(\mathrm{x}\right)=0\end{array}\right.$$

(40)

The solution corresponding to the original problem is always a saddle point of the Lagrangian function [44]. Then, a possible solution of the system (40) $({\mathrm{x}}^{\mathrm{*}},{\mathsf{\lambda }}^{\mathrm{*}})$ allows for the finding of an optimal point of the original problem given by ${\mathrm{f}(\mathrm{x}}^{\mathrm{*}})$.

The method of Lagrange multipliers is generalized by the KKT condition, also allowing for the consideration of inequality constraints; given a minimization non-linear problem (in which the $\mathrm{x}$ variable is a vector) with an objective function $\mathrm{f}\left(\mathrm{x}\right)$ subject to equality constraints ${\mathrm{h}}_{\mathrm{i}}\left(\mathrm{x}\right)=0$ (with $\mathrm{i}=1,\dots ,\mathrm{p}$), and inequality constraints ${\mathrm{g}}_{\mathrm{j}}\left(\mathrm{x}\right)\le 0$ (with $\mathrm{j}=1,\dots ,\mathrm{m}$), the Lagrangian function is set as follows:

$$\mathrm{L}\left(\mathrm{x},\mathsf{\lambda },\mathsf{\mu }\right)\equiv \mathrm{f}\left(\mathrm{x}\right)+{\mathsf{\mu }}^{\mathrm{T}}\mathrm{g}\left(\mathrm{x}\right)+{\mathsf{\lambda }}^{\mathrm{T}}\mathrm{h}(\mathrm{x})$$

(41)

where $\mathsf{\mu }$ and $\mathsf{\lambda }$ are the Lagrange multipliers related to inequalities and equalities, respectively, and they are vectors. If $({\mathrm{x}}^{\mathrm{*}},{\mathsf{\lambda }}^{\mathrm{*}},{\mathsf{\mu }}^{\mathrm{*}})$ is a saddle point of $\mathrm{L}\left(\mathrm{x},\mathsf{\lambda },\mathsf{\mu }\right)$ in the $\mathrm{x}$ domain, and the KKT conditions are satisfied, then ${\mathrm{x}}^{\mathrm{*}}$ identifies an optimal solution. The KKT conditions can be summarized as follows [45]:

$$\begin{array}{l}{\mathrm{g}}_{\mathrm{j}}\left({\mathrm{x}}^{\mathrm{*}}\right)\le 0\\ {\mathrm{h}}_{\mathrm{i}}\left({\mathrm{x}}^{\mathrm{*}}\right)=0\\ \begin{array}{l}{\mathsf{\mu }}_{\mathrm{j}}^{\mathrm{*}}\ge 0\\ {\mathsf{\mu }}_{\mathrm{j}}^{\mathrm{*}}{\mathrm{g}}_{\mathrm{j}}\left({\mathrm{x}}^{\mathrm{*}}\right)=0\\ \nabla \mathrm{f}\left({\mathrm{x}}^{\mathrm{*}}\right)+\sum _{\mathrm{j}}{\mathsf{\mu }}_{\mathrm{j}}^{\mathrm{*}}\nabla {\mathrm{g}}_{\mathrm{j}}({\mathrm{x}}^{\mathrm{*}})+\sum _{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{i}}^{\mathrm{*}}\nabla {\mathrm{h}}_{\mathrm{i}}\left({\mathrm{x}}^{\mathrm{*}}\right)=0\end{array}\end{array}$$

(42)

In general, once the gradient has been obtained, to determine if a constrained local optimum has been found, the KKT conditions are used to verify the necessary conditions for a local optimum; for unconstrained problems, the KKT conditions only require the gradient of the objective function to vanish at the optimum design point.

It is highlighted that the KKT conditions cannot indicate whether a global optimum has been found, and they are useful only for identifying a local optimum.

3.3.2. Lagrange Duality

The duality principle claims that an initial optimization problem called the primal problem can be transformed into a dual problem by applying a specific procedure. In such a way, it allows for the transformation of a minimization problem into a maximization one and vice versa. In some cases, the method might make the solution of the transformed problem easier and more manageable than the primal problem [14].

Mathematically, the Lagrange duality is defined as follows [46]. A not necessarily convex optimization problem is considered and is represented as:

$$\begin{array}{lll}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}:& {\mathrm{f}}_0\left(\mathrm{x}\right)& \\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t} \mathrm{t}\mathrm{o}:& {\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)\le 0& \mathrm{i}=1,\dots ,\mathrm{m}\\ & {\mathrm{h}}_{\mathrm{i}}\left(\mathrm{x}\right)=0& \mathrm{i}=1,\dots ,\mathrm{p}\end{array}$$

(43)

where $\mathrm{x}\in {\mathbb{R}}^{\mathrm{n}}$, with a domain $\mathrm{D}$ and optimal value ${\mathrm{p}}^{\mathrm{*}}$. Then, by using (41), the Lagrangian function of the problem is obtained and given as:

$$\mathrm{L}\left(\mathrm{x},\mathsf{\lambda },\mathsf{\nu }\right)\equiv {\mathrm{f}}_0\left(\mathrm{x}\right)+\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)+\sum _{\mathrm{i}=1}^{\mathrm{p}}{\mathsf{\nu }}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{i}}\left(\mathrm{x}\right)$$

(44)

where $\mathrm{L}\left(\mathrm{x},\mathsf{\lambda },\mathsf{\nu }\right)$ is such that ${\mathbb{R}}^{\mathrm{n}}\times {\mathbb{R}}^{\mathrm{m}}\times {\mathbb{R}}^{\mathrm{p}}\to \mathbb{R}$ with a domain $\mathrm{L}=\mathrm{D}\times {\mathbb{R}}^{\mathrm{m}}\times {\mathbb{R}}^{\mathrm{p}}$, and ${\mathsf{\lambda }}_{\mathrm{i}}$ and ${\mathsf{\nu }}_{\mathrm{i}}$ are the Lagrange multipliers related to inequalities and equalities, respectively. Now, it is possible to define the Lagrangian dual function $\mathrm{g}\left(\mathsf{\lambda },\mathsf{\nu }\right)$ of the problem as the infimum of the Lagrangian function over $\mathrm{x}$ and it is given as:

$$\mathrm{g}\left(\mathsf{\lambda },\mathsf{\nu }\right)=\underset{\mathrm{xϵD}}{\inf } \mathrm{L}\left(\mathrm{x},\mathsf{\lambda },\mathsf{\nu }\right)=\underset{\mathrm{xϵD}}{\inf } \left({\mathrm{f}}_0\left(\mathrm{x}\right)+\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{f}}_{\mathrm{i}}\left(\mathrm{x}\right)+\sum _{\mathrm{i}=1}^{\mathrm{p}}{\mathsf{\nu }}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{i}}\left(\mathrm{x}\right)\right)$$

(45)

where $\mathrm{g}\left(\mathsf{\lambda },\mathsf{\nu }\right)$ is such that ${\mathbb{R}}^{\mathrm{m}}\times {\mathbb{R}}^{\mathrm{p}}\to \mathbb{R}$. It is observed that the infimum is unconstrained (unlike the beginning problem) and $\mathrm{g}\left(\mathsf{\lambda },\mathsf{\nu }\right)$ is concave; moreover, the function $\mathrm{g}$ could be $-\mathrm{\infty }$ for some $\mathsf{\lambda }$, $\mathsf{\nu }$.

At this point, by adopting the lower bound property, if $\mathsf{\lambda }\ge 0$, then $\mathrm{g}\left(\mathsf{\lambda },\mathsf{\nu }\right)\le {\mathrm{p}}^{\mathrm{*}}$; this allows us to define the dual problem by maximizing the obtained Lagrangian dual function in (45). It is given as follows:

$$\begin{array}{l}\underset{\mathsf{\lambda },\mathsf{\nu }}{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}} \mathrm{g}\left(\mathsf{\lambda },\mathsf{\nu }\right)\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t} \mathrm{t}\mathrm{o} \mathsf{\lambda }\ge 0\end{array}$$

(46)

The solution of (46), denoted with ${\mathsf{\lambda }}^{\mathrm{*}}$, ${\mathsf{\nu }}^{\mathrm{*}}$, provides for all convex and non-convex problems a lower bound for the optimum function such that $\mathrm{g}\left({\mathsf{\lambda }}^{\mathrm{*}},{\mathsf{\nu }}^{\mathrm{*}}\right)\le \mathrm{f}({\mathrm{x}}^{\mathrm{*}})$, like it is represented in Figure 6.

Furthermore, Figure 6 clearly represents the weak and strong duality principles. The former provides a bound for the optimal values of the transformed problem, whereas the latter claims that the optimal values of the two problems are equal. The strong duality is usually true for convex optimization problems.

3.3.3. A Localization Non-Convex Problem

The authors in reference [47] propose a method to obtain localization by exploiting only the range measurement data provided by wireless sensors. The robot localization is modelled as a non-convex optimization problem, and to obtain the solution, as a first step, it is reformulated as a convex problem, and then duality theory is applied to find, in an easier way, the optimum.

The localization problem is formulated as follows. The measured range data from anchors are given as follows:

$${\widehat{\mathrm{d}}}_{\mathrm{k}}={‖\mathrm{x}-{\mathrm{y}}_{\mathrm{k}}‖}_2+{\mathsf{\eta }}_{\mathrm{k}}$$

(47)

where $\mathrm{x}$ and ${\mathrm{y}}_{\mathrm{k}}$ are two column vectors to determine the positions of the target and the anchors, respectively, ${\mathsf{\eta }}_{\mathrm{k}}$ is the corresponding measurement noise, $\mathrm{k}=1,\dots ,\mathrm{m}$ identifies the $\mathrm{k}$-th anchor, which are $\mathrm{m}$ in total. The optimal position can be obtained by minimizing the measurement noise between the measured distance ${\widehat{\mathrm{d}}}_{\mathrm{k}}$ and the real distance ${\mathrm{d}}_{\mathrm{k}}={‖\mathrm{x}-{\mathrm{y}}_{\mathrm{k}}‖}_2$; then, the objective function can be expressed as follows:

$${\mathrm{f}}_0\left(\mathrm{x}\right)=\sum _{\mathrm{k}=1}^{\mathrm{m}}{\left({\widehat{\mathrm{d}}}_{\mathrm{k}}-{‖\mathrm{x}-{\mathrm{y}}_{\mathrm{k}}‖}_2\right)}^2$$

(48)

The function in (48) is non-convex, and to solve the problem, the authors proceeded in this way. In the first step, they reformulated the problem into a convex-constrained problem by exploiting a change in the variable; the new problem is expressed as follows:

$$\begin{array}{l}\min \mathrm{f}\left(\mathrm{z}\right)={‖\mathrm{A}\mathrm{z}-\mathrm{b}‖}_2^2\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t} \mathrm{t}\mathrm{o} ({\mathrm{z}}^{\mathrm{T}}\mathrm{C}\mathrm{z}+2{\mathrm{f}}^{\mathrm{T}}\mathrm{z})=0\end{array}$$

(49)

where $\mathrm{z}=\left[\begin{array}{ll}{\mathrm{x}}^{\mathrm{T}}& {\mathrm{x}}^{\mathrm{T}}\mathrm{x}\end{array}\right]={\left[\begin{array}{ll}{\mathrm{x}}^{\mathrm{T}}& \mathrm{t}\end{array}\right]}^{\mathrm{T}}$ and coefficients are defined as $\mathrm{A}=\left[\begin{array}{ll}\begin{array}{l}\begin{array}{l}-2{\mathrm{y}}_1^{\mathrm{T}}\\ -2{\mathrm{y}}_2^{\mathrm{T}}\end{array}\\ \begin{array}{l}⋮\\ -2{\mathrm{y}}_{\mathrm{k}}^{\mathrm{T}}\end{array}\end{array}& \begin{array}{l}\begin{array}{l}1\\ 1\end{array}\\ \begin{array}{l}⋮\\ 1\end{array}\end{array}\end{array}\right]$, $\mathrm{b}=\left[\begin{array}{l}\begin{array}{l}{\mathrm{d}}_1^2-{‖{\mathrm{y}}_1‖}_2^2\\ {\mathrm{d}}_2^2-{‖{\mathrm{y}}_2‖}_2^2\end{array}\\ \begin{array}{l}⋮\\ {\mathrm{d}}_{\mathrm{k}}^2-{‖{\mathrm{y}}_{\mathrm{k}}‖}_2^2\end{array}\end{array}\right]$, $\mathrm{C}=\left[\begin{array}{lll}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]$, and $\mathrm{f}=\left[\begin{array}{l}0\\ 0\\ -1/2\end{array}\right]$. Then, duality theory is adopted, and the obtained Lagrangian function is given as follows:

$$\mathrm{L}\left(\mathrm{z},\mathrm{v}\right)={‖\mathrm{A}\mathrm{z}-\mathrm{b}‖}_2^2+\mathrm{v}\left({\mathrm{z}}^{\mathrm{T}}\mathrm{C}\mathrm{z}+2{\mathrm{f}}^{\mathrm{T}}\mathrm{z}\right)={\mathrm{z}}^{\mathrm{T}}\left({\mathrm{A}}^{\mathrm{T}}\mathrm{A}+\mathrm{v}\mathrm{C}\right)\mathrm{z}-2{\left({\mathrm{A}}^{\mathrm{T}}\mathrm{b}-\mathrm{v}\mathrm{f}\right)}^{\mathrm{T}}\mathrm{z}+{\mathrm{b}}^{\mathrm{T}}\mathrm{b}$$

(50)

where $\mathrm{v}$ is the Lagrange multiplier. Now, by computing the minima of $\mathrm{L}\left(\mathrm{z},\mathrm{v}\right)$ over $\mathrm{z}$, it is possible to obtain the Lagrangian dual function, and the final Lagrange dual problem is given as follows:

$$\begin{array}{ll}\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}& \mathrm{G}(\mathrm{v})={-{\left({\mathrm{A}}^{\mathrm{T}}\mathrm{b}-\mathrm{v}\mathrm{f}\right)}^{\mathrm{T}}\left({\mathrm{A}}^{\mathrm{T}}\mathrm{A}+\mathrm{v}\mathrm{C}\right)\left({\mathrm{A}}^{\mathrm{T}}\mathrm{b}-\mathrm{v}\mathrm{f}\right)+‖\mathrm{b}‖}_2^2\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t} \mathrm{t}\mathrm{o}& {(\mathrm{A}}^{\mathrm{T}}\mathrm{A}+\mathrm{v}\mathrm{C})\ge 0\\ & {(\mathrm{A}}^{\mathrm{T}}\mathrm{b}-\mathrm{v}\mathrm{f})\mathsf{ϵ}\mathrm{R}\left({\mathrm{A}}^{\mathrm{T}}\mathrm{A}+\mathrm{v}\mathrm{C}\right)\end{array}$$

(51)

where $\mathrm{R}$ identifies the domain. The problem in (51) is a convex problem and the solution can be computed.

3.4. Genetic Algorithm

The GA was inspired by Darwin’s principle of survival of the fittest solutions and can solve both constrained and unconstrained optimization problems; it is based on natural selection, the process that drives biological evolution [48]. The algorithm is implemented with the following steps:

• Creation of the initial population.
• Calculation of fitness of each individual.
• New individual generation by adopting the following operators:

  3.1.

  Selection;

  3.2.

  Crossover;

  3.3.

  Mutation;

• Calculation of fitness of each individual.
• Test:

  5.1.

  If the convergence criteria are satisfied, stop.

  5.2.

  Otherwise, repeat from step 3.

It is iterative, and after the first cycle, it is repeated from step 3 until the desired criteria are satisfied. In step 1, the process begins with a set of individuals, which is called the population; each individual is a solution to the problem to be solved. An individual is characterized by a set of parameters (variables) known as genes, which are joined into a string to form a chromosome, i.e., a solution. In the second step, the fitness function of the problem determines how an individual (or solution) fits, and it measures the ability of an individual to compete with other individuals. The output is a fitness score for each individual, and the probability that an individual will be selected for reproduction is based on its score. In step 3, the new individual generation process begins, and it is composed of four operators, i.e., (i) selection, (ii) crossover, and (iii) mutation. In the selection, some individuals are chosen to generate new ones and selected individuals pass some of their genes into the next generation; in other words, some individuals are selected as parents to generate children. In general, individuals are selected based on their fitness scores, and the higher the fitness value, the higher the probability of being selected for reproduction. In the literature, there are many proposals for selection operators, and some of them are the following: roulette wheel selection, sequential selection, tournament selection, dominant selection, hybrid selection, kin selection, back controlled selection, and so on [49]. It is important to highlight that, at every cycle, an elitist population is extracted from the population of the current iteration to directly pass some of the best individuals to the next generation. In some cases, such a procedure might give rise to problems of fast convergence, and to tackle such an issue, an adaptive elitist population strategy which changes the number of the elitist population can be implemented, as proposed in reference [50]. After selection, the crossover operator creates new children by mixing the relative genes. For example, Figure 7 shows the concept of the 1-crossover point operator, where two parents generate two offspring by exchanging portions of their genes. In the literature, there are many proposals for crossover operators and the authors in reference [51] classify them in three categories: standard (such as 1-crossover point, k-crossover point, shuffle-crossover, reduced surrogate-crossover, and so on), binary (random respectful crossover, masked crossover, 1-bit adaption crossover, and so on), and application dependent (or real/tree).

Finally, the mutation operators can introduce random variation of the genes, and for this purpose, the mutation probability value is adopted, which is a random value between 0 and 100%; a value equal to zero means no-mutations are adopted, whereas a 100% value modifies all genes of the child. In general, mutations should not occur so often, but they are important because they allow for the avoidance of being stuck in local minima [52]. For example, some mutation operators in the literature are mirror mutation and binary bit-flipping mutation, mutation based on directed variation techniques, directed mutation, and so on [53].

Once a new population is obtained, step 4 computes the fitness values of the offspring and, finally, step 5 verifies whether the convergence criteria are satisfied or not. If not, the cycle is repeated from step 3 until a satisfactory solution is found or the maximum number of iterations is reached; over successive generations, the population evolves toward an optimal solution [15,48,54].

The GA can solve a variety of optimization problems that are not well-suited for standard optimization algorithms, including problems in which the objective function is discontinuous, nondifferentiable, stochastic, or highly non-linear. Furthermore, GAs are also suitable for problems of mixed-integer programming where solutions have some components that are defined as integers.

Evolutionary Optimization Strategy for Indoor Position Estimation Using Smartphones

The authors in reference [9] propose a smartphone-based indoor localization system adopting a sensor fusion algorithm where localization information from an Inertial Measurement Unit (IMU), Wireless Local Area Network (WLAN), Bluetooth, and magnetic field anomalies is fused to provide the estimated localization.

In the fusion step, the probability distributions of each localization information source are superimposed; in the summation, each distribution is weighted with a coefficient in the interval from 0 to 1.

For determining the weights, optimization methods can be applied; they exploit the result of some tests that were run in buildings to obtain the ground truth points (GTPs). In such a way, it is possible to estimate the weights with a Root-Mean-Square Error (RMSE)-based fitness metric, defined as follows:

$$\mathrm{F}\mathrm{i}\mathrm{t}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s}=\prod _{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t} \mathrm{r}\mathrm{u}\mathrm{n}=1}^{\mathrm{N}}{\displaystyle \frac{10}{\mathrm{R}\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t} \mathrm{r}\mathrm{u}\mathrm{n}}}}$$

(52)

Due to the high ambiguity and the unknown (mixed) stochastic distributions of the input variables, using classical search methods (e.g., gradient-based methods), the optimization quickly converges to a local minimum. Then, a global optimization procedure is needed to maximize the quality metric and avoid local extrema. Because of the building model’s non-linearities, no gradients can be easily derived analytically, and gradient-based optimization methods are not applicable. Therefore, the authors proposed the use of a GA for finding the optimal weights, and in their article, they compared the results with four other local and global optimization strategies: Hill Climbing, the Nelder–Mead method, SA, and PSO. The outcome is that the GA outperforms the other algorithms; however, the drawback is that it needs the most computing time, but it is not a particular issue because the optimization process of the parameters must only be calculated once the end of the offline phase is reached.

3.5. Quantum Genetic Algorithm

The Quantum Genetic Algorithm (QGA) exploits the concept of qubits and the superposition of states in quantum mechanics. The algorithm is explained as follows [11,55]. The information is represented by a quantum bit or qubit which may be in the ‘0’ state, in the ‘1’ state, or any superposition of the two states [55]. A representation of the qubit state is given as follows:

$$\left|\mathsf{\Psi }⟩\right.=\mathsf{\alpha }\left|0⟩\right.+\mathsf{\beta }\left|1⟩\right.$$

(53)

where $\mathsf{\alpha }$ and $\mathsf{\beta }$ are complex numbers that specify the probability amplitudes of the corresponding states, and ${\left|\mathsf{\alpha }\right|}^2$ and ${\left|\mathsf{\beta }\right|}^2$ are the probabilities of the states to be either 0 or 1, respectively, and the condition ${\left|\mathsf{\alpha }\right|}^2+{\left|\mathsf{\beta }\right|}^2=1$ must be satisfied.

A system with $\mathrm{m}$ qubits represents $2^{\mathrm{m}}$ states at the same time, which collapse into a single state during the observing process of a quantum state.

In general, the traditional GA can use binary, numeric, or symbolic representation for a solution encoded with a chromosome. The QGA uses the qubits, where one qubit is represented as $\left[\begin{array}{l}\mathsf{\alpha }\\ \mathsf{\beta }\end{array}\right]$, whilst $\mathrm{m}$ qubits are represented as follows:

$$\left[\left.\begin{array}{l}{\mathsf{\alpha }}_1\\ {\mathsf{\beta }}_1\end{array}\right|\left.\begin{array}{l}{\mathsf{\alpha }}_2\\ {\mathsf{\beta }}_2\end{array}\right|\begin{array}{l}\dots \\ \dots \end{array}\left|\begin{array}{l}{\mathsf{\alpha }}_{\mathrm{m}}\\ {\mathsf{\beta }}_{\mathrm{m}}\end{array}\right.\right]$$

(54)

where ${\left|{\mathsf{\alpha }}_{\mathrm{i}}\right|}^2+{\left|{\mathsf{\beta }}_{\mathrm{i}}\right|}^2=1$ for $\mathrm{i}=1,2,\dots ,\mathrm{m}$. For example, considering $\mathrm{m}$ = 3 and the following chosen amplitudes of probability, a superimposition of states is represented as follows:

$$\left[\begin{array}{lll}{\displaystyle \frac{1}{\sqrt{2}}}& 1& {\displaystyle \frac{1}{2}}\\ {\displaystyle \frac{1}{\sqrt{2}}}& 0& {\displaystyle \frac{\sqrt{3}}{2}}\end{array}\right],$$

(55)

and the system state is expressed as follows:

$${\displaystyle \frac{1}{2\sqrt{2}}}\left|000⟩+\right.{\displaystyle \frac{\sqrt{3}}{2\sqrt{2}}}\left|001⟩+\right.{\displaystyle \frac{1}{2\sqrt{2}}}\left|100⟩+\right.{\displaystyle \frac{\sqrt{3}}{2\sqrt{2}}}\left|101⟩\right..$$

(56)

It has four pieces of information at the same time, i.e., the probabilities of $\left|000⟩\right.$, $\left|001⟩\right.$, $\left|100⟩\right.$, and $\left|101⟩\right.$ are ${\displaystyle \frac{1}{8}}$, ${\displaystyle \frac{3}{8}}$, ${\displaystyle \frac{1}{8}}$, and ${\displaystyle \frac{3}{8}}$, respectively. For a better comprehension of how coefficients are obtained, the coefficient ${\displaystyle \frac{1}{2\sqrt{2}}}$ of the binary number $\left|000⟩\right.$ is obtained by multiplying each ${\mathsf{\alpha }}_{\mathrm{i}}$ in the matrix (55), i.e., ${\mathsf{\alpha }}_1{\mathsf{\alpha }}_2{\mathsf{\alpha }}_3={\displaystyle \frac{1}{\sqrt{2}}}\mathrm{*}1\mathrm{*}{\displaystyle \frac{1}{2}}$, whereas the coefficient ${\displaystyle \frac{\sqrt{3}}{2\sqrt{2}}}$ of the number $\left|101⟩\right.$ is obtained by ${\mathsf{\beta }}_1{\mathsf{\alpha }}_2{\mathsf{\beta }}_3={\displaystyle \frac{1}{\sqrt{2}}}\mathrm{*}1\mathrm{*}{\displaystyle \frac{\sqrt{3}}{2}}$. The square power of a coefficient gives the probability of the related binary number. Moreover, the second digit is always zero, because ${\mathsf{\alpha }}_2=1$.

Because QGA computing can represent a superposition of states, it provides a better characteristic of diversity than classical approaches; Equation (56) is one qubit chromosome that represents four states at the same time, whereas, in a classical representation, at least four chromosomes (000, 001, 100, and 101) would be needed.

During the convergence process, the ${\left|\mathsf{\alpha }\right|}^2$ or ${\left|\mathsf{\beta }\right|}^2$ of the qubit chromosomes can approach 1 or 0, giving rise to a single state.

In the QGA, a population of qubit chromosomes $\mathrm{Q}\left(\mathrm{t}\right)=\left\{{\mathrm{q}}_1^{\mathrm{t}},{\mathrm{q}}_2^{\mathrm{t}},\dots ,{\mathrm{q}}_{\mathrm{n}}^{\mathrm{t}}\right\}$ at generation $\mathrm{t}$ is defined, where $\mathrm{n}$ is the size of the population, and ${\mathrm{q}}_{\mathrm{j}}^{\mathrm{t}}$ is a qubit chromosome with $\mathrm{m}$ qubits defined as follows:

$${\mathrm{q}}_{\mathrm{j}}^{\mathrm{t}}=\left[\left.\begin{array}{l}{\mathsf{\alpha }}_1^{\mathrm{t}}\\ {\mathsf{\beta }}_1^{\mathrm{t}}\end{array}\right|\left.\begin{array}{l}{\mathsf{\alpha }}_2^{\mathrm{t}}\\ {\mathsf{\beta }}_2^{\mathrm{t}}\end{array}\right|\begin{array}{l}\dots \\ \dots \end{array}\left|\begin{array}{l}{\mathsf{\alpha }}_{\mathrm{m}}^{\mathrm{t}}\\ {\mathsf{\beta }}_{\mathrm{m}}^{\mathrm{t}}\end{array}\right.\right].$$

(57)

Considering the main step to implement the QGA, in the initialization step of $\mathrm{Q}\left(\mathrm{t}\right)$, and all ${\mathsf{\alpha }}_{\mathrm{i}}^{\mathrm{t}}$ and ${\mathsf{\beta }}_{\mathrm{i}}^{\mathrm{t}}$ of each ${\mathrm{q}}_{\mathrm{j}}^{\mathrm{t}}$ are set to $1/\sqrt{2}$ and it means that the qubit chromosome at $\mathrm{t}=0$ represents the linear superposition of all possible states with the same probability $\left({\displaystyle \frac{1}{\sqrt{2^{\mathrm{m}}}}}\right)$, and it can be expressed as follows:

$$\left|{\mathsf{\Psi }}_{{\mathbf{q}}_{\mathrm{j}}^0}⟩\right.=\sum _{\mathrm{k}=1}^{2^{\mathrm{m}}}{\displaystyle \frac{1}{\sqrt{2^{\mathrm{m}}}}}\left|{\mathrm{S}}_{\mathrm{k}}⟩\right.$$

(58)

where ${\mathrm{S}}_{\mathrm{k}}$ is the $\mathrm{k}$-th state represented by the binary string ${\mathrm{x}}_1{\mathrm{x}}_2\dots {\mathrm{x}}_{\mathrm{m}}$, where ${\mathrm{x}}_{\mathrm{i}}$ can be equal to 0 or 1.

The next step is to generate a set of solutions from $\mathrm{Q}\left(\mathrm{t}\right)$, and this is called the observing process, which provides $\mathrm{P}\left(\mathrm{t}\right)=\left\{{\mathrm{p}}_1^{\mathrm{t}},{\mathrm{p}}_2^{\mathrm{t}},\dots ,{\mathrm{p}}_{\mathrm{n}}^{\mathrm{t}}\right\}$ at the time $\mathrm{t}$.

One binary solution, ${\mathrm{p}}_{\mathrm{j}}^{\mathrm{t}}$, for $\mathrm{j}$ from 1 to $\mathrm{n}$, is a binary string of length $\mathrm{m}$, and is formed by selecting each bit using the probability of qubit ${\left|{\mathsf{\alpha }}_{\mathrm{i}}\right|}^2$ (or ${\left|{\mathsf{\beta }}_{\mathrm{i}}\right|}^2$), which is compared with an ${\mathrm{r}}_{\mathrm{i}}$ value that is a randomly generated number between 0 and 1; if ${\mathrm{r}}_{\mathrm{i}}\ge {\left|{\mathsf{\alpha }}_{\mathrm{i}}\right|}^2$, then the $\mathrm{i}$-th element ${\mathrm{x}}_{\mathrm{i}}$ of the $\mathrm{j}$-th binary string ${\mathrm{p}}_{\mathrm{j}}^{\mathrm{t}}$ is 1, or, if ${\mathrm{r}}_{\mathrm{i}}<{\left|{\mathsf{\alpha }}_{\mathrm{i}}\right|}^2$, it is 0. Each solution ${\mathrm{p}}_{\mathrm{j}}^{\mathrm{t}}$ is evaluated to give a measure of its fitness. The initial best solution is then selected and stored among the solutions $\mathrm{P}\left(\mathrm{t}\right)$.

After initialization, the iterative part of the algorithm begins, and a set of binary solutions $\mathrm{P}\left(\mathrm{t}\right)$ is formed by observing the $\mathrm{Q}\left(\mathrm{t}-1\right)$ population; each binary solution in $\mathrm{P}\left(\mathrm{t}\right)$ is evaluated to give the fitness value. In the next step, called “update $\mathrm{Q}\left(\mathrm{t}\right)$”, a set of qubit chromosomes ${\mathrm{q}}_{\mathrm{j}}^{\mathrm{t}}$ is updated by applying an appropriate quantum gate $\mathrm{U}\left(\mathrm{t}\right)$. A quantum gate is a reversible gate represented by a unitary operator $\mathrm{U}$ acting on the qubit basis states, and it satisfies the relation ${\mathrm{U}}^{+}\mathrm{U}=\mathrm{U}{\mathrm{U}}^{+}$, where ${\mathrm{U}}^{+}$ is the Hermitian adjoint of $\mathrm{U}$. The quantum gate is selected or designed in agreement with the considered practical problem to be solved; for example, the rotation gates can be used for knapsack problems, and it is given as follows:

$$\mathrm{U}\left(\mathsf{\theta }\right)=\left[\begin{array}{ll}\cos (\mathsf{\theta })& -\mathrm{s}\mathrm{i}\mathrm{n} (\mathsf{\theta })\\ \sin (\mathsf{\theta })& \cos (\mathsf{\theta })\end{array}\right]$$

(59)

where $\mathsf{\theta }$ is the rotation angle. For instance, a qubit amplitude probability ${\left[\begin{array}{ll}{\mathsf{\alpha }}_{\mathrm{i}}& {\mathsf{\beta }}_{\mathrm{i}}\end{array}\right]}^{\mathrm{T}}$ in a chromosome can be updated by a quantum rotation that is given as follows:

$$\left[\begin{array}{l}{\mathsf{\alpha }}_{\mathrm{i}}^{\prime }\\ {\mathsf{\beta }}_{\mathrm{i}}^{\prime }\end{array}\right]=\mathrm{U}\left({\mathsf{\theta }}_{\mathrm{i}}\right)\left[\begin{array}{l}{\mathsf{\alpha }}_{\mathrm{i}}\\ {\mathsf{\beta }}_{\mathrm{i}}\end{array}\right]$$

(60)

where ${\mathsf{\theta }}_{\mathrm{i}}=\mathrm{s}\left({\mathsf{\alpha }}_{\mathrm{i}},{\mathsf{\beta }}_{\mathrm{i}}\right)\Delta {\mathsf{\theta }}_{\mathrm{i}}$, and $\mathrm{s}({\mathsf{\alpha }}_{\mathrm{i}},{\mathsf{\beta }}_{\mathrm{i}})$ is a function to determine the sign of ${\mathsf{\theta }}_{\mathrm{i}}$, and $\Delta {\mathsf{\theta }}_{\mathrm{i}}$ is the variation of ${\mathsf{\theta }}_{\mathrm{i}-1}$. The functions $\mathrm{s}({\mathsf{\alpha }}_{\mathrm{i}},{\mathsf{\beta }}_{\mathrm{i}})$ and $\Delta {\mathsf{\theta }}_{\mathrm{i}}$ are chosen in agreement with the lookup table (as indicated in reference [55]), which is based on experience and experiments; $\Delta {\mathsf{\theta }}_{\mathrm{i}}$ is a key parameter because it influences the computation quality and efficiency of the method; if it is too small, the number of iterations might be too big, whereas when it is too large, the solution might be trapped into divergent locations or give rise to a premature convergence to a local optimum [11].

Thanks to this step, the qubit chromosomes move and can converge to a better solution. In the next step, among the $\mathrm{P}\left(\mathrm{t}\right)$ solutions, the best one is selected and compared with the previous best one; if the new solution is better than the previous one, the older solution is discarded; otherwise, the old one is retained. By the application of some genetic operators, it is also possible to obtain mutations and crossovers of chromosomes; the mutation gives rise to new individuals characterized by a small change, whereas the crossover allows the generation of new individuals by combining parts of two or more different chromosomes. Genetic operators make the probability of the linear superposition of states change; however, mutation and crossover are not needed in the QGA because the quantum representation provides a diversity in chromosomes and the application of genetic gates with a high probability of variation can result in worse performance.

The diversity of the QGA can provide better convergence than conventional GAs, which can deal only with fixed 0 and 1 information. Moreover, the QGA has a speed of calculation $\sqrt{\mathrm{N}}$ times faster than the speed of the traditional genetic method. The main steps of the QGA are summarised as follows:

• Initialize of $\mathrm{Q}\left(\mathrm{t}=0\right)$.
• t = t + 1.
• Observe $\mathrm{Q}\left(\mathrm{t}-1\right)$ and create $\mathrm{P}\left(\mathrm{t}\right)$.
• Evaluation of $\mathrm{P}\left(\mathrm{t}\right).$
• Update $\mathrm{Q}\left(\mathrm{t}\right)$ using quantum gates $\mathrm{U}\left(\mathrm{t}\right)$
• Store the best solution among $\mathrm{P}\left(\mathrm{t}\right)$.
• If max generation is not reached or convergence criteria are not satisfied, repeat from step 2.

An Adaptive Quantum-Inspired Evolutionary Algorithm for Localization of Vehicles

The authors in reference [11] introduce a method based on an Adaptive Quantum-Inspired Evolutionary Algorithm (AQIEA) to estimate the calibration parameters for a real-size intelligent vehicle that can be driven along a non-predefined road path. They use the Differential Global Positioning System (DGPS) as a ground truth and propose an optimization-based calibration method to reduce the systematic errors of wheeled odometers (the differences between the actual wheel diameters and the nominal diameters) and the gyroscope (initial installation error due to misalignment of the sensor with the longitudinal axis of the vehicle). The diameter of the left and right rear wheels and the heading of the gyroscope are considered as the calibration parameters to be optimized to reduce the localization errors in the dead reckoning.

The dead reckoning model can be expressed as follows:

$$\begin{array}{l}{\mathrm{x}}_{\mathrm{t}+\Delta \mathrm{t}}={\mathrm{x}}_{\mathrm{t}}+\Delta \mathrm{d}·\cos ({\mathsf{\vartheta }}_{\mathrm{t}}+\Delta \mathsf{\vartheta }/2)\\ {\mathrm{y}}_{\mathrm{t}+\Delta \mathrm{t}}={\mathrm{y}}_{\mathrm{t}}+\Delta \mathrm{d}·\sin ({\mathsf{\vartheta }}_{\mathrm{t}}+\Delta \mathsf{\vartheta }/2)\\ {\mathsf{\vartheta }}_{\mathrm{t}+\Delta \mathrm{t}}={\mathsf{\vartheta }}_{\mathrm{t}}+\Delta \mathsf{\vartheta }\end{array}$$

(61)

where $\mathrm{x}$, $\mathrm{y}$, and $\mathsf{\vartheta }$ are, respectively, the coordinates of the vehicle along the $\mathrm{x}$ and $\mathrm{y}$-direction in the Cartesian coordinate system and the heading of the vehicle at time $\mathrm{t}$; $\Delta \mathrm{t}$ is the sampling period; $\Delta \mathrm{d}$ is the incremental travel displacement of the vehicle, and it can be expressed as $\Delta \mathrm{d}={\displaystyle \frac{\Delta {\mathrm{d}}_{\mathrm{r}}+\Delta {\mathrm{d}}_{\mathrm{l}}}{2}}$, where $\Delta {\mathrm{d}}_{\mathrm{r}(\mathrm{l})}$ is the displacement of the right (left) rear wheel of the vehicle, and $\Delta \mathsf{\vartheta }$ is the heading change during the sampling period.

The wheel displacement can be expressed as follows:

$$\Delta {\mathrm{d}}_{\mathrm{r}}={\displaystyle \frac{\Delta {\mathrm{N}}_{\mathrm{r}\left(\mathrm{l}\right)}}{{\mathrm{R}}_{\mathrm{r}}\left(\mathrm{l}\right)}}\mathsf{\pi }{\mathrm{D}}_{\mathrm{r}(\mathrm{l})}$$

(62)

where $\Delta {\mathrm{N}}_{\mathrm{r}\left(\mathrm{l}\right)}$ is the increment of the encoder pulses of the right (left) during the sampling period, ${\mathrm{R}}_{\mathrm{r}(\mathrm{l})}$ is the encoder resolution of the right (left) wheel in the unit of pulses per revolution, and ${\mathrm{D}}_{\mathrm{r}(\mathrm{l})}$ is the nominal diameter of the right (left) wheel.

In the real application, this model is influenced by systematic errors and random errors; for example, considering the odometers, the former are due to differences between wheel diameters, different nominal diameters, and the misalignment of wheels; the latter is due to the limited resolution of encoders, the limited sampling rate, the condition of the road (rough surface), the wheel slippage, and so on.

The systematic errors can be compensated and by including them in the model, the displacements and heading variation are expressed as follows:

$$\begin{array}{l}\Delta {\tilde{\mathrm{d}}}_{\mathrm{r}}=\Delta {\mathrm{d}}_{\mathrm{r}}+{\mathrm{c}}_1\\ \Delta {\tilde{\mathrm{d}}}_{\mathrm{l}}=\Delta {\mathrm{d}}_{\mathrm{l}}+{\mathrm{c}}_2\\ {\tilde{\mathsf{\vartheta }}}_{\mathrm{t}}={\mathsf{\vartheta }}_{\mathrm{t}}+{\mathrm{c}}_3\end{array}$$

(63)

where the symbol $~$ denotes the new variables that include systematic errors, and ${\mathrm{c}}_1$, ${\mathrm{c}}_2$, and ${\mathrm{c}}_3$ are the parameters to compensate for the right wheel diameter, the left wheel diameter, and the misalignment of the gyroscope, respectively. By considering these errors, the model represented by Equation (61) is rewritten as follows:

$$\begin{array}{l}{\tilde{\mathrm{x}}}_{\mathrm{t}+\Delta \mathrm{t}}={\tilde{\mathrm{x}}}_{\mathrm{t}}+{\displaystyle \frac{\Delta {\tilde{\mathrm{d}}}_{\mathrm{r}}+\Delta {\tilde{\mathrm{d}}}_{\mathrm{l}}}{2}}·\cos ({\tilde{\mathsf{\vartheta }}}_{\mathrm{t}}+\Delta \tilde{\mathsf{\vartheta }}/2)\\ {\tilde{\mathrm{y}}}_{\mathrm{t}+\Delta \mathrm{t}}={\tilde{\mathrm{y}}}_{\mathrm{t}}+{\displaystyle \frac{\Delta {\tilde{\mathrm{d}}}_{\mathrm{r}}+\Delta {\tilde{\mathrm{d}}}_{\mathrm{l}}}{2}}·\sin ({\tilde{\mathsf{\vartheta }}}_{\mathrm{t}}+\Delta \tilde{\mathsf{\vartheta }}/2)\\ {\tilde{\mathsf{\vartheta }}}_{\mathrm{t}+\Delta \mathrm{t}}={\tilde{\mathsf{\vartheta }}}_{\mathrm{t}}+\Delta \mathrm{t}.\end{array}$$

(64)

The positioning error ${\mathrm{E}}_{\mathrm{i}}$ at the $\mathrm{i}$-th sample is defined as the difference between the data provided by DGPS and the dead reckoning method, and it is computed by the following:

$${\mathrm{E}}_{\mathrm{i}}\left({\mathrm{c}}_1,{\mathrm{c}}_2,{\mathrm{c}}_3,{\mathrm{p}}_{\mathrm{i}}^{\mathrm{G}\mathrm{P}\mathrm{S}},{\mathrm{p}}_{\mathrm{i}}^{\mathrm{D}\mathrm{R}}\right)=\sqrt{{\left({\mathrm{x}}_{\mathrm{i}}^{\mathrm{G}\mathrm{P}\mathrm{S}}-{\tilde{\mathrm{x}}}_{\mathrm{i}}\right)}^2+{\left({\mathrm{y}}_{\mathrm{i}}^{\mathrm{G}\mathrm{P}\mathrm{S}}-{\tilde{\mathrm{y}}}_{\mathrm{i}}\right)}^2}$$

(65)

for $\mathrm{i}=1,2,\dots ,\mathrm{N}$ where $\mathrm{N}$ is the total number of samples.

Then, it is possible to formulate the optimization problem, which is expressed as follows:

$$\begin{array}{l}\begin{array}{ll}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}:& \mathrm{M}\left({\mathrm{c}}_1,{\mathrm{c}}_2,{\mathrm{c}}_3\right) \\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t} \mathrm{t}\mathrm{o}:& \left\{\begin{array}{l}{\mathrm{c}}_1^{\mathrm{L}}\le {\mathrm{c}}_1\le {\mathrm{c}}_1^{\mathrm{U}}\\ {\mathrm{c}}_2^{\mathrm{L}}\le {\mathrm{c}}_2\le {\mathrm{c}}_2^{\mathrm{U}}\\ {\mathrm{c}}_3^{\mathrm{L}}\le {\mathrm{c}}_3\le {\mathrm{c}}_3^{\mathrm{U}}\end{array}\right.\end{array}\\ \end{array}$$

(66)

where $\mathrm{M}$ is the objective function to be minimized with respect to ${\mathrm{c}}_1$, ${\mathrm{c}}_2$, and ${\mathrm{c}}_3$, and the superscript $\mathrm{L}$ and $\mathrm{U}$ denote the low and upper limits, respectively, for the considered variables. Two different optimization objective functions are considered to determine the optimal calibration parameters; the first one is the maximum positioning error, which is defined as follows:

$$\mathrm{M}=\mathrm{m}\mathrm{a}\mathrm{x} {\mathrm{E}}_{\mathrm{i}},$$

(67)

and the second one is the mean positioning error, defined as follows:

$$\mathrm{M}={\displaystyle \frac{1}{\mathrm{N}}}\sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{E}}_{\mathrm{i}}$$

(68)

To solve the optimization problem, the authors use the AQIEA, which is an evolution of the QGA. In the adaptive method, the lookup table is not used because it is substituted by a different approach based on the idea of choosing as input for the quantum gate either a large value for $\Delta {\mathsf{\theta }}_{\mathrm{i}}$ or a small one. The former is used when an individual is away from the optimal solution, and it needs to bring individuals quickly towards the optimal solution; the latter is used to avoid an individual moving away from the optimal solution when it is going to be reached. Then the magnitude of the rotation angle is obtained by the following:

$$\Delta {\mathsf{\theta }}_{\mathrm{i}}=\left\{\begin{array}{ll}{\mathsf{\theta }}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\left({\mathsf{\theta }}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathsf{\theta }}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right){\displaystyle \frac{\mathrm{f}-{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{e}}}{{\mathrm{f}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{e}}}} & \mathrm{f}>{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{e}}\\ {\mathsf{\theta }}_{\mathrm{m}\mathrm{a}\mathrm{x}}& \mathrm{f}\le {\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{e}}\end{array}\right.$$

(69)

where ${\mathsf{\theta }}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and ${\mathsf{\theta }}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ are two positive real numbers ($0.001\mathsf{\pi }\le {\mathsf{\theta }}_{\mathrm{m}\mathrm{i}\mathrm{n}}<{\mathsf{\theta }}_{\mathrm{m}\mathrm{a}\mathrm{x}}<0.5\mathsf{\pi }$), ${\mathrm{f}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum fitness in the current population, ${\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{e}}$ is the average fitness of the individuals in the current population, and $\mathrm{f}$ is the fitness of the selected individual. In such a way, an adaptive change in the rotation angle magnitude of the quantum rotation gate is obtained.

The sign of the rotation angle is determined by the following:

$$\mathrm{s}\left({\mathsf{\alpha }}_{\mathrm{i}},{\mathsf{\beta }}_{\mathrm{i}}\right)=\left\{\begin{array}{ll}\mathrm{s}\mathrm{g}\mathrm{n}\left(\left({\mathrm{o}}_{\mathrm{i}}^{\mathrm{*}}-\mathrm{c}\right){\mathsf{\alpha }}_{\mathrm{i}}{\mathsf{\beta }}_{\mathrm{i}}\right)& {\mathsf{\alpha }}_{\mathrm{i}}{\mathsf{\beta }}_{\mathrm{i}}\ne 0\\ 0& \left({\mathsf{\alpha }}_{\mathrm{i}}=0 \mathrm{a}\mathrm{n}\mathrm{d} {\mathrm{o}}_{\mathrm{i}}^{\mathrm{*}}=1\right) \mathrm{o}\mathrm{r} \left({\mathsf{\beta }}_{\mathrm{i}}=0 \mathrm{a}\mathrm{n}\mathrm{d} {\mathrm{o}}_{\mathrm{i}}^{\mathrm{*}}=0\right)\\ \pm 1& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.$$

(70)

where $\mathrm{s}\mathrm{g}\mathrm{n}$ is a function that provide the sign, ${\mathrm{o}}_{\mathrm{i}}^{\mathrm{*}}$ is the observed state (0 or 1) of the $\mathrm{i}$-th qubit of the best individual, and $\mathrm{c}$ is a positive real number between 0 and 1. The proof of (70) is given in reference [11].

Moreover, to improve global searching and to avoid being trapped in local optimal minima, an adaptive quantum mutation operation is introduced (unlike the traditional QGA, where mutations are not used). Given a selected chromosome in the population, the mutation operation is conducted based on a mutation probability ${\mathrm{p}}_{\mathrm{m}}$:

$${\mathrm{p}}_{\mathrm{m}}=\left\{\begin{array}{ll}{\mathrm{p}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\left({\mathrm{p}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{p}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right){\displaystyle \frac{\mathrm{f}-{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{e}}}{{\mathrm{f}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{e}}}}& \mathrm{f}>{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{e}}\\ {\mathrm{p}}_{\mathrm{m}\mathrm{a}\mathrm{x}}& \mathrm{f}\le {\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{e}}\end{array}\right.$$

(71)

where ${\mathrm{p}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and ${\mathrm{p}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ are two small positive real values between 0 and 1. It is observed, for example, that for an individual with $\mathrm{f}={\mathrm{f}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$, ${\mathrm{p}}_{\mathrm{m}}$ is 0 and no mutations are applied. The mutation is applied by choosing two random positions in the selected individual, and two qubits are mutated by a quantum-not gate that exchanges the values of the amplitude probabilities; the gate is defined as follows:

$$\mathrm{U}\left(\mathsf{\theta }\right)=\left[\begin{array}{ll}0& 1\\ 1& 0\end{array}\right].$$

(72)

3.6. Differential Evolution Algorithm

The DE algorithm is a power population-based and stochastic evolutionary algorithm to solve continuous global optimization problems that can be constrained or unconstrained; it maintains a population of candidate solutions and creates new candidate solutions by combining existing ones according to its formula and then keeping whichever candidate solution has the best score or fitness on the considered optimization problem. In this way, the optimization problem is treated as a black box that provides a measure of quality given by a candidate solution.

The main steps of the DE algorithm are as follows: (i) generate the initial population of the specified size; the population evolves by using (ii) mutation, (iii) crossover, and (iv) selection operators during each generation; and (v) step (ii) is repeated until a termination criterion is obtained [56].

In the initialization step, a number $\mathrm{N}$ of populations is set and each of them has a number $\mathrm{M}$ of $\mathrm{D}$-dimensional vectors ($\mathrm{D}$ is a number related to the nature of the problem), which are potential solutions to the problem. The $\mathrm{i}$-th individual, i.e., the $\mathrm{i}$-th $\mathrm{D}$-dimensional vector, is given as follows:

$${\mathrm{x}}_{\mathrm{i},\mathrm{G}}=\{{\mathrm{x}}_{\mathrm{i},\mathrm{G}}^1,{\mathrm{x}}_{\mathrm{i},\mathrm{G}}^2,\dots ,{\mathrm{x}}_{\mathrm{i},\mathrm{G}}^{\mathrm{D}}\}$$

(73)

where $\mathrm{G}$ identifies the number of the generation. Initial individuals are set with the following formula:

$${\mathrm{x}}_{\mathrm{i},0}={\mathrm{x}}_{\mathrm{m}\mathrm{i}\mathrm{n}}+\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\left(\mathrm{0,1}\right)\ast ({\mathrm{x}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{x}}_{\mathrm{m}\mathrm{i}\mathrm{n}})$$

(74)

where ${\mathrm{x}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ and ${\mathrm{x}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ are vectors which identify the limits of the solution, and $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\left(\mathrm{0,1}\right)$ is a function to generate random individuals with a normal distribution.

In the mutation step, mutation vectors ${\mathrm{v}}_{\mathrm{i},\mathrm{G}}$ are generated for each individual of a population, which is called the target vector. The basic mutation strategies are as follows [57]:

$$(1) \mathrm{DE}/\mathrm{rand}/1: {\mathrm{v}}_{\mathrm{i},\mathrm{G}}={\mathrm{x}}_{\mathrm{r}1,\mathrm{G}}+\mathrm{F}({\mathrm{x}}_{\mathrm{r}2,\mathrm{G}}-{\mathrm{x}}_{\mathrm{r}3,\mathrm{G}})$$

(75)

$$(2) \mathrm{DE}/\mathrm{best}/1: {\mathrm{v}}_{\mathrm{i},\mathrm{G}}={\mathrm{x}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t},\mathrm{G}}+\mathrm{F}({\mathrm{x}}_{\mathrm{r}1,\mathrm{G}}-{\mathrm{x}}_{\mathrm{r}2,\mathrm{G}})$$

(76)

$$(3) \mathrm{DE}/\mathrm{current} \mathrm{to} \mathrm{best}/1: {\mathrm{v}}_{\mathrm{i},\mathrm{G}}={\mathrm{x}}_{\mathrm{i},\mathrm{G}}+\mathrm{F}({\mathrm{x}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t},\mathrm{G}}-{\mathrm{x}}_{\mathrm{i},\mathrm{G}})+\mathrm{F}\left({\mathrm{x}}_{\mathrm{r}1,\mathrm{G}}-{\mathrm{x}}_{\mathrm{r}2,\mathrm{G}}\right)$$

(77)

$$(4) \mathrm{DE}/\mathrm{best}/2: {\mathrm{v}}_{\mathrm{i},\mathrm{G}}={\mathrm{x}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t},\mathrm{G}}+\mathrm{F}\left({\mathrm{x}}_{\mathrm{r}1,\mathrm{G}}-{\mathrm{x}}_{\mathrm{r}2,\mathrm{G}}\right)+\mathrm{F}({\mathrm{x}}_{\mathrm{r}3,\mathrm{G}}-{\mathrm{x}}_{\mathrm{r}4,\mathrm{G}})$$

(78)

$$(5) \mathrm{DE}/\mathrm{rand}/2: {\mathrm{v}}_{\mathrm{i},\mathrm{G}}={\mathrm{x}}_{\mathrm{r}1,\mathrm{G}}+\mathrm{F}\left({\mathrm{x}}_{\mathrm{r}2,\mathrm{G}}-{\mathrm{x}}_{\mathrm{r}3,\mathrm{G}}\right)+\mathrm{F}({\mathrm{x}}_{\mathrm{r}4,\mathrm{G}}-{\mathrm{x}}_{\mathrm{r}5,\mathrm{G}})$$

(79)

where $\mathrm{r}1$, $\mathrm{r}2$, $\mathrm{r}3$, $\mathrm{r}4$, and $\mathrm{r}5$ are integers randomly generated in the range between 1 and $\mathrm{M}$, $\mathrm{F}$ is a scale factor, and ${\mathrm{x}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t},\mathrm{G}}$ is the best individual vector in a population during the $\mathrm{G}$-th generation.

As for the crossover step, trial (or test) vectors ${\mathrm{u}}_{\mathrm{i},\mathrm{G}}$ are created by crossing each individual ${\mathrm{x}}_{\mathrm{i},\mathrm{G}}$ (or target vector) and the corresponding mutation vector ${\mathrm{v}}_{\mathrm{i},\mathrm{G}}$; it needs to form a set of trial vectors used in the selection step; there are two types of crossover methods largely used in the literature: binomial and exponential. In the former, the test vector is given as follows:

$${\mathrm{u}}_{\mathrm{i},\mathrm{G}}^{\mathrm{j}}=\left\{\begin{array}{l}{\mathrm{v}}_{\mathrm{i},\mathrm{G}}^{\mathrm{j}}\hspace{1em}\hspace{1em}\mathrm{i}\mathrm{f} \mathrm{r}\mathrm{a}\mathrm{n}{\mathrm{d}}_{\mathrm{i},\mathrm{j}}(\mathrm{0,1})\le \mathrm{C}\mathrm{R}\\ {\mathrm{x}}_{\mathrm{i},\mathrm{G}}^{\mathrm{j}}\hspace{1em}\hspace{1em}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.$$

(80)

where CR is the crossover factor in the range 0 and 1 and controls the probability that a component of a trial vector is obtained from either a mutant or a previous individual, and j is between 1 and $\mathrm{D}$. In this method, the distribution of the number of parameters inherited from the mutant is (almost) binomial. On the other hand, in the latter, a random integer $\mathrm{r}$ in the range $[1,\mathrm{D}]$ is chosen and it is the component starting point for the crossover process. The components for the trial vector ${\mathrm{u}}_{\mathrm{i},\mathrm{G}}^{\mathrm{j}}$ are donated from the mutated vector until a Bernoulli experiment linked to a probability $\mathrm{C}\mathrm{R}$ is true [58].

$${\mathrm{u}}_{\mathrm{i},\mathrm{G}}^{\mathrm{j}}=\left\{\begin{array}{l}\begin{array}{ll}{\mathrm{v}}_{\mathrm{i},\mathrm{G}}^{\mathrm{j}}& \mathrm{j}=\mathrm{r},\mathrm{r}+1,\dots \mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l} \mathrm{B}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{u}\mathrm{l}\mathrm{l}\mathrm{i} \mathrm{e}\mathrm{x}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t} \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e} \mathrm{o}\mathrm{r} \mathrm{j}=\mathrm{D}\end{array}\\ \begin{array}{ll}{\mathrm{x}}_{\mathrm{i},\mathrm{G}}^{\mathrm{j}}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}.\end{array}\right.$$

(81)

In the end, in the selection step, the objective function of the problem is evaluated for all generated trial vectors and the corresponding target vectors; the best one is selected, and it is a new individual of the next generation. The selection is given as follows:

$${\mathrm{x}}_{\mathrm{i},\mathrm{G}+1}=\left\{\begin{array}{l}{\mathrm{u}}_{\mathrm{i},\mathrm{G}}\hspace{1em}\hspace{1em}\mathrm{i}\mathrm{f} \mathrm{f}({\mathrm{u}}_{\mathrm{i},\mathrm{G}})\le \mathrm{f}({\mathrm{x}}_{\mathrm{i},\mathrm{G}})\\ {\mathrm{x}}_{\mathrm{i},\mathrm{G}}\hspace{1em}\hspace{1em}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.$$

(82)

There are two main differences between the GA and DE algorithm: (i) the GA is conducted by changing some genes in the chromosomes, while in the DE algorithm, the mutant individuals are obtained by adding the difference between two individuals to a third individual, which has a stronger global search capability; (ii) in the traditional GA, offspring individuals replace their parent individuals with a certain probability, while in the DE algorithm, such a condition only happens when the fitness of the offspring is better than the fitness of their parents; in such a way, a greatly increased convergence speed is provided [59].

Visible Light Communication Using a Modified Differential Evolution Algorithm

Authors in reference [59] have developed a Visible Light Communication (VLC) system where the VLC positioning model is transformed into a global optimization problem and a DE algorithm is used to calculate the z-coordinate of the receiver by adopting a modified DE algorithm that is called the self-adaptive DE algorithm; in this work, the 3D problem is reduced to a 1D problem, improving the speed of localization. Because the individual dimension is one (only the z-coordinate), the crossover process is ignored. The 3D positioning strategy exploits some specific equations to calculate the overlapped area of three circles determined by the position of the sources and the target, and is obtained by the projection on the $\mathrm{x}$-y plan; the DE algorithm is used to solve the specific equation and calculate the z-coordinate by finding the optimal fitness value; the $\mathrm{x}$ and $\mathrm{y}$ coordinates are calculated by exploiting three equations for trilateration in the projection plan.

The proposed algorithm has three steps: (i) the generation of the population, (ii) calculation of the fitness, and (iii) mutation operator. The fitness function evaluates the quality of each of the candidate solutions in the population. In the VLC-based 3D positioning system, the individual with the smallest overlapped area can be regarded as the estimated position. The fitness equation can be treated as the overlapped area $\mathrm{S}$ and is represented as follows:

$$\mathrm{f}\mathrm{i}\mathrm{t}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s}\left(\mathrm{z}\right)=\mathrm{S}$$

(83)

If no overlapped area is formed between the three circles, the fitness is zero; the presumed height $\mathrm{z}$ is increased by 0.2 m for each iteration until the overlapped area appears. If the fitness ${\le 10}^{-7}$, then the currently estimated height jumps out of the iteration loop and it is considered as the optimal solution ${\mathrm{z}}_{\mathrm{e}}$; otherwise, the iteration process goes on (step (iii)). In the case that the number of iterations reaches the maximum, and the fitness is still higher than ${10}^{-7}$, the algorithm goes back to step (i).

The authors considered three types of DE algorithms: (i) DE/rand/1/binomial; (ii) DE/current-to-best/1/binomial; and (iii) the proposed self-adaptive DE algorithm. The first and the second are two traditional DE algorithms which adopt the strategy given in (75) and in (77), respectively; furthermore, the DE/current-to-best/1/binomial benefits from its fast convergence because the best solutions in the search space are used and its performance is better than the performance of the DE/rand/1/binomial algorithm. As for the third algorithm, the proposed self-adaptive DE algorithm is derived on the following observation: as the iteration number increases, the difference among individuals is reduced, and therefore a smaller amount of adjustment is needed because an individual has a better fitness and most of its information can be retained; on the contrary, when the number of iterations are low, the solution is still far from the optimal solution and a larger scaling factor might help the exploration process. Then, the authors proposed the following variation:

$${\mathrm{v}}_{\mathrm{i},\mathrm{g}+1}={\mathrm{x}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t},\mathrm{g}}+\mathrm{F}({\mathrm{x}}_{\mathrm{r}2,\mathrm{g}}-{\mathrm{x}}_{\mathrm{r}3,\mathrm{g}})$$

(84)

where $\mathrm{F}$ is defined as follows:

$$\mathrm{F}=\left({\displaystyle \frac{\left|\mathrm{f}\mathrm{i}\mathrm{t}\left({\mathrm{x}}_{\mathrm{i}}\right)-\mathrm{f}\mathrm{i}\mathrm{t}({\mathrm{x}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}})\right|}{\mathsf{\sigma }+\mathrm{f}\mathrm{i}\mathrm{t}\left({\mathrm{x}}_{\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{t}}\right)-\mathrm{f}\mathrm{i}\mathrm{t}({\mathrm{x}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}})}}\right)\times \left({\mathrm{F}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{F}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+{\mathrm{F}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$$

(85)

where ${\mathrm{x}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$ and ${\mathrm{x}}_{\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{t}}$ are the best and worst individuals in the population, namely the vectors with the lowest and highest fitness values in the entire search space, and $\mathsf{\sigma }={10}^{-13}$ is a constant added to avoid the zero division. In such a way, the scaling factor varies between the maximum ${\mathrm{F}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ and minima ${\mathrm{F}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ values according to its difference with the best individual at each generation.

3.7. Mind Evolutionary Algorithm

The Mind Evolutionary Algorithm (MEA) is a method which emulates the learning modes and activities of the human mind. Unlike genetic-based algorithms, which use mutation, crossing, and selection operators, in the MEA, an evolutionary process is also exploited but with the difference that two operations named similartaxis and dissimilation are adopted to overcome the drawback of the bad solutions generated during the mutation process; furthermore, every generation’s evolutionary information is memorized, and it is used to guide the evolution toward the global optimum [60].

The MEA steps are defined as follows:

• Initialization.
• Similartaxis and local competition.
• Dissimilation and global competition.

In the initialization, two main spaces are considered: the superior group space, which is composed of $\mathrm{S}$ superior groups, and the temporary group space with $\mathrm{T}$ temporary groups. Every group is initialized with $\mathrm{m}$ individuals, which are generated by considering an individual as the center and adopting a uniform distribution to generate the other individuals around the center. Furthermore, every individual is scored in agreement with their adaptability to the search space [60,61].

The similartaxis and local competition phase aims to make the subgroups mature by performing a local competition. In every group, individuals compete with each other to find the local optimum. To create a new generation, a winner who has the best score is selected, and a new generation is created around it. Then, individuals are scored again to identify if there is a better new solution or to keep the previous one; in such a way, a new generation can be created, competition starts again, and so on. The information of every new winner is recorded on the local billboard of the group. The generation process in a group continues until the maximal and minimal scores within one group satisfy the mature condition; then, the information related to the mature group is stored in the global billboard. The process to make the group mature is named similartaxis and is performed for both superior and temporary space [60,61].

In the dissimilation and global competition phase, all groups compete against each other by observing the information reported on the global billboard. If, in the temporary space, there is a group with a higher score than any group in the superior space, then this temporary group will replace the bad groups in the superior space. Moreover, if the score of any mature group in the temporary space is lower than the score of any group in the superior one, this temporary group will be discarded, and new ones will be generated in the temporary space with random values. The generation process is called dissimilation, and in such a way, the similartaxis can continue in every group [60,61].

The aforementioned operations are repeated until the scores of the groups in the superior space are so high that they cannot be further improved. At this point, the algorithm is considered convergent, and the winning individual of the group with the best score in the superior space is the global optimum [61]. The architecture of the MEA is reported in Figure 8 [60].

A Mind Evolutionary Algorithm-Based Geomagnetic Positioning System for Indoor Localization

The authors in reference [60] propose a solution to solve the problem in the pedestrian dead reckoning field where a large error is introduced when a walker goes backwards or laterally but the position is still updated forward. To address this problem, they propose a system for the indoor localization of smartphones based on magnetic positioning and the pedestrian dead reckoning (PDR) method, which are fused using an Extended Kalman Filter (EKF).

To improve the performance of magnetic positioning, the MEA-based heuristic searching strategy is adopted to search for the optimal magnetic position.

The geomagnetic features are extracted using a three-axis magnetometer installed on a smartphone. The magnetic vector measured in the device coordinate reference system is transformed into the equivalent vector in the geographic coordinate reference system in order to obtain a measurement that is independent of the attitude of the device. The transformation can be expressed as follows:

$${\mathrm{B}}_{\mathrm{G}}={\left({\mathrm{C}}_{\mathrm{n}}^{\mathrm{b}}\right)}^{\mathrm{T}}{\mathrm{B}}_{\mathrm{D}}$$

(86)

where ${\mathrm{B}}_{\mathrm{G}}$ is the magnetic vector in the geographic coordinate system, ${\mathrm{C}}_{\mathrm{n}}^{\mathrm{b}}$ is the matrix for the transformation, and ${\mathrm{B}}_{\mathrm{D}}$ is the magnetic vector in the device coordinate system.

Then, five magnetic features for matching are extracted from ${\mathrm{B}}_{\mathrm{G}}$ and they are the three magnetic vector components, the horizontal component of the magnetic field, and the magnitude of the vector.

A fine-grained magnetic database is constructed by linearly interpolating the features and the coordinates. Once the fingerprint has been obtained, the MEA can be performed to estimate the magnetic position. The proposed algorithm exploits the large number of samples that are provided by the measurement process, which can be considered to be at a high frequency; over an interval time of one second, these samples generate many temporary positions that are used to generate a population for the MEA; by considering an epoch every second identified by $\mathrm{k}$, the population $\mathrm{M}(\mathrm{k})$ can be represented as follows:

$$\mathrm{M}(\mathrm{k})=\{{\mathrm{g}}_{\mathrm{k}}\left({\mathrm{x}}_1,{\mathrm{y}}_1\right),\dots ,{\mathrm{g}}_{\mathrm{k}}\left({\mathrm{x}}_{\mathrm{i}},{\mathrm{y}}_{\mathrm{i}}\right),\dots ,{\mathrm{g}}_{\mathrm{k}}\left({\mathrm{x}}_{\mathrm{n}},{\mathrm{y}}_{\mathrm{n}}\right)\}$$

(87)

where ${\mathrm{g}}_{\mathrm{k}}\left({\mathrm{x}}_{\mathrm{i}},{\mathrm{y}}_{\mathrm{i}}\right)$ is the $\mathrm{i}$-th individual (i.e., the estimated position) of $\mathrm{M}(\mathrm{k})$, ${\mathrm{x}}_{\mathrm{i}}$ and ${\mathrm{y}}_{\mathrm{i}}$ are the coordinates values, $\mathrm{i}=1,2,\dots ,\mathrm{n}$ is the number of the sample, and $\mathrm{n}$ is the sampling frequency. The score is calculated considering the following formula:

$$\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\{\mathrm{k},\mathrm{i}\}={\displaystyle \frac{1}{{‖\mathrm{G}\left(\mathrm{k}-1\right)-{\mathrm{g}}_{\mathrm{k}}({\mathrm{x}}_{\mathrm{i}},{\mathrm{y}}_{\mathrm{i}})‖}_2}}$$

(88)

where $\mathrm{G}\left(\mathrm{k}-1\right)$ is the previous estimated magnetic position.

As the first step, the MEA initializes the population (using the data generated by measurements), and the temporary individuals (i.e., the geomagnetic positions) are scored using Equation (88). The positions with higher scores are selected as centers of the superior groups, whereas the positions with lower scores are selected as centers of the temporary groups. Around the centers, individuals of the population are generated with a random distribution, and they are scored using Equation (88).

As the second step, the MEA executes the similartaxis where individuals within the same subgroup compete by comparing their scores. The second step is completed when the mature condition $\left|\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right|\le \mathsf{\epsilon }$ is obtained.

In the end, as the third step, the MEA executes dissimilation for the global competition. The MEA executes similartaxis and dissimilation until convergent conditions are satisfied and the individual with the best score in the mature superior groups is selected as the estimation of the final magnetic position.

3.8. Particle Swarm Optimization

PSO algorithms are based on a simplified social model that tries to emulate the behavior of a swarm of animals such as bees or birds that are searching for a food source. It can be explained as follows [15,62].

PSO finds a solution by exploiting information from each individual of the population (or swarm) and also information from the swarm as a whole. At the beginning of the process, the number of particles in the swarm and the number of iterations are set, and an initial population distributed throughout the design space is randomly generated. In the next step, the position of each particle is updated using the following formula:

$${\mathrm{x}}_{\mathrm{i}}^{\mathrm{q}+1}={\mathrm{x}}_{\mathrm{i}}^{\mathrm{q}}+{\mathrm{v}}_{\mathrm{i}}^{\mathrm{q}}\Delta \mathrm{t}$$

(89)

where $\mathrm{i}$ is the $\mathrm{i}$-th individual in the swarm, $\mathrm{q}$ is the $\mathrm{q}$-th iteration, ${\mathrm{v}}_{\mathrm{i}}^{\mathrm{q}}$ is the velocity vector of the $\mathrm{i}$-th individual at the $\mathrm{q}$-th iteration, and $\Delta \mathrm{t}$ is the time increment. At the beginning, the velocity vector of each particle is initialized randomly and is updated at each iteration using the following formula:

$${\mathrm{v}}_{\mathrm{i}}^{\mathrm{q}+1}={\mathsf{\omega }\mathrm{v}}_{\mathrm{i}}^{\mathrm{q}}+{\mathrm{c}}_1{\mathrm{r}}_1{\displaystyle \frac{({\mathrm{p}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}^{\mathrm{q}})}{\Delta \mathrm{t}}}+{\mathrm{c}}_2{\mathrm{r}}_2{\displaystyle \frac{({\mathrm{p}}^{\mathrm{g}}-{\mathrm{x}}_{\mathrm{i}}^{\mathrm{q}})}{\Delta \mathrm{t}}}$$

(90)

where $\mathsf{\omega }$ is the inertial parameter, ${\mathrm{r}}_1$ and ${\mathrm{r}}_2$ are random numbers between 0 and 1, ${\mathrm{c}}_1$ and ${\mathrm{c}}_2$ are the trust parameters, ${\mathrm{p}}_{\mathrm{i}}$ is the best point found so far by the $\mathrm{i}$-th particle, and ${\mathrm{p}}^{\mathrm{g}}$ is the best point found by the swarm. The search behavior of the algorithm is controlled by the inertial parameter $\mathsf{\omega }$; the larger the value, the larger the focus in global searching, whereas the smaller the value, the larger the focus in local searching; reference values for these two extremes are 1.4 and 0.5, for either a more focused global search or for more focused local search, respectively [15]. The ${\mathrm{c}}_1$ and the ${\mathrm{c}}_2$ are two trust parameters; the former, also referred to as cognitive memory, indicates how much the particle trusts itself, whereas the latter, also referred to as social memory, indicates how much the particle trusts the group. The ${\mathrm{p}}^{\mathrm{g}}$ can be defined as either considering some best points among a subset of particles or considering the best point among all particles of the swarm. The values of $\mathsf{\omega }$, ${\mathrm{c}}_1$, and ${\mathrm{c}}_2$ parameters must be tuned in agreement with the managed problem.

The main drawbacks of the method are the tuning parameters and the fact that it is an inherently unconstrained optimization algorithm. The last one is a common characteristic of most of the evolutionary algorithms; researchers have tried to deal with this problem by proposing many constraint-handling techniques for evolution, and among these, the penalty function approach is the most widely used because it has the advantage that it is easy to implement [15].

3.8.1. Particle Swarm Optimization-Based Artificial Neural Network (PSO-ANN)

The authors in reference [63] use an Artificial Neural Network (ANN) to solve the fingerprint localization problem by exploiting the Received Signal Strengths Indicator (RSSI) of the Wi-Fi signal in an indoor environment. The ANN provides the localization of the target, and it is trained by a traditional PSO algorithm; this approach provides a strong learning skill which allows one to better deal with complicated indoor environments. Traditionally ANNs are trained by exploiting a backpropagation algorithm that is characterized by a feedforward propagation of signals and a backpropagation of errors; this approach requires significant time for training and can be easily trapped in a local minimum. To solve these issues, authors adopt the PSO algorithm to find the parameters of the trained ANN; they define the following:

(i)

Each fingerprint is composed by $\mathrm{F}=[{\mathsf{\phi }}_{\mathrm{i}}\left(\mathrm{t}\right),\mathrm{P}]$, where ${\mathsf{\phi }}_{\mathrm{i}}\left(\mathrm{t}\right)$ is the RSSI value at time $\mathrm{t}$ of the $\mathrm{i}$-th Wi-Fi source, and $\mathrm{P}$ is the position of the collected fingerprint.

(ii)

The particle of the PSO is given as $\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{l}{\mathrm{e}}_{\mathrm{j}}=[{\mathrm{W}}_1,\dots ,{\mathrm{W}}_{\mathrm{n}},{\mathrm{B}}_1,\dots ,{\mathrm{B}}_{\mathrm{n}}]$, where ${\mathrm{W}}_{\mathrm{j}}$ and ${\mathrm{B}}_{\mathrm{j}}$ represent the weight and bias of the $\mathrm{j}$-th neuron.

(iii)

The updating equations for the particle position and velocity are defined by the following:

$$\begin{array}{l}{\mathrm{p}}_{\mathrm{k}+1}={\mathrm{p}}_{\mathrm{k}}+{\mathrm{v}}_{\mathrm{k}+1}\\ {\mathrm{v}}_{\mathrm{k}+1}=\mathsf{\omega }{\mathrm{v}}_{\mathrm{k}}+{\mathrm{r}}_1{\mathrm{c}}_1\times \left(\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}{\mathrm{t}}_{\mathrm{k}}-{\mathrm{p}}_{\mathrm{k}}\right)+{\mathrm{r}}_2{\mathrm{c}}_2\times \left(\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}{\mathrm{t}}_{\mathrm{k}}-{\mathrm{p}}_{\mathrm{k}}\right)\end{array}$$

(91)

where ${\mathrm{v}}_{\mathrm{k}}$ is the velocity in the $\mathrm{k}$-th iteration, ${\mathrm{p}}_{\mathrm{k}}$ is the particle in the $\mathrm{k}$-th iteration, ${\mathrm{c}}_1$ and ${\mathrm{c}}_2$ are two acceleration coefficients, and $\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}{\mathrm{t}}_{\mathrm{k}}$ and g$\mathrm{b}\mathrm{e}\mathrm{s}{\mathrm{t}}_{\mathrm{k}}$ are the personal best position and global best position in the $\mathrm{k}$-th iteration, respectively, ${\mathrm{r}}_1$ and ${\mathrm{r}}_2$ are two random constriction coefficients in the range (0, 1) sampled from a uniform distribution to prevent explosion and induce particles to converge on optima, and $\mathsf{\omega }$ is the inertia factor which helps to control the ability of particles to jump out of the local optimum and find the global optimum.

Considering the above definitions, the proposed algorithm follows these steps:

• Set the particle group size M, the max iteration number $\mathrm{m}$, and the expected error e. Then, initialize randomly the initial position and velocity of the particles;
• For all RSSI training sets $[{\mathsf{\phi }}_{\mathrm{i}}\left(\mathrm{t}\right),\mathrm{P}]$, put $\left[{\mathsf{\phi }}_{\mathrm{i}}\left(\mathrm{t}\right)\right]$ in the ANN model and compute the ANN output $\tilde{\mathrm{P}}$;
• Update the pbest and gbest according to the distance between the predicted position $\tilde{\mathrm{P}}$ and the real position $\mathrm{P}$;
• Update the particle position and velocity according to Equation (91);
• Repeat steps 3 and 4 until the distance between $\tilde{\mathrm{P}}$ and the real position $\mathrm{P}$ is less than the expected error or the max iteration number n has been reached.

3.8.2. Particle Swarm Optimization-Based Minimum Residual Algorithm

The authors in reference [64] adopt an RSSI range-based localization method for a mobile robot working in an indoor environment where the localization is estimated by a Minimum Residual (MinRe) localization algorithm based on PSO.

According to the principle of MinRE localization, the sum of the target’s squared residual between the real and the measured distances is given as follows:

$$\mathrm{f}(\mathrm{X})=\sum _{\mathrm{i}}^{\mathrm{N}}{\left\{‖\mathrm{X}-\mathrm{A}{\mathrm{N}}_{\mathrm{i}}‖-{\mathrm{d}}_{\mathrm{i}}^{\mathrm{m}\mathrm{e}\mathrm{s}}\right\}}^2$$

(92)

where $\mathrm{f}(\mathrm{X})$ is the function to be minimized, $\mathrm{X}$ is the real vector position to be found on the target, $\mathrm{A}{\mathrm{N}}_{\mathrm{i}}$ is the vector position of the $\mathrm{i}$-th anchor node, $\mathrm{N}$ is the total number of nodes, and ${\mathrm{d}}_{\mathrm{i}}^{\mathrm{m}\mathrm{e}\mathrm{s}}$ is the measured distance from the RSSI, which can be expressed as ${\mathrm{d}}_{\mathrm{i}}^{\mathrm{m}\mathrm{e}\mathrm{s}}={\mathrm{d}}_{\mathrm{i}}+{\mathsf{\omega }}_{\mathrm{i}}$, where ${\mathrm{d}}_{\mathrm{i}}$ is the unknown real distance and ${\mathsf{\omega }}_{\mathrm{i}}$ is the noise.

The term $‖\mathrm{X}-\mathrm{A}{\mathrm{N}}_{\mathrm{i}}‖$ represents the Euclidean distance between $\mathrm{X}$ and $\mathrm{A}{\mathrm{N}}_{\mathrm{i}}$ and it can be expressed as follows:

$$‖\mathrm{X}-\mathrm{A}{\mathrm{N}}_{\mathrm{i}}‖=\sqrt{{\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{i}}\right)}^2+{\left(\mathrm{y}-{\mathrm{y}}_{\mathrm{i}}\right)}^2{\left(\mathrm{z}-{\mathrm{z}}_{\mathrm{i}}\right)}^2}$$

(93)

where $\mathrm{x}$, $\mathrm{y},$ and $\mathrm{z}$ are the coordinates of the position of the target $\mathrm{X}$, and ${\mathrm{x}}_{\mathrm{i}}$, ${\mathrm{y}}_{\mathrm{i}}$, and ${\mathrm{z}}_{\mathrm{i}}$ are the coordinates of the position of the $\mathrm{A}{\mathrm{N}}_{\mathrm{i}}$. If there is no error, $\mathrm{f}(\mathrm{X})$ would be zero.

The PSO algorithm is applied to estimate the position of the target and the final problem can be written as follows:

$$\mathrm{X}=\arg \mathrm{m}\mathrm{i}\mathrm{n} \mathrm{f}(\mathrm{X})=\sum _{\mathrm{i}}^{\mathrm{N}}{\left\{\sqrt{{\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{i}}\right)}^2+{\left(\mathrm{y}-{\mathrm{y}}_{\mathrm{i}}\right)}^2+{\left(\mathrm{z}-{\mathrm{z}}_{\mathrm{i}}\right)}^2}-{\mathrm{d}}_{\mathrm{i}}^{\mathrm{m}\mathrm{e}\mathrm{s}}\right\}}^2$$

(94)

Setting the number of particles $\mathrm{M}$, the maximum number of iterations $\mathrm{W}$, considering $\mathrm{w}$ as the index of the particle updating iteration for this problem, and $\mathrm{m}$ as the index of the $\mathrm{m}$-th particle, the solution of the problem follows these steps:

• Initialization, where the particles and their initial velocity are generated randomly within the target environment: ${\mathrm{X}}_{\mathrm{m}}^1=[{\mathrm{x}}_{\mathrm{m}}^1,{\mathrm{y}}_{\mathrm{m}}^1,{\mathrm{z}}_{\mathrm{m}}^1]$, and ${\mathrm{V}}_{\mathrm{m}}^1=[\mathrm{v}{\mathrm{x}}_{\mathrm{m}}^1,{\mathrm{v}\mathrm{y}}_{\mathrm{m}}^1,{\mathrm{v}\mathrm{z}}_{\mathrm{m}}^1]$. The fitness $\mathrm{f}({\mathrm{X}}_{\mathrm{m}}^1)$ of the $\mathrm{m}$-th particle is calculated and the initial individual best value is set: ${\mathrm{B}}_{\mathrm{m}}^1={\mathrm{X}}_{\mathrm{m}}^1$. The initial global best fitness of the particle swarm is the particle which achieves the minimum of $\mathrm{f}(\mathrm{X})$ and it is represented as ${\mathrm{B}}_{\mathrm{g}}^1={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n} \mathrm{f}(\mathrm{X}}_{\mathrm{m}}^1)$ where $\mathrm{m}=1:\mathrm{M}$.
• Particle velocities and positions, which are updated in agreement with the individual best and global best values at $\mathsf{\omega }$-th iteration for $2\le \mathrm{w}\le \mathrm{W}$.

  $$\begin{array}{l}{\mathrm{V}}_{\mathrm{m}}^{\mathrm{w}}=\mathsf{\omega }{\mathrm{V}}_{\mathrm{m}}^{\mathrm{w}-1}+{\mathrm{c}}_1\mathsf{\xi }\left({\mathrm{B}}_{\mathrm{m}}^{\mathrm{w}-1}-{\mathrm{X}}_{\mathrm{m}}^{\mathrm{w}-1}\right)+{\mathrm{c}}_2\mathsf{\eta }\left({\mathrm{B}}_{\mathrm{g}}^{\mathrm{w}-1}-{\mathrm{X}}_{\mathrm{m}}^{\mathrm{w}-1}\right) \\ {\mathrm{X}}_{\mathrm{m}}^{\mathrm{w}}={\mathrm{X}}_{\mathrm{m}}^{\mathrm{w}-1}+{\mathrm{V}}_{\mathrm{m}}^{\mathrm{w}-1}\end{array}$$

  (95)

  where $\mathsf{\omega }$ is the inertia weight that controls the exploration search region, ${\mathrm{c}}_1$ and ${\mathrm{c}}_2$ are the acceleration constants between 0 and 2, and $\mathsf{\xi }$ and $\mathsf{\eta }$ are random numbers between 0 and 1.
• The updating of the optimal particle and optimal swarm is initiated:

  $$\begin{array}{l}{\mathrm{B}}_{\mathrm{m}}^{\mathrm{w}}=\left\{\begin{array}{l}{\mathrm{B}}_{\mathrm{m}}^{\mathrm{w}-1}, {\mathrm{f}(\mathrm{X}}_{\mathrm{m}}^{\mathrm{w}})>{\mathrm{f}(\mathrm{B}}_{\mathrm{m}}^{\mathrm{w}-1}) \\ {\mathrm{X}}_{\mathrm{m}}^{\mathrm{w}}, {\mathrm{f}(\mathrm{X}}_{\mathrm{m}}^{\mathrm{w}})\le {\mathrm{f}(\mathrm{B}}_{\mathrm{m}}^{\mathrm{w}-1}) \end{array}\right.\\ {\mathrm{B}}_{\mathrm{g}}^{\mathrm{w}}=\left\{\begin{array}{l}{\mathrm{B}}_{\mathrm{g}}^{\mathrm{w}-1}, {\mathrm{m}\mathrm{i}\mathrm{n}}_{\mathrm{m}=1:\mathrm{M}}{\mathrm{f}(\mathrm{X}}_{\mathrm{m}}^{\mathrm{w}})>{\mathrm{f}(\mathrm{B}}_{\mathrm{g}}^{\mathrm{w}-1}) \\ \arg ({\mathrm{m}\mathrm{i}\mathrm{n}}_{\mathrm{m}=1:\mathrm{M}}{\mathrm{f}(\mathrm{X}}_{\mathrm{m}}^{\mathrm{w}})), {\mathrm{m}\mathrm{i}\mathrm{n}}_{\mathrm{m}=1:\mathrm{M}}{\mathrm{f}(\mathrm{X}}_{\mathrm{m}}^{\mathrm{w}})\le {\mathrm{f}(\mathrm{B}}_{\mathrm{g}}^{\mathrm{w}-1}).\end{array}\right. \end{array}$$

  (96)
• The stopping criterion is checked: if $\mathrm{w}<\mathrm{W}$, then set $\mathrm{w}=\mathrm{w}+1$ and go to step 2; otherwise, the estimated position is ${\tilde{\mathrm{X}}=\mathrm{B}}_{\mathrm{g}}^{\mathrm{w}}$.

The estimated $\tilde{\mathrm{X}}$ is the position at the time $\mathrm{t}$ where the localization is estimated.

3.8.3. Mobile Robot Localization Based on PSO and Simultaneous Perturbation

The authors in reference [65] propose a mobile robot localization based on the PSO estimator where the system exploits the dynamic optimization to find the best robot pose estimate, recursively. The proposed solution is an alternative to using EKF or PF. PSO searches recursively and stochastically along the state space for the best robot pose estimation. Unlike the localization based on PF, in the proposed method, the resampling step (where samples with negligible weights are discarded and substituted by resampling) is not required, and no noise distribution is required either. The robot localization ${\mathrm{x}}_{\mathrm{t}}$ at the $\mathrm{t}$-th epoch is modelled like a Markov process (a random process in which the future is independent of the past given the present), where the initial distribution is ${\mathrm{p}(\mathrm{x}}_{\mathrm{t}})$. Considering a maximum, a posteriori point of view, the localization estimation can be considered as an optimization problem, where the process seeks to estimate the pose that maximizes the a posteriori probability density, defined as follows:

$$\arg \underset{\mathrm{x}}{\max } \mathrm{p}({\mathrm{x}}^{\mathrm{t}}|{\mathrm{z}}^{\mathrm{t}})$$

(97)

where ${\mathrm{x}}^{\mathrm{t}}=\left\{{\mathrm{x}}_0,\dots ,{\mathrm{x}}_{\mathrm{t}}\right\}$ is the history of the estimated positions, and ${\mathrm{z}}^{\mathrm{t}}=\left\{{\mathrm{z}}_0,\dots ,{\mathrm{z}}_{\mathrm{t}}\right\}$ is the history of the observations. By applying Bayes, the posterior over the robot path is given by the following:

$$\mathrm{p}\left({\mathrm{x}}^{\mathrm{t}}|{\mathrm{z}}^{\mathrm{t}}\right)=\mathsf{\eta }\mathrm{p}({\mathrm{z}}_{\mathrm{t}}|{\mathrm{x}}^{\mathrm{t}},{\mathrm{z}}^{\mathrm{t}-1})\mathrm{p}({\mathrm{x}}_{\mathrm{t}}|{\mathrm{x}}^{\mathrm{t}-1},{\mathrm{z}}^{\mathrm{t}-1})\mathrm{p}({\mathrm{x}}^{\mathrm{t}-1}|{\mathrm{z}}^{\mathrm{t}-1})$$

(98)

where $\mathsf{\eta }$ is a normalizing constant. The problem can be reformulated as follows:

$$\arg \underset{\mathrm{x}}{\max } \left(\prod _{\mathrm{i}=1}^{\mathrm{t}}\mathrm{p}\left({\mathrm{z}}_{\mathrm{i}}|{\mathrm{x}}^{\mathrm{i}},{\mathrm{z}}^{\mathrm{i}-1}\right)\prod _{\mathrm{i}=1}^{\mathrm{t}}\mathrm{p}\left({\mathrm{x}}_{\mathrm{i}}|{\mathrm{x}}^{\mathrm{i}-1},{\mathrm{z}}^{\mathrm{i}-1}\right)\mathrm{p}({\mathrm{x}}_0)\right)$$

(99)

where $\mathrm{p}\left({\mathrm{x}}_0\right)=\mathrm{p}({\mathrm{x}}^{\mathrm{t}-1}|{\mathrm{z}}^{\mathrm{t}-1})$ and $\mathrm{i}$ is an index to count the epochs. After some mathematics passages, the objective function can be expressed as follows:

$${\mathrm{f}}_0\left({\mathrm{x}}_{\mathrm{t}}\right)=\log \mathrm{p}\left({\mathrm{z}}_{\mathrm{t}}|{\mathrm{x}}_{\mathrm{t}}\right)+\log \mathrm{p}({\mathrm{x}}_{\mathrm{t}}|{\mathrm{x}}_{\mathrm{t}-1})+{\mathrm{f}}_0({\mathrm{x}}_{\mathrm{t}-1})$$

(100)

where ${\mathrm{f}}_0({\mathrm{x}}_{\mathrm{t}-1})$ considers the previous history for all the $\mathrm{t}-1$ epochs; it is suitable for an iterative application.

Now, the following is considered: (i) the motion model that gives the robot pose at time $\mathrm{t}$. Specifically, this is considered as follows:

$${\mathrm{x}}_{\mathrm{t}}=\mathrm{f}\left({\mathrm{x}}_{\mathrm{t}-1},{\mathrm{u}}_{\mathrm{t}}\right)+{\mathsf{\omega }}_{\mathrm{t}-1}~\mathrm{p}({\mathrm{x}}_{\mathrm{t}}|{\mathrm{x}}_{\mathrm{t}-1},{\mathrm{u}}_{\mathrm{t}})$$

(101)

where $\mathrm{f}$ is a non-linear function, ${\mathsf{\omega }}_{\mathrm{t}}$ is the Gaussian noise distribution with covariance matrix ${\mathrm{Q}}_{\mathrm{t}}$, and ${\mathrm{u}}_{\mathrm{t}}$ is the control; (ii) the measurement model for the observed ${\mathrm{Q}}_{\mathrm{t}}$ feature is defined as follows:

$${\mathrm{z}}_{\mathrm{t}}=\mathrm{h}\left({\mathrm{x}}_{\mathrm{t}},\mathsf{\theta }\right)+{\mathrm{v}}_{\mathrm{t}}~{{\mathrm{p}}_{\mathsf{\theta }}}_{\mathrm{t}}({\mathrm{z}}_{\mathrm{t}}|{\mathrm{x}}_{\mathrm{t}})$$

(102)

where $\mathrm{h}$ is the function of the model, ${\mathrm{v}}_{\mathrm{t}}$ is the measurement Gaussian noise distribution with covariance matrix ${\mathrm{R}}_{\mathrm{t}}$; (iii) the state transition probability density ${{\mathrm{p}}_{\mathsf{\omega }}}_{\mathrm{t}}({\mathrm{x}}_{\mathrm{t}}|{\mathrm{x}}_{\mathrm{t}-1})$ derived from the robot motion model and the observation probability density ${{\mathrm{p}}_{\mathsf{\upsilon }}}_{\mathrm{t}}\left({\mathrm{z}}_{\mathrm{t}}|{\mathrm{x}}_{\mathrm{t}}\right)$ derived from the observation model the problem is rewritten in the following form:

$$\underset{{\mathrm{x}}_{\mathrm{t}}}{\min } (\log {{\mathrm{p}}_{\mathsf{\upsilon }}}_{\mathrm{t}}\left({\mathrm{z}}_{\mathrm{t}}|{\mathrm{x}}_{\mathrm{t}}\right)+\log {{\mathrm{p}}_{\mathsf{\omega }}}_{\mathrm{t}}({\mathrm{x}}_{\mathrm{t}}|{\mathrm{x}}_{\mathrm{t}-1}))$$

(103)

where ${{\mathrm{p}}_{\mathsf{\upsilon }}}_{\mathrm{t}}\left({\mathrm{z}}_{\mathrm{t}}|{\mathrm{x}}_{\mathrm{t}}\right)\propto \exp \left({\displaystyle \frac{1}{2}}\left({\mathrm{z}}_{\mathrm{t}}-{\widehat{\mathrm{z}}}_{\mathrm{t}}\right){\mathrm{R}}_{\mathrm{t}}^{-1}{\left({\mathrm{z}}_{\mathrm{t}}-{\widehat{\mathrm{z}}}_{\mathrm{t}}\right)}^{\mathrm{T}}\right)$ and ${{\mathrm{p}}_{\mathsf{\omega }}}_{\mathrm{t}}({\mathrm{x}}_{\mathrm{t}}|{\mathrm{x}}_{\mathrm{t}-1})\propto \exp \left({\displaystyle \frac{1}{2}}\left({\mathrm{x}}_{\mathrm{t}}-{\widehat{\mathrm{x}}}_{\mathrm{t}}\right){\mathrm{Q}}_{\mathrm{t}}^{-1}{\left({\mathrm{x}}_{\mathrm{t}}-{\widehat{\mathrm{x}}}_{\mathrm{t}}\right)}^{\mathrm{T}}\right)$, where ${\widehat{\mathrm{z}}}_{\mathrm{t}}$ and ${\widehat{\mathrm{x}}}_{\mathrm{t}}$ are the estimated measurement vector and the estimated state vector, respectively; then, the fitness function can be written as follows:

$$\mathrm{F}\mathrm{i}\mathrm{t}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s}=\left({\mathrm{z}}_{\mathrm{t}}-{\widehat{\mathrm{z}}}_{\mathrm{t}}\right){\mathrm{R}}_{\mathrm{t}}^{-1}{\left({\mathrm{z}}_{\mathrm{t}}-{\widehat{\mathrm{z}}}_{\mathrm{t}}\right)}^{\mathrm{T}}+\left({\mathrm{x}}_{\mathrm{t}}-{\widehat{\mathrm{x}}}_{\mathrm{t}}\right){\mathrm{Q}}_{\mathrm{t}}^{-1}{\left({\mathrm{x}}_{\mathrm{t}}-{\widehat{\mathrm{x}}}_{\mathrm{t}}\right)}^{\mathrm{T}}.$$

(104)

It is a quadratic and strictly convex function; the problem is solved using PSO, which can be applied to complex problems such as non-linear and non-differentiable problems.

It is now assumed that the current particles of PSO are given as follows:

$${\mathrm{C}}_{\mathrm{t}}^{\mathrm{G}}=\left\{{\mathrm{x}}_{\mathrm{t},1}^{\mathrm{G}},\dots ,{\mathrm{x}}_{\mathrm{t},{\mathrm{N}}_{\mathrm{p}}}^{\mathrm{G}}\right\}$$

(105)

where ${\mathrm{C}}_{\mathrm{t}}^{\mathrm{G}}$ is an ${\mathrm{N}}_{\mathrm{p}}$ dimensional vector, and ${\mathrm{x}}_{\mathrm{t},\mathrm{i}}^{\mathrm{G}}$ represents each candidate solution $\mathrm{i}$ to the optimization problem at iteration $\mathrm{G}$ at time $\mathrm{t}$. In standard PSO, the local information of an objective function obtained, for example, from the gradient, is not exploited. This means that even if a search point of PSO around the optimal value exists, the particle could go past the optimum point, since local information is not considered; nevertheless, the searching ability of PSO can be improved when gradient information is considered, but with the risk of being trapped in a local minimum and missing the global solution of the problem. In order to overcome this problem, authors combine PSO with the Simultaneous Perturbation (SP) method and gradient. The SP is a method that is given as follows:

$$\begin{array}{l}{\mathrm{x}}_{\mathrm{t}+1}={\mathrm{x}}_{\mathrm{t}}+\mathrm{a}\mathsf{\Delta }{\mathrm{g}}_{\mathrm{t}}\\ \mathsf{\Delta }{\mathrm{g}}_{\mathrm{t}}^{\mathrm{j}}={\displaystyle \frac{\mathrm{f}\left({\mathrm{x}}_{\mathrm{t}}+{\mathrm{c}}_{\mathrm{t}}\right)-\mathrm{f}({\mathrm{x}}_{\mathrm{t}})}{{\mathrm{c}}_{\mathrm{t}}^{\mathrm{j}}}} \end{array}$$

(106)

where $\mathrm{a}$ is a positive constant, ${\mathrm{c}}_{\mathrm{t}}$ and ${\mathrm{c}}_{\mathrm{t}}^{\mathrm{j}}$ are a perturbation vector and its $\mathrm{j}$-th element, which is determined randomly, and $\mathsf{\Delta }{\mathrm{g}}_{\mathrm{t}}^{\mathrm{j}}$ represents the $\mathrm{j}$-th element of $\mathsf{\Delta }{\mathrm{g}}_{\mathrm{t}}$, which estimates the gradient.

As proposed in reference [66], the combination of PSO and SP gives rise to the following equations:

$$\begin{array}{l}{\mathrm{x}}_{\mathrm{t}+1}={\mathrm{x}}_{\mathrm{t}}+\mathsf{\Delta }{\mathrm{x}}_{\mathrm{t}}\\ \mathsf{\Delta }{\mathrm{x}}_{\mathrm{t}+1}=\mathsf{\chi }(-\mathrm{a}\mathsf{\Delta }{\mathrm{g}}_{\mathrm{t}}+{\mathsf{\varphi }}_1\left({\mathrm{p}}_{\mathrm{t}}-{\mathrm{x}}_{\mathrm{t}}\right)+{\mathsf{\varphi }}_2\left({\mathrm{n}}_{\mathrm{t}}-{\mathrm{x}}_{\mathrm{t}}\right))\end{array}$$

(107)

where ${\mathrm{x}}_{\mathrm{t}}$ is the particle, ${\mathrm{p}}_{\mathrm{t}}$ is the best value of the particle, ${\mathrm{n}}_{\mathrm{t}}$ is the best value of the swarm, ${\mathsf{\varphi }}_1$ and ${\mathsf{\varphi }}_2$ are two positive numbers in a predefined range, and $\mathsf{\chi }$ is a gain coefficient. To provide global and local search capabilities at the same time, and keeping in the meantime population diversity, the authors define a combination of PSO and SP, given as follows:

$${\mathrm{x}}_{\mathrm{t}}^{\mathrm{i}}={\mathrm{x}}_{\mathrm{t}-1}^{\mathrm{i}}+{\mathrm{v}}_{\mathrm{t}}^{\mathrm{i}}$$

(108)

$${\mathrm{v}}_{\mathrm{t}}^{\mathrm{i}}=\left\{\begin{array}{ll}{\mathrm{v}}_{\mathrm{m}\mathrm{a}\mathrm{x}}& {\mathrm{v}}_{\mathrm{t}}^{\mathrm{i}}>{\mathrm{v}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\\ \mathrm{w}(-\mathsf{\gamma }\nabla \mathsf{\psi }\left({\mathrm{x}}_{\mathrm{t}-1}^{\mathrm{i}}\right))+{\mathrm{c}}_1{\mathrm{r}}_1\left({\mathrm{P}}_{\mathrm{p}\mathrm{e}\mathrm{s}\mathrm{t}}-{\mathrm{x}}_{\mathrm{t}-1}^{\mathrm{i}}\right)+{\mathrm{c}}_2{\mathrm{r}}_2({\mathrm{P}}_{\mathrm{g}\mathrm{e}\mathrm{s}\mathrm{t}}-{\mathrm{x}}_{\mathrm{t}-1}^{\mathrm{i}})& -{\mathrm{v}}_{\mathrm{m}\mathrm{a}\mathrm{x}}<{\mathrm{v}}_{\mathrm{t}}^{\mathrm{i}}<{\mathrm{v}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\\ {-\mathrm{v}}_{\mathrm{m}\mathrm{a}\mathrm{x}}& {\mathrm{v}}_{\mathrm{t}}^{\mathrm{i}}<-{\mathrm{v}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}\right.$$

(109)

where ${\mathrm{c}}_1=2$, ${\mathrm{c}}_2=2$ are positive coefficients, ${\mathrm{r}}_1$ and ${\mathrm{r}}_2$ are random numbers in the interval (0, 1), ${\mathrm{v}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum velocity limited to 10% of the dynamic range of the variables on each dimension, $\mathrm{w}$ is the inertial weight, $\mathsf{\gamma }$ is a positive constant, and $\nabla \mathsf{\psi }$ is the gradient of the objective function.

3.9. Differential Evolution and Particle Swarm Optimization: DEPSO

The authors in reference [67] propose a hybrid algorithm that combines the DE algorithm and PSO and is composed of two main steps. In the former, mutation and selection operators are used to produce a new population, while in the latter, PSO is adopted to explore local optimums, followed by crossover and selection operations. In this way, it is possible to increase and decrease the extent of searching regions for the population.

The advantages are that the algorithm overcomes degradation defects, enhances the effectiveness of the particles, and improves accuracy and robustness.

Hybrid Algorithm DEPSO

In reference [67], authors use DEPSO in a Particle Filter (PF) algorithm to reduce the particle number and the computational time complexity. The DE algorithm aims to create a new candidate, and to select the best ones among the current individuals and the candidate for the next generation. In this work, the DE steps are developed as follows:

• Mutation. Three individuals from the population are randomly selected and the $\mathrm{i}$-the mutant vector ${\mathrm{v}}_{\mathrm{k}}^{\mathrm{i}}$ is generated with this equation:

  $${\mathrm{v}}_{\mathrm{k}}^{\mathrm{i}}\left(\mathrm{g}+1\right)=\mathrm{F}\left({\mathrm{x}}_{\mathrm{k}}^{{\mathrm{r}}_2}\left(\mathrm{g}\right)-{\mathrm{x}}_{\mathrm{k}}^{{\mathrm{r}}_3}\left(\mathrm{g}\right)\right)+{\mathrm{x}}_{\mathrm{k}}^{{\mathrm{r}}_1}\left(\mathrm{g}\right)$$

  (110)

  where $\mathrm{F}$ is the mutation parameter, ${\mathrm{r}}_1$, ${\mathrm{r}}_2$, and ${\mathrm{r}}_3$ identify the three selected individuals among $\mathrm{N}$, $\mathrm{k}$ identifies the $\mathrm{k}$-th step of the PF, and $\mathrm{g}$ identifies the generation of the DE.
• Crossover. This step performs the crossover operator to generate a candidate for the trial vector and the current individual. The $\mathrm{j}$-th dimension of the candidate is calculated by the following:

  $${\mathrm{u}}_{\mathrm{k}}^{\mathrm{i},\mathrm{j}}\left(\mathrm{g}+1\right)=\left\{\begin{array}{ll}{\mathrm{v}}_{\mathrm{k}}^{\mathrm{i},\mathrm{j}}\left(\mathrm{g}+1\right)& \mathrm{i}\mathrm{f}\left((\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(\mathrm{0,1})\le \mathrm{C}\mathrm{R}\right)\mathrm{o}\mathrm{r}(\mathrm{j}={\mathrm{j}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}))\\ {\mathrm{x}}_{\mathrm{k}}^{\mathrm{i},\mathrm{j}}\left(\mathrm{g}\right)& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.$$

  (111)

  where ${\mathrm{j}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}$ is a random integer in the range of $\left[1,\mathrm{n}\right],$ where $\mathrm{n}$ is the number of dimensions of the decision variables and $\mathrm{C}\mathrm{R}$ is the crossover parameter.
• Selection. In this step, the fitness function of ${\mathrm{u}}_{\mathrm{k}}^{\mathrm{i}}\left(\mathrm{g}+1\right)$ and ${\mathrm{x}}_{\mathrm{k}}^{\mathrm{i}}$ are compared and the best current individual is selected.

  $${\mathrm{x}}_{\mathrm{k}}^{\mathrm{i}}\left(\mathrm{g}+1\right)=\left\{\begin{array}{ll}{\mathrm{u}}_{\mathrm{k}}^{\mathrm{i}}\left(\mathrm{g}+1\right)& \mathrm{i}\mathrm{f}(\mathrm{f}({\mathrm{u}}_{\mathrm{k}}^{\mathrm{i}}\left(\mathrm{g}+1\right))\ge \mathrm{f}({\mathrm{x}}_{\mathrm{k}}^{\mathrm{i}}\left(\mathrm{g}\right)))\\ {\mathrm{x}}_{\mathrm{k}}^{\mathrm{i},\mathrm{j}}\left(\mathrm{g}\right)& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.$$

  (112)
• End. The end condition is reached when a maximum number of iterations are executed or when the population’s average fitness meets requirements. Otherwise, the algorithm is repeated from step 1.

Considering the PSO algorithm, each particle of the swarm has a position $\mathrm{X}$ and a velocity $\mathrm{V}$; it is recalled that their evolution is classically based on the following formulas:

$$\begin{array}{l}{\mathrm{V}}_{\mathrm{i}}=\mathrm{w}·{\mathrm{V}}_{\mathrm{i}}+{\mathrm{c}}_1·\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}·\left(\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}{\mathrm{t}}_{\mathrm{i}}-{\mathrm{X}}_{\mathrm{i}}\right)+{\mathrm{c}}_2·\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}·\left(\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}{\mathrm{t}}_{\mathrm{i}}-{\mathrm{X}}_{\mathrm{i}}\right) \\ {\mathrm{X}}_{\mathrm{i}}={\mathrm{X}}_{\mathrm{i}}+{\mathrm{V}}_{\mathrm{i}}\end{array}$$

(113)

where $\mathrm{w}$ is the inertia weight, $\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}{\mathrm{t}}_{\mathrm{i}}$ is the best position of particle ${\mathrm{X}}_{\mathrm{i}}$ so far; $\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}{\mathrm{t}}_{\mathrm{i}}$ is the best position of the swarm so far; rand is a random number between (0,1); and ${\mathrm{c}}_1$ and ${\mathrm{c}}_2$ are acceleration constants.

The DE algorithm and PSO have both advantages and disadvantages; for example, the differential algorithms have (i) a good global search capability due to the mutation step, which increases the diversity of the population, (ii) good local search capability, which is improved by the crossover step, and (iii) a memory function due to the selection process; however, the convergence speed is slow. On the other hand, PSO is characterized by a faster convergence rate, but it has the drawback that it might be stuck in a local optimum without converging to the global one because the diversity of the population can be lost very easily. Then, to exploit the characteristics of PSO, i.e., the high rate of convergence and the possibility to explore a local minimum quickly, Equation (113) is simplified as follows:

$$\begin{array}{l}{\mathrm{V}}_{\mathrm{i}}=\mathrm{w}·{\mathrm{V}}_{\mathrm{i}}+{\mathrm{c}}_1·\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}·\left(\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}{\mathrm{t}}_{\mathrm{i}}-{\mathrm{X}}_{\mathrm{i}}\right)\\ {\mathrm{X}}_{\mathrm{i}}={\mathrm{X}}_{\mathrm{i}}+{\mathrm{V}}_{\mathrm{i}}\end{array}$$

(114)

After the aforementioned considerations, the hybrid algorithm DEPSO can be represented with the following steps:

• Set $\mathrm{t}=0$, and randomly initialize all individual positions ${\mathrm{X}}_{\mathrm{i}}$ and velocities ${\mathrm{V}}_{\mathrm{i}}$.
• $\mathrm{t}=\mathrm{t}+1$.
• Calculate individual fitness function ${\mathrm{f}(\mathrm{X}}_{\mathrm{i}})$.
• Update the position of particles using the equation for mutation (110) and selection (112).
• Select the individual in the best position of the swarm so far.
• Update the position and velocity of particles according to (114).
• Perform crossover and selection using (111) and (112), respectively.
• If the counter is larger than a threshold, then exit; otherwise, return to step 2.

Authors claim that DEPSO avoids particle degradation, improves positioning accuracy, reduces the number of particles needed for positioning, and enhances the diversity of particles and robustness of the system.

3.10. Bat Algorithm

The bat algorithm (BA) idea is inspired by the echolocation behavior of bats to find food; in nature, they usually emit short pulses, and when they encounter food, they change the frequency to tune the echolocation time and increase the location accuracy of the prey [68]. The parameters of loudness and frequency determine the global search capability, and if a bat is near its scope, the impulse rate increases, whereas loudness decreases. In the standard bat algorithm, each individual $\mathrm{i}$ has a defined position ${\mathrm{x}}_{\mathrm{i}}(\mathrm{t})$ and velocity ${\mathrm{v}}_{\mathrm{i}}(\mathrm{t})$ in the search space, which are updated iteration by iteration with the following set of equations:

$${\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}+1\right)={\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}\right)+{\mathrm{v}}_{\mathrm{i}}(\mathrm{t}+1)$$

(115)

$${\mathrm{v}}_{\mathrm{i}}\left(\mathrm{t}+1\right)={\mathrm{v}}_{\mathrm{i}}\left(\mathrm{t}\right)+{(\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}\right)-\mathrm{p}\left(\mathrm{t}\right){)·\mathrm{f}}_{\mathrm{i}}$$

(116)

$${\mathrm{f}}_{\mathrm{i}}={\mathrm{f}}_{\mathrm{m}\mathrm{i}\mathrm{n}}+({\mathrm{f}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{f}}_{\mathrm{m}\mathrm{i}\mathrm{n}})·\mathsf{\beta }$$

(117)

where $\mathrm{t}$ is the $\mathrm{t}$-th iteration, $\mathrm{p}\left(\mathrm{t}\right)$ is the current global optimal solution, ${\mathrm{f}}_{\mathrm{i}}$ is the frequency of the $\mathrm{i}$-th individual, ${\mathrm{f}}_{\mathrm{m}\mathrm{i}\mathrm{n}}=0$, ${\mathrm{f}}_{\mathrm{m}\mathrm{a}\mathrm{x}}=1$, and $\mathsf{\beta }$ is a random vector in the range [0, 1] with uniform distribution.

To implement the local search strategy, the following equation is considered:

$${\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}+1\right)={\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}\right)+\mathsf{\epsilon }\overline{\mathrm{A}}(\mathrm{t})$$

(118)

where $\mathsf{\epsilon }$ is a random number in the range [−1, 1] and $\overline{\mathrm{A}}(\mathrm{t})$ is the average loudness of the population.

The global search is controlled by loudness ${\mathrm{A}}_{\mathrm{i}}\left(\mathrm{t}+1\right)$ and pulse rate ${\mathrm{r}}_{\mathrm{i}}\left(\mathrm{t}+1\right)$, which are updated with the following equations:

$${\mathrm{A}}_{\mathrm{i}}\left(\mathrm{t}+1\right)=\mathsf{\alpha }{\mathrm{A}}_{\mathrm{i}}(\mathrm{t})$$

(119)

$${\mathrm{r}}_{\mathrm{i}}\left(\mathrm{t}+1\right)={\mathrm{r}}_{\mathrm{i}}\left(\mathrm{t}\right)+(1-\exp (-\mathsf{\gamma }\mathrm{t}))$$

(120)

where $0<\mathsf{\alpha }<1$ is the amplitude decrease rate and is a constant term, and $\mathsf{\gamma }>0$ is a constant. The initial values of loudness ${\mathrm{A}}_{\mathrm{i}}\left(0\right)$ and pulse rate ${\mathrm{r}}_{\mathrm{i}}\left(0\right)$ are set in the initialization.

The steps of the standard algorithm are as follows:

• Initialize, for each bat, the position and velocity parameters and randomly generate the frequency with Equation (117).
• Update the position and velocity of each bat with (115) and (116).
• For each bat, generate a random number $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}1$ in the range [0, 1];
  if $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}1<{\mathrm{r}}_{\mathrm{i}}(\mathrm{t})$, then update the temp position with (118) and calculate the fitness value for the corresponding $\mathrm{i}$-th bat.
• For each bat, generate a random number $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}2$ in the range [0, 1];
  if $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}2<{\mathrm{A}}_{\mathrm{i}}(\mathrm{t})$ and if $\mathrm{F}({\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}\right))<\mathrm{F}({\mathrm{p}}_{\mathrm{i}}\left(\mathrm{t}\right))$, where $\mathrm{F}$ is the fitness function, then update the loudness with (119) and the pulse rate with (120).
• Sort each individual based on fitness values and save the best position.
• If the final condition is met, then exit; otherwise, go to step 2.

Mobile Robot Localization Based on the Leader-Based Bat Algorithm

The authors in reference [69] propose a modified BA to solve the mobile robot global localization problem, which is often characterized by complex multimodal beliefs. In some problems, the standard BA might provide a premature convergence solution. Authors propose the leader-based bat algorithm (LBBA), introducing leaders who can influence the behavior of other individuals. This gives rise to more stable particle movements at multiple local optimum points, providing a convergence to the global optimum in more complex environments.

The LBBA is a nonparametric filter with $\mathrm{N}$ individuals, and each of them approximately corresponds to one region in the search space. The a priori states of individuals are represented by a set of random samples, and this enables the simultaneous representation of multiple optimal locations.

Unlike the standard BA, the LBBA uses both the best global ${\mathrm{g}}^{\mathrm{*}}$ and the best individual ${\mathrm{x}}^{\mathrm{*}}$. The last one enables a higher diversity of solutions in the population, and, in addition, the performance and efficiency in the search for the global optimum are improved thanks to the adoption of leaders, who promote the creation of different colonies, decreasing the randomness of flight.

To determine the target position, the following object function is defined:

$$\mathrm{O}\left({\mathrm{x}}_{\mathrm{i}}\right)={\displaystyle \frac{\left({\mathrm{m}}_{\mathrm{i}}-\mathrm{m}\right)}{\mathrm{m}}}$$

(121)

where ${\mathrm{m}}_{\mathrm{i}}$ is the measurement performed by the simulated sensor of the $\mathrm{i}$-th individual, and $\mathrm{m}$ is the measurement performed by the real robot sensor.

The main steps of a general LBBA are as follows:

• Step 1. Initialization of ${\mathrm{x}}_{\mathrm{i}}$, ${\mathrm{v}}_{\mathrm{i}}$ where $\mathrm{i}=1,\dots ,\mathrm{N}$. Define the number of leaders $\mathrm{L}=1,\dots ,{\mathrm{L}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$.
  Define the number of leaders $\mathrm{L}=1,\dots ,{\mathrm{L}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$.
  Initialization of frequency ${\mathrm{f}}_{\mathrm{i}}$, impulse rate ${\mathrm{r}}_{\mathrm{i}}$, and amplitude ${\mathrm{A}}_{\mathrm{i}}$.
  Assess the fitness of each bat using the fitness function $\mathrm{f}\mathrm{i}\mathrm{t}({\mathrm{x}}_{\mathrm{i}})$.
  Determine leaders ${\mathrm{x}}_{\mathrm{L}}^{\mathrm{*}}$ with the fitness function.
  Define the best global bat and the relative best fitness: ${\mathrm{x}}^{\mathrm{*}}$ and $\mathrm{f}\mathrm{i}\mathrm{t}({\mathrm{x}}^{\mathrm{*}})$.
• Step 2. Compute the following:

  -

  ${\mathrm{f}}_{\mathrm{i}}={\mathrm{f}}_{\mathrm{m}\mathrm{i}\mathrm{n}}+({\mathrm{f}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{f}}_{\mathrm{m}\mathrm{i}\mathrm{n}})·\mathsf{\beta }$;

  -

  ${\mathrm{v}}_{\mathrm{i}}\left(\mathrm{t}+1\right)={\mathrm{v}}_{\mathrm{i}}\left(\mathrm{t}\right)+{({\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}\right)-\mathrm{x}}_{\mathrm{L}}^{\mathrm{*}}{)·\mathrm{f}}_{\mathrm{i}}$;

  -

  ${\widehat{\mathrm{x}}}_{\mathrm{i}}={\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}\right)+{\mathrm{v}}_{\mathrm{i}}\left(\mathrm{t}+1\right)$, where ${\widehat{\mathrm{x}}}_{\mathrm{i}}$ is a temporary variable.

• Step 3. If $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}1>{\mathrm{r}}_{\mathrm{i}}$ then ${\widehat{\mathrm{x}}}_{\mathrm{i}}={\mathrm{x}}^{\mathrm{*}}+\mathsf{\delta }$, where ${\mathrm{x}}^{\mathrm{*}}$ is the best global variable, and $\mathsf{\delta }$ is a random displacement variable.
• Step 4. Estimate the fitness of ${\widehat{\mathrm{x}}}_{\mathrm{i}}$ using $\mathrm{f}\mathrm{i}\mathrm{t}\left({\widehat{\mathrm{x}}}_{\mathrm{i}}\right)$.
• Step 5. Define new leaders ${\widehat{\mathrm{x}}}_{\mathrm{L}}^{\mathrm{*}}$.
• Step 6. If $\mathrm{f}\mathrm{i}\mathrm{t}\left({\widehat{\mathrm{x}}}_{\mathrm{i}}\right)\le \mathrm{f}\mathrm{i}\mathrm{t}\left({\mathrm{x}}_{\mathrm{i}}\right)$ and $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}2<{\mathrm{A}}_{\mathrm{i}}$, then compute the following:

  -

  ${\mathrm{x}}_{\mathrm{i}}={\widehat{\mathrm{x}}}_{\mathrm{i}}$;

  -

  $\mathrm{f}\mathrm{i}\mathrm{t}\left({\mathrm{x}}_{\mathrm{i}}\right)=\mathrm{f}\mathrm{i}\mathrm{t}({\widehat{\mathrm{x}}}_{\mathrm{i}})$;

  -

  ${\mathrm{r}}_{\mathrm{i}}=\mathrm{r}(0)\left(1-{\mathrm{e}}^{-\mathsf{\gamma }\mathrm{t}}\right)$;

  -

  ${\mathrm{A}}_{\mathrm{i}}\left(\mathrm{t}+1\right)=\mathsf{\alpha }{\mathrm{A}}_{\mathrm{i}}(\mathrm{t})$.

• Step 7. If $\mathrm{f}\mathrm{i}\mathrm{t}\left({\widehat{\mathrm{x}}}_{\mathrm{i}}\right)\le \mathrm{f}\mathrm{i}\mathrm{t}\left({\mathrm{x}}^{\mathrm{*}}\right),$ then compute ${\mathrm{x}}^{\mathrm{*}}={\widehat{\mathrm{x}}}_{\mathrm{i}}$, and $\mathrm{f}\mathrm{i}\mathrm{t}({\mathrm{x}}^{\mathrm{*}})={\mathrm{f}\mathrm{i}\mathrm{t}(\widehat{\mathrm{x}}}_{\mathrm{i}})$.
• Step 8. If the condition is met, return ${\mathrm{x}}^{\mathrm{*}}$; otherwise, repeat from step 2.

Authors claim that the LBBA has better performances when compared to other algorithms for robot localization such as PF and PSO, showing the ability to escape from local minima in highly complex environments. Moreover, it also succeeds in solving the robot’s kidnapping problem without premature convergences.

3.11. Artificial Bee Colony Optimization Algorithm

The Artificial Bee Colony (ABC) algorithm is a metaheuristic search algorithm inspired by the intelligent foraging behavior of honeybees in nature [70]. The algorithm is explained as follows [71]. The colony of artificial bees is composed of three groups: employed bees, onlookers, and scouts. The number of employed bees and onlooker bees is equal to the number of food sources. Every bee can only exploit one food source at a time. The scout bee was an employed bee that became a scout bee when it left the food source looking for a new one. A possible solution to the problem corresponds to a food source and the fitness of a solution corresponds to the nectar amount of a source. Each employed bee first explores the associated solution (or food source) according to their memories and explores possible solutions in corresponding neighborhoods. Then, they fly back to the hive and share the information about the food sources with the onlooker bees. The onlooker bees select the food sources, and a source is chosen using a probability criterion based on the value of fitness or amount of nectar: the food sources with more nectar have a higher probability of being selected. If a food source does not provide any better solution after several times that it is explored, the employed bee of that source becomes a scout bee and randomly selects a new source as it continues the work as an employed bee. The main steps of the algorithm are as follows [71]:

• Initialization. A solution is considered a multi-dimensional vector ${\mathrm{x}}_{\mathrm{i}}$ and $\mathrm{S}\mathrm{N}$ solutions are randomly generated at the beginning; the components of the vector are generated by the following equation:

  $${\mathrm{x}}_{\mathrm{i}\mathrm{j}}={\mathrm{x}}_{\mathrm{j}}^{\mathrm{m}\mathrm{i}\mathrm{n}}+\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(\mathrm{0,1})\times ({\mathrm{x}}_{\mathrm{j}}^{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{x}}_{\mathrm{j}}^{\mathrm{m}\mathrm{i}\mathrm{n}})$$

  (122)

  where $\mathrm{i}$ is one of the $\mathrm{S}\mathrm{N}$ solutions, $\mathrm{j}$ is one of the $\mathrm{D}$ dimensions of the vector, the ${\mathrm{x}}_{\mathrm{j}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ and ${\mathrm{x}}_{\mathrm{j}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$ represent the lower and upper bound of the $\mathrm{j}$-th dimension, respectively, and $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(\mathrm{0,1})$ provides a random number in the range [0, 1].
• Employed bee phase. Each employed bee is associated with a solution ${\mathrm{x}}_{\mathrm{i}}$ and then new solutions, ${\mathrm{v}}_{\mathrm{i}}$, are explored (or new food sources are generated) using the following equation:

  $${\mathrm{v}}_{\mathrm{i}\mathrm{j}}={\mathrm{x}}_{\mathrm{i}\mathrm{j}}+\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(-\mathrm{1,1})\times ({\mathrm{x}}_{\mathrm{i}\mathrm{j}}-{\mathrm{x}}_{\mathrm{k}\mathrm{j}})$$

  (123)

  where ${\mathrm{v}}_{\mathrm{i}\mathrm{j}}$ is the $\mathrm{j}$-th component of a new solution, $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(-\mathrm{1,1})$ provides a random number in the range [−1,1], ${\mathrm{x}}_{\mathrm{k}\mathrm{j}}$ is the j-th component of another food source ${\mathrm{x}}_{\mathrm{k}}$ that is selected randomly from the remaining population, and $\mathrm{k}\ne \mathrm{i}$. The old solution, ${\mathrm{x}}_{\mathrm{i}\mathrm{j}}$, and the new one, ${\mathrm{v}}_{\mathrm{i}\mathrm{j}}$, are compared and the employed bee brings back the information of the better one, i.e., the solution with a smaller objective value or the best fitness. It is called the greedy selection strategy. The fitness values, ${\mathrm{f}\mathrm{i}\mathrm{t}}_{\mathrm{i}}$, are calculated with the following equation:

  $${\mathrm{f}\mathrm{i}\mathrm{t}}_{\mathrm{i}}=\left\{\begin{array}{ll}{\displaystyle \frac{1}{1+{\mathrm{f}}_{\mathrm{i}}}}& {\mathrm{f}}_{\mathrm{i}}\ge 0\\ 1+\mathrm{a}\mathrm{b}\mathrm{s}({\mathrm{f}}_{\mathrm{i}})& {\mathrm{f}}_{\mathrm{i}}<0\end{array}\right.$$

  (124)

  where ${\mathrm{f}}_{\mathrm{i}}$ is the output of the objective function calculated using the ${\mathrm{x}}_{\mathrm{i}}$ solution.
• Onlooker bee phase. The employed bees share the fitness information with the onlooker bees. Each onlooker bee chooses a solution for itself using a probability criterion: the higher the probability of a solution, the higher the possibility of being chosen. The selection probability of the $\mathrm{i}$-th solution is called the roulette wheel selection and the probability ${\mathrm{p}}_{\mathrm{i}}$ is calculated as follows:

  $${\mathrm{p}}_{\mathrm{i}}={\displaystyle \frac{\mathrm{f}\mathrm{i}{\mathrm{t}}_{\mathrm{i}}}{\sum _{\mathrm{m}=1}^{\mathrm{S}\mathrm{N}}\mathrm{f}\mathrm{i}{\mathrm{t}}_{\mathrm{m}}}}.$$

  (125)

  After the selection, each onlooker bee generates a new solution in the surroundings of the chosen solution by using (123). Each generated solution is compared with the corresponding previous food source via the greedy selection strategy, and if the new one is better, it is set as ${{\mathrm{x}}_{\mathrm{i}\mathrm{j}}=\mathrm{v}}_{\mathrm{i}\mathrm{j}}$. The difference between the updating process of the employed bees and the onlooker bees is that the employed bees update all solutions, while the onlooker bees only update the selected one.
• Scout bee phase. A solution may be explored several times without providing a better solution in the surroundings. After a pre-set number of explorations in the surroundings of the same solution without better results, the solution is abandoned, and then the corresponding employed bee transforms into a scout bee to randomly generate a new solution by using (122). Then, it continues as an employed bee.

Multistage Localization in Wireless Sensor Networks Using the Artificial Bee Colony Algorithm

The authors in reference [72] propose a method based on the ABC algorithm to solve the sensor localization problem in a wireless sensor network, providing some simulation results. In the considered scenario, $\mathrm{N}$ nodes are deployed in a square plan, and each of them has a communication range of $\mathrm{r}$ units. A few nodes of them are considered as anchor nodes and they are in a known position; the other nodes are deployed in random locations. A node in an unknown position with at least three non-collinear anchors in its range can estimate the distances from anchors by parameters such as RSSI, Time of Arrival (ToA), and so on; if at least three anchor nodes are not available, the sensor cannot be localized. Each localizable node executes the ABC algorithm to estimate the coordinates of its location; once a sensor in an unknown position has been localized, it is used as a new anchor node to identify the other sensors. The aim is to determine the locations of all or as many nodes as possible, and in such a way, the distributed localization problem is modelled as a two-dimensional optimization problem.

The ABC algorithm computes the estimated location $({\widehat{\mathrm{u}}}_{\mathrm{j}\mathrm{x}},{\widehat{\mathrm{u}}}_{\mathrm{j}\mathrm{y}})$, finding the minimum of the mean localization error ${\mathrm{E}}_{\mathrm{j}}$, which is expressed as follows:

$${\mathrm{E}}_{\mathrm{j}}={\displaystyle \frac{1}{3}}\sum _{\mathrm{i}=1}^3\left(\sqrt{{\left({\widehat{\mathrm{u}}}_{\mathrm{j}\mathrm{x}}-{\mathrm{a}}_{\mathrm{i}\mathrm{x}}\right)}^2+{\left({\widehat{\mathrm{u}}}_{\mathrm{j}\mathrm{y}}-{\mathrm{a}}_{\mathrm{i}\mathrm{y}}\right)}^2}-{\widehat{\mathrm{d}}}_{\mathrm{i}}\right)$$

(126)

where $\mathrm{j}$ is the $\mathrm{j}$-th node in an unknown position, $\mathrm{i}$ is the $\mathrm{i}$-th anchor node with known coordinates ${(\mathrm{a}}_{\mathrm{i}\mathrm{x}},{\mathrm{a}}_{\mathrm{i}\mathrm{y}})$, ${\widehat{\mathrm{d}}}_{\mathrm{i}}$ is the distance measured (that includes noise and errors), and ${(\widehat{\mathrm{u}}}_{\mathrm{j}\mathrm{x}},{\widehat{\mathrm{u}}}_{\mathrm{j}\mathrm{y}})$ are the unknown coordinates.

The authors compared the performances of the ABC and PSO algorithms in finding the solution to the specific problem. The ABC algorithm can result in more accurate localization than PSO, but requires a higher computational time.

3.12. Simulated Annealing

The SA algorithm is a random-search technique that is based on the analogy between the process of metal crystalline formation from the liquid phase, and the process of searching for an optimal minimum in a problem. If the cooling process is too rapid, then the metal will solidify into a sub-optimal configuration. If the cooling process takes the proper time, the crystals within the material will solidify optimally into a state of minimum energy, which is the ground state associated with the optimal global solution to the problem. This is a heuristic method that does not need any particular assumptions and can be used for non-linear and non-convex problems; in particular, it results in a suitable technique for solving combinatorial optimization problems.

The main steps of the SA algorithm are as follows [73]:

• Generate a random initial solution vector ${\mathrm{x}=\mathrm{x}}_{\mathrm{o}}$ for the system and then evaluate it using the objective function $\mathrm{f}(\mathrm{x})$.
• Initialize the temperature $\mathrm{T}={\mathrm{T}}_0$.
• Perturbate the previous solution to obtain a neighboring ${\mathrm{x}}_{\mathrm{i}}+\Delta \mathrm{x}$. Evaluate the new solution.
• A new solution $({\mathrm{x}}_{\mathrm{i}}+\Delta \mathrm{x})$ is accepted as a new solution ${{\mathrm{x}}_{\mathrm{i}+1}=\mathrm{x}}_{\mathrm{i}}+\Delta \mathrm{x}$ if ${\mathrm{f}(\mathrm{x}}_{\mathrm{i}}+\Delta \mathrm{x})$ is better than ${\mathrm{f}(\mathrm{x}}_{\mathrm{i}})$.
• Otherwise, accept the new solution ${\mathrm{x}}_{\mathrm{i}}+\Delta \mathrm{x}$ with a probability ${\mathrm{e}}^{-\left({\displaystyle \raisebox{1ex}{$\Delta \mathrm{f}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{T}}_{\mathrm{i}}$}\right.}\right)}$,
  where $\Delta \mathrm{f}={\mathrm{f}(\mathrm{x}}_{\mathrm{i}}+\Delta \mathrm{x}){-\mathrm{f}(\mathrm{x}}_{\mathrm{i}})$.
• Reduce the system temperature according to the cooling schedule.
• Repeat from step 3 until a stopping criterion is met (number of cycles or ${\mathrm{T}}_{\mathrm{i}}=0)$.

In the SA problem, temperature is a parameter that controls the probability in step 5. During the iterations, the lower the temperature, the lower the probability of accepting a bad solution. This mechanism helps the procedure to avoid being stuck in a global minimum at the beginning of the iterations, because it allows the algorithm to jump out of the minima looking for better solutions. To obtain a good result, it is important to properly set the initial value of temperature $\mathrm{T}$. If the value is too high, then it takes more cycles to converge, and vice versa; if the value is too small, the algorithm may provide a local optimum.

In step 6, the choice of an appropriate cooling scheme is important for the algorithm’s success. In the literature, there are several types of schemes, and they can be classified as either monotonic or adaptive. The former consists of a static decrease in the temperature at every iteration independently of the goodness of the solutions and independently of the current explored neighborhood structure, whereas the last one considers the quality of transitions, adopting mechanisms for the temperature decrease that can adopt also reannealing techniques [73].

Simulated Annealing-Based Localization in Wireless Sensor Network

The authors in reference [74] propose a SA localization (SAL) method for wireless sensor networks. It aims to estimate the position of sensors by determining the position of a set of anchors and the measured distance between them and the sensors. However, several factors can influence distance measurements, such as synchronization, multipath, fading, Non-Line of Sight (NLoS) conditions, and other sources of interference. To estimate the real distance of the sensor with a minimum error, the authors propose to solve the optimization problem with the SAL method. The objective function is defined as follows:

$$\mathrm{C}\mathrm{F}=\sum _{\mathrm{i}=1}^{\mathrm{N}}\sum _{\mathrm{j}\in {\mathrm{N}}_{\mathrm{i}}}{\left({\widehat{\mathrm{d}}}_{\mathrm{i}\mathrm{j}}-{\mathrm{d}}_{\mathrm{i}\mathrm{j}}\right)}^2$$

(127)

where $\mathrm{N}$ is number of non-anchor nodes, ${\mathrm{N}}_{\mathrm{i}}$ is the set of neighbors of the $\mathrm{i}$-th node, ${\widehat{\mathrm{d}}}_{\mathrm{i}\mathrm{j}}$ and ${\mathrm{d}}_{\mathrm{i}\mathrm{j}}$ are the estimated distance and the measured distance of node $\mathrm{i}$ with its neighbor $\mathrm{j}$, respectively.

The adopted cooling scheme for temperature is ${\mathrm{T}}_{\mathrm{i}+1}={\mathsf{\alpha }\mathrm{*}\mathrm{T}}_{\mathrm{i}}$, where $\mathsf{\alpha }<1$ is the cooling rate that is determined empirically. The solutions are perturbated as $\Delta {\mathrm{d}}_{\mathrm{i}+1}=\mathsf{\beta }\mathrm{*}\Delta {\mathrm{d}}_{\mathrm{i}}$, where $\Delta \mathrm{d}$ is chosen to span the whole set of allowable moves and $\mathsf{\beta }<1$ is used to bias the generation of random moves at a lower $\mathrm{T}$.

3.13. Ant Colony Optimization

Ant Colony Optimization (ACO) is a probabilistic method, and it was initially proposed to solve the Travelling Salesman Problem (TSP); in general, ACO is suitable for solving problems in which the searched solution is the best path through a graph. It is inspired by the behavior of an ant seeking a path between its colony and food. The ants live in a colony, roam around it looking for food, and communicate with each other by using pheromones which can be released on the ground when insects are moving. Once an ant finds some food, the insect brings it as much as possible back to the colony and, during the travel, it releases pheromones in agreement with the quantity and quality of the found food. Other ants can smell it, and when they come out of the colony looking for food, they choose the path with higher pheromones. The higher the pheromone along a route, the higher the probability that an ant chooses that path. In such a way, more ants follow the path, more pheromones are deposited, and the path has a higher probability of being chosen. After some iterations, the ants will choose the shortest path, which is expected to be the optimum of the problem.

Mathematically, the algorithm is explained as follows [75]. ACO can be understood considering a simple graph, $\mathrm{G}=\left(\mathrm{V},\mathrm{E}\right)$, with two nodes, where $\mathrm{V}$ represents the nodes of the graph (the nest ${\mathrm{v}}_{\mathrm{s}}$ and the food source ${\mathrm{v}}_{\mathrm{d}}$), $\mathrm{E}$ represents two links, ${\mathrm{e}}_1$ and ${\mathrm{e}}_2$, between nodes with length ${\mathrm{l}}_1$ and ${\mathrm{l}}_2$, respectively, such that ${\mathrm{l}}_1<{\mathrm{l}}_2$. The artificial pheromone trails are modelled for each path by ${\mathsf{\tau }}_1$ and ${\mathsf{\tau }}_2$ and they indicate the strength of the pheromone trail on the corresponding path. Every ant chooses with a probability given as follows:

$${\mathrm{p}}_{\mathrm{i}}={\displaystyle \frac{{\mathsf{\tau }}_{\mathrm{i}}}{{\mathsf{\tau }}_1+{\mathsf{\tau }}_2}}$$

(128)

where $\mathrm{i}$ identifies the $\mathrm{i}$-th link; about the lengths, it is ${\mathsf{\tau }}_1>{\mathsf{\tau }}_2$, which means that link $1$ has a higher probability. As for returning, every ant uses the same path chosen in the outward journey, and it updates the artificial pheromone value associated with the used link with the following equation:

$${\mathsf{\tau }}_{\mathrm{i}}={\mathsf{\tau }}_{\mathrm{i}}{\displaystyle \frac{\mathrm{Q}}{{\mathrm{l}}_{\mathrm{i}}}}$$

(129)

where $\mathrm{Q}>0$ is a parameter of the model. It can be inferred that the shorter the path, the higher the amount of added pheromones. At each step (or iteration), all the ants are initially placed in node ${\mathrm{v}}_{\mathrm{s}}$ and each of them moves from the nest to the food source. The evaporation of the pheromone is simulated by the following:

$${\mathsf{\tau }}_{\mathrm{i}}=(1-\mathsf{\rho })$$

(130)

where $\mathsf{\rho }$ is in the range $\left(\mathrm{0,1}\right]$ and regulates the pheromone evaporation. Finally, all ants return to the nest and reinforce their chosen path, as reported above. Iteration by iteration, the shortest path will be the one with the highest pheromones, and then, with the highest probability.

For real combinatory optimization problems with more than two nodes, the solution is constructed as follows [76]. Each of the $\mathrm{m}$ ants is randomly placed on a chosen node and then at each node, a state transition rule is iteratively applied. The unvisited nodes, the pheromone trail strength, and the length between two nodes bias the ant’s choice. The partial tour is stored in a memory called the tabu list which needs to determine at each construction step the nodes to be visited, to guarantee that a feasible solution is built, and to retrace the trail of the related ant once it has been completed. When all ants end up constructing a path, the pheromones are updated using evaporation and new deposition; short paths are updated with a higher quantity of pheromones. The algorithm is summarized as follows [76]:

• Set parameters and initialize pheromone trails.
• Construct solutions.
• Apply local search.
• Update trails.
• If termination conditions are not met, repeat from step 2.

There are different versions of ACO and historically, the first presented version was called the Ant System (AS); some variants of the ACO algorithms (here not investigated) are the Elitist AS (EAS), Rank-based AS (RAS), MAX–MIN Ant System (MMAS), Ant Colony System (ACS), and Hyper-Cube Framework (HCF).

ANN Topology Optimized by ACO

The author in reference [77] deals with the problem of finding a continuous path for a robot from an initial position to a goal position where collisions against static objects must be avoided. The proposed system exploits an ANN that provides the optimum path, and it is trained by adopting the swarm intelligence-based reinforcement learning (SWIRL) method based on PSO and ACO. Because the training process of Multi-Layer Perceptron (MLP) for pattern classification problems can be considered as graph problems where the neurons are vertices and the connections are directed edges (links), PSO is exploited to provide the best ANN connection weights, whilst the ACO is used to optimize the topology of the ANN. In such a way, the SWIRL provides both an appropriate architecture for the problem and the trained connection weights of the network.

At each step of the iteration process, the PSO is applied to adjust ANN connection weights within a given topology structure. Focusing only on the ACO algorithm, it is used to allocate training iterations among a set of candidate network topologies. The desirability $\mathrm{d}(\mathrm{i})$ for the $\mathrm{i}$-th ANN is given as follows:

$$\mathrm{d}(\mathrm{i})={\displaystyle \frac{1}{\mathrm{h}+1}}$$

(131)

where $\mathrm{h}$ is the number of hidden nodes; each $\mathrm{i}$-th ANN represents a topology of a set of ANNs with different numbers of hidden nodes. The pheromone concentration ${\mathsf{\tau }}_{\mathrm{i}}$ is initialized to 0.1, and it is expressed as follows:

$$\mathsf{\tau }\left(\mathrm{i},\mathrm{t}+1\right)=\mathsf{\rho }\mathsf{\tau }\left(\mathrm{i},\mathrm{t}\right)+{\mathrm{n}}_{\mathrm{a}}(\mathrm{i}){\displaystyle \frac{\mathrm{g}(\mathrm{i})}{{\mathrm{g}}_{\mathrm{s}\mathrm{u}\mathrm{m}}}}$$

(132)

where $\mathsf{\rho }$ is the rate of evaporation, ${\mathrm{n}}_{\mathrm{a}}$ is the number of ants returning from the $\mathrm{i}$-th neural network, $\mathrm{g}(\mathrm{i})$ is the global best for $\mathrm{i}$, ${\mathrm{g}}_{\mathrm{s}\mathrm{u}\mathrm{m}}$ is the sum of all the current global bests, and $\mathrm{t}$ is the index iteration of the ACO. Each ant represents one training process for the PSO. During each ACO step, the ants go out into the topology space of ANNs, and they choose an ANN using a probability criteria, given as follows:

$$\mathrm{p}(\mathrm{i})={\displaystyle \frac{{\left[\mathrm{d}\left(\mathrm{i}\right)\right]}^{\mathrm{b}}{\left[\mathsf{\tau }\left(\mathrm{i},\mathrm{t}\right)\right]}^{\mathrm{a}}}{\sum \left\{{\left[\mathrm{d}\left(\mathrm{i}\right)\right]}^{\mathrm{b}}{\left[\mathsf{\tau }\left(\mathrm{i},\mathrm{t}\right)\right]}^{\mathrm{a}}\right\}}}$$

(133)

where $\mathrm{a}$ and $\mathrm{b}$ are constant factors which control the relative influence of pheromones and desirability, respectively [77].

3.14. Whale Optimization Algorithm

The Whale Optimization Algorithm (WOA) is a metaheuristic method inspired by the hunting behavior of humpback whales, and it has been proposed by authors in reference [78]. The algorithm mimics three behaviors: (i) encircling prey or siege predation, (ii) spiral bubble-net feeding maneuver, and (iii) searching for prey. In the first one, a whale is in a location looking for prey to encircle; at the beginning, the locations (or solutions) are random, and the current best solution, defined as the best search agent, is found by exploiting the fitness function; the other search agents or candidate solutions move iteration by iteration towards the current best one. Mathematically, the behavior of each whale can be modelled by the following:

$$\mathrm{D}=\left|\mathrm{C}·{\mathrm{X}}^{\mathrm{*}}\left(\mathrm{t}\right)-\mathrm{X}(\mathrm{t})\right|$$

(134)

$$\mathrm{X}\left(\mathrm{t}+1\right)={\mathrm{X}}^{\mathrm{*}}\left(\mathrm{t}\right)-\mathrm{A}·\mathrm{D}$$

(135)

where $\mathrm{t}$ is the $\mathrm{t}$-th iteration, $\mathrm{A}$ and $\mathrm{C}$ are coefficient vectors, ${\mathrm{X}}^{\mathrm{*}}$ is the position vector of the best solution obtained so far, $\mathrm{X}$ is a candidate position vector, and the dot ($·$) represents the element-by-element multiplication. Furthermore, it is possible to express $\mathrm{A}=2\mathrm{a}·\mathrm{r}-\mathrm{a}$ and $\mathrm{C}=2·\mathrm{r}$, where $\mathrm{a}$ is linearly decreased from 2 to 0 iteration by iteration, and $\mathrm{r}$ is a random vector in the range [0, 1] that allows it to reach any position of the design space. The used vectors in this explanation have the same dimension as the input vector of the problem; however, the algorithm can be applied to a problem with a dimension greater than three even if the algorithm is inspired by the analogy of a moving whale in a three-dimensional space.

In the second behavior (spiral bubble-net feeding manoeuvre), two approaches are adopted: the shrinking encircling mechanism and the spiral updating position. The former needs to push the $\mathrm{X}$ vectors towards ${\mathrm{X}}^{\mathrm{*}}$ and it consists in decreasing the value of $\mathrm{a}$ iteration by iteration; in such a way, considering that $\mathrm{r}$ is random, the result is that $\mathrm{A}$ is composed by random values in the range [$-\mathrm{a},\mathrm{a}$]. As for the second approach, to mimic an encircling movement, a spiral equation is created between the position of a whale and the best position; it is expressed as follows:

$$\mathrm{X}\left(\mathrm{t}+1\right)={\mathrm{D}}_{\mathrm{i}}·{\mathrm{e}}^{\mathrm{b}\mathrm{l}}·\cos \left(2\mathsf{\pi }\mathrm{l}\right)+{\mathrm{X}}^{\mathrm{*}}(\mathrm{t})$$

(136)

where ${\mathrm{D}}_{\mathrm{i}}=\left|{\mathrm{X}}^{\mathrm{*}}\left(\mathrm{t}\right)-\mathrm{X}(\mathrm{t})\right|$ indicates the distance of the $\mathrm{i}$-th candidate solution (or whale) to the best solution obtained so far, $\mathrm{b}$ is a constant for defining the shape of the exponential spiral, and $\mathrm{l}$ is a random number in the range [−1, 1]. Considering the analogy with whales, they are swimming around their candidate location prey within a shrinking circle and along a spiral-shaped path simultaneously. To update the position of whales and to mimic the simultaneous behavior of the shrinking encircling mechanism and the spiral model, in the algorithm they are chosen randomly with a probability of 50% in each iteration. The two approaches are expressed as follows:

$$\mathrm{X}\left(\mathrm{t}+1\right)=\left\{\begin{array}{l}{\mathrm{X}}^{\mathrm{*}}\left(\mathrm{t}\right)-\mathrm{A}·\mathrm{D}\\ {\mathrm{D}}^{\prime }·{\mathrm{e}}^{\mathrm{b}\mathrm{l}}·\cos \left(2\mathsf{\pi }\mathrm{l}\right)+{\mathrm{X}}^{\mathrm{*}}(\mathrm{t})\end{array}\right.\begin{array}{l}\mathrm{i}\mathrm{f} \mathrm{p}<0.5\\ \mathrm{i}\mathrm{f} \mathrm{p}\ge 0.5\end{array}$$

(137)

where $\mathrm{p}$ is a random number in the range [0, 1].

In the search for prey behavior, whales look for a prey location, randomly giving the algorithm the global searching capability. It is expressed as follows:

$$\mathrm{D}=\left|\mathrm{C}·{\mathrm{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}-\mathrm{X}(\mathrm{t})\right|$$

(138)

$$\mathrm{X}\left(\mathrm{t}+1\right)={\mathrm{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}-\mathrm{A}·\mathrm{D}$$

(139)

where ${\mathrm{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}$ is a random position vector chosen from the current whale population.

The WOA pseudo code is composed of four main parts, and they are given as follows [78]:

• Initialize the whale population ${\mathrm{X}}_{\mathrm{i}}$ with $\mathrm{i}=1,2,\dots ,\mathrm{n}$.
  Calculate the fitness of each search agent.
  Define the best search agent ${\mathrm{X}}^{\mathrm{*}}$.
• While $\mathrm{t}<\mathrm{m}\mathrm{a}\mathrm{x} \mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r} \mathrm{o}\mathrm{f} \mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}$:
  • For each search agent:
    Update $\mathrm{a}$, $\mathrm{A}$, $\mathrm{C}$, $\mathrm{l}$, and $\mathrm{p}$.
    If $\mathrm{p}<0.5$: {

    -

    If $\left|\mathrm{A}\right|<0.5$, update the position of the current search agent by (134).

    -

    Elseif $\left|\mathrm{A}\right|\ge 0.5$, select a random search agent ${\mathrm{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}$ and update the position of the current search agent by (139).

    }

    Elseif $\mathrm{p}\ge 0.5$, update the position of the current search agent by (136).
  • End of for loop.
  • Check if any search agent goes beyond the search space and amend it.
  • Calculate the fitness of each search agent.
  • Update ${\mathrm{X}}^{\mathrm{*}}$ whether a better solution has been found.
  • $\mathrm{t}=\mathrm{t}+1$.
• End.
• Return ${\mathrm{X}}^{\mathrm{*}}.$

In this basic version of the proposed WOA, other evolutionary operations are not included, but they might have been included by allowing hybridization with other evolutionary search schemes.

Three Magnetic Targets Using Individual Memory-Based WOA and LM Algorithm

The authors in reference [79] deal with multi-permanent magnet localization, exploiting an array of magnetometers. The aim is to estimate the location of the permanent magnets inside the human body by measuring the magnetic field. To obtain the best estimation of the measurements, the authors propose an improved WOA combined with the LM algorithm.

Each magnetometer in the array measures the magnetic flux that is given by the linear superposition of the flux densities generated by each permanent magnet. The objective function $\mathrm{E}$ of the total error is obtained by considering the sum of the three errors between the measured values by the magnetometers and the theoretical values; it is given as follows:

$$\mathrm{E}=\sum _{\mathrm{l}=1}^{\mathrm{N}}\left[{\left({\mathrm{B}}_{\mathrm{l}\mathrm{x}}^{\prime }-{\mathrm{B}}_{\mathrm{l}\mathrm{x}}\right)}^2+{\left({\mathrm{B}}_{\mathrm{l}\mathrm{y}}^{\prime }-{\mathrm{B}}_{\mathrm{l}\mathrm{y}}\right)}^2+{\left({\mathrm{B}}_{\mathrm{l}\mathrm{z}}^{\prime }-{\mathrm{B}}_{\mathrm{l}\mathrm{z}}\right)}^2\right]$$

(140)

where each ${\mathrm{B}}_{\mathrm{l}}^{\prime }$ represents one of the triaxial components of the magnetic flux density measured by the $\mathrm{l}$-th magnetometer, ${\mathrm{B}}_{\mathrm{l}}$ is the theoretical value of one of the triaxial components of the magnetic flux, and $\mathrm{N}$ is the total number of magnetometers. A generic component ${\mathrm{B}}_{\mathrm{l}}$ is given by the superimposition of the magnetic fields generated by all permanent magnets, and it includes the distance variables between each permanent magnet and each magnetometer.

The name of the proposed algorithm is the Individual Memory-based Whale Optimization and Levenberg–Marquardt (IMWO-LM) algorithm; the IMWO algorithm is used to provide an approximate solution of each magnet pose, which is given as a guess input to the LM algorithm, which is suitable for finding the local minima.

In the proposed IMWO, $\mathrm{S}$ is defined as the number of whales and each of them searches for a solution of dimension $\mathrm{d}$; the vector position for the $\mathrm{i}$-th whale is defined as ${\mathrm{X}}^{\mathrm{i}}=[{\mathrm{x}}_1^{\mathrm{i}},{\mathrm{x}}_2^{\mathrm{i}},\dots ,{\mathrm{x}}_{\mathrm{d}}^{\mathrm{i}}]$. The three behaviors presented in the WOA are here modified to give whale the ability of self-memory and to always remember the previous best position. The IMWO behaviors are defined as follows:

• Search for prey. The current search agent position is updated considering both the position of the randomly selected whale and the historical optimal position. It is given as follows:

  $$\left\{\begin{array}{l}{\mathrm{D}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}^{\mathrm{i}}=\left|\mathrm{C}·{\mathrm{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}^{\mathrm{i}}-{\mathrm{X}}^{\mathrm{t},\mathrm{i}}\right|\\ {\mathrm{D}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\mathrm{i}}=\mathrm{C}{\mathrm{X}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\left(\mathrm{t}-1\right),\mathrm{i}}-{\mathrm{X}}^{\mathrm{t},\mathrm{i}}\\ {\mathrm{X}\prime }^{\mathrm{t},\mathrm{i}}={\mathrm{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}^{\mathrm{i}}-\mathrm{A}{\mathrm{D}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}^{\mathrm{i}}+\mathrm{A}{\mathrm{D}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\mathrm{i}} \end{array}\right.$$

  (141)

  where $\mathrm{t}$ identifies the $\mathrm{t}$-th iteration, $\mathrm{i}$ identifies the search agent, ${\mathrm{X}}^{\mathrm{t},\mathrm{i}}$ is the current position, ${\mathrm{X}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\left(\mathrm{t}-1\right),\mathrm{i}}$ is the historical optimal position, ${\mathrm{D}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\mathrm{i}}$ is the random distance between the current position and its historical optimal position, ${\mathrm{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}^{\mathrm{i}}$ is a randomly selected search agent, ${\mathrm{D}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}^{\mathrm{i}}$ is the random distance between a randomly selected search agent and current search agent, $\mathrm{C}$ and $\mathrm{A}$ are weight vectors defined as in the WOA, and ${\mathrm{X}\prime }^{\mathrm{t},\mathrm{i}}$ is the newly updated location of the $\mathrm{i}$-th search agent.
• Encircling prey or siege predation. The best candidate search agent ${\mathrm{X}}^{\mathrm{*}}$ of the current population is considered the target prey at each iteration; search agents have memory and can remember the previous best position. It is given as follows:

  $$\left\{\begin{array}{l}{\mathrm{D}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{r}}^{\mathrm{i}}=\left|\mathrm{C}{\mathrm{X}}^{\mathrm{*}}-{\mathrm{X}}^{\mathrm{t},\mathrm{i}}\right|\\ {\mathrm{D}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\mathrm{i}}=\mathrm{C}{\mathrm{X}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\left(\mathrm{t}-1\right),\mathrm{i}}-{\mathrm{X}}^{\mathrm{t},\mathrm{i}}\\ {\mathrm{X}\prime }^{\mathrm{t},\mathrm{i}}={\mathrm{X}}^{\mathrm{*}}-\mathrm{A}{\mathrm{D}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{r}}^{\mathrm{i}}+\mathrm{A}{\mathrm{D}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\mathrm{i}}\end{array}\right.$$

  (142)

  where ${\mathrm{D}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{r}}^{\mathrm{i}}$ is the random distance between the optimal candidate search agent and the current search agent, and ${\mathrm{X}}^{\mathrm{*}}$ is the optimal candidate search agent in the whole population so far. In such a way, the current search agent is updated to different positions around ${\mathrm{X}}^{\mathrm{*}}$ and ${\mathrm{X}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\left(\mathrm{t}-1\right),\mathrm{i}}.$
• Spiral bubble-net feeding maneuver or bubble predation. The helix-shaped movement of the whale is obtained using the spiral equation, considering the position of the current search agent, its historical optimal position, and the position of the optimal candidate search agent. The current position is given as follows:

  $${\mathrm{X}}^{\prime \mathrm{t}\mathrm{i}}=\left({\mathrm{r}}_1{\mathrm{D}}^{\mathrm{*}}+{\mathrm{r}}_2{\mathrm{D}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\mathrm{i}}\right)·{\mathrm{e}}^{\mathrm{b}\mathrm{l}}·\cos \left(2\mathsf{\pi }\mathrm{l}\right)+{{\mathrm{r}}_1\mathrm{X}}^{\mathrm{*}}+{\mathrm{r}}_2{\mathrm{X}}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^{\left(\mathrm{t}-1\right),\mathrm{i}}$$

  (143)

  where ${\mathrm{D}}^{\mathrm{*}}=\left|{\mathrm{X}}^{\mathrm{*}}-{\mathrm{X}}^{\mathrm{t},\mathrm{i}}\right|$ is the distance between the current search agent and the optimal candidate search agent, $\mathrm{b}$ is a constant for defining the shape of the exponential spiral, $\mathrm{l}$ is linearly decreased from 1 to −1 throughout iterations, and ${\mathrm{r}}_1$ and ${\mathrm{r}}_2$ are random numbers in the range [0, 1].

The upper and lower boundary of the pose parameters are set, and at each iteration, the algorithm must check if the search agents violate them.

3.15. Summary of the Investigated Methods

In this section,

Table 1. Summary of the investigated methods.

Investigated Methods	Main Characteristics	Significative Advantages (ADV) and/or Disadvantages (DIS)	Applications
Newton’s Method	Iterative method. Fixed step size. Suitable for unconstrained problems. Second-order method. Quadratic rate of convergence. Needing of an initial design point.	High computation cost. (DIS) Needing a twice differentiable function. (DIS) Good initial starting point. (DIS)	No applications investigated for this method.
GD	Iterative method. Variable step size. First-order method. Dealing with convex functions. Unconstrained problems. Convergence time linked to the parameter $\alpha$.	Penalty function can be used. (ADV) The parameter $\alpha$ needs to be set properly to avoid skipping the optimum. (DIS) Function must at least once differentiable. (DIS) It could skip the global optimum and converge to a local one. (DIS)	No applications investigated for this method.
Gauss– Newton	Iterative method. Suitable for convex non-linear least-squares problems. It requires at least once-differentiable function. Assumption: the inverse of the gradient exists.	It is suitable only for a specific class of problems. (DIS)	No applications investigated for this method.
ILS Method/ LAMBDA	Iterative method. Suitable for ILS problems.	It is suitable only for a specific class of problems. (DIS) Increase the resolution speed of the problem. (ADV)	To solve frequency ambiguity problems in GNSS applications [38,39].
LM	Iterative method. Suitable for non-linear optimization problems. It interpolates between the GN and the GD.	A good initial guess is needed. (DIS) If the problem has several local optima, a wrong global optimum could be provided. (DIS)	No applications investigated for this method.
Lagrange Multipliers Method	Analytic method. Suitable for constrained problem (with equality and/or inequality). Differentiable functions are needed.	It guarantees only that an optimum has been found but does not guarantee that is the global one. (DIS)	No applications investigated for this method.
Lagrange Duality	Analytic method. It transforms a minimization or maximization problem into its dual one.	It could provide an easier problem to solve. (ADV)	To estimate the measured range of radio source anchors [47].
GA	Iterative method. Metaheuristic method. Suitable for constrained and unconstrained problems. Inspired by natural selection of animals. It works directly with strings of bits It deals with discontinuous, nondifferentiable, stochastic, or highly non-linear functions. It is suitable for mixed-integer type problems.	It can solve many types of different problems. (ADV) Large solutions space can be spanned. (DIS) Slow solving speed. (DIS) Each considered problem needs a specific tunning. (DIS) Easy to tune. (ADV)	Computation of parameters to superimpose information from different sensors and radio signal [9].
QGA	It is a variation of the GA. Unlike GA, QGA deals with Q-bits and therefore can exploits superimposition of states that allows to reduce the solution time convergence.	The solving speed is $\sqrt{\mathrm{N}}$ times faster than the speed of classic GA. (ADV) High complexity. (DIS)	To calibrate the parameters of a wheeled vehicle to compensate the systematic errors [11]
DE	Iterative method. Metaheuristic method. Suitable for constrained and unconstrained problems. Not inspired by natural selection of animals. It works directly with vectors It deals with continuous problems.	Good global search capability. (ADV) Good local search capability. (ADV) A memory function. (ADV) A slow convergence speed. (DIS) Performance affected by the trial vector. (DIS) Needing a tuning procedure. (DIS)	To calculate the z-coordinate of the receiver in VLC system for localization purposes [59].
MEA	Iterative method. Metaheuristic method. Suitable for constrained and unconstrained problems. Inspired by learning modes and human mind activities. It deals with discontinuous, non-differentiable, stochastic, non-linear, non-convex functions.	Fast rate of convergence. (ADV) It can deal with data from sensors that are used as population generator when work at high frequency. (ADV)	To update the forward or backward movement of a walker in a PDR system [60].
PSO	Iterative method. Metaheuristic method. Suitable for unconstrained and constrained (by using penalization methods) problems. Inspired by natural swarm behaviors of animals. It deals with discontinuous, non-differentiable, stochastic, non-linear, non-convex functions.	Fast rate of convergence. (ADV) Needing a tuning procedure to set the parameters. (DIS) Possibility to be stuck in a local optimum without converging to the global. (DIS)	To train an ANN (used in a fingerprint localization system) [63]. To estimate the RSSI range by minimizing the squared residual between the real and the measured distances [64]. To solve the localization problem modelled as a Markov process in which the position is extracted by maximizing the a posteriori probability density. The approach is an alternative to EKF or PF [65]. To estimate the best ANN connection weights during the training process for an on moving robot that is able to avoid collisions [77]
DEPSO	Hybrid algorithm in which the concepts of DE and PSO are fused.	Better exploration of local optima. (ADV) Better convergence speed than classic DE. (ADV)	To reduce the particle number and the computational time complexity in a PF [67].
BA	Iterative method. Metaheuristic method. Suitable for unconstrained and constrained (by using penalization methods) problems. Inspired by the echolocation behavior It deals with discontinuous, non-differentiable, stochastic, non-linear, non-convex functions.	Fast rate of convergence. (ADV) Good ability to explore the search space and to quit from local optima. (ADV) Standard BA might provide a premature convergence, that it is improved by adopting LBBA. (DIS) Needing a tuning procedure. (DIS)	To localize the robot position by minimizing the error between the real and the measured position [69].
ABC	Iterative method. Metaheuristic method. Suitable for unconstrained and constrained problems. Inspired by the intelligent foraging behavior of honeybees. It deals with discontinuous, non-differentiable, stochastic, non-linear, non-convex functions.	Good ability to explore the search space. (ADV) Slow rate of convergence. (DIS) Needing a tuning procedure. (DIS)	To solve the sensor localization problem in a wireless sensor network to get new anchors; the location is estimated by minimizing the mean localization error [72].
SA	Iterative method. Metaheuristic method. Suitable for unconstrained and constrained problems. Inspired by the process of metal crystalline formation from the liquid phase. It deals with discontinuous, non-differentiable, stochastic, non-linear, non-convex functions.	Good ability to explore the search space. (ADV) Slow rate of convergence. (DIS) Needing a careful tuning procedure. (DIS)	To estimate the real position of sensors by exploiting the anchors positions and the measured distances between them and the sensors [74].
ACO	Iterative method. Metaheuristic method. Inspired by the behavior of ants. Suitable to solve the problems in which the best path needs to be selected.	Solve a specific class of problems. (DIS)	To optimize the best ANN topology during the training process for an on moving robot that is able to avoid collisions [77]
WOA	Iterative method. Metaheuristic method. Suitable for unconstrained and constrained problems. Inspired by the behavior of a whale to capture a prey It deals with discontinuous, non-differentiable, stochastic, non-linear, non-convex functions.	Slow convergence speed. (DIS) In some cases, it is hybridized with other algorithms to look for local optima. (DIS)	To estimate the location of the permanent magnets inside the human body by measuring the magnetic field; the total error is estimated by minimizing the sum of the three errors between the measured values and the theoretical values of the magnetic field [79].

reports a summary of the investigated methods, where the main characteristics, the significative advantages and disadvantages, and the investigated applications are listed to have a comprehensive view of the conducted research.

4. Conclusions

This paper provides an extensive study related to the optimization methods used in the field of localization. In the first part of the article, a light introduction to the optimization problem was reported, providing a general mathematical framework for the discipline. Besides the SO problem, some concepts about the MO problem have been introduced as well, and they aim to provide a clear and complete view of the basis of the subject.

In the second part of the article, extensive research about some interesting algorithms used in localization has been reported, providing for each of them the concept of the algorithm, the related mathematics model, and an example of application. Both types of algorithms for local and global optimization have been described, and it can be noted that metaheuristic algorithms, such as GA, DE, or PSO, are the most used in the field of localization for optimization problems due to their capabilities to manage non-convex and non-linear problems and provide satisfactory results. Moreover, the GA (as well as its modified versions, such as the QGA) is the most versatile optimization tool that allows for different types of problems such as continuous, discrete, or mixed-integer types. However, the standard GA has a slow speed of convergence. In contrast, PSO has a higher speed, and by observing Table 1, it is the most used among the investigated methods. Furthermore, PSO was shown to be the most prone method that allows for hybridization with other types, such as DEPSO or SWIRL; the modified methods improved the performance in solving specific problems.

In general, by comparing metaheuristic algorithms to gradient-based algorithms (that are especially used for local optimum), it is noted that stochastic algorithms do not need any Jacobian (or Hessian) matrix inversion, demonstrating that they are very flexible and suitable optimization tools for many optimization problems. On the other hand, they need to tune the parameters for each specific solved problem.

In the end, it is highlighted that none of the investigated methods is guaranteed to provide a real global optimum.