2.1. Gravitational search algorithm

Gravitational search algorithm, proposed by Rashedi et al. in 2009, is inspired by the law of gravity and motion [Reference Rashedi, Nezamabadi-pour and Saryazdi38, Reference Jiao, Chen, Xin, Li, Zheng and Zhao39]. The solutions are viewed as particles, interacting with each other and moving in $D$ dimensional search space, whose performance is estimated by the mass. The better the particle, the greater mass. Under the action of universal gravitation, the particles will move toward the particle with a large mass, approaching the optimal solution. The flowchart of GSA is presented in Fig. 2.

Figure 2. The flowchart of GSA.

Assume that there are $N$ particles in $D$ dimensional space, whose positions and velocities are described as Eqs. (1) and (2).

(1) \begin{equation} X_{i}(t)=\left(x_{i}^{1}(t),x_{i}^{2}(t),\ldots,x_{i}^{d}(t),\ldots,x_{i}^{D}(t)\right), i=1,2,\ldots,N\end{equation}

(2) \begin{equation} V_{i}(t)=\left(v_{i}^{1}(t),v_{i}^{2}(t),\ldots,v_{i}^{d}(t),\ldots,v_{i}^{D}(t)\right), i=1,2,\ldots,N\end{equation}

The mass $M_{i}(t)$ of particle $X_{i}(t)$ is related to its fitness, for a minimization problem, mass can be calculated by Eqs. (3)-(6):

(3) \begin{equation} M_{i}(t)=\frac{m_{i}(t)}{\sum _{j=1}^{N}m_{j}(t)}\end{equation}

(4) \begin{equation} m_{i}(t)=\frac{fit_{i}(t)-fit_{\textit{worst}}(t)}{fit_{best}(t)-fit_{\textit{worst}}(t)} \end{equation}

(5) \begin{equation} fit_{best}(t)=\mathit{\min }_{i\in \left\{1,\ldots,N\right\}} fit_{i}(t) \end{equation}

(6) \begin{equation} fit_{\textit{worst}}(t)=\mathit{\max }_{i\in \left\{1,\ldots,N\right\}} fit_{i}(t) \end{equation}

where $t$ stands for the number of iterations and $fit_{i}(t)$ represents the fitness of the particle $X_{i}(t)$ .

As the role of gravity, at each generation $t$ , the force $F_{ij}(t)$ between particle $X_{i}(t)$ and $X_{j}(t)$ can be defined as follows:

(7) \begin{equation} F_{ij}^{d}(t)=G(t)\frac{M_{i}(t)\cdot M_{j}(t)}{R_{ij}(t)+\varepsilon }\left(x_{j}^{d}(t)-x_{i}^{d}(t)\right) \end{equation}

where $R_{ij}(t)$ denotes the Euclidean distance between the particle $X_{i}(t)$ and $X_{j}(t) ; \varepsilon$ is a constant, preventing the dominator from being 0. $G(t)$ represents the gravitational constant, which can adjust the gravity and be described as Eq. (8):

(8) \begin{equation} G(t)=G_{0}e^{-\alpha \frac{t}{T}} \end{equation}

where $G_{0}$ refers to the initial gravitational constant, $\alpha$ is the adaptive adjustment coefficient, and $T$ represents the maximum iterations.

To improve the performance of the algorithm, the total force acting on the particle $X_{i}(t)$ from the top $K$ best particles is defined as follows:

(9) \begin{equation} F_{i}^{d}(t)=\sum _{j\in K_{best}, j\neq i}rand_{j}F_{ij}^{d}(t) \end{equation}

where $K_{best}$ is the set storing the top $K$ best particles, decreasing linearly over iterations; $rand_{j}$ stands for a random value between 0 and 1.

According to the law of motion, the acceleration $a$ of the particle $X_{i}(t)$ can be achieved by Eq. (10), using its mass and resultant force.

(10) \begin{equation} a_{i}^{d}(t)=\frac{F_{i}^{d}(t)}{M_{i}(t)} \end{equation}

Finally, the position $x_{i}^{d}(t)$ and velocity $v_{i}^{d}(t)$ of the $i\mathrm{th}$ particle can be updated by Eqs. (11) and (12), respectively.

(11) \begin{equation} v_{i}^{d}\!\left(t+1\right)=rand_{i}\cdot v_{i}^{d}(t)+a_{i}^{d}(t) \end{equation}

