An improved dynamic cooperative random drift particle swarm optimization algorithm based on search history decision

Abstract

A novel dynamic cooperative random drift particle swarm optimization algorithm based on entire search history decision (CRDPSO) is reported. At each iteration, the positions and the fitness values of the evaluated solutions in the algorithm are stored by a binary space partitioning tree structure archive, which leads to a fast fitness function approximation. The mutation is adaptive and parameter-less because of the fitness function approximation enhancing the mutation strategy. The dynamic cooperation between the particles by using the context vector increases the population diversity helps to improve the search ability of the swarm and cooperatively searches for the global optimum. The performance of CRDPSO is tested on standard benchmark problems including multimodal and unimodal functions. The empirical results show that CRDPSO outperforms other well-known stochastic optimization methods.

Introduction

With the development of science and technology, many real-world engineering optimization problems become increasingly complex. Better and more efficient optimization algorithms are urgently needed. In recent years population-based optimization techniques motivated by evolution as seen in nature such as Genetic Algorithm (GA), Evolutionary Programming (EP), Evolution Strategies (ES) and Particle Swarm Optimization (PSO) have been successfully applied to solve various fields of science and technology. Particle Swarm Optimization (PSO) is a metaheuristic optimization algorithm introduced by Kennedy and Eberhart1 to solve optimization problems inspired by the collaborative manner of social animals searching for food, such as birds flocking and fish schooling. In PSO, each volume-less individual in the swarm is called a particle, and exchange of information between particles instead of the crossover and mutation operations used in other evolutionary algorithms. PSO has been widely used to solving a number of applications owing to its simplicity and efficiency.

However, the PSO has been demonstrated that it is not a global optimization algorithm.2 Over the past decades, numerous attempts have been proposed to improve the performance of the standard PSO on both exploitative and exploratory capability. To improve the global search ability of PSO, a variant of the PSO algorithm named random drift particle swarm optimization (RDPSO) was proposed by Sun3 which is inspired by the free electron model of metal conductors in an external electric field.4 RDPSO has been successfully used to obtain effective parameter estimation for complex biochemical dynamic systems and obtained better quality solutions than most global optimization methods for solving inverse problems. RDPSO algorithm can also be used as a general tool for optimization problems that widely exist in the real world. In literature,5 RDPSO is used to obtain the expected optimal cost of the thermal power and ESS operation under different scenarios, showing its better performance than the PSO algorithm in terms of running time, the costs of each scenario and the target expected cost. Like the PSO algorithm, RDPSO also faces premature convergence and falls into local search optimization problems, which lead to huge performance loss and suboptimal solutions. An improved version of the random drift particle swarm optimization (RDPSO) algorithm was proposed for solving the economic dispatch problem.6 The algorithm integrates a self-adaption mechanism and a crossover operation as well as a greedy selection process, with the mean best position of the particles replaced by the personal best position of each particle in the velocity updating equation.

In this work, we integrate RDPSO and the search history with dynamic cooperation shown to be efficient increase the diversity of particles and further improve RDPSO’s performance on complex optimization problems. Evolutionary search history has been used in the form of memory acted as adaptive guide search strategies7–11 However, only partial search history was used in these researches. That means a portion of the information obtained from the search history is retained, and the rest is discarded. The search history includes the operations performed, the location and the fitness values of the evaluated solutions, which are extremely useful information for improving the performance of the Swarm Intelligent Algorithms (SI). Especially, it can maintain the diversity of the swarm, guide the search direction or suggest promising prospective search areas. Most importantly, when the same optimal value reappears again in the search history, it can warn that the search may fall into local optimum. Non-revisiting strategy is integrated with intelligent algorithms12–18 to avoid the re-evaluation of the evaluated solution candidates by using binary space partitioning (BSP). The non-revisiting genetic algorithm (NrGA) is proposed to prevent solution re-evaluation14,15 in the genetic algorithm by using search history and cNrGA16 is another version designed for continuous variables. In cNrGA, a parameter-less mutation operation was introduced to increase the population diversity. Therefore, cNrGA is more robust than GA and has superior performance. However, due to the limitation of sub-regions, the adaptive mutation of cNrGA cannot make the gene fully explore. The HdEA in literature17 improved the disadvantages of the sub-regions portioned by BSP so as to enhance the mutation strategy. Although HdEA is based on the entire previous search history, the global search ability is still not strong enough due to the lack of cooperative operation. Moreover, whatever cNrGA or HdEA has complex computation and slow convergence speed because of the evolutionary operations of the GA algorithm. In literature,18 the ESH-QPSO is proposed to improve the performance of QPSO by overlapped sub-region and One-Particle-Flip mutation. However, fitness values are still not utilized, and the mutation operation is lack of guidance due to randomness. Thus, an improved dynamic cooperative RDPSO algorithm is proposed, which is adaptively mutated with overlapped sub-regions based on evolutionary search history. Since all updated positions are guaranteed to be novel, faster convergence speed is expected.

The rest of the paper is organized as follows. Section 2 introduces the original PSO and RDPSO. Section 3 describes the proposed CRDPSO in detail. Section 4 elaborated the experimental environment. Section 5 presents the simulation results and analysis. Finally, conclusions are given in Section 6.

Random drift particle swarm optimization

Particle swarm optimization

The proposal of the particle swarm optimization algorithm was put forward to solve optimization problems by developing computational simulations of the movement of social animals such as flocks of birds and schools of fish. Supposed in the PSO with a swarm size m, each particle is treated as a volume-less individual in the N-dimensional search space. The position vector and velocity vector of the $\mathit{ith}$ particle at the $\mathit{tth}$ iteration are represented as $X_i\left(t\right)=\left[x_{i,1}\right(t),x_{i,2}(t),\dots ,x_{i,n}(t\left)\right]$ and $V_i\left(t\right)=\left[v_{i,1}\right(t),v_{i,2}(t),\dots ,v_{i,n}(t\left)\right]$ respectively. The particle updates according to the equations:

$$v_{i,j}(t+1)=w\cdot v_{i,j}\left(t\right)+c_1\cdot r_{1,j}\left(t\right)\cdot [\mathit{pbes}{t}_{i,j}-x_{i,j}(t\left)\right]+c_2\cdot r_{2,j}\left(t\right)\cdot \left[\mathit{gbes}{t}_j\right(t)-x_{i,j}(t\left)\right]$$

