Gradient eigendecomposition invariance biogeography-based optimization for mobile robot path planning

The path planning for mobile robots has attracted extensive attention, and evolutionary algorithms have been applied to this problem increasingly. In this paper, we propose a novel gradient eigendecomposition invariance biogeography-based optimization (GEI-BBO) for mobile robot path planning, which has the merits of high rotation invariance and excellent search performance. In GEI-BBO, we design an eigendecomposition mechanism for migration operation, which can reduce the dependency of biogeography-based optimization (BBO) on the coordinate system, improve the rotation invariance and share the information between eigensolutions more effectively. Meanwhile, to find the local optimal solution better, gradient descent is added, and the system search strategy can reduce the occurrence of local trapping phenomenon. In addition, combining the GEI-BBO with cubic spline interpolation will solve the problem of mobile robot path planning through a defined coding method and fitness function. A series of experiments are implemented on benchmark functions, whose results indicated that the optimization performance of GEI-BBO is superior to other algorithms. And the successful application of GEI-BBO for path planning in different environments confirms its effectiveness and practicability.


Introduction
In recent decades, mobile robot path planning (MRPP) has become an indispensable aspect of artificial intelligence in robotics and has been widely studied and discussed (Manikas et al. 2007). The goal of path planning is to search out an optimal or nearly optimal collision free path from the initial state to the target state according to a certain performance index (such as time and distance) (Rasekhipour et al. 2017). According to different environmental information, path planning can be divided into global path planning and local path planning (Sucan et al. 2012). Global path planning is carried out in the known environment. By modeling the current environmental information, mobile robots can use the existing mature path planning algorithm or improved algorithm to plan an optimal path from the start to the end in the established environmental model (Borenstein and Koren 1991). In the local path planning problem, the environmental information is usually unknown, focusing on considering the current local environmental information of the mobile robot to make the robot have good obstacle avoidance ability. The working environment of the mobile robot is detected by sensors to obtain information such as the location and several properties of obstacles (Zeng et al. 2016). However, without complete environmental information, the results of local path planning may not be best or even incorrect. Therefore, global path planning, which can be regarded as an optimization problem, has attracted more attention and became a hot topic in research of path planning (Alam and Rafique 2015).
Researchers have developed many methods for global MRPP. Traditional path planning methods include grid method (Borenstein and Koren 1991), artificial potential field method , visual graph method (Sucan et al. 2012), neural network method (Kim and Lee 2015), etc. However, these methods partly have some problems such as complex calculation, low efficiency and easy to fall into local optimization, which usually leads to certain inconveniences in solving the problem of global MRPP. Meanwhile, evolutionary algorithms for MRPP show a growing number of advantages and attract increasing attention (Lim et al. 2008).
Using the evolutionary algorithm for global MRPP, which can put the problem of path planning down to the problem of finding the optimal path with the minimum cost (Kennedy and Eberhart 1995). In addition, many evolutionary algorithms, such as genetic algorithm (GA) (Bakdi et al. 2017), particle swarm optimization algorithm (PSO) (Tharwat et al. 2019) and ant colony algorithm (ACO) (Kumar et al. 2018), have strong computing capability and robustness, which can achieve optimization results well. At present, many evolutionary algorithms have been applied to the problem of MRPP (Castillo et al. 2006). For example, considering the slow convergence speed of ACO, Liu et al. (2017) proposed an improved ACO to improve its convergence speed and applied it to the MRPP in the grid environment. Marco et al. (2015) combined the artificial bee colony algorithm with the evolutionary planning algorithm to solve the problem of MRPP, through a set of local processes to refine the feasible path. The comparison results of the method on a set of benchmark problems and the experimental results on a real mobile robot show that the method has good performance. Hong et al. (2013) proposed a coevolutionary improved GA for global MRPP, which puts forward an effective fitness function and modifies the genetic operator of traditional GA. For the problem of global smooth MRPP, Song et al. (2017) advanced a new multimodal delayed PSO (MDPSO). The test results based on benchmark function show that the performance of MDPSO is better than other five famous PSO algorithms. Finally, the application of MDPSO in the global smooth MRPP further proves that it has better performance than the global smooth path generated by GA in previous studies. These evolutionary algorithms have been used to solve the problem of MRPP and have achieved remarkable results, but in general, these methods still have some shortcomings in computational complexity, local optimization and adaptability.
Biogeography-based optimization (BBO) (Simon 2008) is an evolutionary algorithm based on the concept of biogeography. Based on a mathematical model, the algorithm describes the migration of species between habitats, that is, from unsuitable habitats to suitable habitats (Simon et al. 2010). BBO is an efficient bionic search algorithm, which has the advantages of relatively simple principle, less adjustment parameters, easier implementation of the algorithm and high operation efficiency. Many researchers have improved BBO and applied it to many fields (Mirjalili et al. 2014), including MRPP. For example, Zhu et al. (2014) proposed a method based on chaotic predator-prey BBO (CPPBBO) to solve the path planning problem of unmanned aerial vehicle, and the simulation results show that CPPBBO is more effective than other algorithms. Mo et al. (2015) advanced a biogeography PSO algorithm to plan the path of mobile robot by combining the BBO and PSO, and then using the BPSO algorithm, the optimal path based on the approved voronoi boundary network is found out. Yang et al. (2017) put forward an improved optimization algorithm based on biogeography to solve the problem of global MRPP in the static environment. It can be seen from these articles that BBO has been effectively applied to MRPP. However, due to the strong dependence of BBO on the coordinate system, the performance of the algorithm is inadequate when solving high-dimensional problems, and the ability of BBO to mine the global optimal solution is also slightly insufficient.
In order to solve these problems, we proposed a novel gradient eigendecomposition invariance biogeography-based optimization algorithm (GEI-BBO). The innovation of this paper is as follows: 1. An eigendecomposition-based migration strategy is proposed, which has the advantages of reducing the dependence on coordinates, improving the rotation invariance of BBO and sharing the information between eigensolutions more effectively. 2. On the basis of the eigendecomposition-based migration strategy, the gradient descent strategy is added to search the neighborhood of the best individual, which can more effectively mine and not only search the optimal solution, so as to improve the local search ability of BBO. 3. The system search strategy is put forward to ensure that the algorithm with gradient descent strategy finds the global optimal solution, which can carry out a range search for each dimension, so as to cover the whole search range well and avoid falling into local optimum. 4. GEI-BBO is combined with the cubic spline interpolation method to solve the problem of MRPP. This paper defines the coding method based on the path node and constructs fitness function which aims at avoiding the obstacle and finding the shortest path. The rest of the paper is organized as follows: Section 2 introduces the basic BBO. Section 3 introduces the gradient eigendecomposition invariance biogeography-based optimization. Section 4 shows the method for path planning, and Section 5 shows the simulation results and analysis.