(12) \begin{equation} x_{i}^{d}\!\left(t+1\right)=x_{i}^{d}(t)+v_{i}^{d}\!\left(t+1\right) \end{equation}

where $rand_{i}$ is a random value within the interval [0, 1].

2.2. Tent map

The tent map is one of the noninvertible and piecewise linear discrete maps [Reference Yang, Liu, Wang and Yang40], which is given by Eq. (13).

(13) \begin{equation} ch_{t+1}=f\!\left(ch_{t}\right)=\begin{cases} \dfrac{ch_{t}}{a}, 0\leq ch_{t}\leq a\\[10pt] \dfrac{1-ch_{t}}{1-a}, a\lt ch_{t}\leq 1 \end{cases} \end{equation}

where $a$ is referred to as the control parameter, and $ch_{t}$ represents the output chaotic sequence. Figure 3 characters the orbit of the tent map during the iteration, and the maximal value is obtained at $ch_{t}=0.7$ .

Figure 3. The tent map, in a case at $a=0.7$ .

3.1. Modeling the environment

Suppose, $S$ and $T$ are the start and terminal points of the task for a UAV, whose corresponding coordinates are $(x_{1}, y_{1},z_{1})$ and $(x_{D}, y_{D},z_{D})$ , respectively. The coordinates of $S$ and $T$ in 3D space are the main objects to be considered, just as shown in Fig. 4, ignoring the obstacles in the environment. $D-1$ equal parts are got by evenly dividing $[x_{1},x_{D}]$ along the $X$ -axis, the vertical planes $(\beta _{1},\beta _{2},\ldots,\beta _{d},\ldots,\beta _{D})$ that are perpendicular to $X$ -axis are attained according to the points $(x_{1},x_{2},\ldots,x_{d},\ldots,x_{D})$ in the $X$ -axis. The point $P_{d}=(x_{d},y_{d},z_{d})$ is taken from the plane $\beta _{d}$ . A path in three-dimensional space can be achieved by connecting these points $(P_{1},\ldots,P_{d},\ldots,P_{D})$ in turn, which can be described as $P=\{S,(x_{2},y_{2},z_{2}),\ldots,(x_{d},y_{d},z_{d}),\ldots,(x_{D-1},y_{D-1},z_{D-1}),T\}$ . Generally, the division method in $X$ -axis can also be revised to fulfill our needs.

Figure 4. The three-dimensional environmental model.

Figure 5. Performance constraints of UAV.

The optimal path problem can be transformed into optimizing a coordinate sequence to satisfy constraints and minimize the objective function in 3D space. Since some coordinates are set in advance and too large variable dimension would magnify the computational complexity, the path variable $P$ can be expressed as $P=(y_{2},z_{2},y_{3},z_{3},\ldots,y_{d},z_{d},\ldots,y_{D-1},z_{D-1})$ , reducing $D-2$ variables.

3.2. Cost function

When performing a task, time and energy consumption are two important factors to be considered, which are bound up with the path length. The greater the path length, the more time and energy required, and vice versa. So, the path length is a very critical factor in the optimal path, which can be described as follows:

(14) \begin{equation} F_{1}=\sum _{d=1}^{D-1}F_{1}^{d} \end{equation}

(15) \begin{equation} F_{1}^{d}=\sqrt{\left(x_{d}-x_{d+1}\right)^{2}+\left(y_{d}-y_{d+1}\right)^{2}+\left(z_{d}-z_{d+1}\right)^{2}} \end{equation}

where $F_{1}^{d}$ denotes the Euclidean distance between the node $d$ and its next node $d+1$ , $(x_{d},y_{d},z_{d})$ and $(x_{d+1},y_{d+1},z_{d+1})$ are the coordinates of node $d$ and its next node $d+1$ , respectively.

Due to the performance constraints of UAV, the UAV will lose control when the yaw angle or pitch angle is too large. Therefore, the smoothness of the path needs to be taken into account, just as shown in Fig. 5, the yaw angle and pitch angle are two primary factors to affect the smoothness and feasibility of the path. $P_{d-1}$ , $P_{d}$ , and $P_{d+1}$ are three adjoining points, whose corresponding projection points are $P_{d-1}^{\prime}$ , $P_{d}^{\prime}$ , and $P_{d+1}^{\prime}$ in the plane $OXY$ , respectively, and the projection of point $P_{d}$ on the line $P_{d+1}P_{d+1}^{\prime}$ is $P_{d+1}^{\prime\prime}$ . The yaw angle $\gamma _{d}$ , indicating the turning of the UAV in the horizontal plane, denotes the angle between $\overrightarrow{P_{d-1}^{\prime}P_{d}^{\prime}}$ and $\overrightarrow{{P_{d}^{\prime}P_{d+1}^{\prime}}}$ , which can be calculated by Eqs. (16) and (17). So the turning cost function could be computed by Eq. (19).

