# Multispecies Coevolution Particle Swarm Optimization Based on Previous Search History

Authors: Danping Wang; Kunyuan Hu; Lianbo Ma; Maowei He; Hanning Chen
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

A hybrid coevolution particle swarm optimization algorithm with dynamic multispecies strategy based onK-means clustering and nonrevisit strategy based on Binary Space Partitioning fitness tree (called MCPSO-PSH) is proposed. Previous search history memorized into the Binary Space Partitioning fitness tree can effectively restrain the individuals’ revisit phenomenon. The whole population is partitioned into several subspecies and cooperative coevolution is realized by an information communication mechanism between subspecies, which can enhance the global search ability of particles and avoid premature convergence to local optimum. To demonstrate the power of the method, comparisons between the proposed algorithm and state-of-the-art algorithms are grouped into two categories: 10 basic benchmark functions (10-dimensional and 30-dimensional), 10 CEC2005 benchmark functions (30-dimensional), and a real-world problem (multilevel image segmentation problems). Experimental results show that MCPSO-PSH displays a competitive performance compared to the other swarm-based or evolutionary algorithms in terms of solution accuracy and statistical tests.


---

## Body

## 1. Introduction

Solving numerical optimization problems is a challenging research endeavor in many scientific areas. Many optimization techniques and search algorithms have drawn their motivation from evolution and social behavior. These include ant colony optimization [1], genetic algorithms [2], differential evolution [3], particle swarm optimization [4], and artificial immune systems [5]. In swarm intelligence (SI) systems, there are many simple individuals who can interact locally with one another and with their environments. Although such systems are decentralized, local interactions between individuals lead to the emergence of global behavior and properties. All of them have many variants, which have excellent performance. These variants are based on various improvement strategies.In the last few years, several heuristics have been developed to improve the performance and set up suitable parameters for the PSO algorithm. van den Bergh [6] analyzed the trajectory of particles under different inertia weights and acceleration coefficients. The original structure of PSO model reflects an intuitive idea that a particle takes any position on which fitness value is better than where it currently is as the reference input. However, many experiments have shown that the basic PSO algorithm easily falls into local optima when solving complex multimodal problems with a huge number of local optima.To overcome this problem, researchers conducted lots of studies. Inspired by the coevolution phenomenon in nature, in Cooperative Particle Swarm Optimization (CPSO) [7], multiple swarms are partitioned uniformly to optimize different components of the solution vector cooperatively. Inspired by this work, a CCPSO integrating the random grouping and adaptive weighting schemes was developed and demonstrated great promise in scaling up PSO on high-dimensional nonseparable problems [8]. In [9], a competitive and cooperative coevolutionary PSO (CCPSO) has considerable potential for solving complex optimization problems by explicitly modeling the coevolution of competing and cooperating species. In Multiswarm Self-Adaptive and Cooperative Particle Swarm Optimization (MSCPSO) [10], particles in each subswarms share the only global historical best optimum to enhance the cooperative capability. In [11], Adaptive Cooperative Particle Swarm Optimizer (ACPSO) facilitates cooperation technique through the usage of the Learning Automata (LA) algorithm. Coevolutionary particle swarm optimization (PSO) algorithm associated with the artificial immune principle (ICPSO) [12] adopts an improved immune-clonal-selection operator for optimizing the elite subpopulation and a migration scheme for the information exchange between elite subpopulation and normal subpopulations.Different from the above works with static-swarms strategy, some cooperative PSO’ variants introduce dynamic multiswarms strategy to enhance the ability of exploring local optima. With the aim of maintaining multiple swarms on different peaks, a clustering PSO (CPSO) proposed in [13] employs a nearest neighbor learning strategy to train particles and a hierarchical clustering method to locate and track multiple optima. In hierarchical cluster-based multispecies particle swarm optimization (HCMSPSO) algorithm [14], a swarm is clustered into multiple species at an upper hierarchical level, and each species is further clustered into multiple subspecies at a lower hierarchical level. And in the lower layer, subspecies within the same species are formed adaptively in each iteration during the particle update. In [15], a dynamic multiswarm particle swarm optimizer (DMS-PSO) merges the harmony search (HS) algorithm into each subswarm. These subswarms are regrouped frequently and information is exchanged among the particles in the whole swarm. In a novel optimizer using Adaptive Heterogeneous Particle Swarms (AHPS2) [16], to best take advantage of the heterogeneous learning strategies, an adaptive competition strategy is designed for dynamically adjusting the population size of an independent swarm with comprehensive learning and another one with subgradient learning based on their group performance. Roles of learning swarm and learnt swarm in particle swarm optimization with interswarm interactive learning strategy (IILPSO) [17] swap to maintain the diversity when interswarm interactive learning (IIL) behavior is triggered to adjust sizes of these two swarms.In order to avoid the different individuals to repetitively exploit the same researched regions or excessive individuals to prematurely search a narrow region, researchers conducted a series of studies. Inspired by biology, niche [18] and speciation [19] techniques are introduced into PSO to prevent the swarm from crowding too closely and to locate as many optimal solutions as possible. Additionally, PSO topological structures are also widely adopted. The ring topology employed in [20] operating as a niching algorithm can drive the particles exploring the search space more broadly. Further, in [21] the multilayer strategy with multiple topologies is adopted to decrease the amount of noneffective exploiting tries. Inspired by ecological cohabitation, chaotic multiswarm particle swarm optimization (CMS-PSO) modifies the generic PSO with the help of the chaotic sequence for multidimension unknown parameter estimation and optimization by forming multiple cooperating swarms [22]. Historical memory strategy in HMPSO [23], which estimates and preserves distribution information of particles’ historical promising, is helpful in preserving the information of optimum solution space and making a comprehensive learning.Although all the above works improve the performance of particle swarm optimization, there is shortage of research into the performance of the hybrid algorithm with the dynamic multiple swarms strategy and nonrepetition search approach.In this paper, we propose a new coevolution particle swarm optimization algorithm with dynamic multispecies strategy based onK-means clustering and nonrevisit strategy based on Binary Space Partitioning fitness tree (called MCPSO-PSH). It is shown that when a Binary Space Partitioning tree archive stores the positions and the fitness values of all evaluated solutions, this archive can be treated as an approximation of the fitness function, which can avoid the individuals’ revisit phenomenon and improve the search efficiency. Moreover, K-means clustering method is introduced to partition whole population into subpopulations frequently. Information communication among subpopulations is implemented in the dynamic repartition process. Therefore, MCPSO-PSH algorithm can accommodate a considerable potential for solving complex problems. To comprehensively evaluate the performance of the proposed algorithm, MCPSO-PSH is compared with other state-of-the-art algorithms on three categories of experiments: (1) ten 10-dimensional and 30-dimensional benchmark problems with various properties, such as unimodal and multimodal functions, are employed to test MCPSO-PSO’s performance for a diverse set of problems; (2) a set of CEC2005 30-dimensional tests including 10 benchmark functions are used to justify MCPSO-PSH’s scalability for complex problems; (3) a real-world problem, the multilevel thresholds based on Otsu method for image segmentation, is employed to benchmark MCPSO-PSH’s applicability for practical problems. The numerical results demonstrate that MCPSO-PSH significantly improves the performance of PSO and outperforms most of the comparison algorithms during the experiments. The main contributions of the proposed methods lie in the following aspects.(1)
A new way of information communication among multiple swarms is designed. Each particle’s position is updated according to not only its personal history information and global information of its locating species but also other global information of other species.(2)
A dynamicK-means cluster is designed. In the process of decreasing the number of clusters and repartitioning the population into species, the search feature general achieves the transition from global search in the earlier evolution period to local search in the later evolution period.(3)
Binary Space Partitioning fitness tree memorizing the previous search history can prevent the individuals from revisiting the unpromising landscapes, which can improve the search efficiency and avoid trapping in local optimum.The rest of this paper is organized as follows: Section2 first gives a review of the related works; Section 3 gives a brief explanation of the proposed coevolution particle swarm optimization algorithm; Section 4 presents the experimental studies of the proposed MCPSO-PSH; in Section 5, the performance of MCPSO-PSH based multilevel thresholding for image segmentation is evaluated. Finally, Section 6 concludes the paper.


## 2. Related Works

### 2.1. Particle Swarm Optimization

Particle swarm optimization (PSO) is an evolutionary computation (EC) algorithm paradigm that emulates the swarm behaviors of birds flocking [24]. It is a population-based search algorithm that exploits a population of individuals to probe promising regions of the search space. In this context, the population is called a swarm, and the individuals are called particles. Each particle moves with an adaptable velocity within the search space and retains the best position it ever encountered in its memory. The standard version of PSO is briefly described here [6].Lets be the swarm size, d be the particle dimension space, and each particle of the swarm have a current position vector Xi, a current velocity vector Vi, and an individual best position vector Pi found by the particle itself so far. The swarm also has the global best position vector Pg found by any particle during all prior iterations in the search space. Assuming that the function f(X) is to be minimized, and describing the following notations in tth iteration, the definitions are as follows:(1)Xit=xi,1t,xi,2t,…,xi,dt,xi,jt∈Rj,i=1,2,…,s,where each dimension of a particle is updated using the following equations:(2)vi,dt+1=w·vi,dt+c1·r1·Pi,dt-xi,dt+c2·r2·Pg,dt-xi,dt,(3)xi,dt+1=xi,dt+vi,dt+1.In (2), c1 and c2 denote constant coefficients and r1 and r2 are elements from random sequences in the range of (0, 1).The inertia weightw plays the role of balancing the global and local searches. It can be a positive constant (e.g., 0.9) or even a positive linear or nonlinear function of time [13, 14]. Although sometimes proper fine-tuning of the parameters c1 and c2 leads to an improved performance, an extended study of the cognitive and social parameters in [6] suggests c1=c2=2 as default values.


### 2.2.K-Means Clustering

K-means clustering is a method of vector quantization, originally from signal processing, that is popular for cluster analysis in data mining. K-means clustering aims to partition n observations into k clusters in which each observation belongs to the cluster with the nearest mean, serving as a prototype of the cluster, as shown in Figure 1. This results in a partitioning of the data space into Voronoi cells [25].Figure 1
K-means clustering result with k=3.The cluster centers are substituted for center positions of food sources and the formula of computing the centers is shown in (4). If the ith cluster contains ni members and the members are denoted as x1i,x2i,…,xnii, then the center (clustericenter) is determined as(4)clustericenter=∑i=1nixini.The radiusI of a cluster is defined as the mean distance (Euclidean) of the members from the center of the cluster. Thus, R can be written as follows:(5)Ri=∑p=1nixp-clustericenterni,where xp is position of the pth member and ni is the number of members in its corresponding cluster.


### 2.3. Fitness Tree

Inspired by Binary Space Partitioning (BSP) tree, fitness tree is proposed by Shiu and Chi [26]. The establishment of a fitness tree is aimed at partitioning the whole search space into several subspaces (subregions) according to the history positions of members.Each node of fitness tree represents a partitioned subspace of the whole search space. And each node has at most two child nodes, left node and right node. The sum of the search spaces of left node and right node is the search space of the parent node. Such fitness tree is constructed by a series of solutions found by the individuals. The process of building fitness tree is random, which means its topology for information communication is constructed stochastically.The individuals of an evolution algorithm can be regarded as a sequence of solutionssq=s1,s2,…,sn, where n is the population size. In each cycle, n parent individuals generate n child individuals. The child individual is inserted into the fitness tree as a leaf node. So fitness tree stores all the search positions SQ={sq1,sq2,…}. This means that there is no revisit for the moment. For normal evolution algorithms, parent individuals with better solutions withstand higher selection pressure, which results in higher probability to revisit the search space. In the extreme case, an individual trapping into a local optimum may be selected over and over, which causes that after several iterations all the individuals are generated from such individual. The phenomenon of revisit leads to a precipitous loss of population diversity.The pseudocode of fitness tree construction is given in Algorithm1. Here, a detailed explanation is given as follows. In the process of inserting a solution into the fitness tree, the current nodeCurr_node represents the survey node. When the survey dimension of current node is equal to the dimension of inserted node Curr_node=s(i), the ith solution s(i) in the sequence sq=s1,s2,… is a revisit.Algorithm 1:Pseudocode of fitness tree construction.
Input: Fitness tree and A new solution
Output: Fitness tree
%% A new solutionz is inserting into the fitness tree
(1)  Initialization: Current node is the root node Curr_node = root, RF = 0
(2)  If (Curr_node has two child nodes x and y)
(3)     Define the comparing dimension j:
(4)         j=argmaxk∈1,D⁡xj-yj
(5)     If xj–zj≤xj–yj
(6)        Curr_node=childnodex;
(7)     Else
(8)        Curr_node=childnodey;
(9)     End If and Goto Step  (2)
(10)  Else
(11)     Insert a “virtual” child node to Curr_node, create a “real” child node that records z “under” the
virtual child node.
(12)  End IfThe core idea of fitness tree is that, under finding out the dimension with max difference, the new individual locates as a leaf node or partitions a parent node (subspace) into two binary nodes (sub-subspace).


## 2.1. Particle Swarm Optimization

Particle swarm optimization (PSO) is an evolutionary computation (EC) algorithm paradigm that emulates the swarm behaviors of birds flocking [24]. It is a population-based search algorithm that exploits a population of individuals to probe promising regions of the search space. In this context, the population is called a swarm, and the individuals are called particles. Each particle moves with an adaptable velocity within the search space and retains the best position it ever encountered in its memory. The standard version of PSO is briefly described here [6].Lets be the swarm size, d be the particle dimension space, and each particle of the swarm have a current position vector Xi, a current velocity vector Vi, and an individual best position vector Pi found by the particle itself so far. The swarm also has the global best position vector Pg found by any particle during all prior iterations in the search space. Assuming that the function f(X) is to be minimized, and describing the following notations in tth iteration, the definitions are as follows:(1)Xit=xi,1t,xi,2t,…,xi,dt,xi,jt∈Rj,i=1,2,…,s,where each dimension of a particle is updated using the following equations:(2)vi,dt+1=w·vi,dt+c1·r1·Pi,dt-xi,dt+c2·r2·Pg,dt-xi,dt,(3)xi,dt+1=xi,dt+vi,dt+1.In (2), c1 and c2 denote constant coefficients and r1 and r2 are elements from random sequences in the range of (0, 1).The inertia weightw plays the role of balancing the global and local searches. It can be a positive constant (e.g., 0.9) or even a positive linear or nonlinear function of time [13, 14]. Although sometimes proper fine-tuning of the parameters c1 and c2 leads to an improved performance, an extended study of the cognitive and social parameters in [6] suggests c1=c2=2 as default values.


## 2.2.K-Means Clustering

K-means clustering is a method of vector quantization, originally from signal processing, that is popular for cluster analysis in data mining. K-means clustering aims to partition n observations into k clusters in which each observation belongs to the cluster with the nearest mean, serving as a prototype of the cluster, as shown in Figure 1. This results in a partitioning of the data space into Voronoi cells [25].Figure 1
K-means clustering result with k=3.The cluster centers are substituted for center positions of food sources and the formula of computing the centers is shown in (4). If the ith cluster contains ni members and the members are denoted as x1i,x2i,…,xnii, then the center (clustericenter) is determined as(4)clustericenter=∑i=1nixini.The radiusI of a cluster is defined as the mean distance (Euclidean) of the members from the center of the cluster. Thus, R can be written as follows:(5)Ri=∑p=1nixp-clustericenterni,where xp is position of the pth member and ni is the number of members in its corresponding cluster.