Biogeography-based optimization
BBO is a new population-based optimization algorithm to solve global optimization problems, based on the concept of biogeography. In biogeography, each habitat is regarded as an individual, and the index to measure its quality of life is called habitat suitability index (HSI) (Mo and Xu 2015). Good habitats for species have high HSI, and bad habitats have low HSI. A habitat has many characteristics, such as area, temperature, rainfall and so on. These characteristics will affect the habitat, so they are referred to as the suitability index variables (SIVs). In BBO, a candidate solution to the problem is considered as a habitat. Each candidate solution is associated with a fitness value that is similar to the HSI of the habitat. High HSI habitats represent better solutions, and low HSI habitats represent worse solutions. Immigration rates and emigration rates for each candidate solution are used to share features probabilistically between habitats.
In BBO, the rate of immigration λ and the rate of emigration μ determine the dynamic movement between habitats, depending on the number of species in the habitat (Simon 2011). The equation for the rate of immigration λ k and the rate of emigration μ k of k species can be written as follows ( Fig. 1): where I and E are the maximum immigration rate and the maximum emigration rate, respectively. n = s max is the maximum number of species in the habitat. If I = E, then λ k + μ k = I . BBO is different from other evolutionary algorithms because of its migration strategy. The migration strategy is based on the migration model, that is, the immigration rate and the emigration rate equation. The migration model generally adopts the linear model. When a habitat X i needs to be relocated according to the migration rate, some methods, such as roulette, will be used to select a source habitat X j probabilistically according to the migration rate, and then a SIV will be selected randomly in the source habitat to replace directly or modify X i . BBO can improve the solution through the migration operation.
The mutation operation of BBO is to randomly mutate the habitat according to the mutation rate m(s), that is, changing the SIVs of the habitat. m(s) is determined by the following equation: where m(s) is the mutation rate of habitat, m max is the maximum mutation rate, and P max is the maximum probability. BBO can increase the diversity of species through mutation operation (Du et al. 2009).