(16) \begin{equation} \gamma _{d}=\mathit{\arccos }\!\left(\frac{\overrightarrow{P_{d-1}^{\prime}P_{d}^{\prime}}\cdot \overrightarrow{P_{d}^{\prime}P_{d+1}^{\prime}}}{\left\| \overrightarrow{P_{d-1}^{\prime}P_{d}^{\prime}}\right\| \cdot \left\| \overrightarrow{P_{d}^{\prime}P_{d+1}^{\prime}}\right\| }\right) \end{equation}

(17) \begin{equation} \gamma _{d}=\mathit{\arccos }\\ \left(\frac{\left(x_{d}-x_{d-1},y_{d}-y_{d-1}\right)\cdot \left(x_{d+1}-x_{d},y_{d+1}-y_{d}\right)^{T}}{\left\| x_{d}-x_{d-1},y_{d}-y_{d-1}\right\| \cdot \left\| x_{d+1}-x_{d},y_{d+1}-y_{d}\right\| }\right) \end{equation}

(18) \begin{equation} F_{2}^{d}=\begin{cases} 0, if \gamma _{d}\leq \gamma _{max} \\[5pt] k_{t}\cdot \gamma _{d}, \textit{otherwise} \end{cases} \end{equation}

(19) \begin{equation} F_{2}=\sum _{d=2}^{D-1}F_{2}^{d} \end{equation}

where $\| \cdot \|$ is the Euclidean distance, $k_{t}$ and $\gamma _{max}$ are the coefficient and threshold of yaw angle, respectively, and $(x_{d-1},y_{d-1},z_{d-1})$ denotes the coordinates of the point $P_{d-1}$ .

The pitch angle, denoting the climbing angle in the vertical plane, represents the angle between the projection $\overrightarrow{P_{d}P_{d+1}^{\prime\prime}}$ of $\overrightarrow{P_{d}P_{d+1}}$ on the horizontal plane and $\overrightarrow{P_{d}P_{d+1}}$ , which can be calculated using Eqs. (20) and (21). The corresponding climbing cost function could be described by Eq. (23).

(20) \begin{equation} \theta _{d}=\mathit{\arctan }\left(\frac{\left| z_{d+1}-z_{d}\right| }{\left\| \overrightarrow{P_{d}^{\prime}P_{d+1}^{\prime\prime}}\right\| }\right) \end{equation}

(21) \begin{equation} \theta _{d}=\mathit{\arctan }\!\left(\frac{\left| z_{d+1}-z_{d}\right| }{\left\| \left(x_{d}-x_{d+1},y_{d}-y_{d+1}\right)\right\| }\right) \end{equation}

(22) \begin{equation} F_{3}^{d}=\begin{cases} 0, if \theta _{d}\leq \theta _{max}\\[5pt] k_{c}\cdot \theta _{d}, \textit{otherwise} \end{cases} \end{equation}

(23) \begin{equation} F_{3}=\sum _{d=1}^{D-1}F_{3}^{d} \end{equation}

where $| \cdot |$ is the absolute value operation, $k_{c}$ and $\theta _{max}$ stand for the coefficient and threshold of pitch angle, respectively.

A UAV cannot touch buildings or break through the maximum altitude of flight, and the higher the flight altitude, the greater the possibility of being monitored by radar. Considering the constraints of geography, radar, and UAV performance, it is significant for a UAV to fly only within the feasible altitude range. Assume that $h_{min}$ and $h_{max}$ are the maximum and minimum feasible flight altitudes, respectively. Then the height cost function can be described by Eq. (25):

(24) \begin{equation} F_{4}^{d}=\begin{cases} 0, \quad if h_{min}^{\prime}\leq z_{d}\leq h_{max}\\[5pt] k_{h}\cdot max\!\left(z_{d}-h_{max}, h_{min}^{\prime}-z_{d}\right), \textit{otherwise} \end{cases} \end{equation}

\begin{equation*} h_{min}^{\prime}=h_{min}+h\!\left(x_{i},y_{i}\right) \end{equation*}