## 2.3. Fitness Tree

Inspired by Binary Space Partitioning (BSP) tree, fitness tree is proposed by Shiu and Chi [26]. The establishment of a fitness tree is aimed at partitioning the whole search space into several subspaces (subregions) according to the history positions of members.Each node of fitness tree represents a partitioned subspace of the whole search space. And each node has at most two child nodes, left node and right node. The sum of the search spaces of left node and right node is the search space of the parent node. Such fitness tree is constructed by a series of solutions found by the individuals. The process of building fitness tree is random, which means its topology for information communication is constructed stochastically.The individuals of an evolution algorithm can be regarded as a sequence of solutionssq=s1,s2,…,sn, where n is the population size. In each cycle, n parent individuals generate n child individuals. The child individual is inserted into the fitness tree as a leaf node. So fitness tree stores all the search positions SQ={sq1,sq2,…}. This means that there is no revisit for the moment. For normal evolution algorithms, parent individuals with better solutions withstand higher selection pressure, which results in higher probability to revisit the search space. In the extreme case, an individual trapping into a local optimum may be selected over and over, which causes that after several iterations all the individuals are generated from such individual. The phenomenon of revisit leads to a precipitous loss of population diversity.The pseudocode of fitness tree construction is given in Algorithm1. Here, a detailed explanation is given as follows. In the process of inserting a solution into the fitness tree, the current nodeCurr_node represents the survey node. When the survey dimension of current node is equal to the dimension of inserted node Curr_node=s(i), the ith solution s(i) in the sequence sq=s1,s2,… is a revisit.Algorithm 1:Pseudocode of fitness tree construction.
Input: Fitness tree and A new solution
Output: Fitness tree
%% A new solutionz is inserting into the fitness tree
(1)  Initialization: Current node is the root node Curr_node = root, RF = 0
(2)  If (Curr_node has two child nodes x and y)
(3)     Define the comparing dimension j:
(4)         j=argmaxk∈1,D⁡xj-yj
(5)     If xj–zj≤xj–yj
(6)        Curr_node=childnodex;
(7)     Else
(8)        Curr_node=childnodey;
(9)     End If and Goto Step  (2)
(10)  Else
(11)     Insert a “virtual” child node to Curr_node, create a “real” child node that records z “under” the
virtual child node.
(12)  End IfThe core idea of fitness tree is that, under finding out the dimension with max difference, the new individual locates as a leaf node or partitions a parent node (subspace) into two binary nodes (sub-subspace).


## 3. Algorithm

### 3.1. Multispecies Coevolution Model for Optimization

Inspired by the mutualism phenomenon, the single population PSO is extended into the interacting multispecies model by enhancing particle dynamics.In such coevolution model, the information for the particle’s position updating includes not only its personal history information and global information of its locating species but also other global information of other species. Information communication simultaneously occurs in and between the local species and neighboring species.In mathematical terms, the original particle’s fly velocity equation and position update equation (please see (2) and (3)) are improved as follows:(6)vi,dkk+1=w·vi,dkt+c1·r1·Pi,dt-xi,dkt+c2·r2·Pg,dkt-xi,dkt+c3·r3·Pg,dθt-xi,dkt,(7)xi,dkt+1=xi,dkt+vi,dkt+1,where the superscript k represents the kth species, xi,dk denotes the position of the ith particle of the kth species, Pi,d is the current particle’s history best position, Pg,dk is the history best position found by its locating species, and Pg,dθ is the history best position found by the whole population. In (6), c1, c2, and c3 are the constant coefficients. And in this paper, c1, c2, and c3 are set to have the same value, c1=c2=c3=1.6667; r1, r2, and r3 are the random values in the range of (0, 1).In (6), c1·r1·Pi,dt–xi,dkt takes into account the individual’s own experiences, c2·r2·Pg,dkt–xi,dkt means the internal communication within the kth species, and c3·r3·Pg,dθt–xi,dkt indicates the symbiotic coevolution between dissimilar species.


### 3.2. DynamicallyK-Means Cluster

In order to implement the multispecies coevolution model, the stochastically generated population should be partitioned inton subspecies. The widely adopted K-means cluster method (please see Section 2.2) is employed for achieving the population division. The main drawback of K-means cluster is that the effect of population division depends on the initial distribution of stochastically generated individuals and the selection of initial division points. A relatively effective means is to set a large number of clusters.In our proposed algorithm, the number of clusters is predefined asC=c1,c2,…,cm, where c1>c2>⋯>cm. During optimization, each individual in every cluster operates independently as the multispecies coevolution model (please see Section 3.1). After several optimization iterations Ti, the whole population is repartitioned into ci+1 clusters; ci and ci+1 are orderly got from C=c1,c2,…,cm. Here, T1+T2+⋯+Tm=T, whereT is the total number of iterations. As ci+1 is smaller than ci, the individuals in small clusters may be redistributed into the relatively larger clusters when the number of the clusters changes.To balance the exploration and exploitation, the value ofTi is not constant. To enhance the ability of local search especially in early stages of searching, Ti should be larger than Ti+1. Therefore, in our proposed algorithm, Ti is calculated as follows:(8)Ti=m-i+1∑j=1mj·T.Moreover, in the process of optimization, the centers of two or more clusters move too close or overlap each other. In order to avoid this phenomenon and to control every cluster to locally search a part of landscape, the relatively small cluster in two or more adjacent clusters decomposes and the individuals in such cluster redistribute into other adjacent clusters. Evaluation criterion is listed as follows:(9)Disclusters-Rcluster1-Rcluster2<0.2·min⁡Rcluster1,Rcluster2,where Disclusters is the distance between two adjacent clusters and Rcluster1 and Rcluster2 are the radiuses of two adjacent clusters (please see (5)); the left term (Disclusters-Rcluster1-Rcluster2) represents the gap between two adjacent clusters.


### 3.3. Adaptive Normalization of Population Members

In the proposed multispecies particle swarm optimization, the idea of fitness tree is also employed to overcome the revisit phenomenon. Here, a modified adaptive normalized fitness tree is proposed.In the process of inserting a new individual into fitness tree (please see Algorithm1), there is a potential drawback: when determining the comparing dimension, it depends on the maximum difference among the dimensions of two child nodes.For overcoming this drawback, the normalization method is employed to determine which dimension has the largest distinguishable difference. First, the extreme points of the new generated populationSnew and constructed fitness tree Stree are determined by identifying the minimum and maximum values zjmin and zjmax for each dimension j=1,2,…,n and by constructing the extreme points zmin=(z1min,z2min,…,znmin) and zmax=(z1max,z2max,…,znmax). Each dimension value of nodes in fitness tree snode is then translated by subtracting the corresponding dimension value of extreme points and dividing the distribution range of corresponding dimension as follows:(10)znsi,j=snode,j-zjminzjmax-zjmin.The pseudocode of adaptive normalization of population members is listed in Algorithm2.Algorithm 2:Pseudocode of adaptive normalization of population members.
Input: Fitness tree and new generated individuals
Output: Fitness tree
(1)  All nodes of Fitness tree Stree=snode1,snode2,… and new generated individuals Snew=s1,s2,…
(2)  S=Stree∪Snew
(3)  For j=1:n
(4)    Compute minimum extreme point of each dimension:
(5)        znewjmin=mins∈Snew⁡si,j
zjmin=min⁡zjmin,znewjmin
(6)    Compute maximum extreme point of each dimension:
(7)        zjmax=maxs∈Snew⁡si,j
zjmax=max⁡zjmax,znewjmax
(8)    For sj=1:num(S)
(9)      Normalize each dimension of each individual znsi,j = snode,j–zjminzjmax–zjmin
(10)   End For
(11) End For


### 3.4. Algorithmic Framework

In this paper, a new multispecies coevolution particle swarm optimization based on previous search history is proposed. TheK-means clustering is applied to dynamically partition the whole population into multiple species. Fitness tree is employed to avoid that the individuals revisit the previous search regions.The pseudocode of the proposed algorithm is listed in Algorithm3.Algorithm 3:Pseudocode of multispecies coevolution particle swarm optimization based on previous search history.
Initialization:
(1)    Stochastically generate a population of N individuals in search space;
2    Set the numbers of clusters C=c1,c2,…,cm and the total number of iteration T;
3    Calculate fitness values of individuals;
(4)    Normalize the dimensions of each individual as Algorithm 2;
(5)    Initialize the fitness tree as Algorithm 1;
(6)    t=1; I=1;
Loop:
(7)  While (t<T)
(8)    If (t>∑ti)
(9)        Partition the whole population into ci multiple species based on K-means clustering;
I=I+1;
(10)      End If
(11)      Update the position of individuals according to Eq. (6) and (7);
(12)      Normalize the dimensions of each individual as Algorithm 2;
(13)      Insert new individuals into fitness tree as Algorithm 1;
(14)      Update some individuals from fitness tree;
(15)      Calculate fitness values of individuals;
(16)      t=t+1
(17) End


## 3.1. Multispecies Coevolution Model for Optimization

Inspired by the mutualism phenomenon, the single population PSO is extended into the interacting multispecies model by enhancing particle dynamics.In such coevolution model, the information for the particle’s position updating includes not only its personal history information and global information of its locating species but also other global information of other species. Information communication simultaneously occurs in and between the local species and neighboring species.In mathematical terms, the original particle’s fly velocity equation and position update equation (please see (2) and (3)) are improved as follows:(6)vi,dkk+1=w·vi,dkt+c1·r1·Pi,dt-xi,dkt+c2·r2·Pg,dkt-xi,dkt+c3·r3·Pg,dθt-xi,dkt,(7)xi,dkt+1=xi,dkt+vi,dkt+1,where the superscript k represents the kth species, xi,dk denotes the position of the ith particle of the kth species, Pi,d is the current particle’s history best position, Pg,dk is the history best position found by its locating species, and Pg,dθ is the history best position found by the whole population. In (6), c1, c2, and c3 are the constant coefficients. And in this paper, c1, c2, and c3 are set to have the same value, c1=c2=c3=1.6667; r1, r2, and r3 are the random values in the range of (0, 1).In (6), c1·r1·Pi,dt–xi,dkt takes into account the individual’s own experiences, c2·r2·Pg,dkt–xi,dkt means the internal communication within the kth species, and c3·r3·Pg,dθt–xi,dkt indicates the symbiotic coevolution between dissimilar species.


## 3.2. DynamicallyK-Means Cluster

In order to implement the multispecies coevolution model, the stochastically generated population should be partitioned inton subspecies. The widely adopted K-means cluster method (please see Section 2.2) is employed for achieving the population division. The main drawback of K-means cluster is that the effect of population division depends on the initial distribution of stochastically generated individuals and the selection of initial division points. A relatively effective means is to set a large number of clusters.In our proposed algorithm, the number of clusters is predefined asC=c1,c2,…,cm, where c1>c2>⋯>cm. During optimization, each individual in every cluster operates independently as the multispecies coevolution model (please see Section 3.1). After several optimization iterations Ti, the whole population is repartitioned into ci+1 clusters; ci and ci+1 are orderly got from C=c1,c2,…,cm. Here, T1+T2+⋯+Tm=T, whereT is the total number of iterations. As ci+1 is smaller than ci, the individuals in small clusters may be redistributed into the relatively larger clusters when the number of the clusters changes.To balance the exploration and exploitation, the value ofTi is not constant. To enhance the ability of local search especially in early stages of searching, Ti should be larger than Ti+1. Therefore, in our proposed algorithm, Ti is calculated as follows:(8)Ti=m-i+1∑j=1mj·T.Moreover, in the process of optimization, the centers of two or more clusters move too close or overlap each other. In order to avoid this phenomenon and to control every cluster to locally search a part of landscape, the relatively small cluster in two or more adjacent clusters decomposes and the individuals in such cluster redistribute into other adjacent clusters. Evaluation criterion is listed as follows:(9)Disclusters-Rcluster1-Rcluster2<0.2·min⁡Rcluster1,Rcluster2,where Disclusters is the distance between two adjacent clusters and Rcluster1 and Rcluster2 are the radiuses of two adjacent clusters (please see (5)); the left term (Disclusters-Rcluster1-Rcluster2) represents the gap between two adjacent clusters.


## 3.3. Adaptive Normalization of Population Members

In the proposed multispecies particle swarm optimization, the idea of fitness tree is also employed to overcome the revisit phenomenon. Here, a modified adaptive normalized fitness tree is proposed.In the process of inserting a new individual into fitness tree (please see Algorithm1), there is a potential drawback: when determining the comparing dimension, it depends on the maximum difference among the dimensions of two child nodes.For overcoming this drawback, the normalization method is employed to determine which dimension has the largest distinguishable difference. First, the extreme points of the new generated populationSnew and constructed fitness tree Stree are determined by identifying the minimum and maximum values zjmin and zjmax for each dimension j=1,2,…,n and by constructing the extreme points zmin=(z1min,z2min,…,znmin) and zmax=(z1max,z2max,…,znmax). Each dimension value of nodes in fitness tree snode is then translated by subtracting the corresponding dimension value of extreme points and dividing the distribution range of corresponding dimension as follows:(10)znsi,j=snode,j-zjminzjmax-zjmin.The pseudocode of adaptive normalization of population members is listed in Algorithm2.Algorithm 2:Pseudocode of adaptive normalization of population members.
Input: Fitness tree and new generated individuals
Output: Fitness tree
(1)  All nodes of Fitness tree Stree=snode1,snode2,… and new generated individuals Snew=s1,s2,…
(2)  S=Stree∪Snew
(3)  For j=1:n
(4)    Compute minimum extreme point of each dimension:
(5)        znewjmin=mins∈Snew⁡si,j
zjmin=min⁡zjmin,znewjmin
(6)    Compute maximum extreme point of each dimension:
(7)        zjmax=maxs∈Snew⁡si,j
zjmax=max⁡zjmax,znewjmax
(8)    For sj=1:num(S)
(9)      Normalize each dimension of each individual znsi,j = snode,j–zjminzjmax–zjmin
(10)   End For
(11) End For


## 3.4. Algorithmic Framework

In this paper, a new multispecies coevolution particle swarm optimization based on previous search history is proposed. TheK-means clustering is applied to dynamically partition the whole population into multiple species. Fitness tree is employed to avoid that the individuals revisit the previous search regions.The pseudocode of the proposed algorithm is listed in Algorithm3.Algorithm 3:Pseudocode of multispecies coevolution particle swarm optimization based on previous search history.
Initialization:
(1)    Stochastically generate a population of N individuals in search space;
2    Set the numbers of clusters C=c1,c2,…,cm and the total number of iteration T;
3    Calculate fitness values of individuals;
(4)    Normalize the dimensions of each individual as Algorithm 2;
(5)    Initialize the fitness tree as Algorithm 1;
(6)    t=1; I=1;
Loop:
(7)  While (t<T)
(8)    If (t>∑ti)
(9)        Partition the whole population into ci multiple species based on K-means clustering;
I=I+1;
(10)      End If
(11)      Update the position of individuals according to Eq. (6) and (7);
(12)      Normalize the dimensions of each individual as Algorithm 2;
(13)      Insert new individuals into fitness tree as Algorithm 1;
(14)      Update some individuals from fitness tree;
(15)      Calculate fitness values of individuals;
(16)      t=t+1
(17) End