(1)

$$x_{i,j}(t+1)=x_{i,j}\left(t\right)+v_{i,j}(t+1)$$

(2)

for $i=1\text{,}2,\dots ,m;j=1\text{,}2,\dots ,N$ , where $c_1$ and $c_2$ are acceleration coefficients indicated cognitive and social factors, respectively, and $w$ is the inertia weight used to balance the global and local search abilities.19 The vector $\mathit{pbes}{t}_i\left(t\right)=\left[\mathit{pbes}{t}_{i,1}\right(t),\mathit{pbes}{t}_{i,2}(t),\dots ,\mathit{pbes}{t}_{i,n}(t\left)\right]$ is the best historical position found by particle $i$ called personal best position (the position giving the best objective function value), and the vector $\mathit{gbest}\left(t\right)=\left[\mathit{gbes}{t}_1\right(t),\mathit{gbes}{t}_2(t),\dots ,\mathit{gbes}{t}_n(t\left)\right]$ is the position of the best particle in the swarm and called global best position. Without loss of generality, taking a class of minimization problems into account, $\mathit{pbes}{t}_i\left(t\right)$ can be updated by:

$$\mathit{pbes}{t}_i(t)=\{\begin{array}{ll}{x}_i(t)& \text{ if }f\left[x_i(t)\right]<f\left[\mathit{pbes}{t}_i(t-1)\right]\\ \mathit{pbes}{t}_i(t-1)& \text{ if }f\left[x_i(t)\right]\ge f\left[\mathit{pbes}{t}_i(t-1)\right]\end{array}$$

(3)

where $f$ is the objective function. Consequently, the $\mathrm{ }\mathit{gbest}\left(t\right)$ can be found by:

$$\mathit{gbest}\left(t\right)=\mathit{pbes}{t}_g\left(t\right),\mathrm{ }\mathrm{where} g=\mathit{arg}\underset{1\le i\le m}{\mathit{min}}\left\{f\left[\mathit{pbes}{t}_i\left(t\right)\right]\right\}$$

(4)

The parameters $r_{1,j}\left(t\right)$ and $r_{2,j}\left(t\right)$ are two sequences of random numbers distributed uniformly within the range $\left[0\text{,}1\right]$ , that is $\frac{r_{1,j}\left(t\right),r_{2,j}\left(t\right)}{U\left(0\text{,}1\right)}$ . Generally, the value of $v_{i,j}\left(t\right)$ is restricted in the interval [–V _max , V _max ].

The random drift particle swarm optimization

Trajectory analysis20 illustrated the fact the PSO algorithm's convergence may be attained supposing that each particle converges to its local attractor, $p_i^k=(p_{i,1}^k,p_{i,2}^k,\cdots p_{i,N}^k)$ defined at the coordinates

$$p_{i,j}^k=\frac{c_1{r}_{i,j}^k{P}_{i,j}^k+c_2{R}_{i,j}^k{G}_j^k}{c_1{r}_{i,j}^k+c_2{R}_{i,j}^k},1\le j\le N$$

(5)

which can be restated as

$$p_{i,j}^k={\varphi }_{i,j}^k{P}_{i,j}^k+\left(1-{\varphi }_{i,j}^k\right),1\le j\le N, \mathrm{where}\mathrm{ }{\varphi }_{i,j}^k=\frac{c_1{r}_{i,j}^k}{c_1{r}_{i,j}^k+c_2{R}_{i,j}^k}$$

(6)

In general, the acceleration coefficients c₁ and c₂ are set to be the same in the PSO algorithm so that ${\varphi }_{i,j}^k$ is a uniformly distributed random number over [0, 1], i.e. ${\varphi }_{i,j}^k\sim U\left(0\text{,}1\right)$ . In effect, as the particles are convergent to their own local attractors, their personal best positions, current position, the global best positions and local attractors all converge to one region, which leads to the convergence of PSO algorithm. The particle's directional motion toward to local attractor is similar to an election's drift movement in metal conductors in an external electric field.4 As reported by the free electron model, in addition to the directional movement originated from the electric field, the electron is also in thermal movement which seems to be random migration. The whole consequence of motion of the electron is that the electron careens towards the position of minimum potential energy. Obviously, if the location of an electron is regarded as a candidate solution and the potential energy function as the objective function of the problem, the electron's movement is very similar to the procedure of obtaining a minimal solution of the minimization problem.

Illuminated by the above facts, the particle in the swarm of the PSO behaving like an electron that moves in a metal conductor in an external electric field. Thus, the motion of the particles is superposed by the directional movement and the thermal movement toward $p_i^k$ . Consequently, the velocity of the particles can be depicted by $V_{i,j}^{k+1}=V{1}_{i,j}^{k+1}+V{2}_{i,j}^{k+1}$ , where $V{1}_{i,j}^{k+1}$ and $V{2}_{i,j}^{k+1}$ represents the velocities of the directional motion and thermal motion to $p_i^k$ , respectively. In the drift particle swarm optimization,3 $V{1}_{i,j}^{k+1}$ is the superposition of the social and the cognitive parts in equation (1) and the $V{2}_{i,j}^{k+1}$ is assumed to obey the Maxwell distribution, that is, Gaussian probability distribution. This velocity update equation seems to have some effects on the search performance of PSO, but there are some defects. First of all, the tail of the Gaussian distribution is very thin, so there is less chance of generating outliers. In consequence, the thermal movement of the particle owns less randomness which cannot significantly improve the global search ability of the particle. Secondly, the update equation of $V{1}_{i,j}^{k+1}$ can be represented as:

$$V{1}_{i,j}^{k+1}=c_1{r}_{i,j}^k(P_{i,j}^k-X_{i,j}^k)+c_2{r}_{i,j}^k(G_j^k-X_{i,j}^k)$$

(7)

The two random scaling factors $r_{i,j}^k$ add randomness to the movement of a particle in equation (7), which means that the position of the particle is sampled uniformly and distributed random within a super-rectangular shape around its personal best position and the global optimal position. Although this can guarantee the particle converge to its local attractor, it is unable to improve the overall search ability of particles because of the limited range of hyperrectangle, but may weaken the local search which is the responsibility of directional motion. Thus, the particle velocity decided by the combination of $V{1}_{i,j}^{k+1}$ and $V{2}_{i,j}^{k+1}$ may not able to achieve a good balance between the local search and the global search of the particle. A new method to determine $V{1}_{i,j}^{k+1}$ and $V{2}_{i,j}^{k+1}$ is proposed in literature.3

The velocity of directional motion $V{1}_{i,j}^{k+1}$ can be simply expressed as

$$V{1}_{i,j}^{k+1}=\beta (p_{i,j}^k-X_{i,j}^k)$$

(8)

where $p_{i,j}^k$ is determined by

$$\begin{array}{l}{p}_{i,j}^k={\varphi }_{i,j}^k{\mathit{pbest}}_{i,j}^k+(1-{\varphi }_{i,j}^k){\mathit{gbest}}^k,\hspace{1em}{\varphi }_{i,j}^k\sim U\left(0\text{,}1\right),\\ 1\le j\le N\end{array}$$

(9)

where ${\mathit{pbest}}_{i,j}^k$ is the personal best (pbest) of a particle and ${\mathit{gbest}}^k$ is the global best (gbest).

As for the velocity of thermal motion $V{2}_{i,j}^{k+1}$ obeys double exponential distribution whose probability density function and probability distribution function are

$$f_{V{2}_{i,j}^{k+1}}\left(v\right)=\frac{1}{{\sigma }_{i,j}^k}{e}^{\frac{-2\left|v\right|}{{\sigma }_{i,j}^k}}$$

(10)

and

$$F_{V{2}_{i,j}^{k+1}}\left(v\right)=1-e^{\frac{-2\left|v\right|}{{\sigma }_{i,j}^k}}$$

(11)

respectively, where v represents the value of the random variable $V{2}_{i,j}^{k+1}$ , ${\sigma }_{i,j}^k$ is the standard deviation of the distribution. $V{2}_{i,j}^{k+1}$ can be expressed as fellows by stochastic simulation method

$$V{2}_{i,j}^{k+1}=\frac{{\sigma }_{i,j}^k}{2}{\varphi }_{i,j}^k$$

(12)

where

$${\varphi }_{i,j}^k=\left\{\begin{array}{l}\begin{array}{ll}+\mathit{ln}\left(1/u_{i,j}^k\right)& \begin{array}{ll}\text{i}\text{f}& s>0.5\end{array}\end{array}\\ \begin{array}{lll}-\mathit{ln}\left(1/u_{i,j}^k\right)& \text{i}\text{f}& s\le 0.5\end{array}\end{array}\right.$$

(13)

In equation (13) s and $u_{i,j}^k$ are two different random numbers uniformly distributed on (0,1), i.e. $s,k\sim U\left(0\text{,}1\right)$ . It can be directly proved that $V_{i,j}^{k+1}=V{2}_{i,j}^{k+1}$ , when $k\to \infty$ , $X_{i,j}^k\to p_{i,j}^k$ . Therefore, the expression of $V{1}_{i,j}^{k+1}$ in equation (8) does ensure that the particle points to $p_i^k$ as an overall result. Then, the value of ${\sigma }_{i,j}^k$ is determined by an adaptive strategy, i.e. ${\sigma }_{i,j}^k=2\alpha |C_j^k-X_{i,j}^k|$ , where $C^k=(C_1^k,C_2^k,\dots ,C_D^k)$ is defined as the mean of the pbest positions of all the particles named mean best (mbest) position as followed

$$\begin{array}{ll}{C}_j^k=(1/m)\sum _{i=1}^m{P}_{i,j}^k& (1\le j\le N)\end{array}$$

(14)

Finally, the particles in RDPSO are updated equations as follows:

$$V_{i,j}^{k+1}=\alpha |C_j^k-X_{i,j}^k|{\varphi }_{i,j}^k+\beta (p_{i,j}^k-X_{i,j}^k)$$

(15)

$$X_{i,j}^{k+1}=X_{i,j}^k+V_{i,j}^{k+1}$$

(16)

where $\alpha$ means compression-expansion coefficient and $\beta$ is named acceleration coefficient. The pseudo-code of the standard RDPSO algorithm is showed as A1.

Algorithm A1: Random Drift Particle Swarm
Input: 1) a D-dimension minimization problem F(.) 2) search space S⊂RD 3) a particle swarm with population size m
Output：the optimal solution gbest
1. Begin
2. Initialize the current positions and velocities of a swarm with m particles randomly, evaluate the fitness of each particle. Set the pbest position of each particle to be its current position and generate the gbest of the swarm according to formula (4);
3. for t = 1 to maximum number of iterations
4.  Compute the mbest position Ck according to formula (14)
5.  for each particle i in the swarm
6.   Update the velocity and current position of particle i according to equations (13), (15) and (16);
7.   Evaluate the objective function value $f\left(X_i^k\right)$ and update pbest $P_i^k$ ;
8.  end for
9.  Update the gbest $G^k$ of the swarm;
10. end for
11. end