(25) \begin{equation} F_{4}=\sum _{d=1}^{D}F_{4}^{d} \end{equation}

where $z_{d}$ is the current flight altitude or the coordinate of UAV in the $Z$ -axis, and $k_{h}$ denotes the flight height coefficients of the UAV.

In conclusion, the overall cost function can be represented by Eq. (26), considering the path length, UAV performance, and environmental factors for the path $P$ .

(26) \begin{equation} F\!\left(P\right)=\sum _{t=1}^{4}F_{t}\!\left(P\right) \end{equation}

3.3. Enhanced gravitational search algorithm for UAV path planning

As we know, exploration denotes the process of accessing the entire search space, while exploitation is the process that visiting the neighborhood of the most promising points, and these two processes have a great influence on the population-based algorithm with evolutionary behavior [Reference Chen, Xin, Peng, Dou and Zhang41–Reference Vikas and Parhi43]. The particles in particle swarm optimization, having memory function, can employ the social information to move toward the current optimal particle, which can promote optimization. Since the current optimal solution is uncertain to be the global optimal solution, perturbation is necessary, aiding the particle to jump out of the local optimal solution. Therefore, an EGSA is proposed, where each particle has the memory ability and random disturbance ability to improve the performance of the algorithm.

Levy flight is a non Gaussian random process, whose essence is a random walk derived from Levy stable distribution. The step length is vital for Levy flight, which determines whether a particle can jump out of the local solution. Since the proper step length is intractable to set, chaos is a good idea, considering randomness, ergodicity, and sensitivity to the initial value for chaos. Therefore, the combination of chaos and Levy flight is applied as a random disturbance, boosting the balance ability of exploration and exploitation. The velocity Eq. (11) can be rewritten as Eq. (27).

(27) \begin{equation} v_{i}^{d}\!\left(t+1\right)=rand_{i}\cdot v_{i}^{d}(t)+c_{1}(t)\cdot a_{i}^{d}(t)+\\ c_{2}(t)\cdot \left(\textit{gbest}^{d}-x_{i}^{d}(t)\right)+c_{3}(t)\cdot ch_{t}\cdot Levy\!\left(\beta \right) \end{equation}

(28) \begin{equation} c_{1}(t)=-5\cdot \left(\frac{t}{T}\right)^{3}+5 \end{equation}

(29) \begin{equation} c_{2}(t)=2.6\cdot \left(\frac{t}{T}\right)^{2} \end{equation}

(30) \begin{equation} c_{3}(t)=2.58\cdot \left(\frac{t}{T}\right)^{10} \end{equation}

(31) \begin{equation} Levy\!\left(\beta \right)=\frac{\mu }{\left| \nu \right| ^{1/\beta }} \end{equation}

where $\textit{gbest}^{d}$ denotes the $d\mathrm{th}$ dimension of the current optimal solution; $c_{1}$ , $c_{2}$ , and $c_{3}$ are named as the coefficients of acceleration, memory, and disturbance, respectively; $\mu$ and $\nu$ follow the normal distribution as follows:

\begin{equation*} \mu \sim N\!\left(0,1\right) \end{equation*}

\begin{equation*} \nu \sim N\!\left(0,\sigma _{\varsigma }^{2}\right) \end{equation*}

(32) \begin{equation} \sigma _{\varsigma }=\left(\frac{\Gamma \!\left(1+\beta \right)\mathit{\sin } \!\left(\pi \beta /2\right)}{\Gamma \!\left(\left(1+\beta \right)/2\right)\beta 2^{\left(\beta -1\right)/2}}\right)^{1/\beta }\end{equation}

where $\Gamma$ is referred to as the standard gamma function.

In summary, the pseudocode of our proposed EGSA for UAV path planning in 3D space is shown in Algorithm 1.

3.4. Computational complexity

Here, for the convenience of complexity analysis, let $D$ and $N$ be the dimension of decision variables and the population size. For EGSA, the mass calculation, acceleration computing, update of velocity and position, and evaluation consume the main calculation, whose maximum computational complexity is $O(ND)$ in one loop. Therefore, the overall computational complexity of EGSA is $O(ND)$ for one main cycle.

Algorithm 1. EGSA for UAV path planning

4.1. Parameter setting