## 4. Experimental Study

### 4.1. Experimental Settings

The proposed multispecies coevolution particle swarm optimization based on previous search history (MCPSO-PSH) was executed with four state-of-the-art EA and SI algorithms for comparisons: (1) canonical particle swarm optimization (PSO) [24]; (2) differential evolution (DE) [27]; (3) Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [28]; (4) nonrevisiting genetic algorithm (NrGA) [26].In all experiments in this section, the values of the common parameters used in each algorithm such as population size and total generation number are chosen to be the same. Population size is set as 100 and the maximum evaluation number is set as 100000. For the continuous testing functions used in this paper, the dimensions are all set as 10D and 30D.All the control parameters for the EA and SI algorithms are set to be default of their original literatures.For canonical PSO, the learning ratesc1 and c2 are both set as 2.0.For canonical DE, the zoom factor is set asF=0.9 and crossover rate is set as CR=0.9.For CMA-ES, the population sizeλ is chosen by the suggested setting in [28] (i.e., λ=4+3∗ln⁡D).For NrGA, the crossover rate of 0.8 and the mutation rate of 0.01 are adopted [26].For our proposed algorithm, the number of clusters is set asC={10,5,4,2,1}, and so T is calculated as T=T1,T2,T3,T4,T5=333,267,200,133,67. The learning rates c1, c2, and c3 in (6) are all set as c1=c2=c3=1.3667.


### 4.2. Test and Benchmark Problems

According to the No Free Lunch theorem, “for any algorithm, any elevated performance over one class of problems is exactly paid for in performance over another class” [29]. To fully evaluate the performance of the MCPSO-PSH without a biased conclusion towards some chosen problems, a set of ten basic benchmark functions and ten CEC2005 benchmark functions [30] is employed. The formulas of these functions are presented below.


#### 4.2.1. Basic Benchmark Functions

(1) Rosenbrock function(11)f1x=∑i=1n100×xi+1-xi22+1-xi2.(2) Quadric function(12)f2x=∑i=1n∑j=1ixj2.(3) Sum of different powers(13)f3x=∑i=1nxii+1.(4) Axis parallel hyperellipsoid function(14)f4x=∑i=1ni·xi2.(5) Sphere function(15)f5x=∑i=1nxi2.(6) Ackley’s function(16)f6x=-20exp⁡-0.21n∑i=1nxi2-exp⁡1n∑i=1ncos⁡2πxi+20+e.(7) Griewank’s function(17)f7x=14000∑i=1nxi2-∏i=1ncos⁡xii+1.(8) Michalewicz function(18)f8x=-∑i=1nsin⁡xi·sin⁡i·xi2π2·m,m=10.(9) Rastrigin’s function(19)f9x=∑i=1nxi2-10cos⁡2πxi+10.(10) Schwefel’s function(20)f10x=∑i=0nxi+∏i=1nxi.These ten benchmark functions can be grouped as unimodal continuous functionsf1–f5 and multimodal continuous functionsf6–f10.


#### 4.2.2. CEC2005 Benchmark Functions

(11) Shifted sphere function(21)f11x=∑i=1Dzi2+f_bias11,z=x-o,x=x1,x2,…,xD.(12) Shifted Schwefel’s problem 1.2(22)f12x=∑i=1D∑j=1izj2+f_bias12,z=x-o,x=x1,x2,…,xD.(13) Shifted rotated high conditioned elliptic function(23)f13x=∑i=1D106i-1/D-1zi2+f_bias13,z=x-o∗M,x=x1,x2,…,xD.(14) Shifted Schwefel’s problem 1.2 with noise in fitness(24)f14x=∑i=1D100zi2-zi+12+zi-12+f_bias14,z=x-o+1,x=x1,x2,…,xD.(15) Schwefel’s Problem 2.6 on bounds(25)f15x=max⁡Aix-Bi+f_bias15,i=1,…,D,x=x1,x2,…,xD.(16) Shifted Rosenbrock’s function(26)f16x=∑i=1D∑j=1izj2∗1+0.4N0,1+f_bias16,z=x-o,x=x1,x2,…,xD.(17) Shifted rotated Griewank’s function without bounds(27)f17x=∑i=1Dzi24000-∏i=1Dcos⁡zii+1+f_bias17,z=x-o∗M,x=x1,x2,…,xD.(18) Shifted rotated Ackley’s function on bounds(28)f18x=-20exp⁡-0.21D∑i=1Dzi2-exp⁡1D∑i=1Dcos⁡2πzi+20+e+f_bias18,z=x-o∗M.(19) Shifted Rastrigin’s function(29)f19x=∑i=1Dzi2-10cos⁡2πzi+10+f_bias19,z=x-o∗M,x=x1,x2,…,xD.(20) Schwefel’s Problem 2.13(30)f20x=∑i=1DAi-Bix+f_bias20,x=x1,x2,…,xD,Ai=∑j=1Daijsin⁡αj+bijcos⁡αj,Bi=∑j=1Daijsin⁡xj+bijcos⁡xj.These ten benchmark functions can be grouped as unimodal continuous functionsf11–f15 and multimodal continuous functionsf16–f20.


### 4.3. Result and Statistical Analysis

#### 4.3.1. Results for Basic Benchmark Functions

In this experiment, the proposed MCPSO-PSH algorithm is compared with nonrevisiting genetic algorithm (NrGA), Covariance Matrix Adaptation Evolution Strategy (CMA-ES), the canonical particle swarm optimization (PSO) algorithm, and differential evolution (DE) algorithm.For testing 10D benchmark test functionsf1–f10 in 10 runs, the basic statistical results including the mean and standard deviations of the function values are listed in Table 1. From Table 1, it is obvious that the canonical PSO and DE algorithms demonstrate a worse performance than MCPSO-PSH, NrGA, and CMA-ES. For the low-dimensional (10D) unimodal functionsf1–f5, the uniformity of CMA-ES algorithm is relatively good. For functionf1, CMA-ES gets the best results. And, for functionsf2–f5, the performance of CMA-ES is rather competitive. Compared with the canonical PSO and other state-of-the-art algorithms, MCPSO-PSH brings in the dynamical multispecies strategy based onK-means to enhance the ability of local search and fitness tree to overcome the revisit phenomenon, which results in a huge improvement of performance. The MCPSO-PSH and NrGA algorithms with nonrevisit strategy have nearly the same results for low-dimensional problems. Benefitting from the dynamicK-means clustering and the adaptive normalization of population members, the improved multispecies PSO has a slightly better performance than NrGA on functionsf1–f5. For the low-dimensional (10D) multimodal functionsf6–f10, intuitively, the proposed MCPSO-PSH algorithm has rather smaller mean values and standard deviations than other algorithms for functionf6. For functionsf7 andf8, the MCPSO-PSH and NrGA obtain the global optimum with a high efficiency. It is due to the fact that the nonrevisit strategy based on BSP fitness tree effectively enhances the global search especially when a large number of local optima exist in the search space. For functionsf9 andf10, the CMA-ES algorithm also gets the global optimum because the MCPSO-PSH, NrGA, and CMA-ES own their special mechanisms for jumping out of the local optimum. On the whole, regardless of the unimodal functions or the multimodal functions, the performance of the proposed MCPSO-PSH is uniform.Table 1
Performance of all algorithms on 10D test functionsf1–f10.
Problem
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
f
1
Mean
4.3136
E
-
02
3.9407
E
-
01
9.4574E − 05
7.3831
E
+
04
3.0851
E
+
13
Std
7.3787
E
-
07
3.1066
E
-
05
1.4783E − 05
6.3480
E
+
03
1.0345
E
+
12
f
2
Mean
8.0378E − 33
3.7721
E
-
27
8.4991
E
-
18
1.0386
E
-
03
3.9501
E
-
01
Std
1.0378E − 30
7.0550
E
-
30
3.8731
E
-
21
2.7801
E
-
03
4.0404
E
+
00
f
3
Mean
9.4931E − 38
5.7707
E
-
33
3.4959
E
-
14
1.4301
E
+
47
3.5430
E
-
03
Std
3.9093E − 41
7.3900
E
-
37
8.0455
E
-
20
1.3784
E
+
47
5.7791
E
-
04
f
4
Mean
2.9043E − 73
3.8844
E
-
36
8.3995
E
-
18
9.4034
E
+
08
3.1893
E
-
43
Std
4.7898E − 99
1.3894
E
-
45
2.9393
E
-
20
1.0373
E
+
09
7.8440
E
-
53
f
5
Mean
4.9550E − 64
3.9090
E
-
39
4.4139
E
-
43
8.0041
E
-
17
6.0903
E
-
02
Std
4.8414E − 89
3.4568
E
-
40
6.3941
E
-
46
9.4389
E
-
18
5.9031
E
-
02
f
6
Mean
2.2137E − 14
7.3678
E
-
13
3.8931
E
-
13
7.3491
E
+
01
1.3940
E
+
02
Std
3.7890E − 17
2.4654
E
-
14
3.9404
E
-
14
3.4791
E
-
03
4.9448
E
+
01
f
7
Mean
0
0
7.4327
E
-
03
3.9051
E
+
00
2.6461
E
-
2
Std
0
0
4.7801
E
-
03
9.4803
E
-
01
8.4780
E
-
3
f
8
Mean
−9.8452E − 01
−9.8452E − 01
9.0377
E
-
01
2.0773
E
-
01
2.4743
E
+
00
Std
0
0
6.9010
E
-
01
1.3434
E
-
01
9.0370
E
-
1
f
9
Mean
0
0
0
3.8944
E
+
02
1.9640
E
+
03
Std
0
0
0
9.0378
E
+
01
8.9901
E
+
01
f
10
Mean
0
0
0
1.3355
E
+
04
7.9140
E
+
04
Std
0
0
0
8.9440
E
+
03
4.4708
E
+
04For testing 30D benchmark test functionsf1–f10 in 10 runs, the basic statistical results including the mean and standard deviations of the function values are listed in Table 2. With the increase of the dimensions, the convergence rate becomes slow. The results of Table 2 show great disparities in convergence accuracy. For the high-dimensional unimodal problemsf1–f5, the performance of MCPSO-PSH algorithm is more outstanding than the other four algorithms because MCPSO-PSH with dynamic clustering strategy can accelerate the convergence rate especially in the later evolution period. The performance of NrGA algorithm becomes slightly poor although it also has the nonrevisit strategy for improving the search efficiency. Only for functionsf1–f5 does CMA-ES outperform NrGA on convergence rate and accuracy. For solving the high-dimensional multimodal problemsf6–f10, the MCPSO-PSH has a good performance. With the dynamic multispecies strategy based onK-means clustering and nonrevisit strategy based on BSP fitness tree, the MCSPSO-PSH algorithm can well avoid premature convergence to local optimum in the process of solving the high-dimensional multimodal problems. From the basic statistical results, NrGA also owns the most outstanding performance on testing functionsf6 andf10, which means the nonrevisit strategy based on BSP fitness tree plays an important role in guiding search direction. On the whole, the performance of all algorithms tested here is ordered as MCPSO-PSH > NrGA > CMA-ES > PSO > DE.Table 2
Performance of all algorithms on 30D test functionsf1–f10.
Problem
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
f
1
Mean
4.8910E − 01
1.3047
E
+
01
3.1794
E
+
00
3.4903
E
+
12
2.4941
E
+
21
Std
3.3890E − 01
9.0360
E
+
00
2
.
0935
E
+
00
8.0410
E
+
11
9.3892
E
+
19
f
2
Mean
1.6851E − 14
1.8440
E
-
08
5.3882
E
+
02
3.7840
E
+
04
9.9038
E
8
Std
9.0376E − 16
9.5040
E
-
10
3.0318
E
+
02
9.9011
E
+
03
1.3931
E
7
f
3
Mean
5.9990E − 37
6.3890
E
+
01
1.4890
E
-
12
3.8964
E
+
308
7.3280
E
+
223
Std
9.5490E − 41
1.4378
E
+
01
2.1809
E
-
13
1.9430
E
+
310
1.9653
E
+
224
f
4
Mean
1.0099E − 16
8.3701
E
-
08
1.0454
E
-
23
7.9016
E
+
13
3.7890
E
+
08
Std
8.4098E − 18
1.0871
E
-
07
9.5487
E
-
25
1.0346
E
+
14
1.3944
E
+
08
f
5
Mean
1.3933E − 61
2.0031
E
-
11
2.9890
E
-
16
2.3498
E
+
02
5.3903
E
+
05
Std
9.0809E − 75
9.6036
E
-
13
7.3892
E
-
23
7.3954
E
+
01
1.2309
E
+
05
f
6
Mean
1.3343E − 14
1.3408
E
-
12
6.3378
E
-
06
8.4723
E
+
02
2.5044
E
+
02
Std
6.6066E − 15
8.0651
E
-
14
1.6999
E
-
06
3.1704
E
+
02
7.3909
E
+
01
f
7
Mean
0
0
3.0966
E
-
2
2.0330
E
-
01
8.9036
E
-
01
Std
0
0
7.8012
E
-
3
3.3891
E
-
01
6.7009
E
-
01
f
8
Mean
−2.8743E + 01
-
2.7945
E
+
01
1.7640
E
+
01
6.2307
E
-
01
-
3.3883
E
+
00
Std
1.664E + 00
2.0874
E
+
00
7.3370
E
+
00
3.1700
E
-
01
6.0310
E
-
01
f
9
Mean
3.3880E − 12
8.5801
E
+
00
2.1388
E
+
00
4.3880
E
+
04
8.4891
E
+
03
Std
7.3038E − 13
5.8817
E
+
00
2.9863
E
+
00
8.3909
E
+
03
2.7450
E
+
03
f
10
Mean
7.9063
E
-
01
3.6421E − 04
2.4578
E
+
03
2.3841
E
+
04
4.4093
E
+
05
Std
4.1939
E
-
01
1.3940E − 04
2.7309
E
+
03
1.7301
E
+
04
1.0938
E
+
06The convergence curves of the MCPSO-PSH algorithms and other EAs against test functionsf1–f10 in 10D and 30D are shown in Figure 2. It can be observed that, with the increase of the dimensions, the convergence rate of the proposed MCPSO-PSH algorithm becomes faster. Predictably, with a further increase of the dimensions of the test functions, the MCPSO-PSH algorithm will have an overwhelming superiority in solving the high-dimensional problems.Figure 2
Convergence curves of MCPSO-PSH, NrGA, CMA-ES, PSO, and DE against test functionsf1–f10 in low dimension, 10D, and high dimension, 30D, respectively.


#### 4.3.2. Results for CEC2005 Benchmark Functions