Gradient eigendecomposition invariance biogeography-based optimization
In this section, we introduce a method to solve the problem of two-dimensional static global MRPP. We introduce three innovative points of the proposed algorithm, namely the eigendecomposition-based migration, gradient descent and system search strategy. Finally, the algorithm is described.

Eigendecomposition-based migration
In order to improve the rotation invariance of BBO algorithm, the eigendecomposition-based migration is proposed in this paper. The core of the eigendecomposition-based migration is to rotate the original coordinate system into the eigenvector-based coordinate system, in which habitats can share their information more effectively. The proposed method illustrates that the migration of BBO can be carried out more effectively in the eigenvector-based coordinate system. The eigendecomposition-based migration is to carry out the eigendecomposition of population and migrate the decomposed population, which will make BBO run more effectively (Fig. 2).
where i = 1, 2, ..., n, j = 1, 2, ..., D, n is the population size, D is the number of independent variables, G is the number of iterations, and H G i is the ith population in the Gth iteration.
H is not a square matrix of order n; it cannot be eigendecomposed directly, so we introduce the concept of covariance matrix. Make: Then, the covariance matrix is where cov(c i , c j ) is the covariance between the ith independent variable and the jth independent variable in Gth generation, which is defined as follows: where ξ i and ξ j represent the mean values of the ith and jth independent variables, respectively. Thus, it can be seen that covariance matrix cov(H ) is a symmetric matrix of order n, and the following standard form can be obtained by eigendecomposition: where Q H is a D × D eigenmatrix composed of feature vector of cov(H ) and H is a diagonal matrix composed of all the eigenvalues of cov(H ). After eigendecomposition, the eigenvector is obtained. We also need to rotate the matrix H with the eigenvectors. That is: where i = 1, 2, ..., n and j = 1, 2, ..., D. Migration of the eigH G i obtained by the rotation operation results in the generation of eigH G+1 i . Because of the eigendecomposition and rotation operation, the migration can run more efficiently. After the migration operation, H should be rotated back and used as the population of the new generation. The following operation is enough.
In summary, the steps of the migration algorithm based on eigendecomposition are as follows:

Algorithm 1 Eigendecomposition-based migration
Compute the matrix cov(H ) 3: Apply eigendecomposition according to Equation 8 4: Rotate according to Equation ?? 5: for k = 1 to D do 6: if rand < λ i then 7: select eig H G j with probability μ j 8: 10: end if 11: end for 12: Get H G+1 i according to Equation 10 13: end for

Gradient decent strategy
BBO algorithm is effective in searching the global optimal region, but it is difficult in mining the global optimal region. The local search performance of BBO can be improved by adding gradient descent into BBO to search the best individual domain.
The following conditions are required to activate the local search of gradient descent: If both of these conditions are met, then we apply gradient descent strategy to N G optimal individuals. G flag starts at zero and increases by 1 when the following conditions are met.
where R I is the improvement ratio, which is the relative improvement in the minimum generation value from the i(t) iteration to the i(t + 1) iteration. ε 1 is a predefined threshold. Moreover, by adding G A through G INC , the incremental G A delays the next call to gradient descent local search. The reason why the gradient is called infrequently is that when the algorithm falls into a local optimal state, multiple local searches may not be helpful, and the function evaluation times will be wasted. In short, gradient descent is activated when R I is continuously smaller than ε 1 . We use the f mincon function in MATLAB to achieve the gradient descent of the best individual. The gradient descent algorithm is shown as follows: Apply gradient descent to the N G best individuals in the population 9: G f lag = 0 10:

System search strategy
In order to ensure that the algorithm with gradient descent strategy can find the global optimal solution, we add the system search strategy, which can cover the whole search range well and avoid falling into the local optimal.
The conditions for activating the system search strategy are similar to those for gradient descent in 3.2, and the improved ratio R I is also used. When the following formula is satisfied, the system search strategy will be activated.
where ε is the computer precision, that is, when there is no improvement in the minimum cost value from the i(t) iteration to the i(t +1) iteration, a global search will be conducted for N s optimal individuals to execute the system search strategy. The steps of the system search strategy are given as Algorithm 3.