To verify the performance of the proposed algorithm, let EGSA compare with moth flame optimization algorithm (MFO) [Reference Mirjalili44], GSA [Reference Rashedi, Nezamabadi-pour and Saryazdi38], chaotic kbest gravitational search algorithm (CKGSA) [Reference Mittal, Pal, Kulhari and Saraswat45], gbest guided gravitational search algorithm (GGSA) [Reference Mirjalili and Lewis46], gaussian particle swarm optimization and gravitational search algorithm (GPSOGSA) [Reference Kumar and John47], hybrid gravitational search particle swarm optimization algorithm (HGSPSO) [Reference Khan and Ling48], and logarithmic kbest gravitational search algorithm (LKGSA) [Reference Mittal and Saraswat49], and the parameters of these eight algorithms are listed in Table I, and the population size $N$ and the maximum number of iterations $T$ are 50 and 3000, respectively. The coefficients $k_{t}$ , $k_{c}$ , and $k_{h}$ of yaw angle, pitch angle, and flight height are defined to 20, 20, and 30, respectively. The thresholds $\gamma _{max}$ , $\theta _{max}$ , $h_{min}$ , and $h_{max}$ of yaw angle, pitch angle, minimum flight height, and maximum flight height are set to 60, 45, 5, and 200, respectively. To be fair, all the experiments are run on the computer with an Intel Core i7 2.6 GHz and Windows 10 operating system.

Table I. Parameter settings for eight algorithms.

4.2. Experimental results on CEC 2020 benchmark functions

Here, the performance of EGSA is verified on CEC 2020 with dimension $D=20$ , compared with MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, and LKGSA [Reference Yue, Price, Suganthan, Liang, Ali, Qu, Awad and Biswas50]. These eight algorithms were run 30 times independently, and several indicators, including the best results, mean results, and standard deviation (Std) of the optimal values, are applied to evaluate the performance of these algorithms. Furthermore, the Wilcoxon rank sum test at a 0.05 significance level is employed between EGSA and the other seven algorithms, and the symbols ‘+’, ‘−’, and ‘∼’ denote that the performance of the corresponding algorithm is significantly better, worse, and similar than that obtained by EGSA on the benchmark function, respectively [Reference Li, Guo, Lerch, Li, Wang and Wu51, Reference Derrac, García, Molina and Herrera52]. The experimental results on CEC 2020 obtained by these eight algorithms are recorded in Table II and Fig. 6.

Figure 6. The convergence curves of mean values obtained by EGSA and compared algorithms on CEC 2020 benchmark functions.

Table II. Best values, mean values, standard deviations, and p-values of the minimum values obtained by EGSA and compared algorithms on CEC 2020 benchmark functions.

From the view of the mean value in Table II, EGSA reveals a better performance than the other seven algorithms on most of these benchmark functions, winning the best place except CEC1 and CEC8. For function CEC1, the best performance is obtained by GSA, in fact, EGSA behaves similarly to GSA from the perspective of the Wilcoxon rank sum test. For function CEC8, GGSA performs best, but the performance of EGSA is not bad and ranks second among that displayed by eight algorithms, maybe the random disturbance cannot improve the performance of EGSA on CEC8. At the same time, the convergence curves in Fig. 6 verified the excellent performance of EGSA in CEC 2020 again.

As denoted by the symbols of “+/−/∼” at the bottom of Table II, EGSA performs better than MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, and LKGSA on 10, 8, 8, 7, 10, 10, and 7 benchmark functions and gets similar results with GSA, CKGSA, GGSA, and LKGSA on 2, 2, 2, and 3 out of the 10 test instances. As shown by the average ranking, the average ranking of EGSA is better than that of comparison algorithms, demonstrating that EGSA owns a better balance between exploitation and exploration and has superiority over the peer algorithms.

4.3. The effect of acceleration, memory ability of particle, and random disturbance

To investigate the influence of acceleration, memory ability of particle, and random disturbance, EGSA is compared with six variants on CEC 2020, meanwhile, explaining the rationality of parameter setting. The coefficients $c_{1}$ , $c_{2}$ , and $c_{3}$ are set to the minimum and maximum in the range of variation, respectively. The changed parameters and functions of six variants about EGSA are listed in Table III, and all the other parameters remain unchanged. Table IV and Fig. 7 reveal the experimental results attained by six variants and EGSA.

Table III. The changed parameter and function of six variants about EGSA.

Table IV. Best values, mean values, standard deviations, and p-values of the minimum values obtained by six variants and EGSA on CEC 2020 benchmark functions.