The statistical results in Table3 show that MCPSO-PSH has significantly outstanding performance for all these cases except f14 and f20. For f14, although all the algorithms provide similar results (except DE), the solutions of MCPSO-PSH are still rather competitive and accessible. Forf20 MCPSO-PSH achieved the second best solution just behind NrGA. The remarkable performance of MCPSO-PSH for f11–f15 indicates that MCPSO-PSH also applies to unimodal functions though it is primarily designed for multimodal functions. For the multimodal functions f16–f19, MCPSO-PSH outperforms NrGA, CMA-ES, PSO, and DE. In terms of average obtained results, MCPSO-PSH is better than NrGA on f12,f13,f15,f16,f17, and f18, and is somehow equivalent to it on f11 and f19. MCPSO-PSH and NrGA outperform other three algorithms on most cases, which demonstrates great promise of nonrevisiting strategy for maintaining population diversity and improving search efficiency synchronously. Obviously, MCPSO-PSH is superior to NrGA which benefits from its unique dynamic multispecies strategy based on K-means clustering.Table 3
Basic statistics off11–f20 CEC2005 benchmark test functions obtained by MCPSO-PSH, NrGA, CMA-ES, PSO, and DE.
Problem
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
f
11
Mean
−4.5000E + 02
−4.5000E + 02
-
3.1857
E
+
02
-
4.3912
E
+
02
-
4.4313
E
+
02
Std
0
0
1.8576
E
+
02
2.3
441
E
+
01
9.7454
E
+
00
f
12
Mean
−4.5000E + 02
3.1554
E
+
03
-
2.2058
E
+
02
1.2360
E
+
03
-
4.4506
E
+
02
Std
0
1.9571
E
+
03
8.3398
E
+
01
3.2960
E
+
02
3.6262
E
-
09
f
13
Mean
−4.5000E + 02
6.9857
E
+
06
2.4164
E
+
05
2.5967
E
+
07
3.5243
E
+
07
Std
0
2.5378
E
+
06
3.0876
E
+
04
4.3459
E
+
06
6.5900
E
+
07
f
14
Mean
3.8043
E
+
02
2.3818E + 02
5.3723
E
+
02
2.7801
E
+
02
3.5976
E
+
03
Std
3.2991
E
+
00
6.8326E + 00
6.3167
E
+
01
8.3767
E
+
00
3.2048
E
+
03
f
15
Mean
2.3174E + 03
1.6341
E
+
04
5.2245
E
+
03
3.3347
E
+
03
3.8727
E
+
03
Std
5.0134E + 02
1.0450
E
+
03
2.1529
E
+
03
5.7811
E
+
02
4.3668
E
+
03
f
16
Mean
−4.5000E + 02
2.1850
E
+
04
-
1.8444
E
+
02
6.5597
E
+
03
-
4.4761
E
+
02
Std
0
3.1742
E
+
03
1.2313
E
+
02
1.3545
E
+
03
3.4378
E
-
12
f
17
Mean
−1.8000E + 02
1.9705
E
+
03
1.8865
E
+
04
4.6964
E
+
03
−1.8000E + 02
Std
0
1.9839
E
+
02
3.9520
E
+
02
8.4681
E
+
02
0
f
18
Mean
−1.1997E + 02
-
1.1911
E
+
02
-
1.1958
E
+
02
-
1.1909
E
+
02
-
1.1995
E
+
02
Std
8.9147E − 03
7.8361
E
-
02
3.7314
E
-
01
4.2087
E
-
02
6.0499
E
-
01
f
19
Mean
−3.3000E + 02
−3.300e + 02
-
1.9906
E
+
02
-
3.1776
E
+
02
-
3.2842
E
+
02
Std
0
0
1.6472
E
+
01
7.5984
E
+
01
8.5989
E
+
01
f
20
Mean
-
4.3294
E
+
02
−4.3933E + 02
-
4.2932
E
+
02
-
3.9982
E
+
02
-
4.0363
E
+
02
Std
3.1803
E
+
01
1.7580E + 01
2.0839
E
+
01
6.1074
E
+
01
3.5741
E
+
01In order to intuitively depict the search effect, the box plots are employed to demonstrate the statistical performance of five different algorithms, as shown in Figure3. Although nonrevisiting strategy is implemented in both NrGA and MCPSO-PSH, the robustness of NrGA is still undesirable. MCPSO-PSH with dynamic multispecies strategy has great enhancement of exploring ability.Figure 3
Box plots of results obtained by all algorithms. Here, 1 to 5 in the horizontal axis are indices of MCPSO-PSH, NrGA, CMA-ES, PSO, and DE, respectively (the box has lines at the lower quartile, median, and upper quartile values. The whiskers extend to the most extreme data points not considered outliers. Outliers (denoted by +) are data with values beyond 100 units of interquartile range).From the above observations, the superiority of MCPSO-PSH over other peer algorithms can be easily presumed.


## 4.1. Experimental Settings

The proposed multispecies coevolution particle swarm optimization based on previous search history (MCPSO-PSH) was executed with four state-of-the-art EA and SI algorithms for comparisons: (1) canonical particle swarm optimization (PSO) [24]; (2) differential evolution (DE) [27]; (3) Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [28]; (4) nonrevisiting genetic algorithm (NrGA) [26].In all experiments in this section, the values of the common parameters used in each algorithm such as population size and total generation number are chosen to be the same. Population size is set as 100 and the maximum evaluation number is set as 100000. For the continuous testing functions used in this paper, the dimensions are all set as 10D and 30D.All the control parameters for the EA and SI algorithms are set to be default of their original literatures.For canonical PSO, the learning ratesc1 and c2 are both set as 2.0.For canonical DE, the zoom factor is set asF=0.9 and crossover rate is set as CR=0.9.For CMA-ES, the population sizeλ is chosen by the suggested setting in [28] (i.e., λ=4+3∗ln⁡D).For NrGA, the crossover rate of 0.8 and the mutation rate of 0.01 are adopted [26].For our proposed algorithm, the number of clusters is set asC={10,5,4,2,1}, and so T is calculated as T=T1,T2,T3,T4,T5=333,267,200,133,67. The learning rates c1, c2, and c3 in (6) are all set as c1=c2=c3=1.3667.


## 4.2. Test and Benchmark Problems

According to the No Free Lunch theorem, “for any algorithm, any elevated performance over one class of problems is exactly paid for in performance over another class” [29]. To fully evaluate the performance of the MCPSO-PSH without a biased conclusion towards some chosen problems, a set of ten basic benchmark functions and ten CEC2005 benchmark functions [30] is employed. The formulas of these functions are presented below.


### 4.2.1. Basic Benchmark Functions

(1) Rosenbrock function(11)f1x=∑i=1n100×xi+1-xi22+1-xi2.(2) Quadric function(12)f2x=∑i=1n∑j=1ixj2.(3) Sum of different powers(13)f3x=∑i=1nxii+1.(4) Axis parallel hyperellipsoid function(14)f4x=∑i=1ni·xi2.(5) Sphere function(15)f5x=∑i=1nxi2.(6) Ackley’s function(16)f6x=-20exp⁡-0.21n∑i=1nxi2-exp⁡1n∑i=1ncos⁡2πxi+20+e.(7) Griewank’s function(17)f7x=14000∑i=1nxi2-∏i=1ncos⁡xii+1.(8) Michalewicz function(18)f8x=-∑i=1nsin⁡xi·sin⁡i·xi2π2·m,m=10.(9) Rastrigin’s function(19)f9x=∑i=1nxi2-10cos⁡2πxi+10.(10) Schwefel’s function(20)f10x=∑i=0nxi+∏i=1nxi.These ten benchmark functions can be grouped as unimodal continuous functionsf1–f5 and multimodal continuous functionsf6–f10.


### 4.2.2. CEC2005 Benchmark Functions

(11) Shifted sphere function(21)f11x=∑i=1Dzi2+f_bias11,z=x-o,x=x1,x2,…,xD.(12) Shifted Schwefel’s problem 1.2(22)f12x=∑i=1D∑j=1izj2+f_bias12,z=x-o,x=x1,x2,…,xD.(13) Shifted rotated high conditioned elliptic function(23)f13x=∑i=1D106i-1/D-1zi2+f_bias13,z=x-o∗M,x=x1,x2,…,xD.(14) Shifted Schwefel’s problem 1.2 with noise in fitness(24)f14x=∑i=1D100zi2-zi+12+zi-12+f_bias14,z=x-o+1,x=x1,x2,…,xD.(15) Schwefel’s Problem 2.6 on bounds(25)f15x=max⁡Aix-Bi+f_bias15,i=1,…,D,x=x1,x2,…,xD.(16) Shifted Rosenbrock’s function(26)f16x=∑i=1D∑j=1izj2∗1+0.4N0,1+f_bias16,z=x-o,x=x1,x2,…,xD.(17) Shifted rotated Griewank’s function without bounds(27)f17x=∑i=1Dzi24000-∏i=1Dcos⁡zii+1+f_bias17,z=x-o∗M,x=x1,x2,…,xD.(18) Shifted rotated Ackley’s function on bounds(28)f18x=-20exp⁡-0.21D∑i=1Dzi2-exp⁡1D∑i=1Dcos⁡2πzi+20+e+f_bias18,z=x-o∗M.(19) Shifted Rastrigin’s function(29)f19x=∑i=1Dzi2-10cos⁡2πzi+10+f_bias19,z=x-o∗M,x=x1,x2,…,xD.(20) Schwefel’s Problem 2.13(30)f20x=∑i=1DAi-Bix+f_bias20,x=x1,x2,…,xD,Ai=∑j=1Daijsin⁡αj+bijcos⁡αj,Bi=∑j=1Daijsin⁡xj+bijcos⁡xj.These ten benchmark functions can be grouped as unimodal continuous functionsf11–f15 and multimodal continuous functionsf16–f20.


## 4.2.1. Basic Benchmark Functions

(1) Rosenbrock function(11)f1x=∑i=1n100×xi+1-xi22+1-xi2.(2) Quadric function(12)f2x=∑i=1n∑j=1ixj2.(3) Sum of different powers(13)f3x=∑i=1nxii+1.(4) Axis parallel hyperellipsoid function(14)f4x=∑i=1ni·xi2.(5) Sphere function(15)f5x=∑i=1nxi2.(6) Ackley’s function(16)f6x=-20exp⁡-0.21n∑i=1nxi2-exp⁡1n∑i=1ncos⁡2πxi+20+e.(7) Griewank’s function(17)f7x=14000∑i=1nxi2-∏i=1ncos⁡xii+1.(8) Michalewicz function(18)f8x=-∑i=1nsin⁡xi·sin⁡i·xi2π2·m,m=10.(9) Rastrigin’s function(19)f9x=∑i=1nxi2-10cos⁡2πxi+10.(10) Schwefel’s function(20)f10x=∑i=0nxi+∏i=1nxi.These ten benchmark functions can be grouped as unimodal continuous functionsf1–f5 and multimodal continuous functionsf6–f10.


## 4.2.2. CEC2005 Benchmark Functions

(11) Shifted sphere function(21)f11x=∑i=1Dzi2+f_bias11,z=x-o,x=x1,x2,…,xD.(12) Shifted Schwefel’s problem 1.2(22)f12x=∑i=1D∑j=1izj2+f_bias12,z=x-o,x=x1,x2,…,xD.(13) Shifted rotated high conditioned elliptic function(23)f13x=∑i=1D106i-1/D-1zi2+f_bias13,z=x-o∗M,x=x1,x2,…,xD.(14) Shifted Schwefel’s problem 1.2 with noise in fitness(24)f14x=∑i=1D100zi2-zi+12+zi-12+f_bias14,z=x-o+1,x=x1,x2,…,xD.(15) Schwefel’s Problem 2.6 on bounds(25)f15x=max⁡Aix-Bi+f_bias15,i=1,…,D,x=x1,x2,…,xD.(16) Shifted Rosenbrock’s function(26)f16x=∑i=1D∑j=1izj2∗1+0.4N0,1+f_bias16,z=x-o,x=x1,x2,…,xD.(17) Shifted rotated Griewank’s function without bounds(27)f17x=∑i=1Dzi24000-∏i=1Dcos⁡zii+1+f_bias17,z=x-o∗M,x=x1,x2,…,xD.(18) Shifted rotated Ackley’s function on bounds(28)f18x=-20exp⁡-0.21D∑i=1Dzi2-exp⁡1D∑i=1Dcos⁡2πzi+20+e+f_bias18,z=x-o∗M.(19) Shifted Rastrigin’s function(29)f19x=∑i=1Dzi2-10cos⁡2πzi+10+f_bias19,z=x-o∗M,x=x1,x2,…,xD.(20) Schwefel’s Problem 2.13(30)f20x=∑i=1DAi-Bix+f_bias20,x=x1,x2,…,xD,Ai=∑j=1Daijsin⁡αj+bijcos⁡αj,Bi=∑j=1Daijsin⁡xj+bijcos⁡xj.These ten benchmark functions can be grouped as unimodal continuous functionsf11–f15 and multimodal continuous functionsf16–f20.


## 4.3. Result and Statistical Analysis

### 4.3.1. Results for Basic Benchmark Functions