Algorithm 3 System search strategy
Require: As shown in the table, we increment or decrement the independent variables in accordance with in the search space, where α 0 is a specific score that can be evaluated as α 0 = 0.1. The system search strategy decreases the value of the given dimension po p (i,k) by an increment , equal to 10% of the size of the search space, one increment at a time, until the value reaches the lower bound of the search space. Then, the system search strategy increases the value of the given dimension, one increment at a time, until it reaches the upper limit of the search space. By performing this process

Description of the GEI-BBO algorithm
The GEI-BBO algorithm proposed in this paper includes the contents mentioned in Sects. 3.1, 3.2 and 3.3. Before the number of iterations reaches the maximum number of iterations m, n habitats are selected to migrate according to the proportional parameter P in each iteration. If the random number is less than P, the migration is based on the eigendecomposition, and if the random number is greater than P, the migration is based on the standard BBO. After migration, the mutation was carried out. Then, the gradient descent strategy is implemented on the N G best individuals and the system search strategy is implemented on the N s best individuals. N G and N s are preset parameters. The specific algorithm flow is as follows:

Method for path planning
Combining the improved algorithm with cubic spline interpolation method, the coding method based on the path node is defined, and the method and fitness function aiming at Algorithm 4 GEI-BBO 1: Initialize a population of N habitats 2: for it = 1 to M do 3: for t = 1 to N do 4: if thenrand < P 5: Apply eigendecomposition-based migration (3.1) 6: else 7: Apply standard BBO migration 8: Apply mutation operation according to mutation rate m(s) 9: end if 10: Apply gradient descent to the N G optimal individuals in the population (3.2) 11: Apply system search strategy to the N S optimal individuals in the population (3.3) 12: end for 13: end for solving the obstacle avoidance and shortest path of mobile robot are constructed to solve the problem of MRPP. Figure 3 shows a conceptual map of habitat migration in MRPP using GEI-BBO to understand how the proposed approach works.
In the figure, habitat 1 has the highest HSI, representing the most suitable habitat, followed by habitat 2, habitat 3 and habitat 4. The higher the HSI, the more suitable for species growth, the less need to be changed, so the lower the immigration rate, the higher the emigration rate. It can be seen that habitat 1 has the highest emigration rate, while habitat 4 has the highest emigration rate. Residents represent the path nodes. So the fourth habitat accepts many residents (path nodes) from other habitats, as shown in different colors. Purple nodes and links also describe mutations that occur in all habitats, regardless of their HSI values. This example shows how path planning evolved using the proposed approach. The specific coding method and fitness function are shown below.

Cubic spline interpolation
Cubic spline interpolation is a kind of piecewise interpolation method to form a smooth curve through a train of interpolation point intervals based on cubic polynomials. The curve of moving path of mobile robot fitted by cubic spline interpolation method is smoother than the curve fitted by straight line and circular arc. In this paper, the cubic spline interpolation method is integrated into the improved GEI-BBO algorithm to solve the problem of static global MRPP.
On the interval [a, b], taking n + 1 nodes a = x 0 < x 1 < · · · < x n = b, if s(x) satisfies the following conditions: The cubic spline interpolation function is a piecewise cubic polynomial, on each small interval [x i , x i+1 ], which can be written as: is the undetermined coefficient, so s(x) has 4n undetermined coefficients. To solve for s of x, we need 4n conditions. From s(x i ) = f i , i = 0, 1, ..., n, n + 1 interpolation conditions can be obtained. And from s(x) ∈ C 2 [a, b], s(x) is second differentiable on the interval [a, b], so the first order is differentiable and s(x) is continuous. Therefore, the following conditions can be obtained: So far, there are 4n−2 conditions. In the actual calculation, two boundary conditions need to be introduced to calculate s(x). Commonly used boundary conditions are: (1) The value of the first derivative at the two endpoints is given.
(2) The value of the second derivative at the two endpoints is given.