Figure 7. The convergence curves of mean values obtained by six variants and EGSA on CEC 2020 benchmark functions.

Our proposed algorithm EGSA achieves the top two on eight benchmark functions and wins the first place on six benchmark functions. For the functions CEC1 and CEC6, the EGSA-Variant5 and EGSA-Variant4 win the championship, respectively, however, these two algorithms do not perform well in the other five functions, which implies that the performance of EGSA cannot be promoted by overstrengthening the memory ability of particle or abandoning the random disturbance. For the algorithms EGSA-Variant2, EGSA-Variant3, and EGSA-Variant6, the performance of these three algorithms seem to be deteriorating, which implies that acceleration and memory ability of particle can evoke the performance, but large random disturbance is not desirable. For EGAS-Variant1, boosting the effect of acceleration, the performance is a little worse, which illustrates that overemphasized acceleration cannot facilitate performance.

The results of the Wilcoxon test, shown in the last row of Table IV, indicate that EGSA is superior to six variants in 1, 10, 7, 6, 6, and 8 out of the 10 test functions and performs similarly to six variants in 9, 0, 3, 3, and 2 out of the 10 test functions. The average ranking, put on the penultimate row of Table IV, explains that the performance of EGSA surpasses that of six variants, and the convergence curves in Fig. 7 declare the outstanding performance of EGSA again. Therefore, the acceleration, memory ability of particle, and random disturbance can affect the performance, which should not be too weakened or strengthened.

4.4. Experimental results in the simple environment

Assume that there is no more than one obstacle in an environment, which is named as a simple environment. The coordinates of the starting point and end point for a UAV are defined as $(1,5,10)$ and $(201,190,70)$ , respectively. The paths achieved by eight algorithms are shown in Fig. 8, and the relevant convergence curves are represented in Fig. 9. Intuitively, in Fig. 8, the path got by EGSA is much better than that obtained by MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, and LKGSA in terms of smoothness and length. From the perspective of convergence curves, just as presented in Fig. 9, EGSA converges faster than the other seven algorithms and wins the lowest overall cost value in the end, which means that EGSA has a good balance between exploration and exploitation for this question. So EGSA can perform better than MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, and LKGSA in this simple environment. As a matter of fact, we discover that $F_{2}$ , $F_{3}$ , and $F_{4}$ in the overall cost function are all equal to 0 for the last obtained paths, which implies that the proposed method can obtain the path satisfying the three constraints and the path length becomes the main factor affecting the optimality at last.

Figure 8. Path planning results obtained by eight algorithms in the simple environment; blue triangle dotted line, blue pentagram solid line, green triangle dotted line, green pentagram solid line, purple triangle dotted line, purple pentagram solid line, red triangle dotted line, and red pentagram solid line are acquired by MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, LKGSA, and EGSA, respectively.

Figure 9. Convergence curves of eight algorithms in the simple environment.

Figure 10. Path planning results obtained by eight algorithms in the complex environment 1; blue triangle dotted line, blue pentagram solid line, green triangle dotted line, green pentagram solid line, purple triangle dotted line, purple pentagram solid line, red triangle dotted line, and red pentagram solid line are acquired by MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, LKGSA, and EGSA, respectively.

Figure 11. Convergence curves of eight algorithms in the complex environment 1.

Figure 12. Path planning results obtained by eight algorithms in the complex environment 2; blue triangle dotted line, blue pentagram solid line, green triangle dotted line, green pentagram solid line, purple triangle dotted line, purple pentagram solid line, red triangle dotted line, and red pentagram solid line are acquired by MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, LKGSA, and EGSA, respectively.

Figure 13. Convergence curves of eight algorithms in the complex environment 2.

Figure 14. Path planning results obtained by eight algorithms in the complex environment 3; blue triangle dotted line, blue pentagram solid line, green triangle dotted line, green pentagram solid line, purple triangle dotted line, purple pentagram solid line, red triangle dotted line, and red pentagram solid line are acquired by MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, LKGSA, and EGSA, respectively.

Figure 15. Convergence curves of eight algorithms in the complex environment 3.

Figure 16. Path planning results obtained by eight algorithms in the complex environment 4; blue triangle dotted line, blue pentagram solid line, green triangle dotted line, green pentagram solid line, purple triangle dotted line, purple pentagram solid line, red triangle dotted line, and red pentagram solid line are acquired by MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, LKGSA, and EGSA, respectively.

Figure 17. Convergence curves of eight algorithms in the complex environment 4.