In this experiment, the proposed MCPSO-PSH algorithm is compared with nonrevisiting genetic algorithm (NrGA), Covariance Matrix Adaptation Evolution Strategy (CMA-ES), the canonical particle swarm optimization (PSO) algorithm, and differential evolution (DE) algorithm.For testing 10D benchmark test functionsf1–f10 in 10 runs, the basic statistical results including the mean and standard deviations of the function values are listed in Table 1. From Table 1, it is obvious that the canonical PSO and DE algorithms demonstrate a worse performance than MCPSO-PSH, NrGA, and CMA-ES. For the low-dimensional (10D) unimodal functionsf1–f5, the uniformity of CMA-ES algorithm is relatively good. For functionf1, CMA-ES gets the best results. And, for functionsf2–f5, the performance of CMA-ES is rather competitive. Compared with the canonical PSO and other state-of-the-art algorithms, MCPSO-PSH brings in the dynamical multispecies strategy based onK-means to enhance the ability of local search and fitness tree to overcome the revisit phenomenon, which results in a huge improvement of performance. The MCPSO-PSH and NrGA algorithms with nonrevisit strategy have nearly the same results for low-dimensional problems. Benefitting from the dynamicK-means clustering and the adaptive normalization of population members, the improved multispecies PSO has a slightly better performance than NrGA on functionsf1–f5. For the low-dimensional (10D) multimodal functionsf6–f10, intuitively, the proposed MCPSO-PSH algorithm has rather smaller mean values and standard deviations than other algorithms for functionf6. For functionsf7 andf8, the MCPSO-PSH and NrGA obtain the global optimum with a high efficiency. It is due to the fact that the nonrevisit strategy based on BSP fitness tree effectively enhances the global search especially when a large number of local optima exist in the search space. For functionsf9 andf10, the CMA-ES algorithm also gets the global optimum because the MCPSO-PSH, NrGA, and CMA-ES own their special mechanisms for jumping out of the local optimum. On the whole, regardless of the unimodal functions or the multimodal functions, the performance of the proposed MCPSO-PSH is uniform.Table 1
Performance of all algorithms on 10D test functionsf1–f10.
Problem
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
f
1
Mean
4.3136
E
-
02
3.9407
E
-
01
9.4574E − 05
7.3831
E
+
04
3.0851
E
+
13
Std
7.3787
E
-
07
3.1066
E
-
05
1.4783E − 05
6.3480
E
+
03
1.0345
E
+
12
f
2
Mean
8.0378E − 33
3.7721
E
-
27
8.4991
E
-
18
1.0386
E
-
03
3.9501
E
-
01
Std
1.0378E − 30
7.0550
E
-
30
3.8731
E
-
21
2.7801
E
-
03
4.0404
E
+
00
f
3
Mean
9.4931E − 38
5.7707
E
-
33
3.4959
E
-
14
1.4301
E
+
47
3.5430
E
-
03
Std
3.9093E − 41
7.3900
E
-
37
8.0455
E
-
20
1.3784
E
+
47
5.7791
E
-
04
f
4
Mean
2.9043E − 73
3.8844
E
-
36
8.3995
E
-
18
9.4034
E
+
08
3.1893
E
-
43
Std
4.7898E − 99
1.3894
E
-
45
2.9393
E
-
20
1.0373
E
+
09
7.8440
E
-
53
f
5
Mean
4.9550E − 64
3.9090
E
-
39
4.4139
E
-
43
8.0041
E
-
17
6.0903
E
-
02
Std
4.8414E − 89
3.4568
E
-
40
6.3941
E
-
46
9.4389
E
-
18
5.9031
E
-
02
f
6
Mean
2.2137E − 14
7.3678
E
-
13
3.8931
E
-
13
7.3491
E
+
01
1.3940
E
+
02
Std
3.7890E − 17
2.4654
E
-
14
3.9404
E
-
14
3.4791
E
-
03
4.9448
E
+
01
f
7
Mean
0
0
7.4327
E
-
03
3.9051
E
+
00
2.6461
E
-
2
Std
0
0
4.7801
E
-
03
9.4803
E
-
01
8.4780
E
-
3
f
8
Mean
−9.8452E − 01
−9.8452E − 01
9.0377
E
-
01
2.0773
E
-
01
2.4743
E
+
00
Std
0
0
6.9010
E
-
01
1.3434
E
-
01
9.0370
E
-
1
f
9
Mean
0
0
0
3.8944
E
+
02
1.9640
E
+
03
Std
0
0
0
9.0378
E
+
01
8.9901
E
+
01
f
10
Mean
0
0
0
1.3355
E
+
04
7.9140
E
+
04
Std
0
0
0
8.9440
E
+
03
4.4708
E
+
04For testing 30D benchmark test functionsf1–f10 in 10 runs, the basic statistical results including the mean and standard deviations of the function values are listed in Table 2. With the increase of the dimensions, the convergence rate becomes slow. The results of Table 2 show great disparities in convergence accuracy. For the high-dimensional unimodal problemsf1–f5, the performance of MCPSO-PSH algorithm is more outstanding than the other four algorithms because MCPSO-PSH with dynamic clustering strategy can accelerate the convergence rate especially in the later evolution period. The performance of NrGA algorithm becomes slightly poor although it also has the nonrevisit strategy for improving the search efficiency. Only for functionsf1–f5 does CMA-ES outperform NrGA on convergence rate and accuracy. For solving the high-dimensional multimodal problemsf6–f10, the MCPSO-PSH has a good performance. With the dynamic multispecies strategy based onK-means clustering and nonrevisit strategy based on BSP fitness tree, the MCSPSO-PSH algorithm can well avoid premature convergence to local optimum in the process of solving the high-dimensional multimodal problems. From the basic statistical results, NrGA also owns the most outstanding performance on testing functionsf6 andf10, which means the nonrevisit strategy based on BSP fitness tree plays an important role in guiding search direction. On the whole, the performance of all algorithms tested here is ordered as MCPSO-PSH > NrGA > CMA-ES > PSO > DE.Table 2
Performance of all algorithms on 30D test functionsf1–f10.
Problem
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
f
1
Mean
4.8910E − 01
1.3047
E
+
01
3.1794
E
+
00
3.4903
E
+
12
2.4941
E
+
21
Std
3.3890E − 01
9.0360
E
+
00
2
.
0935
E
+
00
8.0410
E
+
11
9.3892
E
+
19
f
2
Mean
1.6851E − 14
1.8440
E
-
08
5.3882
E
+
02
3.7840
E
+
04
9.9038
E
8
Std
9.0376E − 16
9.5040
E
-
10
3.0318
E
+
02
9.9011
E
+
03
1.3931
E
7
f
3
Mean
5.9990E − 37
6.3890
E
+
01
1.4890
E
-
12
3.8964
E
+
308
7.3280
E
+
223
Std
9.5490E − 41
1.4378
E
+
01
2.1809
E
-
13
1.9430
E
+
310
1.9653
E
+
224
f
4
Mean
1.0099E − 16
8.3701
E
-
08
1.0454
E
-
23
7.9016
E
+
13
3.7890
E
+
08
Std
8.4098E − 18
1.0871
E
-
07
9.5487
E
-
25
1.0346
E
+
14
1.3944
E
+
08
f
5
Mean
1.3933E − 61
2.0031
E
-
11
2.9890
E
-
16
2.3498
E
+
02
5.3903
E
+
05
Std
9.0809E − 75
9.6036
E
-
13
7.3892
E
-
23
7.3954
E
+
01
1.2309
E
+
05
f
6
Mean
1.3343E − 14
1.3408
E
-
12
6.3378
E
-
06
8.4723
E
+
02
2.5044
E
+
02
Std
6.6066E − 15
8.0651
E
-
14
1.6999
E
-
06
3.1704
E
+
02
7.3909
E
+
01
f
7
Mean
0
0
3.0966
E
-
2
2.0330
E
-
01
8.9036
E
-
01
Std
0
0
7.8012
E
-
3
3.3891
E
-
01
6.7009
E
-
01
f
8
Mean
−2.8743E + 01
-
2.7945
E
+
01
1.7640
E
+
01
6.2307
E
-
01
-
3.3883
E
+
00
Std
1.664E + 00
2.0874
E
+
00
7.3370
E
+
00
3.1700
E
-
01
6.0310
E
-
01
f
9
Mean
3.3880E − 12
8.5801
E
+
00
2.1388
E
+
00
4.3880
E
+
04
8.4891
E
+
03
Std
7.3038E − 13
5.8817
E
+
00
2.9863
E
+
00
8.3909
E
+
03
2.7450
E
+
03
f
10
Mean
7.9063
E
-
01
3.6421E − 04
2.4578
E
+
03
2.3841
E
+
04
4.4093
E
+
05
Std
4.1939
E
-
01
1.3940E − 04
2.7309
E
+
03
1.7301
E
+
04
1.0938
E
+
06The convergence curves of the MCPSO-PSH algorithms and other EAs against test functionsf1–f10 in 10D and 30D are shown in Figure 2. It can be observed that, with the increase of the dimensions, the convergence rate of the proposed MCPSO-PSH algorithm becomes faster. Predictably, with a further increase of the dimensions of the test functions, the MCPSO-PSH algorithm will have an overwhelming superiority in solving the high-dimensional problems.Figure 2
Convergence curves of MCPSO-PSH, NrGA, CMA-ES, PSO, and DE against test functionsf1–f10 in low dimension, 10D, and high dimension, 30D, respectively.


### 4.3.2. Results for CEC2005 Benchmark Functions

The statistical results in Table3 show that MCPSO-PSH has significantly outstanding performance for all these cases except f14 and f20. For f14, although all the algorithms provide similar results (except DE), the solutions of MCPSO-PSH are still rather competitive and accessible. Forf20 MCPSO-PSH achieved the second best solution just behind NrGA. The remarkable performance of MCPSO-PSH for f11–f15 indicates that MCPSO-PSH also applies to unimodal functions though it is primarily designed for multimodal functions. For the multimodal functions f16–f19, MCPSO-PSH outperforms NrGA, CMA-ES, PSO, and DE. In terms of average obtained results, MCPSO-PSH is better than NrGA on f12,f13,f15,f16,f17, and f18, and is somehow equivalent to it on f11 and f19. MCPSO-PSH and NrGA outperform other three algorithms on most cases, which demonstrates great promise of nonrevisiting strategy for maintaining population diversity and improving search efficiency synchronously. Obviously, MCPSO-PSH is superior to NrGA which benefits from its unique dynamic multispecies strategy based on K-means clustering.Table 3
Basic statistics off11–f20 CEC2005 benchmark test functions obtained by MCPSO-PSH, NrGA, CMA-ES, PSO, and DE.
Problem
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
f
11
Mean
−4.5000E + 02
−4.5000E + 02
-
3.1857
E
+
02
-
4.3912
E
+
02
-
4.4313
E
+
02
Std
0
0
1.8576
E
+
02
2.3
441
E
+
01
9.7454
E
+
00
f
12
Mean
−4.5000E + 02
3.1554
E
+
03
-
2.2058
E
+
02
1.2360
E
+
03
-
4.4506
E
+
02
Std
0
1.9571
E
+
03
8.3398
E
+
01
3.2960
E
+
02
3.6262
E
-
09
f
13
Mean
−4.5000E + 02
6.9857
E
+
06
2.4164
E
+
05
2.5967
E
+
07
3.5243
E
+
07
Std
0
2.5378
E
+
06
3.0876
E
+
04
4.3459
E
+
06
6.5900
E
+
07
f
14
Mean
3.8043
E
+
02
2.3818E + 02
5.3723
E
+
02
2.7801
E
+
02
3.5976
E
+
03
Std
3.2991
E
+
00
6.8326E + 00
6.3167
E
+
01
8.3767
E
+
00
3.2048
E
+
03
f
15
Mean
2.3174E + 03
1.6341
E
+
04
5.2245
E
+
03
3.3347
E
+
03
3.8727
E
+
03
Std
5.0134E + 02
1.0450
E
+
03
2.1529
E
+
03
5.7811
E
+
02
4.3668
E
+
03
f
16
Mean
−4.5000E + 02
2.1850
E
+
04
-
1.8444
E
+
02
6.5597
E
+
03
-
4.4761
E
+
02
Std
0
3.1742
E
+
03
1.2313
E
+
02
1.3545
E
+
03
3.4378
E
-
12
f
17
Mean
−1.8000E + 02
1.9705
E
+
03
1.8865
E
+
04
4.6964
E
+
03
−1.8000E + 02
Std
0
1.9839
E
+
02
3.9520
E
+
02
8.4681
E
+
02
0
f
18
Mean
−1.1997E + 02
-
1.1911
E
+
02
-
1.1958
E
+
02
-
1.1909
E
+
02
-
1.1995
E
+
02
Std
8.9147E − 03
7.8361
E
-
02
3.7314
E
-
01
4.2087
E
-
02
6.0499
E
-
01
f
19
Mean
−3.3000E + 02
−3.300e + 02
-
1.9906
E
+
02
-
3.1776
E
+
02
-
3.2842
E
+
02
Std
0
0
1.6472
E
+
01
7.5984
E
+
01
8.5989
E
+
01
f
20
Mean
-
4.3294
E
+
02
−4.3933E + 02
-
4.2932
E
+
02
-
3.9982
E
+
02
-
4.0363
E
+
02
Std
3.1803
E
+
01
1.7580E + 01
2.0839
E
+
01
6.1074
E
+
01
3.5741
E
+
01In order to intuitively depict the search effect, the box plots are employed to demonstrate the statistical performance of five different algorithms, as shown in Figure3. Although nonrevisiting strategy is implemented in both NrGA and MCPSO-PSH, the robustness of NrGA is still undesirable. MCPSO-PSH with dynamic multispecies strategy has great enhancement of exploring ability.Figure 3
Box plots of results obtained by all algorithms. Here, 1 to 5 in the horizontal axis are indices of MCPSO-PSH, NrGA, CMA-ES, PSO, and DE, respectively (the box has lines at the lower quartile, median, and upper quartile values. The whiskers extend to the most extreme data points not considered outliers. Outliers (denoted by +) are data with values beyond 100 units of interquartile range).From the above observations, the superiority of MCPSO-PSH over other peer algorithms can be easily presumed.


## 4.3.1. Results for Basic Benchmark Functions

