Hybridizing Particle Swarm Optimization and Differential Evolution for the Mobile Robot Global Path Planning

Abstract

Global path planning is a challenging issue in the filed of mobile robotics due to its complexity and the nature of non-deterministic polynomial-time hard (NP-hard). Particle swarm optimization (PSO) has gained increasing popularity in global path planning due to its simplicity and high convergence speed. However, since the basic PSO has difficulties balancing exploration and exploitation, and suffers from stagnation, its efficiency in solving global path planning may be restricted. Aiming at overcoming these drawbacks and solving the global path planning problem efficiently, this paper proposes a hybrid PSO algorithm that hybridizes PSO and differential evolution (DE) algorithms. To dynamically adjust the exploration and exploitation abilities of the hybrid PSO, a novel PSO, the nonlinear time-varying PSO (NTVPSO), is proposed for updating the velocities and positions of particles in the hybrid PSO. In an attempt to avoid stagnation, a modified DE, the ranking-based self-adaptive DE (RBSADE), is developed to evolve the personal best experience of particles in the hybrid PSO. The proposed algorithm is compared with four state-of-the-art evolutionary algorithms. Simulation results show that the proposed algorithm is highly competitive in terms of path optimality and can be considered as a vital alternative for solving global path planning.

1 Introduction

Over the past few decades, mobile robotics has been successfully applied in industry, military and security environments to perform crucial unmanned missions such as planet exploration, surveillance and landmine detection [39]. A fundamental issue of mobile robotics is path planning, which aims to find an optimal or sub-optimal obstacle-free path from a starting location to a destination location, while optimizing some performance criteria. Due to its widespread application in different areas such as computer games, animation, assembly planning and computational biology, path planning has drawn increasing research interest since the mid-1960s [36].

Generally, path planning can be categorized into global path planning (GPP) and local path planning (LPP), with respect to whether full knowledge of the environment is available or not [39]. GPP, executed offline, determines an obstacle-free path while being aware of complete information pertaining to the environment [35, 39]. LPP, also called online path planning, generates an efficient path by using no or partial information of the environment [35, 39]. This paper focuses on the GPP problem, which is good at generating an efficient path offline. For simplicity, the GPP problem is referred to as the path planning problem throughout the rest of this paper.

To date, path planning has been intensively studied in the literature and numerous path planning methods have been proposed. However, owing to its non-deterministic polynomial-time hard (NP-hard) nature, efficiently solving the path planning problem remains a challenge [5, 6, 14, 21, 35, 36, 39]. Thanks to their population-based nature and good search ability for producing high-quality solutions within a tractable time, even for complex problems [39], evolutionary algorithms (EAs) such as the genetic algorithm [1], simulated annealing [28] and artificial neural network [9] have been proposed in the field of path planning in recent years.

As one of the most powerful EAs, due to its simplicity and high convergence speed, the particle swarm optimization (PSO) algorithm has been widely implemented to solve path planning [12, 19, 32, 36]. Nevertheless, since the basic PSO algorithm cannot balance exploration and exploitation sufficiently, it is easily trapped in a local optimum [2]. Furthermore, the basic PSO may also suffer from stagnation when no particle can find a position that is better than its previous best positions [16]. In order to enhance the performance of PSO, it is essential to overcome these two shortcomings.

The differential evolution (DE) algorithm is another promising evolutionary algorithm that emerged as a highly competitive algorithm more than a decade ago [4]. Due to its simplicity, reliability, high performance and easy implementation, a considerable amount of research has been dedicated to applying DE to solve different optimization problems [8, 20, 38]. It is well-known that different generation strategies (i.e., mutation and crossover operators) and control parameters (i.e., scaling factor and crossover rate) significantly influence the performance of DE [8, 20, 38]. However, how to set generation strategies and control parameters are problem-dependent, which results in difficulties when designing an efficient DE algorithm [8, 38].

Since both PSO and DE are EAs that deal with population evolution, it is natural to mix these two algorithms together to leverage their advantages when developing an integrated method for efficiently solving optimization problems. Motivated by this idea, some research such as [11, 16, 17, 24, 33, 34] has been devoted to hybridizing PSO and DE in order to improve the performance of hybrid algorithms.

In this study, we also combine PSO with DE and propose a hybrid PSO-DE algorithm, the HNTVPSO-RBSADE, for mobile robot global path planning. In the proposed hybrid PSO algorithm, a novel PSO, called NTVPSO, is first developed to update the velocities and positions of particles. Then, after each iteration, a modified DE algorithm, called RBSADE, is proposed to evolve the personal best positions of particles in order to force particles to jump out of stagnation. In order to sufficiently balance the exploration and exploitation capabilities of particles, a new self-adaptive mechanism is proposed to adaptively adjust the three control parameters (inertia weight and cognitive and social acceleration parameters) of particles in NTVPSO. To improve the performance and easily adjust the control parameters of RBSADE, a ranking-based mutation operator and a self-adaptive strategy are developed in RBSADE.

To verify the proposed algorithm for solving mobile robot global path planning, it is compared to JADE [35], the time-varying particle swarm optimization (TVPSO) method [15], the gravitational search (GS) method [18] and the modified genetic algorithm (mGA) [25] under different simulation scenarios. For rigorous verification, path optimality and computation time are examined and compared for each algorithm. The simulation results confirm that the proposed algorithm outperforms the algorithms it is compared to in terms of path optimality. Moreover, the computation time of the proposed algorithm is also comparable with that of the other algorithms.

The remainder of the paper is organized as follows. Section 2 formulates the path planning problem and introduces the scheme for encoding the path, as well as the constraint handling technique. Section 3 presents the HNTVPSO-RBSADE algorithm and analytical investigations pertaining to this algorithm. Section 4 presents the framework for applying HNTVPSO-RBSADE to solve mobile robot global path planning. Simulations and comparisons are performed in section 5. Section 6 summarizes the paper by way of a conclusion and options for future work.

2 Statement of the Path Planning Problem

2.1 Modelling of the workspace

This paper focuses on finding an obstacle-free path for a mobile robot working in a 2-D complex workspace. The model presented in [13] is adopted to establish the workspace, due to its simplicity and small sensitivity to the shapes of obstacles. As shown in Fig. 1, a global coordinate system o - xy is established, where st and ta denote the start location and the destination location of a robot, respectively. The x -axis of the global coordinate system coincides with line st - ta and the y -axis of the global coordinate system is perpendicular to line st - ta. In the global coordinate system, line st - ta is equally divided into n + 1 subsections by n points, where n is a predefined parameter. After drawing n vertical lines l₁, l₂,…, l _n , as shown in Fig. 1, a point-to-point path, denoted as ph = [st, p₁, p₂,…,p _n , ta], is constructed by random sampling on vertical lines l₁, l₂,…, l _n . All obstacles in the workspace are assumed to be static and are represented by polygons. The robot is considered a mass point and obstacles are enlarged according to the size of the robot to compensate [9, 39].