The RDPSO algorithm is quite different from the drift particle swarm optimization algorithm (DPSO)21 because it can make a better balance between the local search and the global search of the particle population leading to better performance than DPSO. In the RDPSO algorithm, ɑ and ß are two very important user specified parameters that can be adjusted to balance the local and global search of the particles. How to choose these parameter values to prevent particle explosion is an unresolved issue. Setting a larger value for ɑ and ß means that the particle has better global search ability, and vice versa means that it has better local search ability. Thus, a good balance between the local search and the global search of the RDPSO algorithm is vital to the algorithmic performance when it is applied to solving a problem. It is recommended in literature3 that ɑ should be set to be no greater than 1.0 and ß to be no greater than 1.5. Furthermore, the performance of the population-based technique can also be affected by the maximum number of iterations (MaxIteration) and the population size. As well as other PSO variants, the swarm size of RDPSO is suggested that should be greater than 20. Configuration of the MaxIteration depends on the complexity of the problem. In general, the more complex the problem, the larger MaxIteration is required. And the velocity of a particle in RDPSO is still restricted in the interval [-V _max , V _max ].

Dyanmic cooperative RDPSO based on search history

Search history scheme

Definition 1: The Sub-region of a solution Suppose a solution x is in the search space D, i.e. $x\in D$ . D is partitioned into partition set $\mathrm{ }D={\cup }_i{d}_i$ by a binary space portioning (BSP) tree T and each partition d_i is recorded as a leaf node of T. The partition $d_i\subseteq D$ is defined as the ‘the sub-region of x’ if $x\in d_i$ .