Encoding
It can be seen from the above that cubic spline interpolation is a piecewise interpolation method, and the junction between segments is called the path node. The splines between segments are different, and the whole spline curve is continuous in the first order and is continuous in the second order at the path node. Therefore, the directions between segmented splines may be different, and the path node represents the maximum turning times of the entire path. According to this, this paper takes all nodes on a path as the code of a habitat individual, that is, a habitat individual represents all nodes on the corresponding path. The x-coordinate set of all m path nodes on a path constitutes the m-dimensional xcoordinate of habitat individuals, and correspondingly, the y-coordinate set of all m path nodes on a path constitutes the m-dimensional y-coordinate of habitat individuals.
Suppose we know the coordinates of m path node (x m1 , y m1 ), (x m2 , y m2 ), ..., (x mm , y mm ), starting point coordinates (x s , y s ) and terminal coordinates (x t , y t ), on the interval (x s , x m1 , x m2 , ..., x mm , x t ) and (y s , y m1 , y m2 , ..., y mm , y t ), the abscissa and ordinate of n interpolation points are obtained by cubic spline interpolation. In this way, we get n interpolation points, and the line between the path node and the interpolation point and the starting point and the ending point is the path of the mobile robot that we want.

Fitness function
There are two conditions for global MRPP: (1) Do not collide with obstacles; (2) The path length should be as short as possible. In this paper, the evaluation criterion of fitness function is the shortest path length satisfying the above conditions. In this paper, the constructed fitness function is where β is a very large number, and it is used to impose a very large cost value to exclude paths that have collisions, and it can take on a value of 1000. L is the distance from the starting point to the terminal point; the calculation formula is where v is a flag variable. v = 0 means the path is a noncollision path; otherwise, there is a collision path. The initial value of v is set to 0. In the rectangular obstacle environment, the following equation is used to determine whether the obstacle is encountered: < y y j &&y y j < (yobs k + w_obs k +space) where j is the number of nodes, including path nodes and interpolation points. x x j and y y j represent the abscissa and ordinate of the node. (xobs k , yobs k ) is the starting point of the kth rectangle obstacle, that is the bottom left corner of the rectangle. l_obs k , w_obs k is the length and width of the rectangular obstacle. space is a small positive number to avoid hitting the edge of the rectangle.
At each node, judge whether it collides with the obstacle according to formula 18. When it encounters the obstacle, d will generate a nonzero change value. When d changes, v will be added with one. If the path does not pass through any obstacle, then v = 0; otherwise, v is equal to a number that is not zero.

Path planning process
We combine GEI-BBO with cubic spline interpolation to design the path planning steps; the flowchart of path planning is shown in Fig. 6.
Step 1 Determine the number of path nodes m and initialize the algorithm.
Step 2 The coordinates of n interpolation points are obtained by cubic spline interpolation.
Step 3 The path length L and the flag variable v are calculated, respectively, to obtain the value of fitness function.
Step 4 Using GEI-BBO to update the ordered habitat individuals.
Step 5 Determine whether the maximum number of iterations is reached. If yes, output the optimal result directly. If not, return to the step 2 for circulation (Fig. 4).

Simulation results and analysis
In order to prove the effectiveness and practicability of the proposed method, we first use 23 benchmark test functions to compare and test the GEI-BBO and then use it in the simulation experiment of path planning and conduct in-depth analysis and summary of the experimental results.

Experiments on benchmark test functions
We use the benchmark test function proposed by Yao et al. (1999) for comparative test to verify the effectiveness of the algorithm proposed in this paper. There are 23 test functions in total, and the basic information is shown in Table 1. It can be seen from Table 1 that function F1-F13 is a highdimensional problem, function F1-F5 is unimodal, function F6 is a step function, function F7 is a quartic function with noise, function F8-F13 is a multimodal function, the number of local minimum values increases exponentially with the problem dimension, and function F14-F23 is a lowdimensional function with only a few local minimum values (Yao et al. 1999). This group of test functions includes not only low-dimensional and single-mode functions, but also many high-dimensional and multi-mode functions. It can detect the convergence speed of the algorithm well and also reflect the ability of the algorithm to get rid of the poor local optimization and find the better near global optimization. We selected some of the more typical test functions, whose three-dimensional diagram is shown in Fig. 5.
In Tables 2 and 3, MEAN is the average of ten test results, BEST represents the best result of 10, and the bold part is the best one among 11 algorithms. Our proposed GEI-BBO obtains 16 optimal BEST and 14 optimal MEAN solutions in the optimization process of 23 functions, and its performance is better than the other 10 comparison algorithms. The experimental results verify the effectiveness of the optimization process of GEI-BBO on the benchmark function.