4.5. Experimental results in the complex environment

Suppose that plenty of obstacles exist in the environment, like hills and buildings, which is considered as a complex environment. In complex environment 1, the obstacles are scattered, while in complex environments 2, 3, and 4, the obstacles are relatively concentrated. The starting point and terminal point settings are identical to those in the simple environment. For the complex environment 1, the paths and convergence curves obtained by eight algorithms are drawn in Fig. 10 and Fig. 11. The path obtained by EGSA is shorter and smoother than those obtained by the other seven algorithms. The convergence rate of EGSA is faster than those of comparison algorithms and its corresponding cost value is the smallest among the eight algorithms. For the complex environments 2, 3, and 4, Fig. 12, Fig. 13, Fig. 14, Fig. 15, Fig. 16, and Fig. 17 represent the paths and convergence curves obtained by eight algorithms, respectively. Intuitively, the paths obtained by MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, and LKGSA are longer than those obtained by EGSA. At the same time, from the views of convergence speed and cost value, the performance of EGSA is much better than that of the other seven algorithms, and EGSA stands out as the top performer. Therefore, the performance of EGSA is the best among those eight algorithms in these four complex environments. The path length, like that in the simple environment, also becomes the primary contributor affecting the optimality for the last obtained path.

Table V. Comparison of simulation results of eight algorithms.

Figure 18. Boxplots of the results obtained by eight algorithms in different environments.

Figure 19. Average convergence curves of EGSA with different population size in the simple environment.

To objectively and accurately verify the performance of the proposed algorithm, the eight algorithms have been run 30 times in four kinds of environments, and Table V presents the mean, Std, and p-vlaue of the overall cost values. Obviously, the mean and Std of EGSA are much smaller than the results of MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, and LKGSA. From the Wilcoxon test, our proposed EGSA has significant advantages over other comparison algorithms. Meanwhile, it can be seen that chaotic levy flight could regulate exploration and exploitation compared with GGSA. The data distribution characteristics are presented in Fig. 18. The boxplots of the proposed algorithm, with the lowest values, are very narrow compared to the results of the other seven algorithms. Hence, the proposed EGSA has much better effectiveness and robustness than MFO, GSA, CKGSA, GGSA, GPSOGSA, HGSPSO, and LKGSA for path planning in simple and complex environments.

Figure 20. Boxplots of the results obtained by EGSA with different population size in the simple environment.

Figure 21. Average convergence curves of EGSA with different population size in the complex environment 1.

Figure 22. Boxplots of the results obtained by EGSA with different population size in the complex environment 1.

Figure 23. Average convergence curves of EGSA with different population size in the complex environment 2.

Figure 24. Boxplots of the results obtained by EGSA with different population size in the complex environment 2.

Figure 25. Average convergence curves of EGSA with different population size in the complex environment 3.

Figure 26. Boxplots of the results obtained by EGSA with different population size in the complex environment 3.

Figure 27. Average convergence curves of EGSA with different population size in the complex environment 4.

Figure 28. Boxplots of the results obtained by EGSA with different population size in the complex environment 3.

Figure 29. The mean running time of EGSA with different population size $N$ .

4.6. Parameter analysis

As we know, the population size $N$ can affect the performance of algorithm. If $N$ is too small, it will harm the convergence of the population to the optimal solution. On the contrary, too large $N$ will consume more computing resource. Here, to observe the influence of parameter $N$ on our proposed algorithm, all the parameters in EGSA remain the same except for the population size $N$ , the population $N$ is set to 10, 30, 50, and 70 in simple and complex environments, respectively. Figure 19, Fig. 21, Fig. 23, Fig. 25, and Fig. 27 show the average convergence with iteration in different environments, respectively. The convergence rates are very similar when the populations are 30, 50, and 70, respectively. The boxplots with diverse population sizes for three environments are revealed in Fig. 20, Fig. 22, Fig. 24, Fig. 26, and Fig. 28, respectively. Figure 29 describes the running time for our proposed path planning method in three environments, showing that the running time increases with the population size. Table VI represents the mean and Std of overall cost values for various environments with different population sizes. It is obvious that the performance of EGSA becomes better with the increase of population size $N$ , but in both cases of $N=50$ and $N=70$ , the performance of EGSA is not very different, the performance improvement is very small for $N$ increasing from 50 to 70. So, a compromise between the obtained path length and the running time is made, and 50 is set as the size of the population.