In the entire search history recorder scheme, the binary space (BSP) tree is used as a tree-structure archive to memorize all visited solutions and its fitness {[x_i, f(x_i)]}. During the iterative process, the search space is partitioned into a set of regions $D={\cup }_i{d}_i$ according to the distribution of the evaluated solutions {x_i}. A region $d_i$ of the search space is regarded as a leaf node of the BSP tree. Supposing a parent node N has two child nodes n₁ and n₂. The sub-region of N is linearly partitioned into two overlapped sub-regions by the child nodes. The corresponding partitioning of the search space is cut along the j^th dimension where $\mathrm{ }j=\arg \mathit{max}\left|n_1\left(j\right)-n_2\left(j\right)\right|$ . With this approach, each previous solution generated by the RDPSO with its fitness value can be stored in a node of the tree. Therefore, throughout the solution process, the BSP tree records all the solutions and analyzes the current solution by using the evolutionary search history of previous solutions. Since the construction of the tree depends on the order of solution set generated and evaluated by RDPSO, the BSP tree is a random tree whose topology varies from iteration to iteration.

The nodes in the BSP tree linearly divide adjacent sub-regions into a certain overlap with each other, that is, overlapped sub-regions. The method creates a channel in which particles could migrate from one sub-region to another with better fitness, as shown in Figure 1. The search for the overlapped sub-region d of a solution x is executed as tree nodes search. The BSP tree T is initialized to consist of the root node D only representing the entire search space. The search starts with checking the root node. Every iteration the search moves downward from the root node and d is shrunk along a specified direction until the search reaches a leaf node. The procedure to find the overlapped sub-region of a solution x is summarized in algorithm A2. Assume that there are two evaluated solutions in the two-dimensional search space indicated as m = [0.25, 0.5] and n = [0.75, 0.5]. The squares shown in Figure 1(a) and (b) represent the same sub-region of the parent node P. Since m and n have maximum distance along x axis, the partitioning cuts along x axis. The overlapped sub-regions of m and n, namely p_a and p_b, are defined as the gray-filled regions shown in literature 1(a) and (b) respectively. It is worth noting that the union of r_a and r_b creates a path where particles in r_a can move to r_b by mutation, and vice versa. In addition, the sub-region contraction scheme of algorithm A2 ensures that the resultant sub-region overlaps with all its adjacent sub-regions. Consequently, the opinion of overlapping sub-regions integrated with mutation allows the particles to gradually approach their optimal solution nearby. Figure 2 gives an example of the BSP tree memorized four evaluated solutions: x₁, x₂, x₃ and x₄ in two dimensional. The distribution of the evaluated solutions in search space S = [0,1]² is shown in Figure 2(a). Figure 2(b) shows the corresponding BSP tree. The archive tree consists of seven nodes. The solutions x₁ and x₂ are stored in node C and D respectively whilst the solutions x₃ and x₄ are stored in node E and F respectively.

Figure 1.

The example of overlapping sub-regions.

Figure 2.

The BSP Tree of solutions.

Algorithm A2: Overlapped Sub-region Search
Input: 1) BSP Tree T, 2) solution x $\in$ RD (D is dimension) and 3) search space S is partitioned into a set of regions $h={\prod }_i[l_i,u_i]=S$
Output: the sub-region r of solution x
1. Current_node = root node of T
2. While (Current _node has two child nodes n1 and n2 represented left and right respectively)
3. $i=\arg \mathit{max}\left|n_1\left(i\right)-n_2\left(i\right)\right|$ , i = 1, 2, …, D
4. if ( $\left|n_1\left(i\right)-x\left(i\right)\right|\le \left|n_2\left(i\right)-x\left(i\right)\right|$ )
5. Curren_tnode =  $n_1$
6. ${u_i=n}_2\left(i\right)$
7. else
8. Current_node =  $n_2$
9. ${l_i=n}_1\left(i\right)$
10. End if
11. Loop

Considering that the fitness values of the evaluated solutions are also stored in the memory archive by the BSP tree; this archive can also be regarded as the approximation $\tilde{f}\left(x\right)$ of a fitness function f(x). If the sub-region of an evaluated solution x_i is d_i, the fitness value of an unseen solution y inside the sub-region d_i can be approximated as y_i i.e. $f\left(y\right)\approx \tilde{f}\left(y\right)=y_i=f\left(x_i\right)$ . Because the fitness of all solutions in the sub-region of a BSP tree node is approximated to the same value, i.e., $\tilde{f}\left(a\right)=\tilde{f}\left(b\right)=y_i=f\left({\mathrm{x}}_i\right)$ for all $a,b,x_i\in d_i$ , $\tilde{f}\left(x\right)$ is a step-wise function. In the course of evolution, the approximation error of $\tilde{f}\left(x\right)$ decreases monotonously because estimated solutions gradually increase and recorded one by one. Hence the approximation process can be thought of as a simple incremental learning method of machine learning.

An adaptive parameter-less mutation mechanism

Definition 2: The neighborhood of a sub-region The sub-region N of node n is a neighborhood of the sub-region M of the leaf node m if M ⊆ N.

Each sub-region can be viewed as optimal in a certain amount neighborhood named neighborhood size that refers to the number of neighbor points (evaluated solutions) focused. The neighborhood size of a node decreases with its node depth due to the topology of the BSP tree. Thus, the neighborhood size is constituted by the tree node depth difference. If the sub-region M of node m is the neighborhood of the sub-region N of node n, the neighborhood size(denoted by L) of M related to N is defined as the depth difference of m and n. M is assigned as an optimal sub-region in case the fitness value of node n is the minimum within all evaluated solutions in N. Since L is selected to be the depth of the BSP tree, it is presumed that the approximated fitness functions consist of only one optimum. All particles are imposed to move toward the best-found optimal value. If L equals to zero, each sub-region is considered as a local optimum. Then the particles are randomly mutated within its own sub-region. To make the approximated fitness function locally guide the search with adaptive mutation, L is set to be 2 in this study because L should be slightly larger than zero but much smaller than the depth of the BSP tree.