Simulation on path planning
To verify the effectiveness and feasibility of the proposed algorithm in solving the problem of MRPP, we compare the proposed GEI-BBO algorithm with GA (Bakdi et al. 2017), PSO (Kennedy and Eberhart 1995), ACO (Mavrovouniotis and Yang 2011), BBO (Simon 2011) and IWO/BBO (Khademi et al. 2017) algorithms, so as to comprehensively test the optimization performance of different algorithms. Specifically, we try to verify the path planning performance of different algorithms in three environments with different complexity. Circular obstacles, rectangular obstacles and labyrinth obstacles are designed in turn in the three environments, and the best path planning results of GEI-BBO in the three environments are shown in Fig. 6.  It can be seen intuitively from Fig. 6 that the path of the mobile robot planned by the proposed GEI-BBO algorithm rarely has unnecessary twists and turns. Specifically, in environment 1, the starting point coordinate is (0,0) and the ending point coordinate is (4,6). The obstacle consists of three circular obstacles, and the optimal path length is 7.5465. There are 12 rectangular obstacles in environment 2, the starting point coordinate is (0,0), the ending point coordinate is (20,20), and the optimal path length is 28.8721. Environment 3 is composed of more complex rectangular obstacles, the starting point coordinate is (0,0), the ending point coordinate is (45,45), and the optimal path length is 64.1673. Table  4 shows the path results of each algorithm running independently for 10 times, including best, worst, mean and standard deviation (SD) of path length. Meanwhile, the success rate of each algorithm in the 10 times algorithm is shown in Table 5.
Combining Tables 4 and 5, we can see that GEI-BBO has good performance in finding the optimal path. In environment 1, each algorithm can find its own optimal path, but compared with GA, PSO, ACO, BBO and IWO/BBO, GEI-BBO has the shortest path length, and the standard deviation is 0. In environment 2, except GEI-BBO, the success rate of other algorithms is less than 100%, and the optimal value and mean value obtained by GEI-BBO are optimal. In environment 3, GEI-BBO has the shortest path length of 64.1673, and the standard deviation is 0. Although other algorithms can also find the approximate optimal path with GEI-BBO, they are not very stable, and the success rate is low.
To make the experimental results more intuitive, we further add box plots as shown in Fig. 7. In environment 1, it can be seen that the variance of GA and PSO data is relatively large. Although the variance of ACO, BBO and IWO/BBO is 0, GEI-BBO has a smaller path length. In environment 2, the variance of GA and PSO is still larger, followed by BBO and IWO/BBO. The variance of ACO is similar to GEI-BBO, but its optimal path length is larger. In environment 3, except for ACO, the optimal path length of the other four algorithms is similar to GEI-BBO, but the variance is greater than GEI-BBO, especially GA. Generally, the path planned by the proposed GEI-BBO algorithm is better than the other five comparison algorithms. This is because the strategies of eigendecomposition-based migration, gradient descent and system search are introduced, which balances the capabilities of global search and local mining. Simulation experiments verify the feasibility and reliability of the combination of GEI-BBO and cubic spline interpolation method to solve the problem of MRPP.

Conclusion
In this paper, we proposed a new improved algorithm of BBO for MRPP, which combines the eigendecompositionbased migration, gradient descent and system search strategy. This method can effectively reduce the dependence of BBO on coordinate system and improve the local search ability. The comparison of 23 benchmark functions has proven the effectiveness of GEI-BBO. Combining GEI-BBO with cubic spline interpolation function, a method and fitness function for solving the obstacle avoidance and shortest path of mobile robot are constructed to solve the problem of MRPP. The simulation results also show that the proposed method has higher accuracy and success rate, which can prove the feasibility of the algorithm.
In the future, we will study how to use the new approaches to further enhance the performance of BBO, and apply the proposed algorithm to the more complex global path planning for mobile robots.
Funding This work was supported by the State Key Laboratory of Robotics (2019-O18) and the Fundamental Research Funds for the Central Universities (DUT20LAB114, DUT2018TB06).
Data availability Enquiries about data availability should be directed to the authors.

Conflict of interest
The authors have not disclosed any competing interests.