Figure 1.

The illustration of the 2-D workspace

2.2 Formulation of the path planning problem

As searching for the best path is often associated with finding the shortest and safest path, as in [9, 13, 18, 36, 37], path length and path safety are considered as being two performance criteria of the path planning problem in this paper. Let p₀ and p_n+1 denote the start location and the destination location of a robot, respectively. The path planning problem can be mathematically formulated as [18, 37]:

$$\{\begin{array}{c}\text{finding: }ph=[st=p_0\mathrm{,}{p}_1\mathrm{,}{p}_2\mathrm{,\dots ,}{p}_n\mathrm{,}ta=p_{n+1}]\\ \text{minimize: }{J}_{cost}=w_1\cdot J_L+w_2\cdot \frac{1}{J_S}\\ \text{s.t. }{p}_i{p}_{i+1}\in \text{semi}-\text{free workspace}\mathrm{,}0\le i\le n\end{array}$$

(1)

where J _L and J _S denote path length and path safety, respectively, and w₁ and w₂ are two weighting parameters indicating the relative importance of J _L and J _S (w₁ = w₂ = 1 in this paper).

To generate a point-to-point path, J _L is calculated as [6]:

$$J_L={\displaystyle \sum _{i=0}^n}dis(p_i\mathrm{,}{p}_{i+1})$$

(2)

where dis(p _i , p_i+1) denotes the Euclidean distance between waypoint p _i and waypoint p_i+1.

J _L is the summation of the minimum distance of each path segment from the nearest obstacle along the path, which is calculated as [6]:

$$J_S={\displaystyle \sum _{i=0}^n}\underset{1\le j\le Nob}{{\displaystyle \min }}\{Mindis(p_i{p}_{i+1}\mathrm{,}{O}_j)\}$$

(3)

where Nob denotes the number of obstacles. Mindis(p _i p_i+1,O _j ) represents the minimum distance between path segment p _i p_i+1 and obstacle O _j .

2.3 Path encoding scheme

As stated in section 2.1, a path is constructed by a set of waypoints p₁, p₂,…, p _n , while a set of lines l₁, l₂,…, l _n determines the x -axis values of these points. Since the set of lines l₁, l₂,…, l _n is given beforehand during the construction of the workspace, values of p₁,…, p _n lie only with the y -axis values decided by l₁, l₂,…, l _n . These y -axis values, denoted as (y _{p 1} ,y _{p 2} ,…,y _{p n} ), are applied in order to encode the path [13, 36]. The following saturation strategy is used to modify y _{p i} (i = 1,2,…,n) when y _{p i} is outside of the workspace as follows:

$$y_{p_i}=\{\begin{array}{c}\frac{width}{2}\mathrm{,}\text{ if }{y}_{p_i}>\frac{width}{2}\\ \frac{-width}{2}\mathrm{,}\text{ if }{y}_{p_i}<\frac{-width}{2}\\ y_{p_i}\mathrm{,}\text{ otherwise}\end{array}$$

(4)

where width represents the width of the workspace as shown in Fig. 1.

2.4 Handling constraint and evaluating path

From (1), it can be seen that the path planning problem is a constrained optimization problem. In order to solve the constraint problem easily and efficiently, the task of how to handle the constraint must be addressed. As the constraint of the path planning problem is to generate an obstacle-free path, it is reasonable to evaluate the constraint violation degree of each candidate path by counting its collision times with obstacles [36]. Given Nob obstacles, the total constraint violation degree of candidate path m₁ is calculated as [36]:

$$\begin{array}{l}violatio{n}_{m_1}=\frac{1}{Nob}{\displaystyle \sum _{j=1}^{Nob}}vio{l}_{m_{1}j}\hfill \\ vio{l}_{m_{1}j}=\{\begin{array}{cc}1& \text{if }{m}_1\text{ collides with obstacle }j\\ 0& \text{otherwise}\end{array}\hfill \end{array}$$

(5)

After calculating the constraint violation degree and the fitness value of each candidate path, the feasibility-based rule [7, 27] is used to evaluate and select the elite path between any two candidate paths. The feasibility-based rule is described as: (1) for any two paths with the same constraint violation degree, the path with better fitness is preferred; (2) for any two paths with different constraint violation degrees, the path with a smaller constraint violation degree is preferred.

Since the fitness and constraint violation information of each candidate path are considered separately in the feasibility-based rule, no additional parameter is needed to transform the constrained path planning problem into an unconstrained one, which decreases optimization difficulty. Additionally, although the non-feasible paths collide with some obstacles, they may also contain high-quality path segments. When these high-quality path segments are considered, the diversity of the paths and the possibility of finding high-quality paths are increased. This is one main reason why non-feasible solutions are considered in the feasibility-based rule.

3 The Optimization Method

3.1 Review of the basic PSO

Inspired by birds flocking and fish schooling, Kennedy and Eberhart first proposed PSO in 1995 [10]. Each particle in PSO represents a potential solution to an optimization problem and is associated with a velocity that is dynamically adjusted according to its own flight experience, as well as those of its companions. Therefore, each particle is attracted to a stochastically weighted average of its personal best position and the global best position of the swarm. In the basic PSO algorithm, from iteration k to iteration k + 1, each particle updates its velocity and position as follows [10]:

$${\mathrm{V}}_m^{k+1}=w{\mathrm{V}}_m^k+c_1{r}_1\left(\mathrm{p}bes{\mathrm{t}}_m^k-{\mathrm{X}}_m^k\right)+c_2{r}_2\left(\mathrm{g}bes{\mathrm{t}}^k-{\mathrm{X}}_m^k\right)$$

(6)

$${\mathrm{X}}_m^{k+1}={\mathrm{X}}_m^k+{\mathrm{V}}_m^{k+1}$$

(7)

where V ^k _m and X ^k _m denote the velocity and position vectors of particle m at iteration k; w is a real coefficient denoting the inertia weight; c₁ and c₂ are two positive real coefficients representing the cognitive and social acceleration parameters; r₁ and r₂ are two random numbers uniformly distributed in [0,1]; pbest ^k _m denotes the personal best position of particle m at iteration k and gbest ^k denotes the global best position of the swarm at iteration k.

3.2 Review of the conventional DE

DE is a population-based algorithm that initializes a population of SP individuals in a D-dimensional search space. Each individual, known as either a genome or chromosome, denotes a potential solution to an optimization problem. After mutation, crossover and selection operators, an offspring is allowed to the next generation only if it improves on the fitness of the parent [23].

The mutation operator modifies an individual via random changes in order to generate a new offspring, which aims to add some diversification into the population to avoid the local optimum [23]. At iteration k, for each parent x _i (k), a trial vector v _i (k) is generated by mutating three (or more) individuals as follows:

$${\mathrm{v}}_i\left(k\right)={\mathrm{x}}_{i_3}\left(k\right)+{\mathrm{F}}_i\left({\mathrm{x}}_{i_1}\left(k\right)-{\mathrm{x}}_{i_2}\left(k\right)\right)$$

(8)

where x _{i 1} (k), x _{i 2} (k) and x _{i 3} (k) are randomly selected from the current swarm with i ≠ i₁ ≠ i₂ ≠ i₃ and F _i denotes the scaling vector that is applied to control the amplification of the differential variation.

The crossover operator follows a discrete recombination approach, where elements from the parent vector x _i (k) are combined with those from the trial vector v _i (k) to produce a new offspring u _i (k). The aim of the crossover operator is to explore new areas of the search space [23]. Using the binomial recombination, each element in u _i (k) is obtained as follows:

$$u_{ij}(k)=\{\begin{array}{l}{v}_{ij}(k)\mathrm{,}\text{ if }rn{d}_{ij}[\mathrm{0,1}]\le CR\text{ or }j=j_{rnd}\hfill \\ x_{ij}(k)\mathrm{,}\text{ otherwise}\hfill \end{array}$$

(9)

where j(j = 1,2,…,D) refers to a specific dimension; x _ij , v _ij (k) and u _ij (k) are the jth elements of x _i (k), v _i (k) and u _i (k), respectively; j _rnd indicates a random integer in [1,D]; rnd _ij is a random number uniformly distributed in 0,1 and CR represents the crossover rate.

The selection operator is a one-to-one spawning strategy in which the new offspring u _i (k) competes against the parent x _i (k) and the better one will be selected for the next generation. In a minimization optimization problem, the selection operator is given as follows:

$${\mathrm{x}}_i(k+1)=\{\begin{array}{l}{\mathrm{u}}_i(k)\mathrm{,}\text{ if }f({\mathrm{u}}_i(k))\le f({\mathrm{x}}_i(k))\hfill \\ {\mathrm{x}}_i(k)\mathrm{,}\text{ otherwise}\hfill \end{array}$$

(10)

where f(u _i (k)) and f(x _i (k)) denote the fitness values of u _i (k) and x _i (k), respectively.

3.3 The proposed HNTVPSO-RBSADE

Since the proposed HNTVPSO-RBSADE algorithm incorporates NTVPSO and RBSADE, NTVPSO and RBSADE are first presented separately below. Then, the HNTVPSO-RBSADE-based framework is presented in this subsection.

3.3.1 Nonlinear time-varying PSO (NTVPSO)

When using PSO to solve optimization problems, it is necessary to properly control the exploration and exploitation abilities of PSO in order to efficiently find an optimal solution [2, 22]. Ideally, on the one hand, in the early stage of the evolution, the exploration ability of PSO needs to be strengthened so that particles can wander through the entire search space, rather than clustering around the current population-best solution [2, 22]. On the other hand, in the later stage of the evolution, the exploitation capability of PSO must be promoted so that particles can search carefully in a local region to find optimal solutions efficiently [2, 22].

Proverbially, the exploration and exploitation capabilities of PSO heavily depend on its three control parameters. The basic philosophies concerning how the three control parameters influence such abilities of PSO can be summarized as follows: (1) a large inertia weight enhances exploration, while a small inertia weight facilitates exploitation [22, 26]; (2) a large cognitive component, compared to the social component, results in the wandering of particles through the entire search space, which strengthens exploration [22, 26]; (3) a large social component, compared with the cognitive component, leads particles to a local search, which strengthens exploitation [22, 26].

Considering all of the concerns noted above, this paper proposes a novel PSO called NTVPSO. The main purpose of this PSO is to adaptively balance the exploration and exploitation capabilities of NTVPSO, so that the performance of the algorithm can be improved. In order to achieve this goal, a novel self-adaptive mechanism that is used to update the three control parameters of each particle in NTVPSO is proposed as follows:

$$w=\frac{w_{max}-w_{min}}{\exp ({\delta }_m\cdot {\beta }_m\cdot k)}+w_{min}$$

(11)

$$c_1=\frac{c_{1s}-c_{1f}}{\exp ({\delta }_{c1}\cdot {\beta }_m\cdot k)}+c_{1f}$$

(12)

$$c_2=\frac{c_{2s}-c_{2f}}{\exp (-{\delta }_{c2}\cdot {\beta }_m\cdot k)}+c_{2f}$$

(13)

$${\delta }_m=\frac{w_{max}-w_{min}}{k_{max}}$$

(14)

$${\delta }_{c1}=\frac{c_{1s}-c_{1f}}{k_{max}}$$

(15)

$${\delta }_{c2}=\frac{c_{2s}-c_{2f}}{k_{max}}$$

(16)

$${\beta }_m=\frac{\Vert {\mathrm{V}}_m^{k-1}\Vert +\Vert {\mathrm{V}}_m^k\Vert }{2\Vert {\mathrm{V}}_m^{k-1}\Vert +\mathrm{\Delta }}$$

(17)

where w _max and w _min denote the upper and lower bounds of w; c_1s and c_1f denote the initial and final values of c₁; c_2s and c_2f represent the initial and final values of c₂; k _max represents the maximum iteration number; ‖V^k-1 _m ‖ and ‖V ^k _m ‖ denote the L₂-norm of the velocity vector of particle m at iteration (k - 1) and iteration k. Δ is a very small positive value in case the denominator of β _m becomes zero (Δ = 1e - 25 in this paper). Note that w _max > w _min , c_1s > c_1f and c_2s<c_2f in NTVPSO. Furthermore, note that particles in NTVPSO update their velocities and positions based on (6) and (7).

3.3.2 Parametric analysis for NTVPSO

From (11) to (13), with increasing iteration number k, it is clear that w and c₁ decrease, while c₂ increases in NTVPSO. Therefore, according to the aforementioned basic philosophies, NTVPSO is likely to start the search with a high exploration tendency, which will be reduced over time, so that exploitation may be favoured in the late stage of the evolution. Note that following the update rule of a fixed β _m , the balance between exploration and exploitation varies only with respect to the iteration number k.

We also adapt the balance of the search in NTVPSO using an additional parameter β _m . From (11) and (12), it is trivial that w and c₁ decrease as β _m increases. On the other hand, the variation in c₂ becomes larger as β _m increases, as shown in (13). This implies that, for large β _m , the exploration capability of NTVPSO tends to be retained, compared to the case where there is small β _m . In contrast, the exploitation ability of the algorithm takes over the exploration ability more rapidly as β _m decreases.

Briefly, by utilizing the proposed self-adaptive strategy defined from (11) to (17) in NTVPSO, the three control parameters of particles are adaptively adjusted, complying with the basic philosophies of PSO improvement. Therefore, the NTVPSO algorithm is expected to improve the ability for finding high-quality solutions. Fig. 2 demonstrates the tendency of these changes in the three control parameters with respect to different values of β _m . Note that w _max = 0.9, w _min = 0.1, c_1s = c_2f = 2.5, c_1f = c_2s = 0.5 and k _max = 100 in Fig. 2.

3.3.3 Ranking-based self-adaptive DE (RBSADE)

In the conventional DE algorithm, the three different individuals used for the mutation operator are randomly selected from the current swarm. However, it is natural that good species always contain good genetic information, making them more likely to be utilized for mutating better offspring. Based on this consideration, this paper develops a modified DE called RBSADE, in which a ranking-based strategy is developed to select the good species. In the ranking-based strategy, all parent individuals at the current iteration are first sorted in ascending order according to their fitness values and constraint violation degrees. Then, at the current iteration, the ranking value of the i th parent individual is assigned as follows:

Figure 2.

Changes of the three control parameters under different β _m in NTVPSO

$$R{V}_i=SP-i$$

(18)

where i =1,2,…,SP and SP denotes the size of the swarm.

After calculating the ranking value of each parent individual, the selection probability that the ith individual is allowed to participate in the mutation operator is calculated as follows:

$$p_i=\frac{R{V}_i}{SP}$$

(19)

It is clear from the ranking-based strategy that the better genetic information a parent contains, the higher the probability that it will be selected for participating in the mutation operator and the more likely a better offspring will be produced. After calculating the selection probability of each parent individual at the current iteration, the roulette-wheel mechanism is applied to select three different individuals from the current swarm to mutate a new offspring according to the mutation operator given in (8).

In order to decrease the difficulties in controlling the crossover rate, the self-adaptive strategy proposed in [31] is applied to dynamically update the crossover rate of each individual in RBSADE as follows:

$$C{R}_i(k+1)=\{\begin{array}{l}C{R}_i(k)\mathrm{,}\text{ if }f({\mathrm{u}}_i(k))\le f({\mathrm{x}}_i(k))\hfill \\ N(\mathrm{0.5,0.1})\mathrm{,}\text{ otherwise}\hfill \end{array}$$

(20)

where N(0.5,0.1) is a random number generated by a normal distribution of mean 0.5 and standard deviation 0.1. When f(u _i (k)) ⩽ f(x _i (k)) for individual i, it means the current crossover rate of the individual has a higher possibility to improve the quality of the candidate solutions at the next iteration. Hence, keeping the current crossover rate may obtain better solutions [31]. When the trial vector u _i (k) is worse than the current vector x _i (k), it implies that the current crossover rate of the individual has a lower chance of generating better candidate solutions. Therefore, changing the current crossover rate may be more suitable for the evolution process [31].

To adaptively adjust the scaling vector F _i , a new self-adaptive strategy is developed to update each element of F _i by using a diversity measure as follows:

$$F_{i\mathrm{,}j}=1-di{v}_{i\mathrm{,}j}(k)$$

(21)

$$di{v}_{i\mathrm{,}j}(k)=\frac{{\displaystyle \sum _{i=1}^{SP}}\sqrt{E}}{SP\cdot max(\sqrt{E})}$$

(22)

$$E={\displaystyle \sum _{j=1}^D}{\left[x_{i\mathrm{,}j}\left(k\right)-{\bar{x}}_j\left(k\right)\right]}^2$$

(23)

$${\bar{x}}_j\left(k\right)=\frac{1}{SP}{\displaystyle \sum _{i=1}^{SP}}{x}_{i\mathrm{,}j}\left(k\right)$$

(24)

where F _i,j is the j th element of scaling vector F _i ; SP denotes the size of the swarm; div _i,j (k) represents the normalized diversity for the j th dimension of the i th parent individual at iteration k; ${\bar{x}}_j(k)$ is the average of the j th dimension over all parent individuals. One advantage of the developed self-adaptive strategy is that the mutation step-size F _i (x _{i 1} (k)-x _{i 2} (k)) can be adaptively controlled based on the diversity measure, without adding any additional parameters, which consequently releases the burden of controlling scaling vector F _i .

3.3.4 The framework of the HNTVPSO-RBSADE algorithm

Let NP denote the size of the swarm and BP = {pbest ^k ₁,…,pbest ^k _NP } define the personal best experience set of particles at iteration k. The algorithmic scheme of HNTVPSO-RBSADE is illustrated in Table 1.

Table 1.

The algorithmic scheme of HNTVPSO-RBSADE

1. Randomly generate an initial swarm 2. Gain gbest and BP for the initial swarm 3. while not exit condition do 4. Evaluate and sort BP 5. for i = 1: NP do 6. Calculate selection probability p i for pbest i in BP by (18) and (19) 7. end for 8. for i = 1: NP do 9. Select three different individuals for pbest i from BP by the roulette-wheel 10. Mutation crossover to generate v i by (8) 11. Crossover operator to generate u i by (9) 12. Selection operator to obtain better pbest i by (10) 13. Update CR i by (20) 14. Update F i by (21)–(24) 15. Update velocity of particle i by (6) 16. Update position of particle i by (7) 17. Update w, c1 and c2 of particle i by (11)–(17) 18. end for 20. Update gbest and BP for the current swarm 21. end while

3.4 Convergence analysis of HNTVPSO-RBSADE

When designing a PSO algorithm, the convergence of the algorithm is of great importance [3, 29]. Since each dimension of the velocity and position vectors of each particle is updated independently from one another in (6) and (7), HNTVPSO-RBSADE can be simplified into a one-dimensional case in order to analyse its convergence. Without loss of generality, by omitting the subscript m in (6) and (7) for simplicity, the one-dimensional HNTVPSO-RBSADE can be rewritten into matrix form as follows:

$$\left[\begin{array}{c}X(k+1)\\ V(k+1)\end{array}\right]=\mathrm{A}\left[\begin{array}{c}X(k)\\ V(k)\end{array}\right]+\mathrm{B}P$$

(25)

where:

$$\mathrm{A}=\left[\begin{array}{cc}1-\phi & w\\ -\phi & w\end{array}\right]$$

(26)

$$B={\left[\phi \phi \right]}^{\mathrm{T}}$$

(27)

$${\phi }_1=c_1{r}_1$$

(28)

$${\phi }_2=c_2{r}_2$$

(29)

$$\phi ={\phi }_1+{\phi }_2$$

(30)

$$P=\frac{{\phi }_1\cdot pbest+{\phi }_2\cdot gbest}{{\phi }_1+{\phi }_2}$$

(31)

Solving |λE - A| = 0 where E is the identity matrix with the same size as A, the characteristic equation of the dynamic system (25) is derived as follows:

$${\lambda }^2-\left(w+1-\phi \right)\cdot \lambda +w=0$$

(32)

where two roots, denoted as λ_1,2, are obtained as follows:

$${\lambda }_{\mathrm{1,2}}=\frac{1+w-\phi \pm \sqrt{{\left(1+w-\phi \right)}^2-4w}}{2}$$

(33)

In the context of dynamic system theory, the necessary and sufficient condition for the convergence of system (25) is that magnitudes of λ₁ and λ₂ should be less than 1 [30]. Thus, the system (25) converges, if and only if:

$$\text{Max}\left\{\left|{\lambda }_1\right|\mathrm{,}\left|{\lambda }_2\right|\right\}<1$$

(34)

Since it appears from (33) that λ₁ and λ₂ are two real or complex numbers, we will discuss two cases separately where λ_1,2 are two real and complex numbers in order to analyse the convergence of system (25).

1. The case where λ_1,2 are two complex numbers, denoted as λ_1,2 ∈ ℂ.

   Lemma 1. For the system (25), λ_1,2 ∈ ℂ, if and only if:

   $$\{\begin{array}{l}w-2\sqrt{w}+1<\phi <w+2\sqrt{w}+1\hfill \\ w\ge 0\hfill \end{array}$$

   (35)

   Proof. For system (25), it is obvious that:

   $${\lambda }_{\mathrm{1,2}}\in C\iff {\left(1+w-\phi \right)}^2-4w<0$$

   (36)

   Solving the right-hand side of (36) using the classic approach, Lemma 1 can be easily proven. □

   Now, let us find conditions on φ and w guaranteeing the convergence of system (25) in the case where λ_1,2 ∈ ℂ. It is trivial that the system converges if and only if Max{|λ₁|,|λ₂|}<1.

   Lemma 2. In the case where λ_1,2 ∈ ℂ, the system (25) converges, if and only if:

   $$\{\begin{array}{l}w-2\sqrt{w}+1<\phi <w+2\sqrt{w}+1\hfill \\ 0\le w<1\hfill \end{array}$$

   (37)

   Proof. Note that the magnitude of a complex number Z can be calculated as |Z| = √Z² _r + Z² _c , where Z _r and Z _c denote the real and imaginary parts of Z. Thus, for λ_1,2 ∈ ℂ, it is clear that:

   $$\begin{array}{l}\text{Max}\left\{\left|{\lambda }_1\right|\mathrm{,}\left|{\lambda }_2\right|\right\}=\left|{\lambda }_1\right|=\left|{\lambda }_2\right|=\sqrt{w}\hfill \end{array}$$

   (38)

   Therefore:

   $$\text{Max}\left\{\left|{\lambda }_1\right|\mathrm{,}\left|{\lambda }_2\right|\right\}<1\iff \sqrt{w}<1$$

   (39)

   In the case where λ_1,2 ∈ ℂ, according to Lemma 1, (35) must be satisfied. Hence, when considering the conditions λ_1,2 ∈ ℂ and Max{|λ₁|,|λ₂|}<1, for λ_1,2 ∈ ℂ, the system (25) converges, if and only if:

   $$\{\begin{array}{l}w-2\sqrt{w}+1<\phi <w+2\sqrt{w}+1\hfill \\ 0\le w<1\hfill \end{array}$$

   (40)

   Fig. 3 shows the convergent region of the system (25), i.e., HNTVPSO-RBSADE, in the case where λ_1,2 ∈ ℂ.

   Figure 3.

   Convergent region of HNTVPSO-RBSADE in the case where λ_1,2 ∈ ℂ

2. The case where λ₁ and λ₂ are two real numbers, denoted as λ_1,2 ∈ ℝ.

Lemma 3. For the system (25), λ_1,2 ∈ ℝ, if and only if:

$$\{\begin{array}{ll}\phi \in \mathrm{R}\mathrm{,}w<0\hfill & \hfill \\ \phi \le w-2\sqrt{w}+1\text{ or }\phi \ge w+2\sqrt{w}+\mathrm{1,}w\ge 0\hfill & \hfill \end{array}$$

(41)

Proof. For system (25), it is clear that:

$${\lambda }_{\mathrm{1,2}}\in \mathrm{R}\iff {\left(1+w-\phi \right)}^2-4w\ge 0$$

(42)

Solving the right-hand side of (42) using the classic approach, Lemma 3 can be easily proven. □

Now, let us find conditions on φ and w guaranteeing the convergence of the system (25) in the case where λ_1,2 ∈ ℝ. According to (33) and (34), in the case where λ_1,2 ∈ ℝ, Max{|λ₁|, |λ₂|}<1 holds, if and only if:

$$-1<\frac{1+w-\phi \pm \sqrt{{\left(1+w-\phi \right)}^2-4w}}{2}<1$$

(43)

Hence:

$$-3-w+\phi <\pm \sqrt{{\left(1+w-\phi \right)}^2-4w}<1-w+\phi$$

(44)

As λ_1,2 ∈ ℝ, it is clear that:

$$(44)\iff \{\begin{array}{l}-3-w+\phi <-\sqrt{{\left(1+w-\phi \right)}^2-4w}\hfill \\ \sqrt{{\left(1+w-\phi \right)}^2-4w}<1-w+\phi \hfill \end{array}$$

(45)

Solving the right-hand inequalities in (45), yields:

$$(44)\iff \{\begin{array}{l}2w+2-\phi >0\hfill \\ \phi >0\hfill \end{array}$$

(46)

According to Lemma 3, (41) must be satisfied in the case where λ_1,2 ∈ ℝ. Considering both the conditions that λ_1,2 ∈ ℝ and Max{|λ₁|, |λ₂|}<1, in the case where λ_1,2 ∈ ℝ, the system (25) converges, if and only if:

$$\{\begin{array}{l}0<\phi <2w+\mathrm{2,}-1<w<0\hfill \\ 0<\phi \le w-2\sqrt{w}+1\text{ or }w+2\sqrt{w}+1\le \phi <2w+\mathrm{2,}0\le w<1\hfill \end{array}$$

(47)

Considering both cases where λ_1,2 ∈ ℂ and λ_1,2 ∈ ℂ together, the system (25), i.e., HNTVPSO-RBSADE, converges, if and only if:

$$\{\begin{array}{ll}0<\phi <2w+2\hfill & \hfill \\ -1<w<1\hfill & \hfill \end{array}$$

(48)

Fig. 4 shows the convergent region of HNVPSO-RBSADE in both cases where λ_1,2 ∈ ℝ and λ_1,2 ∈ ℂ, which is a triangle area. Only if a parameter selection of w and φ is located in this area does HNTVPSO-RBSADE converge. Fig. 5 illustrates the 3-D representation of the value of Max{|λ₁|,|λ₂|}.

Figure 4.

Convergent region of HNTVPSO-RBSADE

Now, let us find the equilibrium point of HNTVPSO-RBSADE. Calculating limits on both sides of (25) yields the following:

$$\{\begin{array}{l}\underset{k\to \infty }{{\displaystyle \lim }}X\left(k+1\right)=w\underset{k\to \infty }{{\displaystyle \lim }}V(k)+\phi \underset{k\to \infty }{{\displaystyle \lim }}\left(P-X\left(k\right)\right)\hfill \\ \underset{k\to \infty }{{\displaystyle \lim }}V\left(k+1\right)=\underset{k\to \infty }{{\displaystyle \lim }}X(k)+\underset{k\to \infty }{{\displaystyle \lim }}V(k)\hfill \end{array}$$

(49)

When HNTVPSO-RBSADE converges, $\underset{k\to \infty }{{\displaystyle \lim }}X(k+1)=\underset{k\to \infty }{{\displaystyle \lim }}X(k)$ and $\underset{k\to \infty }{{\displaystyle \lim }}V(k+1)=\underset{k\to \infty }{{\displaystyle \lim }}V(k)$ . Therefore, substituting these two equations into (49) yields:

$$\{\begin{array}{l}\underset{k\to \infty }{{\displaystyle \lim }}X(k)=P=\frac{{\phi }_1\cdot pbest+{\phi }_2\cdot gbest}{{\phi }_1+{\phi }_2}\hfill \\ \underset{k\to \infty }{{\displaystyle \lim }}V(k)=0\hfill \end{array}$$

(50)

where φ = φ₁ + φ₂, φ₁ = c₁r₁ and φ₂ = c₂r₂. pbest denotes the personal best position of the particle and gbest represents the global best position of the swarm.

Figure 5.

3-D representation of the value of Max{|λ₁|, |λ₂|}

3.5 Convergence behaviour of particles in HNTVPSO-RBSADE

Prior to particles converging to the equilibrium point given in (50), they may oscillate in different ways around the equilibrium point as a result of the different values of w and φ. Since different convergence oscillations may influence the quality of the final solution found by particles [29], it is necessary to investigate the different convergence oscillations of particles. Four typical oscillations of particles in HNTVPSO-RBSADE are shown in Fig. 6.

Non-oscillatory behaviour, as shown in Fig. 6(a), leads particles to only search on one side of the equilibrium point, which will be useful when the search space is bounded. Particles exhibit non-oscillatory convergence behaviour when λ₁ and λ₂ are two real roots and at least one of them is positive, which is equivalent to 0⩽(1 + w - φ)²-4w and 0 < 1 + w - φ. Harmonic oscillation behaviour, as shown in Fig. 6(b), may be beneficial for the exploitation stage, since particles smoothly oscillate around the equilibrium point. Harmonic oscillation behaviour occurs when two roots, λ₁ and λ₂ are complex, that is, (1 + w - φ)²-4w < 0. Zigzagging convergence, as shown in Fig. 6(c), may also facilitate exploitation as particles zigzag around the equilibrium point, which may be useful when the search space is rugged. Particles exhibit zigzagging convergence behaviour when at least one of λ₁ and λ₂ has a negative real part, i.e., w < 0 or 1 + w - φ < 0. The combined harmonic with zigzagging behaviour, as shown in Fig. 6(d), can be beneficial for the transition from exploration to exploitation, due to its mixed nature, which emerges when at least one of the two complex roots λ₁ and λ₂ has a negative real part, i.e., (1 + w - φ)² - 4w < 0 ∩ w < 0 ∪ (1 + w - φ)² - 4w < 0 ∩ 1 + w - φ < 0.

If the boundaries of coefficients associated with these oscillations are known beforehand, one may easily design an adaptive method to change the values of these coefficients, so that the convergence of PSO can be guaranteed and the quality of the final solution searched by particles can be improved.

Figure 6.

Convergence behaviour of particles in HNTVPSO-RBSADE

3.6 Convergence of HNTVPSO-RBSADE without considering its stochastic nature

Due to the stochastic nature of HNTVPSO-RBSADE it is difficult to rigorously establish the exact relationship between its stochastic nature and its convergence property. Thus it is necessary to first analyse the convergence of HNTVPSO-RBSADE without considering its stochastic nature Note that the stochastic nature of HNTVPSO-RBSADE is attributed to the existence of the two random numbers r₁ and r₂.

Lemma 4 Without considering its stochastic nature HNTVPSO-RBSADE converges only if:

$$\{\begin{array}{l}2w+2>c_1+c_2\hfill \\ -1<w<1\hfill \\ c_1\mathrm{,}{c}_2>0\hfill \end{array}$$

(51)

Proof Since φ = φ₁ + φ₂ = c₁r₁ + c₂r₂ the necessary and sufficient condition given by (48) can be rewritten as:

$$\{\begin{array}{l}0<c_1{r}_1+c_2{r}_2<2w+2\hfill \\ -1<w<1\hfill \end{array}$$

(52)

Since c₁ and c₂ are two positive parameters, and r₁ and r₂ are two random numbers uniformly distributed in [0,1], we have c₁ ⩾ c₁r₁ and c₂ ⩾ c₂r₂. Therefore:

$$\{\begin{array}{l}2w+2>c_1+c_2\hfill \\ -1<w<1\hfill \\ c_1\mathrm{,}{c}_2>0\hfill \end{array}\Rightarrow \{\begin{array}{l}0<r_1{c}_1+r_2{c}_2<2w+2\hfill \\ -1<w<1\hfill \end{array}$$

(53)

Since the right-hand side inequality in (53) is the necessary and sufficient condition for the convergence of HNTVPSO-RBSADE, it is trivial that Lemma 4 holds. □

3.7 The parameter selection principle of HNTVPSO-RBSADE

After analytically investigating the convergence of HNTVPSO-RBSADE, we still need to discover how to set the initial and final values of w, c₁ and c₂ to guarantee the convergence of HNTVPSO-RBSADE. In the following lemma, an alternative parameter selection principle for how to set the initial and final values of w, c₁ and c₂ to guarantee the convergence of HNTVPSO-RBSADE is provided.

Figure 7.

Convergent position and velocity trajectories of the particle in HNTVPSO-RBSADE under the suggested parameter selection

Lemma 5. The HNTVPSO-RBSADE algorithm converges only if the initial and final values of its three control parameters satisfy the following:

$$\{\begin{array}{l}2w_{min}+2>c_{1s}+c_{1f}\hfill \\ 1>w_{max}>w_{min}>-1\hfill \\ c_{1s}=c_{2f}>c_{1f}=c_{2s}>0\hfill \end{array}$$

(54)

Proof. When c_1s = c_2f and c_1f = c_2s, it appears from (12), (13), (15) and (16) that c₁ + c₂ = c_1s + c_1f for any particle at any iteration in HNTVPSO-RBSADE. From (11) and (13), it is clear that w _m ⩽ w ⩽ w _max , c_1f ⩽ c₁ ⩽ c_1s and c_2s ⩽ c₂ ⩽ c_2f for any particle at any iteration in HNTVPSO-RBSADE. Therefore:

$$\{\begin{array}{l}2w_{min}+2>c_{1s}+c_{1f}\hfill \\ 1>w_{max}>w_{min}>-1\hfill \\ c_{1s}=c_{2f}>c_{1f}=c_{2s}>0\hfill \end{array}\Rightarrow \{\begin{array}{l}2w+2>c_1+c_2\hfill \\ -1<w<1\hfill \\ c_1\mathrm{,}{c}_2>0\hfill \end{array}$$

(55)

The right-hand side inequality in (55) is the sufficient condition for the convergence of HNTVPSO-RBSADE, which is proved in Lemma 4. This completes the proof for Lemma 5. □

Since w _max , w _min , c_1s, c_1f, c_2s and c_2f are predefined parameters, the convergent condition given by (54) can be easily satisfied by setting proper values of these parameters. Fig. 7 shows the convergent position and velocity trajectories of the particle in HNTVPSO-RBSADE under the suggested parameter selection: w _max = 0.9, w _min = 0.4, c_1s = c_2f = 2 and c_1f = c_2s = 0.1.

4 Apply HNTVPSO-RBSADE to Solve Global Path Planning

The algorithmic scheme of applying HNTVPSO-RBSADE to solve mobile robot global path planning is illustrated in Table 2. Note that the evolution will not exit until the iteration number reaches the maximum iteration number. In Table 2, BP = {pbest ^k ₁,…,pbest ^k _NP } denotes the personal best experience set of particles at iteration k.

Table 2.

The HNTVPSO-RBSADE-based framework for global path planning

Randomly generate an initial swarm Update gbest for the initial swarm Update BP for the initial swarm while not exit condition do Evaluate and sort BP by the feasibility-based rule Calculate selection probability for each individual in BP by (18) and (19) for m = 1: NP do Select three different individuals for pbest m from BP Mutation crossover to generate v m by (8) Crossover operator to generate u m by (9) Selection operator to get better pbest m by (10) Update CR m by (20) Update F m by (21)–(24) Update velocity vector of particle m by (6) Update position vector of particle m by (7) Modify position vector of particle m by (4) Calculate fitness value of particle m by (1)–(3) Calculate constraint violation of particle m by (5) Update w, c1 and c2 of particle m by (11)–(17) end for Update gbest and BP for the current swarm end while Output gbest to navigate the robot

5 Simulations and Analysis

In order to validate the proposed method, it is compared with JADE [35], TVPSO [15], GS [18] and mGA [25]. The simulation parameters for each method are shown in Table 3. To fully investigate the performance of the proposed method, four numerical simulations and a Monte-Carlo experiment with 50 independent runs are conducted. In each simulation, the optimal path of each method is obtained after 200 iterations of 40 particles using MATLAB 2012B software on a Windows 8 personal computer with i3-2350@2.30GHz and 2GB RAM. For rigorous comparison, the path optimality and the computation time of each method are examined and compared. Path optimality is measured by the fitness value of the optimal path. Note that this paper focuses on the numerical value of the path fitness, which is calculated by J _L + 1/J _S , as given in (1), where J _L and J _S denote the path length and path safety.

Table 3.

The simulation parameters for each method

Method	Simulation parameters
Proposed	w max = 0.9	w min = 0.4	c1s = c2f = 2	c1f = c2s = 0.1
TVPSO [15]	w max = 0.7	w min = 0.4	c1s = 5c1f = 2.5	c2f = 5c2s = 2.5
GS [18]	c1 = c2 = 0.5	c max = 5	c min = 1
mGA [25]	cr = 0.5
JADE [35]	c1 = 0.1

5.1 Numerical simulations

This subsection shows four numerical simulations that are all executed in a 100m x 100m workspace. From the first simulation to the last, there are, respectively, 8, 5, 17 and 11 obstacles in the workspace. Line st - ta is equally divided into 8, 12, 12 and 13 subsections from the first simulation to the last. The obtained paths of the five methods for the four simulations are demonstrated in Fig. 8. Fig. 9 displays the evolution curves of the five methods for the four simulations. The simulation results of the five methods for the four simulations are shown in Tables 4–7.

From Fig. 8, it is clear that each method is efficient at finding an obstacle-free path in each simulation. From Table 4 through to Table 7, it can be seen that Proposed method ≻ JADE ≻ GS ≻ mGA ≻ TVPSO in terms of path optimality and TVPSO ≻ Proposed method ≻ mGA ≻ GS ≻ JADE with respect to computation time. Here, “≻” denotes “dominate” or “outperform”. Although TVPSO is the fastest algorithm, it exhibits the worst performance at finding an optimal path. Though the proposed method is the second fastest algorithm, it is the best one at finding an optimal path. It is worth noting that the difference in computation time between TVPSO and the proposed method is around 2s in each simulation. Considering all of the concerns raised, it can be concluded that the proposed method performs superior to the other four methods in terms of establishing path optimality. Furthermore, the computation time of the proposed method is comparable with those of the other four methods.

From Table 4 through to Table 7, it is also interesting to note that the computation time of each method changes significantly from the first simulation to the last. For example, the computation time of the proposed method varies from 24.59s, 35.44s, 79.33s to 72.25s from the first simulation to the last. This is because the number of obstacles and the number of subsections of line st - ta are different in each simulation. Based on the computation time of a specific method from the first simulation to the last, as shown in Tables 4–7, it is found that when the workspace contains more obstacles and st - ta is divided into more subsections, the path planner will take a longer time to generate an efficient path. This is because when the workspace contains more obstacles and line st - ta is divided into more subsections, the complexity of the workspace and the number of waypoints along a candidate path are increased. Thus, the path planner will take longer to calculate the cost and the constraint violation degree of a candidate path, in addition to updating the position information of waypoints of a candidate path.

Table 4.

Simulation results for the first numerical simulation

Method	Optimal path length (m)	Optimal path safety (m)	Optimal path fitness	Computation time (s)
TVPSO	129.42	8.73	129.53	22.72
mGA	129.37	9.10	129.49	26.11
GS	127.10	11.47	127.18	28.02
JADE	125.27	10.81	125.36	30.12
Proposed	121.97	10.57	122.07	24.59

Table 5.

Simulation results for the second numerical simulation

Method	Optimal path length (m)	Optimal path safety (m)	Optimal path fitness	Computation time (s)
TVPSO	141.80	23.06	141.84	34.50
mGA	129.85	23.17	129.89	36.46
GS	127.46	21.27	127.51	38.55
JADE	125.68	20.86	125.73	40.41
Proposed	122.09	22.56	122.14	35.44

From Fig. 9, it can be observed that the cost curves of some methods for some numerical simulations initially rise in the early part of the evolution. This happens because, in this paper, the feasibility-based rule [7, 27] is used to evaluate and select the elite path between any two candidate paths. Since the initial swarm of each method is randomly generated, the initial global best path of a method may collide with some obstacles, which implies that the initial global best path of the method is not feasible. In this case, with the evolution continuing, if any particle in the method searches for a new solution that has a smaller constraint violation degree but a larger path fitness value than those of the initial global best path, the initial global best path will be replaced by the newly-found solution, based on the feasibility-based rule. Therefore, in this case, the cost curve of the method will initially rise in the early phase of the evolution. However, with the evolution continuing, the cost curves of all methods will either keep dropping or will remain unchanged in the late phase of the evolution, as shown in Fig. 9. This is because each method can search for a feasible global best path in the late phase of the evolution. The feasible global best path can only be replaced by feasible solutions that have smaller fitness values than those of the feasible global best path, according to the feasibility-based rule.

Figure 8.

Optimal paths of all methods for the four numerical simulations

Table 6.

Simulation results for the third numerical simulation

Method	Optimal path length (m)	Optimal path safety (m)	Optimal path fitness	Computation time (s)
TVPSO	143.07	21.04	143.12	77.35
mGA	139.44	17.92	139.49	81.24
GS	124.25	17.94	124.31	82.78
JADE	122.09	20.66	122.14	84.34
Proposed	118.61	15.58	118.68	79.73

Table 7.

Simulation results for the fourth numerical simulation

Method	Optimal path length (m)	Optimal path safety (m)	Optimal path fitness	Computation time (s)
TVPSO	137.98	17.08	138.03	71.23
mGA	131.56	15.50	131.63	74.50
GS	117.36	22.98	117.40	76.36
JADE	116.75	20.06	116.80	79.14
Proposed	113.91	17.75	113.97	72.25

5.2 Monte-Carlo experiment

A Monte-Carlo experiment with 50 independent runs is carried out to further study the performance of the proposed method in this subsection. In each run of the Monte-Carlo experiment, the workspace is a 100mx100m area that contains 18 obstacles. Each obstacle has four vertexes that are randomly and uniformly distributed in the workspace. Line st - ta is equally divided into 15 subsections in each run of the Monte-Carlo experiment. After performing 50 independent runs, the statistical results, with respect to the fitness value of the optimal path and the computation time of each method, are reported see Table 8 and Table 9. The results corresponding to the statistical results shown in Table 8 and Table 9 are visualized in Fig. 10.

Figure 9.

Evolution curves of all methods for the four numerical simulations

From Table 8, it is clear that Proposed method ≻ JADE ≻ GS ≻ mGA ≻ TVPSO in terms of the best, worst and average fitness value of the optimal path. Compared with TVPSO, the proposed method improves 9.75% on the average fitness value of the optimal path. From Table 9, it is obvious that TVPSO ≻ Proposed method ≻ mGA ≻ GS ≻ JADE in terms of the best, worst and average computation time. Compared to the proposed method, TVPSO improves 1.38% on the average computation time. However, despite the best performance of TVPSO concerning computation time, TVPSO has the worst performance in terms of path optimality. Although the proposed method ranks second with respect to computation time, it is the most efficient for producing optimal path. Since the difference in average computation time between the proposed method and TVPSO is less than 2s in the Monte-Carlo experiment, it can be concluded that the proposed method generally outperforms the other four methods in the context of finding the optimal path. In addition, the computation time of the proposed method is comparable with those of the other methods compared.

Table 8.

Statistical results of the fitness value for the optimal paths by different methods in the Monte-Carlo experiment

	Proposed	JADE	GS	mGA	TVPSO
Best	121.63	126.60	127.71	132.54	137.78
Wost	136.25	141.64	143.77	145.66	147.87
Average	128.86	134.35	135.85	139.33	142.78
Std.Dev.	4.3630	4.5799	4.7780	4.5306	3.2588

Table 9.

Statistical results of computation time by different methods in the Monte-Carlo experiment

	TVPSO	Proposed	mGA	GS	JADE
Best	102.31	104.09	106.46	107.24	109.28
Wost	105.34	106.10	109.28	110.56	112.30
Average	103.72	105.17	107.79	108.75	110.71
Std.Dev.	0.8610	0.5897	0.7786	0.9736	0.8851

Figure 10.

Simulation results for the fitness value of the optimal paths and the computation time of different methods in the Monte-Carlo experiment

6 Conclusion and Future Work

In this study, a hybrid PSO algorithm called HNTVPSO-RBSADE, which combines NTVPSO with RBSADE, is proposed to solve mobile robot global path planning. In the proposed algorithm, particles first follow the moving rules defined in NTVPSO to update their velocities and positions. Then, the RBSADE algorithm is developed to evolve the personal best positions of particles in order to avoid particles becoming stagnant. By fine-tuning the three control parameters, the exploration and exploitation capabilities of NTVPSO can be well-balanced. Consequently, NTVPSO promotes particles to search for high-quality paths. In order to enhance the performance and easily adjust the control parameters of RBSADE, a rank-based mutation operator and a self-adaptive strategy are developed in RBSADE. Since the convergence of HNTVPSO-RBSADE remains important, this paper first analytically investigates the convergence of HNTVPSO-RBSADE through standard results yielded by dynamic system theory. Then, a parameter selection principle is proposed to guarantee the convergence of HNTVPSO-RBSADE.

The proposed algorithm is tested using four numerical simulations and a Monte-Carlo experiment with multiple runs against four evolutionary algorithms. The simulation results reveal that the proposed algorithm outperforms the other four algorithms in terms of path optimality; furthermore, the computation time of the proposed algorithm is comparable with those of the other algorithms compared.

There are some issues that deserve future study. The first is how to theoretically analyse the optimality of the equilibrium point of HNTVPSO-RBSADE. Since the three control parameters also influence the convergence speed of the PSO algorithm, the second issue concerns the analysis of the sensitivity of the convergence speed of the proposed PSO algorithm to the three control parameters. Furthermore, since there does not exist a widely-accepted criterion for selecting an optimum number of waypoints along a candidate path [35], our work has to predefine the number of subsections that line st - ta is divided into, which remains problem-dependent. Thus, developing a strategy that can autonomously choose the number of dividing subsections for line st - ta can be considered as a future study. In order to generate a higher quality global path, more performance criteria, such as path length, path safety and the path smoothness must be considered when addressing the path planning problem in future research. Moreover, for further verifying the proposed method, we are considering comparing the proposed method with other non-evolutionary methods under more complex environments in the near future. Finally, we are also considering the possibility of extending the proposed method into 3-D path planning, since such an approach is more practical for real-world applications.