In this experiment, the proposed MCPSO-PSH algorithm is compared with nonrevisiting genetic algorithm (NrGA), Covariance Matrix Adaptation Evolution Strategy (CMA-ES), the canonical particle swarm optimization (PSO) algorithm, and differential evolution (DE) algorithm.For testing 10D benchmark test functionsf1–f10 in 10 runs, the basic statistical results including the mean and standard deviations of the function values are listed in Table 1. From Table 1, it is obvious that the canonical PSO and DE algorithms demonstrate a worse performance than MCPSO-PSH, NrGA, and CMA-ES. For the low-dimensional (10D) unimodal functionsf1–f5, the uniformity of CMA-ES algorithm is relatively good. For functionf1, CMA-ES gets the best results. And, for functionsf2–f5, the performance of CMA-ES is rather competitive. Compared with the canonical PSO and other state-of-the-art algorithms, MCPSO-PSH brings in the dynamical multispecies strategy based onK-means to enhance the ability of local search and fitness tree to overcome the revisit phenomenon, which results in a huge improvement of performance. The MCPSO-PSH and NrGA algorithms with nonrevisit strategy have nearly the same results for low-dimensional problems. Benefitting from the dynamicK-means clustering and the adaptive normalization of population members, the improved multispecies PSO has a slightly better performance than NrGA on functionsf1–f5. For the low-dimensional (10D) multimodal functionsf6–f10, intuitively, the proposed MCPSO-PSH algorithm has rather smaller mean values and standard deviations than other algorithms for functionf6. For functionsf7 andf8, the MCPSO-PSH and NrGA obtain the global optimum with a high efficiency. It is due to the fact that the nonrevisit strategy based on BSP fitness tree effectively enhances the global search especially when a large number of local optima exist in the search space. For functionsf9 andf10, the CMA-ES algorithm also gets the global optimum because the MCPSO-PSH, NrGA, and CMA-ES own their special mechanisms for jumping out of the local optimum. On the whole, regardless of the unimodal functions or the multimodal functions, the performance of the proposed MCPSO-PSH is uniform.Table 1
Performance of all algorithms on 10D test functionsf1–f10.
Problem
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
f
1
Mean
4.3136
E
-
02
3.9407
E
-
01
9.4574E − 05
7.3831
E
+
04
3.0851
E
+
13
Std
7.3787
E
-
07
3.1066
E
-
05
1.4783E − 05
6.3480
E
+
03
1.0345
E
+
12
f
2
Mean
8.0378E − 33
3.7721
E
-
27
8.4991
E
-
18
1.0386
E
-
03
3.9501
E
-
01
Std
1.0378E − 30
7.0550
E
-
30
3.8731
E
-
21
2.7801
E
-
03
4.0404
E
+
00
f
3
Mean
9.4931E − 38
5.7707
E
-
33
3.4959
E
-
14
1.4301
E
+
47
3.5430
E
-
03
Std
3.9093E − 41
7.3900
E
-
37
8.0455
E
-
20
1.3784
E
+
47
5.7791
E
-
04
f
4
Mean
2.9043E − 73
3.8844
E
-
36
8.3995
E
-
18
9.4034
E
+
08
3.1893
E
-
43
Std
4.7898E − 99
1.3894
E
-
45
2.9393
E
-
20
1.0373
E
+
09
7.8440
E
-
53
f
5
Mean
4.9550E − 64
3.9090
E
-
39
4.4139
E
-
43
8.0041
E
-
17
6.0903
E
-
02
Std
4.8414E − 89
3.4568
E
-
40
6.3941
E
-
46
9.4389
E
-
18
5.9031
E
-
02
f
6
Mean
2.2137E − 14
7.3678
E
-
13
3.8931
E
-
13
7.3491
E
+
01
1.3940
E
+
02
Std
3.7890E − 17
2.4654
E
-
14
3.9404
E
-
14
3.4791
E
-
03
4.9448
E
+
01
f
7
Mean
0
0
7.4327
E
-
03
3.9051
E
+
00
2.6461
E
-
2
Std
0
0
4.7801
E
-
03
9.4803
E
-
01
8.4780
E
-
3
f
8
Mean
−9.8452E − 01
−9.8452E − 01
9.0377
E
-
01
2.0773
E
-
01
2.4743
E
+
00
Std
0
0
6.9010
E
-
01
1.3434
E
-
01
9.0370
E
-
1
f
9
Mean
0
0
0
3.8944
E
+
02
1.9640
E
+
03
Std
0
0
0
9.0378
E
+
01
8.9901
E
+
01
f
10
Mean
0
0
0
1.3355
E
+
04
7.9140
E
+
04
Std
0
0
0
8.9440
E
+
03
4.4708
E
+
04For testing 30D benchmark test functionsf1–f10 in 10 runs, the basic statistical results including the mean and standard deviations of the function values are listed in Table 2. With the increase of the dimensions, the convergence rate becomes slow. The results of Table 2 show great disparities in convergence accuracy. For the high-dimensional unimodal problemsf1–f5, the performance of MCPSO-PSH algorithm is more outstanding than the other four algorithms because MCPSO-PSH with dynamic clustering strategy can accelerate the convergence rate especially in the later evolution period. The performance of NrGA algorithm becomes slightly poor although it also has the nonrevisit strategy for improving the search efficiency. Only for functionsf1–f5 does CMA-ES outperform NrGA on convergence rate and accuracy. For solving the high-dimensional multimodal problemsf6–f10, the MCPSO-PSH has a good performance. With the dynamic multispecies strategy based onK-means clustering and nonrevisit strategy based on BSP fitness tree, the MCSPSO-PSH algorithm can well avoid premature convergence to local optimum in the process of solving the high-dimensional multimodal problems. From the basic statistical results, NrGA also owns the most outstanding performance on testing functionsf6 andf10, which means the nonrevisit strategy based on BSP fitness tree plays an important role in guiding search direction. On the whole, the performance of all algorithms tested here is ordered as MCPSO-PSH > NrGA > CMA-ES > PSO > DE.Table 2
Performance of all algorithms on 30D test functionsf1–f10.
Problem
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
f
1
Mean
4.8910E − 01
1.3047
E
+
01
3.1794
E
+
00
3.4903
E
+
12
2.4941
E
+
21
Std
3.3890E − 01
9.0360
E
+
00
2
.
0935
E
+
00
8.0410
E
+
11
9.3892
E
+
19
f
2
Mean
1.6851E − 14
1.8440
E
-
08
5.3882
E
+
02
3.7840
E
+
04
9.9038
E
8
Std
9.0376E − 16
9.5040
E
-
10
3.0318
E
+
02
9.9011
E
+
03
1.3931
E
7
f
3
Mean
5.9990E − 37
6.3890
E
+
01
1.4890
E
-
12
3.8964
E
+
308
7.3280
E
+
223
Std
9.5490E − 41
1.4378
E
+
01
2.1809
E
-
13
1.9430
E
+
310
1.9653
E
+
224
f
4
Mean
1.0099E − 16
8.3701
E
-
08
1.0454
E
-
23
7.9016
E
+
13
3.7890
E
+
08
Std
8.4098E − 18
1.0871
E
-
07
9.5487
E
-
25
1.0346
E
+
14
1.3944
E
+
08
f
5
Mean
1.3933E − 61
2.0031
E
-
11
2.9890
E
-
16
2.3498
E
+
02
5.3903
E
+
05
Std
9.0809E − 75
9.6036
E
-
13
7.3892
E
-
23
7.3954
E
+
01
1.2309
E
+
05
f
6
Mean
1.3343E − 14
1.3408
E
-
12
6.3378
E
-
06
8.4723
E
+
02
2.5044
E
+
02
Std
6.6066E − 15
8.0651
E
-
14
1.6999
E
-
06
3.1704
E
+
02
7.3909
E
+
01
f
7
Mean
0
0
3.0966
E
-
2
2.0330
E
-
01
8.9036
E
-
01
Std
0
0
7.8012
E
-
3
3.3891
E
-
01
6.7009
E
-
01
f
8
Mean
−2.8743E + 01
-
2.7945
E
+
01
1.7640
E
+
01
6.2307
E
-
01
-
3.3883
E
+
00
Std
1.664E + 00
2.0874
E
+
00
7.3370
E
+
00
3.1700
E
-
01
6.0310
E
-
01
f
9
Mean
3.3880E − 12
8.5801
E
+
00
2.1388
E
+
00
4.3880
E
+
04
8.4891
E
+
03
Std
7.3038E − 13
5.8817
E
+
00
2.9863
E
+
00
8.3909
E
+
03
2.7450
E
+
03
f
10
Mean
7.9063
E
-
01
3.6421E − 04
2.4578
E
+
03
2.3841
E
+
04
4.4093
E
+
05
Std
4.1939
E
-
01
1.3940E − 04
2.7309
E
+
03
1.7301
E
+
04
1.0938
E
+
06The convergence curves of the MCPSO-PSH algorithms and other EAs against test functionsf1–f10 in 10D and 30D are shown in Figure 2. It can be observed that, with the increase of the dimensions, the convergence rate of the proposed MCPSO-PSH algorithm becomes faster. Predictably, with a further increase of the dimensions of the test functions, the MCPSO-PSH algorithm will have an overwhelming superiority in solving the high-dimensional problems.Figure 2
Convergence curves of MCPSO-PSH, NrGA, CMA-ES, PSO, and DE against test functionsf1–f10 in low dimension, 10D, and high dimension, 30D, respectively.


## 4.3.2. Results for CEC2005 Benchmark Functions

The statistical results in Table3 show that MCPSO-PSH has significantly outstanding performance for all these cases except f14 and f20. For f14, although all the algorithms provide similar results (except DE), the solutions of MCPSO-PSH are still rather competitive and accessible. Forf20 MCPSO-PSH achieved the second best solution just behind NrGA. The remarkable performance of MCPSO-PSH for f11–f15 indicates that MCPSO-PSH also applies to unimodal functions though it is primarily designed for multimodal functions. For the multimodal functions f16–f19, MCPSO-PSH outperforms NrGA, CMA-ES, PSO, and DE. In terms of average obtained results, MCPSO-PSH is better than NrGA on f12,f13,f15,f16,f17, and f18, and is somehow equivalent to it on f11 and f19. MCPSO-PSH and NrGA outperform other three algorithms on most cases, which demonstrates great promise of nonrevisiting strategy for maintaining population diversity and improving search efficiency synchronously. Obviously, MCPSO-PSH is superior to NrGA which benefits from its unique dynamic multispecies strategy based on K-means clustering.Table 3
Basic statistics off11–f20 CEC2005 benchmark test functions obtained by MCPSO-PSH, NrGA, CMA-ES, PSO, and DE.
Problem
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
f
11
Mean
−4.5000E + 02
−4.5000E + 02
-
3.1857
E
+
02
-
4.3912
E
+
02
-
4.4313
E
+
02
Std
0
0
1.8576
E
+
02
2.3
441
E
+
01
9.7454
E
+
00
f
12
Mean
−4.5000E + 02
3.1554
E
+
03
-
2.2058
E
+
02
1.2360
E
+
03
-
4.4506
E
+
02
Std
0
1.9571
E
+
03
8.3398
E
+
01
3.2960
E
+
02
3.6262
E
-
09
f
13
Mean
−4.5000E + 02
6.9857
E
+
06
2.4164
E
+
05
2.5967
E
+
07
3.5243
E
+
07
Std
0
2.5378
E
+
06
3.0876
E
+
04
4.3459
E
+
06
6.5900
E
+
07
f
14
Mean
3.8043
E
+
02
2.3818E + 02
5.3723
E
+
02
2.7801
E
+
02
3.5976
E
+
03
Std
3.2991
E
+
00
6.8326E + 00
6.3167
E
+
01
8.3767
E
+
00
3.2048
E
+
03
f
15
Mean
2.3174E + 03
1.6341
E
+
04
5.2245
E
+
03
3.3347
E
+
03
3.8727
E
+
03
Std
5.0134E + 02
1.0450
E
+
03
2.1529
E
+
03
5.7811
E
+
02
4.3668
E
+
03
f
16
Mean
−4.5000E + 02
2.1850
E
+
04
-
1.8444
E
+
02
6.5597
E
+
03
-
4.4761
E
+
02
Std
0
3.1742
E
+
03
1.2313
E
+
02
1.3545
E
+
03
3.4378
E
-
12
f
17
Mean
−1.8000E + 02
1.9705
E
+
03
1.8865
E
+
04
4.6964
E
+
03
−1.8000E + 02
Std
0
1.9839
E
+
02
3.9520
E
+
02
8.4681
E
+
02
0
f
18
Mean
−1.1997E + 02
-
1.1911
E
+
02
-
1.1958
E
+
02
-
1.1909
E
+
02
-
1.1995
E
+
02
Std
8.9147E − 03
7.8361
E
-
02
3.7314
E
-
01
4.2087
E
-
02
6.0499
E
-
01
f
19
Mean
−3.3000E + 02
−3.300e + 02
-
1.9906
E
+
02
-
3.1776
E
+
02
-
3.2842
E
+
02
Std
0
0
1.6472
E
+
01
7.5984
E
+
01
8.5989
E
+
01
f
20
Mean
-
4.3294
E
+
02
−4.3933E + 02
-
4.2932
E
+
02
-
3.9982
E
+
02
-
4.0363
E
+
02
Std
3.1803
E
+
01
1.7580E + 01
2.0839
E
+
01
6.1074
E
+
01
3.5741
E
+
01In order to intuitively depict the search effect, the box plots are employed to demonstrate the statistical performance of five different algorithms, as shown in Figure3. Although nonrevisiting strategy is implemented in both NrGA and MCPSO-PSH, the robustness of NrGA is still undesirable. MCPSO-PSH with dynamic multispecies strategy has great enhancement of exploring ability.Figure 3
Box plots of results obtained by all algorithms. Here, 1 to 5 in the horizontal axis are indices of MCPSO-PSH, NrGA, CMA-ES, PSO, and DE, respectively (the box has lines at the lower quartile, median, and upper quartile values. The whiskers extend to the most extreme data points not considered outliers. Outliers (denoted by +) are data with values beyond 100 units of interquartile range).From the above observations, the superiority of MCPSO-PSH over other peer algorithms can be easily presumed.


## 5. Multilevel Threshold for Image Segmentation

### 5.1. Otsu Method

The Otsu multithreshold measure [31] proposed by Otsu has been popularly employed in determining whether the optimal threshold method can provide image segmentation with satisfactory results. The traditional Otsu method can be described as follows.Let the grey levels of a given image withN pixels range over [0, L-1] and h(i) denote the number of the ith grey level pixels. P(i) is the probability of i.(31)N=∑i=0L-1hi,Pi=hiN,for0≤i≤L-1.Otsu defined the between-class variance as the sum of sigma functions of each class by (22)(32)ft=δ0+δ1,δ0=w0·u0-uT2,δ1=w1u1-uT2,δ1=w1·u1-uT2,where uT is the mean intensity of the original image. In bilevel thresholding, the mean level of each class, (ui), can be calculated as follows: (33)u0=∑i=0i=t-1i×Piw0,u1=∑i=ti=Li×Piw0,where(34)w0=∑i=0i=t-1Pi,w1=∑i=ti=LPi.The optimal threshold can be derived by maximizing the between-class variance function(35)t∗=arg⁡max⁡ft.Bilevel thresholding based on the between-class variance can be extended to multilevel thresholding. The extended between-class variance function is calculated as follows:(36)ft=∑i=0Mδi,where M is the number of thresholds. And (37)δ0=w0u0-uT2,δj=wjuj-uT2,δM=wMuM-uT2,u0=∑i=0ti-1i×Piwi,uj=∑i=tjtj+1-1i×Piwi,uM=∑i=tML-1i×Piwi.The optimal thresholds are found by maximizing the between-class variance function:(38)t1,t2,…,tM=arg⁡max⁡ft1,t2,…,tM.Equation (38) is used as the objective function which is to be optimized (maximized). A close look into this equation will show that it is very similar to the expression for uniformity measure.


### 5.2. Experiment Setup

For evaluating the performance of the proposed algorithm (MCPSO-PSH), we implemented it on a wide variety of images. These images include the popular tested images used in previous studies [32, 33], named avion.ppm, lena.ppm, and hunter.pgm (available in [34]) and 41033.jpg, 62096.jpg, and 145086.jpg provided by the Berkeley Segmentation Database (BSDS) (available in [35]). The sizes of avion.ppm, lena.ppm, and hunter.pgm are 512∗512 and the sizes of 41033.jpg, 62096.jpg, and 145086.jpg are 321∗481. Each image has a unique grey-level histogram. Most of the images are difficult to segment due to multimodality of the histograms.The performance of the proposed algorithm is demonstrated by comparing with those of other popular algorithms developed in the literature so far, for example, canonical particle swarm optimization (PSO), differential evolution (DE), Covariance Matrix Adaptation Evolution Strategy (CMA-ES), and nonrevisiting genetic algorithm (NrGA). A comparison between the proposed algorithm and other methods is evaluated based on Otsu, which means (38) is regarded as fitness function to evaluate all involved algorithms. The parameters settings are identical to those introduced in Section 4.1.The numbers of thresholdsM investigated in the experiments were 2, 3, 4, 5, 7, and 9, while all the experiments were repeated 30 times for each image. The population size is 100 and the maximum number of FEs is 2000. The numbers of thresholds M-1=2,3,4 can be realized as low-dimensional problems, which also can be solved through the traditional search method although suffering from the long CPU time. The numbers of thresholds M-1=5,7,9 are regarded as medium levels for image segmentation, which means the high-dimensional problem for image segmentation cannot be solved by the traditional search method.In all experiments, to ensure that the initial values of each random algorithm are equal, we used the MATLAB commandrand⁡X,Y∗max⁡L-min⁡L+min⁡L. max⁡L and min⁡L are the maximum and minimum greys of the tested images. Figure 4 presents six original images and their histograms.Figure 4
Test images and their histograms.


### 5.3. Experimental Results of Multilevel Threshold