Definition 3: Sub-region Distance Suppose the sub-regions of nodes m and n are M and N respectively, let the sub-region X of node x is the smallest sub-region which contains both M and N, i.e. M, N ⊆ X, the definition of the sub-region distances from M to N is the depth difference between m and x.

Once the optimal sub-region of the approximate fitness function is found, the direction of the sub-region M pointing to its nearest optimal sub-region is taken as its mutation direction. Suppose H is the optimal sub-region set in search space D and the sub-regions of nodes m and n are M and N respectively. If the distance from M to N is the minimum amongst the distance from M to all optimal sub-regions, i.e. $N=\underset{A\subseteq H}{\mathrm{argmin}}\mathit{Distance}(M,A)$ , the sub-region N ⊆ H is the nearest optimal sub-region of the sub-region M. The more details about Definition 2 and 3 can refer to literature.15

The mutation direction of the nodes in sub-region point to its nearest optimal neighborhood sub-region. In the adaptive parameter-less mutation mechanism, a particle m moves along the direction of the local maximum increasement of approximation fitness $\tilde{f}\left(x\right)$ at m. Suppose v as the mutation direction of node n pointing to the nearest historical optimum m i.e., v = m – n. The best form of approximation fitness $\tilde{f}\left(x\right)$ is a D-dimensional sub-region rather than a D-dimensional point because it can be regarded as a step function. Furthermore, the topological structure of the BSP tree describes the adjacencies of sub-regions amongst the steps. Therefore, the process to find the nearest optimum of node m equates to find the nearest optimal sub-region of m in $\tilde{f}\left(x\right)$ . A parameter α selected in the interval (0, 1) randomly is set as the mutation step size to balance the exploitative effect of this decreasing gradient like direction allocation, hence the mutant m′ of m is a linear combination of the m and n, m′ = m+αv = (1-α) m+αn.

Algorithm A3 shows the procedure of the adaptive mutation. Seen from A3, the mutation step size is randomly assigned according to the search history within an adaptive adjustment range, whilst its search direction is determined by the approximation fitness function from the BSP tree. Therefore, the mutation mechanism is an adaptive, parameter-less and guided mutation operator. Moreover, it naturally avoids the creation of any out-of-bound particles, which is superior to some other methods that might produce out-of-bound solutions and an additional repair operator is needed to change the solutions back to valid ones.

Algorithm A3: Adaptive Parameter-less Mutate in CRDPSO
Input: 1) the particle p needs to mutate, 2) BSP tree T stored the evaluated solution set $\left\{p_i\right\}$
Output: p′
1. The partitioned sub-region set $H=\cup h_i$ where $h_i$ is the sub-region of $p_i$
2. Search for the sub-region h⊆H of p
3. Find the nearest optimal sub-region of h
4. y = The evaluated solution inside h
5. α = Random(0, 1)
6. p′= αp + (1 −α)y

Dynamic cooperative evolution

Concerning high-dimensional problems, stochastic search algorithms (including random drift particle swarm optimization and genetic algorithms) suffer from the “curse of dimensionality” of the scaling problems. Since the overall fitness of a particle will be affected by every dimension, even if the dimension values of some particles are very close to the corresponding dimensions of the global optimal solution, the particles still get lower fitness. The significance of cooperation operation is that it preserves “good” dimensional information within the particle swarm and can prevent potentially useful information from being discarded unnecessarily. By performing cooperation one dimension at a time, each dimension can be evaluated separately and the most useful information saved to speed up convergence. The algorithm performance based on cooperative evolution is proved to be efficient in the high dimension problems.22,23 In the dynamic cooperative evolution, each particle p_i generates four particles {p_i1, … ,p_i4} by equation (14) to equation (16). Therefore, five individuals are obtained and then select the best individual p_{i_b} according to their fitness value among five individuals. Then p_{i_b} is assigned as the current local context vector for the particle p_i and cooperates with the other four individuals. Each dimension of the context vector is replaced with the corresponding dimension of each of the remaining four individuals successively, and the alternative dimension is tested whether it improves the particle fitness value. The replacement value will be adopted if it does, and the resulting vector is used to substitute for the context vector. The process is repeated until all dimensions of all five measurement vectors have been evaluated, and the best arrangement is selected to form a new particle for the new generation. In order to reduce the fitness function calculations for each dimension replacement, the dynamic cooperative evolution is improved by selecting the optimal vector and collaborating with other particles. The procedure of the dynamic cooperative evolution is shown in algorithm A4.

Algorithm A4: Dynamic Cooperative Evolution in CRDPSO
Input: Particle swarm P = {p1,p2, … ,pn}
Output: Particle swarm
1. for each particle
2. generate four particles using equations (14) to (16)
3. pi = pi_b = $\arg \min \{f\left(p_i\right),f\left(p_{i1}\right),\hspace{0.17em}f\left(p_{i2}\right),f\left(p_{i3}\right)\hspace{0.17em},f\left(p_{i4}\right)\}$
4. pc = pi;
5. f = f(pi_b);
6. for each particle pk, k = 1, … ,5, $k\ne$ i_b
7. for each dimension j
8. if f(pc(j,pkj))< f
9. pij = pkj
10. End if
11. End for
12. End for
13. End for

Dynamic cooperative RDPSO algorithm based on search history decision