Case  I: Multilevel Threshold Results with
M
-
1
=
2,3
,
4. We employ (38) as the fitness function to provide a comparison of performances. Since the classical Otsu method is an exhaustive searching algorithm, we regarded its results as the “standard.” Table 4 shows the fitness function values (with M-1=2,3,4) attained by Otsu methods. In Table 4, the mean computation time and corresponding optimal thresholds are listed. For M-1>4, owing to the unbearable consumption of CPU time, we do not list the correlative values in our experiment. These “standard” results will be used to compare with the results attained by optimization methods.Table 4
Objective values and thresholds by the Otsu method.
Image
M
=
2
M
=
3
M
=
4
Objective values
Optimal thresholds
Objective values
Optimal thresholds
Objective values
Optimal thresholds
Avion
3.4413
E
4
113, 173
3.4489
E
4
93, 145, 191
3.4535
E
4
84, 129, 172, 203
Lena
1.7597
E
4
93, 151
1.7764
E
4
81, 127, 171
1.7827
E
4
75, 114, 145, 180
Hunter
9.2768
E
3
51, 116
9.4261
E
3
36, 86, 135
9.4822
E
3
30, 72, 111, 146
41033.jpg
1.2753
E
4
64, 131
1.2885
E
4
57, 107, 165
1.2945
E
4
49, 79, 122, 171
62096.jpg
2.7469
E
4
102, 196
2.7585
E
4
86, 156, 218
2.7631
E
4
71, 121, 176, 225
145086.jpg
1.5754
E
4
89, 175
1.5982
E
4
69, 122, 186
1.6085
E
4
41, 85, 130, 189
Mean CPU time (s)
1.4427
61.421
2644.543It is well known that the EA-based Otsu method for multilevel thresholding segmentation can only accelerate the computation velocity. The attained results will just get incredibly close to the optimal results but fail to improve the optimal results. Using the proposed MCPSO-PSH and other evolution algorithms, the mean fitness function values and their standard deviations obtained by the test algorithms (withM-1=2,3,4) are shown in Table 5. The bigger mean values and smaller standard deviation listed in Table 5 mean that the algorithm gets the better results.Table 5
Objective value and standard deviation by the compared EA-based methods on Otsu algorithm.
Image
M
Objective values (standard deviation)
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
Avion
2
3.4413
E
4
3.4413
E
4
3.4054
E
4
3.4411
E
4
3.4054
E
4
0.1390
7.66953
E
-
12
6.633
E
-
2
2.9091
7.6695
E
-
12
3
3.4489
E
4
3.4489
E
4
3.4129
E
4
3.4485
E
4
3.4130
E
4
0.2708
1.9659
E
-
2
0.2735
4.4573
0.1892
4
3.4534
E
4
3.4535
E
4
3.4172
E
4
3.4507
E
4
3.4174
E
4
0.5044
0
1.5801
28.4693
0.4677
Lena
2
1.7597
E
4
1.7597
E
4
1.7348
E
4
1.7590
E
4
1.7348
E
4
3.8348
E
-
12
3.8347
E
-
12
0.1153
10.9993
3.8348
E
-
12
3
1.7764
E
4
1.7764
E
4
1.7513
E
4
1.7723
E
4
1.7515
E
4
6.3873
E
-
3
3.83477
E
-
1
1.1303
65.6575
5.6065
E
-
2
4
1.7827
E
4
1.7828
E
4
1.7573
E
4
1.7797
E
4
1.7578
E
4
0.4322
1.19275
E
-
4
1.9239
31.2048
0.7919
Hunter
2
9.2768
E
3
9.2769
E
3
9.1196
E
3
9.1641
E
3
9.1200
E
3
0
0
0.4918
1.4411
E
2
0
3
9.4261
E
3
9.4261
E
3
9.2670
E
3
9.3724
E
3
9.2686
E
3
1.3909
E
-
2
1.9174
E
-
12
1.2988
65.2063
0.6674
4
9.4821
E
3
9.4821
E
3
9.3215
E
3
9.4147
E
3
9.3214
E
3
0.1261
0.1262
1.7285
47.8078
3.2285
41033
2
1.2753
E
4
1.2753
E
4
1.2753
E
4
1.2746
E
4
1.2753
E
4
0
4.0595
E
-
3
0.1169
14.7680
0
3
1.2885
E
4
1.2885
E
4
1.2884
E
4
1.2866
E
4
1.2885
E
4
0
0
1.0243
42.8287
0.1079
4
1.2945
E
4
1.2945
E
5
1.2940
E
4
1.2911
E
4
1.2944
E
4
1.398
E
-
3
1.3245
E
-
3
2.0414
31.5746
2.2078
62096
2
2.7469
E
4
2.7469
E
4
2.7469
E
4
2.7457
E
4
2.7469
E
4
0
3.8348
E
-
12
0.1693
15.2460
0
3
2.7585
E
4
2.7585
E
4
2.7584
E
4
2.7568
E
4
2.7585
E
4
0.5133
0
0.8584
35.4894
0.1933
4
2.7631
E
4
2.7631
E
4
2.7629
E
4
2.7596
E
4
2.7629
E
4
0.2255
15222
E
-
2
1.1273
33.2765
1.0023
145086
2
1.5754
E
4
1.5754
E
4
1.5753
E
4
1.5745
E
4
1.5754
E
4
0
3.8348
E
-
12
0.5892
24.8969
3.8348
E
-
12
3
1.5982
E
4
1.5982
E
4
1.5980
E
4
1.5955
E
4
1.5982
E
4
2.3331
E
-
2
1.9174
E
-
12
1.3726
31.9229
0.1242
4
1.6085
E
4
1.6085
E
4
1.6078
E
4
1.6054
E
4
1.6083
E
4
0
1.9174
E
-
12
3.2000
19.4620
1.8526According to the mean CPU times shown in Table6, there is no obvious difference between the EA-based Otsu methods. All the EA-based Otsu methods are suitably efficient and practical in terms of time complexity for high-dimensional image segmentation problems. Therefore, in comparison with the traditional Otsu algorithm, all the EA-based Otsu methods effectively shorten the computation time.Table 6
The mean CPU time of the compared EA-based methods on Otsu algorithm.
Dim
Algorithms
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
2 dim
0.3398
0.3088
0.2802
0.2553
0.2478
3 dim
0.3588
0.3193
0.2206
0.2497
0.3170
4 dim
0.3221
0.3897
0.2264
0.2520
0.3418
5 dim
0.3087
0.3903
0.1991
0.2355
0.3481
7 dim
0.3021
0.3818
0.2068
0.2311
0.3583
9 dim
0.3406
0.3949
0.1932
0.2308
0.3390Among the EA-based Otsu methods tested listed in Table5, MCPSO-PSH obtains results equal to or close to that obtained in the traditional methods when M-1=2,3,4. Moreover, the standard deviations obtained by MCPSO-PSH are considerably small among the results obtained by the EA-based Otsu methods. This is attributed to the fact that the MCPSO-PSH has a good balance between the global and local search abilities. Therefore, the MCPSO-PSH-based Otsu method is suitable for the low-dimensional segmentation problems (M-1=2,3,4).Case II: Multilevel Threshold Results with
M
-
1
=
5,7
,
9. To further investigate the search performance of each population-based method on high-dimensional segmentation problems, we conduct these algorithms on image segment with M-1=5,7,9. Due to the unacceptable computation time for the traditional Otsu method, it does not need to be assessed in this test. The average fitness value and standard deviation obtained by each of the EA-based Otsu methods are listed in Table 7, where the correlative results with the larger values and smaller standard deviations indicate the better achievements.Table 7
Objective value and standard deviation by the compared EA-based methods on Otsu algorithm.
Image
M
Objective values (standard deviation)
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
Avion
5
3.4561
E
4
3.4561
E
4
3.4197
E
4
3.4547
E
4
3.4198
E
4
2.2681
E
-
2
6.6013
E
-
2
2.2062
10.5393
1.2606
7
3.4585
E
4
3.4585
E
4
3.4222
E
4
3.4576
E
4
3.4218
E
4
0.1227
0.4133
1.3292
7.2486
3.1975
9
3.4598
E
4
3.4597
E
4
3.4233
E
4
3.4580
E
4
3.4230
E
4
0.2757
1.0760
1.1419
11.0739
2.1899
Lena
5
1.7853
E
4
1.7853
E
4
1.7598
E
4
1.7829
E
4
1.7603
E
4
1.1762
E
-
2
0.5437
1.3751
26.0899
1.9589
7
1.7885
E
4
1.7885
E
4
1.7626
E
4
1.7857
E
4
1.7631
E
4
0.5759
0.1317
1.9112
23.0333
4.2468
9
1.7898
E
4
1.7898
E
4
1.7641
E
4
1.7882
E
4
1.7642
E
4
0.7920
0.1250
2.3463
6.4354
3.2478
Hunter
5
9.5208
E
3
9.3860
E
3
9.3567
E
3
9.4591
E
3
9.3581
E
3
1.3509
0.9022
3.0022
46.7651
3.9539
7
9.5559
E
3
9.4292
E
3
9.3892
E
3
9.5000
E
3
9.3895
E
3
0.132607
1.6085
3.6985
61.6729
3.6627
9
9.5736
E
3
9.4400
E
3
9.4060
E
3
9.5248
E
3
9.4041
E
3
0.661602
0.8640
2.6400
24.2906
3.0252
41033
5
1.2974
E
4
1.2974
E
4
1.2967
E
4
1.2938
E
4
1.2970
E
4
1.7191
E
-
2
0.1151
2.5816
19.4046
2.1793
7
1.3006
E
4
1.3006
E
4
1.2995
E
4
1.2968
E
4
1.2995
E
4
0.220026
0.1466
1.9642
12.1857
3.9443
9
13021.1
1.3020
E
4
1.3013
E
4
1.3004
E
4
1.3012
E
4
0.883857
1.2077
2.4069
5.8058
2.4286
62096
5
2.7658
E
4
2.7658
E
4
2.7654
E
4
2.7650
E
4
2.7654
E
4
3.5718
E
-
2
0.2520
1.0799
5.7383
2.2228
7
2.7682
E
4
2.7682
E
4
2.7677
E
4
2.7669
E
4
2.7677
E
4
0.1319
6.7590
E
-
2
1.1261
5.7592
2.1289
9
2.7695
E
4
2.7695
E
4
2.7689
E
4
2.7682
E
4
2.7689
E
4
0.545201
0.4215
1.6069
9.1974
1.1693
145086
5
1.6143
E
4
1.6143
E
4
1.6134
E
4
1.6108
E
4
1.6139
E
4
2.7021
E
-
2
3.0505
E
-
2
4.0455
43.2172
3.4949
7
1.6195
E
4
1.6195
E
4
1.6184
E
4
1.6160
E
4
1.6184
E
4
0.1298
5.3333
E
-
2
3.5967
23.9921
2.7644
9
1.6217
E
4
1.6217
E
4
1.6205
E
4
1.6192
E
4
1.6205
E
4
0.4759
0.2663
2.7680
15.9752
2.4381From Table7, it is obvious that, with the growth in dimension, there is statistically significant gap among experiments using these EA-based Otsu methods, in terms of both efficiency (fitness values) and stability (standard deviation). Depending on the nonrevisit strategy based on Binary Space Partitioning fitness tree and dynamic multispecies strategy based on K-means clustering, MCPSO-PSH owns the best performance and stability on a high dimension, which is more efficient than the canonical PSO and other classical EAs for image segmentation with Otsu method. Moreover, NrGA, which used nonrevisit strategy, also gets more acceptable results than the other algorithms, except MCPSO-PSH. With the increase of dimensions, the fitness values and standard deviation demonstrate that the performance of CMA-ES-based Otsu method has greatly improved. The results shown in Table 7 prove that the MCPSO-PSH-based Otsu method is more suitable for resolving multilevel image segmentation problems.As mentioned above, the EA-based Otsu methods adopt different strategies to conduct global and local searches. An elaborate optimizer should have a better balance between exploration and exploitation search abilities. The above tests show that the MCPSO-PSH-based Otsu method outperforms the other EA-based Otsu method for multilevel image segmentation, especially on high-dimensional thresholding (withM-1=5,7,9). Moreover, compared with the traditional Otsu method, this method shortens the computational time.


## 5.1. Otsu Method

The Otsu multithreshold measure [31] proposed by Otsu has been popularly employed in determining whether the optimal threshold method can provide image segmentation with satisfactory results. The traditional Otsu method can be described as follows.Let the grey levels of a given image withN pixels range over [0, L-1] and h(i) denote the number of the ith grey level pixels. P(i) is the probability of i.(31)N=∑i=0L-1hi,Pi=hiN,for0≤i≤L-1.Otsu defined the between-class variance as the sum of sigma functions of each class by (22)(32)ft=δ0+δ1,δ0=w0·u0-uT2,δ1=w1u1-uT2,δ1=w1·u1-uT2,where uT is the mean intensity of the original image. In bilevel thresholding, the mean level of each class, (ui), can be calculated as follows: (33)u0=∑i=0i=t-1i×Piw0,u1=∑i=ti=Li×Piw0,where(34)w0=∑i=0i=t-1Pi,w1=∑i=ti=LPi.The optimal threshold can be derived by maximizing the between-class variance function(35)t∗=arg⁡max⁡ft.Bilevel thresholding based on the between-class variance can be extended to multilevel thresholding. The extended between-class variance function is calculated as follows:(36)ft=∑i=0Mδi,where M is the number of thresholds. And (37)δ0=w0u0-uT2,δj=wjuj-uT2,δM=wMuM-uT2,u0=∑i=0ti-1i×Piwi,uj=∑i=tjtj+1-1i×Piwi,uM=∑i=tML-1i×Piwi.The optimal thresholds are found by maximizing the between-class variance function:(38)t1,t2,…,tM=arg⁡max⁡ft1,t2,…,tM.Equation (38) is used as the objective function which is to be optimized (maximized). A close look into this equation will show that it is very similar to the expression for uniformity measure.


## 5.2. Experiment Setup

For evaluating the performance of the proposed algorithm (MCPSO-PSH), we implemented it on a wide variety of images. These images include the popular tested images used in previous studies [32, 33], named avion.ppm, lena.ppm, and hunter.pgm (available in [34]) and 41033.jpg, 62096.jpg, and 145086.jpg provided by the Berkeley Segmentation Database (BSDS) (available in [35]). The sizes of avion.ppm, lena.ppm, and hunter.pgm are 512∗512 and the sizes of 41033.jpg, 62096.jpg, and 145086.jpg are 321∗481. Each image has a unique grey-level histogram. Most of the images are difficult to segment due to multimodality of the histograms.The performance of the proposed algorithm is demonstrated by comparing with those of other popular algorithms developed in the literature so far, for example, canonical particle swarm optimization (PSO), differential evolution (DE), Covariance Matrix Adaptation Evolution Strategy (CMA-ES), and nonrevisiting genetic algorithm (NrGA). A comparison between the proposed algorithm and other methods is evaluated based on Otsu, which means (38) is regarded as fitness function to evaluate all involved algorithms. The parameters settings are identical to those introduced in Section 4.1.The numbers of thresholdsM investigated in the experiments were 2, 3, 4, 5, 7, and 9, while all the experiments were repeated 30 times for each image. The population size is 100 and the maximum number of FEs is 2000. The numbers of thresholds M-1=2,3,4 can be realized as low-dimensional problems, which also can be solved through the traditional search method although suffering from the long CPU time. The numbers of thresholds M-1=5,7,9 are regarded as medium levels for image segmentation, which means the high-dimensional problem for image segmentation cannot be solved by the traditional search method.In all experiments, to ensure that the initial values of each random algorithm are equal, we used the MATLAB commandrand⁡X,Y∗max⁡L-min⁡L+min⁡L. max⁡L and min⁡L are the maximum and minimum greys of the tested images. Figure 4 presents six original images and their histograms.Figure 4
Test images and their histograms.


## 5.3. Experimental Results of Multilevel Threshold

Case  I: Multilevel Threshold Results with
M
-
1
=
2,3
,
4. We employ (38) as the fitness function to provide a comparison of performances. Since the classical Otsu method is an exhaustive searching algorithm, we regarded its results as the “standard.” Table 4 shows the fitness function values (with M-1=2,3,4) attained by Otsu methods. In Table 4, the mean computation time and corresponding optimal thresholds are listed. For M-1>4, owing to the unbearable consumption of CPU time, we do not list the correlative values in our experiment. These “standard” results will be used to compare with the results attained by optimization methods.Table 4
Objective values and thresholds by the Otsu method.
Image
M
=
2
M
=
3
M
=
4
Objective values
Optimal thresholds
Objective values
Optimal thresholds
Objective values
Optimal thresholds
Avion
3.4413
E
4
113, 173
3.4489
E
4
93, 145, 191
3.4535
E
4
84, 129, 172, 203
Lena
1.7597
E
4
93, 151
1.7764
E
4
81, 127, 171
1.7827
E
4
75, 114, 145, 180
Hunter
9.2768
E
3
51, 116
9.4261
E
3
36, 86, 135
9.4822
E
3
30, 72, 111, 146
41033.jpg
1.2753
E
4
64, 131
1.2885
E
4
57, 107, 165
1.2945
E
4
49, 79, 122, 171
62096.jpg
2.7469
E
4
102, 196
2.7585
E
4
86, 156, 218
2.7631
E
4
71, 121, 176, 225
145086.jpg
1.5754
E
4
89, 175
1.5982
E
4
69, 122, 186
1.6085
E
4
41, 85, 130, 189
Mean CPU time (s)
1.4427
61.421
2644.543It is well known that the EA-based Otsu method for multilevel thresholding segmentation can only accelerate the computation velocity. The attained results will just get incredibly close to the optimal results but fail to improve the optimal results. Using the proposed MCPSO-PSH and other evolution algorithms, the mean fitness function values and their standard deviations obtained by the test algorithms (withM-1=2,3,4) are shown in Table 5. The bigger mean values and smaller standard deviation listed in Table 5 mean that the algorithm gets the better results.Table 5
Objective value and standard deviation by the compared EA-based methods on Otsu algorithm.
Image
M
Objective values (standard deviation)
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
Avion
2
3.4413
E
4
3.4413
E
4
3.4054
E
4
3.4411
E
4
3.4054
E
4
0.1390
7.66953
E
-
12
6.633
E
-
2
2.9091
7.6695
E
-
12
3
3.4489
E
4
3.4489
E
4
3.4129
E
4
3.4485
E
4
3.4130
E
4
0.2708
1.9659
E
-
2
0.2735
4.4573
0.1892
4
3.4534
E
4
3.4535
E
4
3.4172
E
4
3.4507
E
4
3.4174
E
4
0.5044
0
1.5801
28.4693
0.4677
Lena
2
1.7597
E
4
1.7597
E
4
1.7348
E
4
1.7590
E
4
1.7348
E
4
3.8348
E
-
12
3.8347
E
-
12
0.1153
10.9993
3.8348
E
-
12
3
1.7764
E
4
1.7764
E
4
1.7513
E
4
1.7723
E
4
1.7515
E
4
6.3873
E
-
3
3.83477
E
-
1
1.1303
65.6575
5.6065
E
-
2
4
1.7827
E
4
1.7828
E
4
1.7573
E
4
1.7797
E
4
1.7578
E
4
0.4322
1.19275
E
-
4
1.9239
31.2048
0.7919
Hunter
2
9.2768
E
3
9.2769
E
3
9.1196
E
3
9.1641
E
3
9.1200
E
3
0
0
0.4918
1.4411
E
2
0
3
9.4261
E
3
9.4261
E
3
9.2670
E
3
9.3724
E
3
9.2686
E
3
1.3909
E
-
2
1.9174
E
-
12
1.2988
65.2063
0.6674
4
9.4821
E
3
9.4821
E
3
9.3215
E
3
9.4147
E
3
9.3214
E
3
0.1261
0.1262
1.7285
47.8078
3.2285
41033
2
1.2753
E
4
1.2753
E
4
1.2753
E
4
1.2746
E
4
1.2753
E
4
0
4.0595
E
-
3
0.1169
14.7680
0
3
1.2885
E
4
1.2885
E
4
1.2884
E
4
1.2866
E
4
1.2885
E
4
0
0
1.0243
42.8287
0.1079
4
1.2945
E
4
1.2945
E
5
1.2940
E
4
1.2911
E
4
1.2944
E
4
1.398
E
-
3
1.3245
E
-
3
2.0414
31.5746
2.2078
62096
2
2.7469
E
4
2.7469
E
4
2.7469
E
4
2.7457
E
4
2.7469
E
4
0
3.8348
E
-
12
0.1693
15.2460
0
3
2.7585
E
4
2.7585
E
4
2.7584
E
4
2.7568
E
4
2.7585
E
4
0.5133
0
0.8584
35.4894
0.1933
4
2.7631
E
4
2.7631
E
4
2.7629
E
4
2.7596
E
4
2.7629
E
4
0.2255
15222
E
-
2
1.1273
33.2765
1.0023
145086
2
1.5754
E
4
1.5754
E
4
1.5753
E
4
1.5745
E
4
1.5754
E
4
0
3.8348
E
-
12
0.5892
24.8969
3.8348
E
-
12
3
1.5982
E
4
1.5982
E
4
1.5980
E
4
1.5955
E
4
1.5982
E
4
2.3331
E
-
2
1.9174
E
-
12
1.3726
31.9229
0.1242
4
1.6085
E
4
1.6085
E
4
1.6078
E
4
1.6054
E
4
1.6083
E
4
0
1.9174
E
-
12
3.2000
19.4620
1.8526According to the mean CPU times shown in Table6, there is no obvious difference between the EA-based Otsu methods. All the EA-based Otsu methods are suitably efficient and practical in terms of time complexity for high-dimensional image segmentation problems. Therefore, in comparison with the traditional Otsu algorithm, all the EA-based Otsu methods effectively shorten the computation time.Table 6
The mean CPU time of the compared EA-based methods on Otsu algorithm.
Dim
Algorithms
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
2 dim
0.3398
0.3088
0.2802
0.2553
0.2478
3 dim
0.3588
0.3193
0.2206
0.2497
0.3170
4 dim
0.3221
0.3897
0.2264
0.2520
0.3418
5 dim
0.3087
0.3903
0.1991
0.2355
0.3481
7 dim
0.3021
0.3818
0.2068
0.2311
0.3583
9 dim
0.3406
0.3949
0.1932
0.2308
0.3390Among the EA-based Otsu methods tested listed in Table5, MCPSO-PSH obtains results equal to or close to that obtained in the traditional methods when M-1=2,3,4. Moreover, the standard deviations obtained by MCPSO-PSH are considerably small among the results obtained by the EA-based Otsu methods. This is attributed to the fact that the MCPSO-PSH has a good balance between the global and local search abilities. Therefore, the MCPSO-PSH-based Otsu method is suitable for the low-dimensional segmentation problems (M-1=2,3,4).Case II: Multilevel Threshold Results with
M
-
1
=
5,7
,
9. To further investigate the search performance of each population-based method on high-dimensional segmentation problems, we conduct these algorithms on image segment with M-1=5,7,9. Due to the unacceptable computation time for the traditional Otsu method, it does not need to be assessed in this test. The average fitness value and standard deviation obtained by each of the EA-based Otsu methods are listed in Table 7, where the correlative results with the larger values and smaller standard deviations indicate the better achievements.Table 7
Objective value and standard deviation by the compared EA-based methods on Otsu algorithm.
Image
M
Objective values (standard deviation)
MCPSO-PSH
NrGA
CMA-ES
PSO
DE
Avion
5
3.4561
E
4
3.4561
E
4
3.4197
E
4
3.4547
E
4
3.4198
E
4
2.2681
E
-
2
6.6013
E
-
2
2.2062
10.5393
1.2606
7
3.4585
E
4
3.4585
E
4
3.4222
E
4
3.4576
E
4
3.4218
E
4
0.1227
0.4133
1.3292
7.2486
3.1975
9
3.4598
E
4
3.4597
E
4
3.4233
E
4
3.4580
E
4
3.4230
E
4
0.2757
1.0760
1.1419
11.0739
2.1899
Lena
5
1.7853
E
4
1.7853
E
4
1.7598
E
4
1.7829
E
4
1.7603
E
4
1.1762
E
-
2
0.5437
1.3751
26.0899
1.9589
7
1.7885
E
4
1.7885
E
4
1.7626
E
4
1.7857
E
4
1.7631
E
4
0.5759
0.1317
1.9112
23.0333
4.2468
9
1.7898
E
4
1.7898
E
4
1.7641
E
4
1.7882
E
4
1.7642
E
4
0.7920
0.1250
2.3463
6.4354
3.2478
Hunter
5
9.5208
E
3
9.3860
E
3
9.3567
E
3
9.4591
E
3
9.3581
E
3
1.3509
0.9022
3.0022
46.7651
3.9539
7
9.5559
E
3
9.4292
E
3
9.3892
E
3
9.5000
E
3
9.3895
E
3
0.132607
1.6085
3.6985
61.6729
3.6627
9
9.5736
E
3
9.4400
E
3
9.4060
E
3
9.5248
E
3
9.4041
E
3
0.661602
0.8640
2.6400
24.2906
3.0252
41033
5
1.2974
E
4
1.2974
E
4
1.2967
E
4
1.2938
E
4
1.2970
E
4
1.7191
E
-
2
0.1151
2.5816
19.4046
2.1793
7
1.3006
E
4
1.3006
E
4
1.2995
E
4
1.2968
E
4
1.2995
E
4
0.220026
0.1466
1.9642
12.1857
3.9443
9
13021.1
1.3020
E
4
1.3013
E
4
1.3004
E
4
1.3012
E
4
0.883857
1.2077
2.4069
5.8058
2.4286
62096
5
2.7658
E
4
2.7658
E
4
2.7654
E
4
2.7650
E
4
2.7654
E
4
3.5718
E
-
2
0.2520
1.0799
5.7383
2.2228
7
2.7682
E
4
2.7682
E
4
2.7677
E
4
2.7669
E
4
2.7677
E
4
0.1319
6.7590
E
-
2
1.1261
5.7592
2.1289
9
2.7695
E
4
2.7695
E
4
2.7689
E
4
2.7682
E
4
2.7689
E
4
0.545201
0.4215
1.6069
9.1974
1.1693
145086
5
1.6143
E
4
1.6143
E
4
1.6134
E
4
1.6108
E
4
1.6139
E
4
2.7021
E
-
2
3.0505
E
-
2
4.0455
43.2172
3.4949
7
1.6195
E
4
1.6195
E
4
1.6184
E
4
1.6160
E
4
1.6184
E
4
0.1298
5.3333
E
-
2
3.5967
23.9921
2.7644
9
1.6217
E
4
1.6217
E
4
1.6205
E
4
1.6192
E
4
1.6205
E
4
0.4759
0.2663
2.7680
15.9752
2.4381From Table7, it is obvious that, with the growth in dimension, there is statistically significant gap among experiments using these EA-based Otsu methods, in terms of both efficiency (fitness values) and stability (standard deviation). Depending on the nonrevisit strategy based on Binary Space Partitioning fitness tree and dynamic multispecies strategy based on K-means clustering, MCPSO-PSH owns the best performance and stability on a high dimension, which is more efficient than the canonical PSO and other classical EAs for image segmentation with Otsu method. Moreover, NrGA, which used nonrevisit strategy, also gets more acceptable results than the other algorithms, except MCPSO-PSH. With the increase of dimensions, the fitness values and standard deviation demonstrate that the performance of CMA-ES-based Otsu method has greatly improved. The results shown in Table 7 prove that the MCPSO-PSH-based Otsu method is more suitable for resolving multilevel image segmentation problems.As mentioned above, the EA-based Otsu methods adopt different strategies to conduct global and local searches. An elaborate optimizer should have a better balance between exploration and exploitation search abilities. The above tests show that the MCPSO-PSH-based Otsu method outperforms the other EA-based Otsu method for multilevel image segmentation, especially on high-dimensional thresholding (withM-1=5,7,9). Moreover, compared with the traditional Otsu method, this method shortens the computational time.


## 6. Conclusions

In this paper, we have presented a new coevolution particle swarm optimization algorithm with dynamic multispecies strategy based onK-means clustering and nonrevisit strategy based on Binary Space Partitioning fitness tree (called MCPSO-PSH). With the help of storing and analyzing the positions and the fitness values of all evaluated solutions in Binary Space Partitioning tree archive, MCPSO-PSH can prevent its individuals from revisiting the same researched regions and improve the search efficiency. The proposed algorithm adopts K-means clustering method to dynamically partition the population into many clusters and the number of clusters is changing frequently to implement information exchange among the different clusters. The proposed algorithm keeps a good balance between exploration and exploitation search abilities. As a result, the proposed algorithm showed promising results among the compared algorithms. Finally, we applied the proposed algorithm in the multilevel image segmentation problem. The proposed algorithm gets favorable results compared to the other compared methods. The use of the MCPSO-PSH algorithm in other real-world problems merits further research in our future work.While this research is promising, there is a large margin of improvement for future research:(1)
The number of clustersC and the repartitioned cycle T are preset in this work. However, according to our preliminary experiments, different values of C and T will impact the results of several functions. C and T determine the range and rate of information exchange among species. For getting more suitable results of most test cases, we preset these two values. Therefore, how to adaptively adjust the values of C and T with regard to various functions is to be explored in the next stage.(2)
A nonrevisit strategy based on Binary Space Partitioning fitness tree is employed to partition the whole search space into several subregions according to the history positions of members. According to our preliminary observation of particles’ trajectory, most of the global best positions are local under a parent node or a small patch of leaf nodes in the fitness tree, especially in the later evolution period. Most of the nodes in fitness tree have a lack of vitality. Therefore, in the future work, we may modify the static fitness tree into a dynamic one with a node age or a lifecycle feature.(3)
We would also like to introduce the proposed dynamic multiswarm and nonrevisit strategies to other swarm-based algorithms. At the same time, applying MCPSO-PSH to address more real-world complex problems will also be studied in the future work.



---

*Source: 5193013-2017-06-27.xml*