Dynamic cooperative RDPSO algorithm based on entire search history is a real coded stochastic swarm intelligence algorithm. The whole algorithm searches for the optimum by main four parts that are swarm initialization, evolution, mutation and dynamic cooperative operation. CRDPSO algorithm focuses on selecting the optimal particle through the dynamic coordination mechanism and enhancing adaptive mutation based on the evolutionary search history. The procedure of CRDPSO is outlined in algorithm A5. For a D-dimensional minimization problem F(.) with the search space S⊂R^D, a particle p in the swarm is a 1 by D real-valued vector, i.e., $p\in R^D$ . The algorithm starts from initializing the current population of n particles P = {p₁,p₂, … ,p_n} whilst the BSP tree T is initialized to be composed of a root node only. Afterward the population is evaluated and its fitness is both stored by T. Each particle in P performs mutation adaptively, i.e. $p_i\stackrel{\mathrm{mutate}}{\to }\mathrm{\hspace{0.5em}}{p}_i^{\mathrm{\prime }}$ . Afterwards each particle p_i generates four particles {p_i1,…,p_i4} by RDPSO and dynamically cooperate with the corresponding mutant $p_i^{\mathrm{\prime }}$ on each dimension, that is, the particle group of the ith particle for dynamic cooperative evolution is $\left\{{p_i^{\mathrm{\prime }},p}_{i1},\dots ,p_{i4}\right\}$ . The best permutation is selected to form a new particle for the new offspring. The fitness value of the new offspring particles needs to be recalculated and inserted into the BSP tree T. The processes are repeated until the termination criterion is satisfied.

Algorithm A5: CRDPSO
Input: 1) a D-dimension minimization problem F(.) 2) search space S⊂RD 3) a particle swarm P with population size n
Output：the optimal solution gbest
1. Initialize the particle swarm P = {p1,p2,…,pn} and evaluate the fitness of each particle $F_i=F(p_i$ )
2. Set the pbest position of each particle to be its current position and generate the gbest
3. Initialize the BSP tree T which consists of root node only
4. for i = 1 to n
5. $[p_i,F(p_i$ )] insert to T
6. end for
7. While stop criteria is not met
8. for i = 1 to n
9. ${ p}_i$ adaptively mutate to be $p_i^{\mathrm{\text{'}}}$
10. $\left\{p_{i1},\dots ,p_{i4}\right\}$ updated from equations (14) to (16)
11. $\left\{{pi\prime \mathrm{`},p}_{i1},\dots ,p_{i5}\right\}$ create new next generation according to the dynamic coordination mechanism
12. end for
13. for i = 1 to n
14. $[p_i,F(p_i)$ ] insert to T
15. end for
16. for i = 1 to n
17. if $F\left(p_i\right)<F\left({\mathit{pbest}}_i\right)$
18. ${\mathit{pbest}}_i=p_i$
19. end if
20. end for
21. Update global best fitness gbest
22. Loop

Simulation setup

Test function

A set of multimodal or unimodal benchmark functions F = {f₁(x), f₂(x),…, f₁₀(x)}2425 is adopted to test for the proposed algorithm and verify its performance. The ten test functions are shown in Table 1. All the empirical studies in this paper were carried out in the MATLAB.

Table 1.

Test functions.

		X	Global fmin
f1 Spherical model	$\sum _{i=1}^D{x}_i^2$	[−100, 100] D	0
f2 Schwefel’s Problem 2.22	$\sum _{i=1}^D\left|x_i\right|+\prod _{i>1}^D\left|x_i\right|$	[−10, 10] D	0
f3 Schwefel’s Problem 1.2	$\sum _{i=1}^D{\left\{\sum _{j=1}^i{x}_j\right\}}^2$	[−100, 100] D	0
f4 Schwefel’s Problem 2.21	$\underset{i\in [1,D]}{\max }\left|x_i\right|$	[−10, 10] D	0
f5 Restrigin’s function	$\sum _{i>1}^D[x_i^2-10\cos \left(2\pi x_i\right)+10]$	[−5.12,5.12] D	0
f6 Ackley function	$-20exp\left(-0.2\sqrt{\frac{1}{D}\sum _{i=1}^D{x}_i^2}\right)-\mathit{exp}\left(\frac{1}{D}\sum _{i=1}^D\cos 2\pi x_i\right)$ +20+e	[−32, 32] D	0
f7 Rosenbrock’sfuncton	$\sum _{i=1}^{D-1}[100{\left(x_{i+1}-x_i^2\right)}^2+{(x_i-1)}^2$	[−29, 31] D	0
f8 Griewank function	$\frac{1}{4000}\sum _{i=1}^D{x}_i^2-{\prod }_{i=1}^D\cos \frac{x_i}{\sqrt{i}}+1$	[−600, 600] D	0
f9 Quaric function	$\sum _{i=1}^{D}i{x}_1^4+\mathit{random}\left(0\text{,}1\right)$	[−1.28,1.28] D	0
f10 Pathological	$\sum _{i=1}^D(0.5+\frac{{\sin \sqrt{100x_i^2+x_{i+1}^2}}^2-0.5}{1+{0.001\left(x_i^2-2x_i{x}_{i+1}+x_{i+1}^2\right)}^2})$	[−100, 100] D	0

Algorithms setting

Besides RDPSO and CRDPSO, Differential Evolution (DE),26 Covariance Matrix Adaptation Evolutionary Strategies(CMA-ES),27 cNrGA and two versions of QPSO, i.e. CCQPSO23 and CLQPSO,28 were also tested on the test functions for performance comparison. The parameters of the CRDPSO and RDPSO in our experiments are set the same as those recommended in literature.3 To CLQPSO the β decreases linearly from 1 down to 0.5 and other factors are set as literature.28 To cNrGA, the cross rate is uniform crossover and sets as r_x  = 0.5. For all compared algorithms the population size is set to 100 and the maximum number of iterations is set to 2000. The dimension of all the test functions is set as D = 30. Each algorithm independently runs 30 times in experiments, and the mean value of optimal solutions and standard deviations are achieved after 30 independent trials.

Simulation results

For all the test functions, the statistical values from 30 search runs of 2000 iterations by each algorithm are listed in Table 2 where “Best”, “Worst” and “Ave.” stand for the best, the worst and the mean values of optimal solutions obtain in 30 independent runs, respectively. “Std.” represents the standard deviations of the achieved optimal solutions. The value of the boldface means that the corresponding algorithm gets the optimal value amongst the test algorithm for a particular benchmark function. It can clearly be seen from Table 2, the proposed CRDPSO is the most superior to other six algorithms for all the test functions. Although RDPSO and CCQPSO have both shown outstanding performance on f₁, f₂ and f₆, CRDPSO has outperformed them on all the rest functions. It shows the standard deviation of CRDPSO is the best for almost the test functions except f_9. This demonstrates that CRDPSO has the best stability for the most test functions. Moreover, the detailed simulation results (Best, Worst, Ave. and Std.) show that the best performance of CRDPSO ranks first in total of ten cases. The analysis of RDPSO29 shows that the RDPSO algorithm is more stable and better for most test environments. The algorithm can maintain the diversity of the swarm and improve the global search ability of particles by using the average position waiting mechanism. Therefore, the diversity in the algorithm declines in a smooth way in the later stage of the search process so that there is good balance between the global search and the local search of the algorithm, which results in an overall better performance than the other variants of PSO. Moreover, the simulation results of the empirical studies listed in Table 2 show that the performance of CRDPSO is better than the other compared algorithms. Using the search history mechanism, CRDPSO can ensure non-repeating access during the search process and rapid convergence to the local optimal. Meanwhile, it employs adaptive mutation to guide particles to search better neighborhoods and further improve the local search ability. In general, CRDPSO obtained the highest average accuracy and it also is the most stable.

Table 2.

The optimize results of the test function.

	RDPSO	CRDPSO	CMA-ES	cNrGA	CCQPSO	DE	CLQPSO
f 1
Best	0	0	1.75897E-25	1.1139E-08	0	4.18904E-07	7.40207E-07
Worst	0	0	4.3217E-25	9.7340E-01	0	4.29766E-06	2.13162E-06
Ave.	0	0	2.84821E-25	3.5001E-02	0	1.34E-06	1.40809E-06
Std.	0	0	6.32276E-26	1.7732E-01	0	8.67663E-07	3.58137E-07
f 2
Best	0	0	1.3856E-13	3.1145E-06	0	0.00137362	2.9643E-05
Worst	0	0	2.8599E-13	1.0248E-01	0	0.00579298	6.8636E-05
Ave.	0	0	2.3057E-13	1.4545E-02	0	0.00273726	4.6353E-05
Std.	0	0	3.4456E-14	2.6594E-02	0	0.00114531	9.9336E-06
f 3
Best	9.3627E-04	0	8.2294E-24	4.8387E + 02	9.2864E-12	39.3725608	5608.60016
Worst	5.6610E-02	0	3.4376E-23	2.7991E + 03	0.00042061	162.940307	13287.5458
Ave.	1.1162E-02	0	1.6878E-23	1.7262E + 03	1.7261E-05	88.3684554	9.8525E + 03
Std.	1.2250E-02	0	6.1177E-24	6.1711E + 02	7.6411E-05	35.5411464	1.9354E + 03
f 4
Best	6.9065E-12	0	1.7053E-13	1.4496E-01	0	1.2243E-01	1.6449E + 01
Worst	6.5872E-10	0	4.1211E-13	6.9400E-01	5.8236E-11	1.7228E + 00	2.2702E + 01
Ave.	1.5596E-10	0	2.4727E-13	3.5697E-01	1.9668E-12	4.3307E-01	1.9765E + 01
Std.	1.6442E-10	0	5.2988E-14	1.4813E-01	1.0628E-11	3.4863E-01	1.4157E + 00
f 5
Best	9.0126E + 00	0	2.6864E + 01	3.1406E-08	2.9849E + 00	1.8476E + 02	1.3495E-06
Worst	6.5233E + 01	0	8.4571E + 01	2.0674E + 00	1.4924E + 01	2.2501E + 02	8.2318E-06
Ave.	1.7756E + 01	0	4.6498E + 01	6.7054E-01	9.3899E + 00	2.0197E + 02	3.1481E-06
Std.	9.8049E + 00	0	1.2827E + 01	8.4448E-01	2.7248E + 00	9.7150E + 00	1.5211E-06
f 6
Best	8.8818E-16	8.8818E-16	2.0000E + 01	5.6498E-05	8.8818E-16	1.7317E-04	3.5003E-04
Worst	8.8818E-16	8.8818E-16	2.1593E + 01	2.3660E-01	8.8818E-16	5.6041E-04	9.2455E-04
Ave.	8.8818E-16	8.8818E-16	2.1436E + 01	1.8326E-02	8.8818E-16	3.5209E-04	6.1720E-04
Std.	1.0029E-31	1.0029E-31	2.8602E-01	4.4425E-02	1.0029E-31	1.0143E-04	1.2824E-04
f 7
Best	2.0862E-03	0	4.8427E + 00	1.5554E + 00	1.3355E + 01	2.2035E + 01	2.6395E + 01
Worst	9.2749E + 01	1.4777E-08	1.4112E + 02	4.3660E + 02	9.4946E + 01	2.6257E + 01	6.7362E + 01
Ave.	3.7207E + 01	5.1687E-10	1.9915E + 01	8.6977E + 01	3.3428E + 01	2.3721E + 01	4.1579E + 01
Std.	3.0608E + 01	2.6966E-09	2.9182E + 01	8.4565E + 01	2.8903E + 01	9.1871E-01	1.0091E + 01
f 8
Best	0	0	0	7.4061E-04	0	9.4913E-07	5.8648E-06
Worst	2.9459E-02	0	1.2321E-02	3.0395E-01	3.6770E-02	2.6043E-05	3.4202E-05
Ave.	6.4840E-03	0	1.9714E-03	9.6984E-02	1.0330E-02	3.9483E-06	2.0352E-05
Std.	7.7636E-03	0	4.1634E-03	8.5924E-02	1.1312E-02	4.5272E-06	7.6931E-06
f 9
Best	1.2928E-03	6.7247E-07	5.0404E-02	6.6570E-03	8.4983E-06	1.1918E-01	4.7632E-03
Worst	5.5234E-03	3.4890E-05	5.7714E-01	3.1475E-02	3.3232E-05	3.6087E-01	1.3281E-02
Ave.	2.9361E-03	9.5058E-06	2.5975E-01	1.8056E-02	1.8129E-05	2.1710E-01	2.1328E-02
Std.	1.0048E-03	9.0303E-06	1.2435E-01	7.4958E-03	5.4642E-06	5.9533E-02	3.1031E-02
f 10
Best	9.2795E + 00	0	1.3231E + 01	2.5737E + 00	5.0896E + 00	1.1207E + 01	4.7307E + 00
Worst	1.1470E + 01	1.0629E-12	1.4500E + 01	4.5666E + 00	7.3723E + 00	1.2348E + 01	6.0035E + 00
Ave.	1.0537E + 01	5.3144E-14	1.3970E + 01	3.5592E + 00	6.1297E + 00	1.1974E + 01	5.5151E + 00
Std.	5.8348E-01	2.3767E-13	2.5855E-01	6.3615E-01	5.2576E-01	2.2182E-01	3.0625E-01

Convergence curves versus iteration for the different algorithms are shown in Figure 3. The results in Table 2 show that for f₁, f₂ and f₆, all of RDPSO, CRDPSO, and CCQPSO finally achieve the optimal values. However, as can be seen from the convergence curves of Figure 3 (1)(2)(6) respectively, the convergence speed of CRDPSO is much better than the other two algorithms. It has converged to the global optimal value within less than 200 iterations of the search. For f₃- f₅ and f₇- f₁₀, the convergence speed and convergence accuracy of CRDPSO are far beyond the other six algorithms. The global optimal values of each benchmark were found within a short number of iterations. CRDPSO algorithm can search the optimal value much faster because the particles can adaptive mutate once their own nearest optimal neighborhood sub-region is found by using the entire evolutionary history information. Moreover, the coordination mechanisms can be used to keep the optimal information, so the algorithm can quickly converge to the optimal value. In Figure 3, it is shown that for f₁ to f₅ and f₈ the CRDPSO was able to reached 0 that the Matalab can recognize and the logarithmic scale of 0 is minus infinity, so that the fitness value cannot be shown anymore in the figure. This also indicates that the simulation results for the objective function values obtained by CRDPSO for these test functions are infinitely close to the theoretical values, i.e. 0. For f₇ in Figure 3(7), even though the algorithm converged, but it was able to further find a new better solution in the later stage (i.e. after 1500 iterations), implying that CRDPSO algorithm still has a better search ability in the final stage of the search process and has a large chance to escape local optimum to get a better solution. In conclusion, CRDPSO algorithm compared with other algorithms has the best convergence performance and optimization accuracy with a strong search ability.

Figure 3.

Convergence curve of the algorithms.

In order to determine whether there was a significant difference of the performance between CRDPSO algorithm with other six algorithms in addition to the precision and convergence. Table 3 lists Wilcoxon rank sum test results for a significance level of 0.05 between different algorithms. Seen from Table 3 for f₁, f₂ and f₆, CRDPSO, RDPSO and CCQPSO have no significant difference. Nevertheless, for the vast majority of the test functions, CRDPSO algorithm can get h = 1, shows that compared with other algorithms, the optimized performance of CRDPSO has improved significantly. It is observed that the optimal performance of CRDPSO algorithm improved prominently with better convergence and higher precision whether in the unimodal functions and multi-peak functions.

Table 3.

The Wilcoxon rank sum test.

	CRDPSO vs CCQPSO		CRDPSO vs RDPSO		CRDPSO vs CLQPSO		CRDPSO vs cNrGA		CRDPSO vs DE		CRDPSO vsCMA_ES
	p	h	p	h	p	h	p	h	p	h	p	h
f 1	NaN	0	NaN	0	1.2118E-12	1	1.2118E-12	1	1.2108E-12	1	1.2108E-12	1
f 2	NaN	0	NaN	0	1.2118E-12	1	1.2118E-12	1	1.2118E-12	1	1.2059E-12	1
f 3	1.2118E-12	1	1.2118E-12	1	1.2118E-12	1	1.2118E-12	1	1.2118E-12	1	1.2118E-12	1
f 4	3.6097E-13	1	1.2118E-12	1	1.2118E-12	1	1.2118E-12	1	1.1001E-12	1	1.1001E-12	1
f 5	1.1613E-12	1	1.2108E-12	1	1.2118E-12	1	1.2118E-12	1	1.1001E-12	1	1.2088E-12	1
f 6	NaN	0	NaN	0	1.2118E-12	1	1.2118E-12	1	1.1001E-12	1	1.2118E-12	1
f7	9.9025E-09	1	4.0279E-09	1	1.5647E-11	1	1.6027E-10	1	5.6013E-08	1	2.8306E-08	1
f 8	2.2016E-06	1	5.2014E-06	1	1.2118E-12	1	1.2118E-12	1	1.2118E-12	1	2.8000E-03	1
f 9	7.8252E-05	1	4.4618E-11	1	4.4618E-11	1	6.5095E-11	1	4.4618E-11	1	6.5095E-11	1
f 10	4.6908E-07	1	3.1575E-12	1	3.1575E-12	1	3.1575E-12	1	3.1549E-12	1	3.1549E-12	1

Conclusion

Premature convergence and diversity are the most need to solve two problems with the RDPSO algorithm. Therefore, an improved dynamic cooperative RDPSO algorithm based on search history assisted (CRDPSO) is proposed that integrates the search history scheme and a standard random drift particle swarm optimization (RDPSO). The proposed algorithm utilizes a binary space partitioning tree (BSP) structure to store the positions and the fitness values of the evaluated solutions. The dynamic cooperation mechanism between the solutions in the swarm is introduced for better use of context information of each dimension to take advantage of any new information, which improves the capability of the variation of the particles and enhances search ability helping to prevent premature convergence. Benefiting from the spatial partitioning scheme, a fast fitness function for improving the mutation strategy in CRDPSO is obtained using the archive structure. The resultant mutation is a kind of adaptive and parameter-less. The overlap between adjacent sub-regions also provides a channel for the movement of particles in adjacent areas, making it easier for particles to move from the area to a better area. In comparison to other well-known algorithms, the experimental results on ten standard testing functions indicate that the proposed algorithm has the best global optimization ability, with enhancement in both convergence speed and precision that demonstrate the effectiveness and efficiency of the CRDPSO. The improved algorithm can be proved in theory to be superiority; however, it also needs to be further tested in practical problems. Hence the CRDPSO algorithm applied to practical engineering problems becoming the focus research in the future. The next research job is to verify the CRDPSO performance of the engineering applications such as image segmentation, data clustering and parameter identification of biological pathways.