<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Enhanced Opposition Differential Evolution Algorithm for Multimodal
Optimization

Shatendra Singh  Aruna Tiwari  [  [  Computer Science Engineering, IIT Indore,
Indore, Madhya Pradesh, India

###### Abstract

Most of the real-world problems are multimodal in nature that consists of
multiple optimum values. Multimodal optimization is defined as the process of
finding multiple global and local optima (as opposed to a single solution) of
a function. It enables a user to switch between different solutions as per the
need while still maintaining the optimal system performance. Classical
gradient-based methods fail for optimization problems in which the objective
functions are either discontinuous or non-differentiable. Evolutionary
Algorithms (EAs) are able to find multiple solutions within a population in a
single algorithmic run as compared to classical optimization techniques that
need multiple restarts and multiple runs to find different solutions. Hence,
several EAs have been proposed to solve such kinds of problems. However,
Differential Evolution (DE) algorithm is a population-based heuristic method
that can solve such optimization problems, and it is simple to implement. The
potential challenge in Multi-Modal Optimization Problems (MMOPs) is to search
the function space efficiently to locate most of the peaks accurately. The
optimization problem could be to minimize or maximize a given objective
function and we aim to solve the maximization problems on multimodal functions
in this study. Hence, we have proposed an algorithm known as Enhanced
Opposition Differential Evolution (EODE) algorithm to solve the MMOPs. The
proposed algorithm has been tested on IEEE Congress on Evolutionary
Computation (CEC) 2013 benchmark functions, and it achieves competitive
results compared to the existing state-of-the-art approaches.

###### keywords:

Multimodal Optimization, Evolutionary Algorithms, Niching, Differential
Evolution

##  1 Introduction

Evolutionary algorithms (EAs) are population-based approaches that have been
widely used for solving optimization problems having a single global optimum.
It is better to find several optima so that if an optimal solution cannot be
used due to resource constraints, other equally good solutions can be used
enabling the system performance stay at peak. Also, the knowledge of multiple
optimal solutions in the search space enables a user to get useful insights
into the properties of optimal solutions of the function. Many real-world
problems consist of multiple optima, such as holographic designapp1 , electro-
magnetic designapp2 , protein structure prediction app3 , and data mining app4
. Therefore, it is essential to find as many global optimizers as possible as
it gives users the flexibility to choose alternate solutions if one solution
cannot be chosen due to limited resource constraints. However, it is more
challenging to find multiple global optima than to find a single global
optimum. Many researchers have used EAsEAs  and swarm intelligence
algorithmsSwarm  to solve Multi-Modal Optimization Problems (MMOPs), such as
the Genetic Algorithm (GA) GA , Ant Colony Optimization (ACO) ACO , Estimation
of Distribution Algorithm (EDA)EDA , Particle Swarm Optimization (PSO)PSO ,
and Differential Evolution (DE)DE . However, they tend to lose the effective
balance between exploration and exploitation of the search space. Exploration
here means increasing the diversity of the population in the hope to find
better individuals while the term Exploitation here means refining the
existing candidate solutions to reach the best one in a local region(sub-
population). While solving the exploration vs exploitation dilemma, two major
problems arise. First, if diversity is increased too much then it may lead to
unnecessary fitness computations as there may be no optima in those areas of
search space where the search is being performed. Second, if exploitation is
performed in the undesired region of the search space, then it may get trapped
in local optima. Therefore, a different mechanism, based on classical EAs, is
required to locate multiple optima simultaneously. Niching niching  has been
widely used in the literature to help an EA maintain population diversity in
multimodal optimization. Some well-known niching techniques include crowding
niching1 , clearingniching2 , fitness sharing niching3 , and speciation
niching4 . Niching divides the population into several sub-populations and
each subpopulation is responsible for finding the optima within that sub-
population. Each sub-population can be considered as a species. Since the
fitness landscape may be uneven and complex, the niches could be of different
shapes and sizes. Therefore, it becomes a great challenge to form multiple
sub-populations out of the global population in such a way that a maximum
number of peaks are located by the sub-populations. Classical niching methods
are sensitive to the parameters they use; for example, crowding is sensitive
to the crowding size, and speciation is sensitive to the niching radius.
Hence, parameter-free or parameter-insensitive techniques have been developed
to improve niching, such as history-based topological speciation topological ,
and clustering clustering . Other adaptive learning strategies have been
developed to enhance the diversity for EAs, such as GAada1 , DEada2 , and
PSOada3 .

##  2 Preliminaries

###  2.1 Differential Evolution

DE is a stochastic method, and it was first introduced in 1997 by Storn and
PriceDE . It has been used for solving a wide variety of problems, including
multi-objective, constrained, dynamic, and multimodal optimization problems.
DE consists of several operations such as initialization, mutation, crossover,
and selection, similar to operations in traditional evolutionary algorithms.
The population is initialized randomly within the bounds of the dimensions of
the problem. Mutation is the process of changing the genetic sequence of
individuals. DE uses vector difference-based mutation operation, which
perturbs the current generation population members with a scaled difference of
randomly selected and distinct population members.

$\vec{V}_{i,G}=\vec{X}_{r_{1},G}+F*(\vec{X}_{r_{2},G}-\vec{X}_{r_{3}^{j},G})$
(1)

$\vec{V}_{i,G}$ is the offspring generated by adding parent
$\vec{X}_{r_{1},G}$ by scaled vector difference of parents $\vec{X}_{r_{2},G}$
and $\vec{X}_{r_{3}^{j},G}$. After mutation, DE generates a trial vector
$\vec{u}_{i,G}=u_{1,i,G},u2,i,G,\ldots,u_{D,i,G}$ by crossover operation to
enhance the potential diversity of the population. In the basic version, a
binomial crossover operation is employed by DE as

$u_{j,i,G}=\left\\{\begin{array}[]{c}v_{j,i,G},\text{ if
}\left(\operatorname{rand}_{i,j}[0,1]\leq\mathrm{CR}\right)\text{ or
}\left(j=j_{\text{rand }}\right)\\\ x_{j,i,G,}\indent\text{otherwise\hskip
113.81102pt}\end{array}\right.$ (2)

where i=1,2, $\ldots$, NP, j=1,2, $\ldots$, D, and the crossover rate (CR) is
a user-specified constant between 0 and 1 to control the fraction of parameter
values copied from the donor vector. $j_{\text{rand }}$ is an integer randomly
chosen from 1 to D, which ensures that $\vec{u}_{i,G}$ gets at least one
component from $\vec{v}_{i,G}$. Selection compares the target vector
$\vec{x}_{i,G}$ with the trial vector $\vec{u}_{i,G}$ in terms of the fitness
values to determine which vector may survive to the next generation

$\vec{x}_{i,G+1}=\left\\{\begin{array}[]{ll}\vec{u}_{i,G},&\text{ if
}f\left(\vec{u}_{i,G}\right)\leq f\left(\vec{x}_{i,G}\right)\\\
\vec{x}_{i,G},&\text{ otherwise }\end{array}\right.$ (3)

###  2.2 ODE Framework

ODE is a combination of the concept of opposition based learning and the DE
algorithm. Let $k\in[a,b]$ then opposite number (o) is defined as :
$o=a+b-k$b1 . ODE framework mainly consists of two sections, i.e. opposition-
based population initialization and opposition-based generation jumping.

####  2.2.1 Opposition-Based Population Initialization

The randomly generated initial population is used to produce an opposition-
based population. The opposite population is calculated as :
$OP_{i,j}=a_{j}+b_{j}-P_{i,j},i=1,2,\ldots,N_{p}\mathchar
24635\relax\;j=1,2,\ldots,D,$ where $P_{i,j},$ and $OP_{i,j}$ denote the
$j^{th}$ variable of the $i^{th}$ vector of the population and opposite
population, respectively b1 .

####  2.2.2  Opposition Based Generation Jumping

This concept helps in the generation of new individuals that are likely to be
more fitter than their corresponding parents. Jumping Rate($J_{r}$) is used to
control the degree of opposition of the offspring candidates. $J_{r}$ also
indirectly controls the diversity of the population. The higher the value of
the jumping rate, the more the search space shrinks. The value of $J_{r}$ is
chosen to be 0.3 based on b1 . The search space is reduced in each generation
and the bounds of the dimensions are updated according to the generation.
These bounds of each generation are used to compute the opposite population
using

$OP_{i,j}=\mathrm{MIN}_{j}^{p}+\mathrm{MAX}_{j}^{p}-P_{i,j}$ (4)

Here $\mathrm{MIN}_{j}^{p},\mathrm{MAX}_{j}^{p}$ represents the minimum and
maximum values of the variable in $j^{th}$ dimension in the current generation
b1 .

###  2.3 NBC-Minsize

Nearest Better Clustering (NBC)NBC  is a tree-based technique used in MMOPs to
divide the population into several species and each species tries to find an
optimal solution. NBC-MinsizeNBC-Minsize  is an enhanced version of NBC. The
parameter $\varphi$, in NBC, is responsible for controlling the number of
species in the population. Smaller values of $\varphi$ result in formation of
increasingly multiple species whereas higher values result in fewer species.
It is important to note that species with few members (e.g. 1 or 2) find it
difficult to evolve with the mutation operators of DE. Therefore, the
parameter $minsize$ is used to limit the minimal number of individuals in the
species. It is important to note that too small value of $minsize$ makes the
process ineffective whereas too large value would result in merging two
species as one. To be able to locate multiple peaks, it is wise to allocate
small values in the early stages and relatively higher values in the later
stages of the evolution. In the early stages, the focus of evolution must be
to explore the function space well, and hence multiple species formation is
encouraged by using lower values of $minsize$. In the later stages, the focus
should be on exploiting the search space and hence a larger $minsize$ value is
preferred. Thus, the species converge to the global optima more quickly.
Algorithm 1 refers to the NBC-Minsize.

Algorithm 1 NBC-Minsize

1:Set minsize by equation (5), (6);

2:Construct the spanning tree T ;

3:Calculate the mean distance $\mu_{dist}$ ;

4:Calculate the follow vector;

5:Sort the edges in T from the longest to the shortest;

6:for each $e\in T$ do

7: if $dist(e)$ $\textgreater$ $\varphi*\mu_{dist}$ then

8: Set $e_{f}$ to the follower individual of e;

9: Set $e_{r}$ to the root of the subtree containing $e_{f}$ ;

10: if $follow(e_{f})\geq minsize$ and

11: $follow(e_{r})-follow(e_{f}$) $\geq$ $minsize$ and

12: f($\frac{e_{f}+e_{r}}{2}$)<$f(e_{r}$) and
f($\frac{e_{f}+e_{r}}{2}$)<$f(e_{f}$) then

13: Cut off $e$;

14: Set $e_{l}$ to the leader individual of e;

15: for each $x$ on the path from $e_{l}$ to $e_{r}$ do

16: $follow(x)=follow(x)-follow(e_{f})$

17: end for

18: end if

19: end if

20:end for

The value of $minsize$ of species is set by equations (5) and (6) which are as
follows:

$minsize(g)=5+g/2$ (5)

where g represents the generation and $minsize(g)$ represents $minsize$ at
generation g. During the later stages, the values of g becomes increasingly
large, and hence to bound it, below equation (6) is used.

$bound=max(10,3*D)$ (6)

where D represents the dimensionality of the problem. A spanning tree is
constructed similar to the standard NBC and the mean distance $\mu_{dist}$ is
calculated. Next, we calculate the follow vector where every element of follow
represents the number of nodes in the subtree rooted at the corresponding
individual. The follow value is initially set to 1. The edges are sorted in
descending order according to the fitness values of the follower individuals.
Finally, for each edge, the follow value of each follower individual is added
to that of a leader individual. Also, the edges in T are sorted by their
Euclidean distance in descending order. Subsequently, edges satisfying the
conditions are cut off. Finally, the species are obtained after cutting the
edges.

##  3 Proposed Methodology

In this section, we present the opposition DE-based method for solving MMOPs.
The proposed method consists of five different components which are presented
in the form of Algorithms 2-8 where Algorithm 2 represents the generic
framework used for solving MMOPs. Algorithms 1, and 3-8 are further components
of Algorithm 2. The proposed method is tested on the IEEE CEC 2013 benchmark
functions.

In the next sub-section, we describe the multi-species framework used for
solving MMOPs.

###  3.1 Multi-Species Framework

In a multi-species framework, the idea is to divide the randomly distributed
population into sub-populations called species and perform the evolutionary
process on each species independentlymulti-pop . The block diagram represented
by Figure 1 broadly depicts the components of the EODE Algorithm.

Figure 1: Block Diagram of EODE Algorithm

In Figure 1, initialization is the process of randomly assigning values to the
population within bounds and it corresponds to step 1 of Algorithm 2. Two-
level speciation is used to divide the global population into local sub-
populations or species. Two-level speciation corresponds to section 3.1.1 and
step 4 of Algorithm 2. Species balance strategy is used to balance the
oversized and undersized species. It corresponds to section 3.1.2 and step 5
of Algorithm 2. An evolutionary process is employed to evolve the hybrid
population that has been generated by DE and ODE operations. The evolutionary
process corresponds to section 3.1.3 and step 7 of Algorithm 2 while local
search corresponds to section 3.1.5 and step 8 of Algorithm 2. Merge Archive
procedure checks for the peaks belonging to the same species to locate
multiple different peaks. It corresponds to section 3.1.6 and step 9 of
Algorithm 2. Fes represents the current count of fitness evaluations while
MaxFes represents the count of maximum fitness evaluations allowed. The
procedures such as two-level speciation, species balance strategy,
evolutionary process, local search and merge archive keep on running until the
current fitness count becomes greater than or equal to the maximum count of
fitness evaluations allowed. Algorithm 2 represents the generic framework
which further consists of various other algorithms.

Algorithm 2 EODE

1:: Function ($f(\vec{X})$),Population Size (NP), $MaxFes$

2:Initialize the population randomly within the bounds of the dimensions

3:$Gen=0$,$Fes=0$,$archive=[\hskip 2.84544pt]$

4:while $Fes\leq MaxFes$ do

5: Obtain multiple species by two-level application of Algorithm 2

6: balanceSpecies($multi$-$species$)

7: for $species\in multi$-$species$ do

8: $localbest$=Modified Opposition DE($species,speciesfitness$)

9: $bestfit$=localSearch($localbest$,$species$)

10: mergeArchive($archive$,$bestfit$)

11: end for

12: $Gen$+=1

13:end while

In algorithm 2, the population is initialized randomly within the bounds
specified for each dimension. The parameters $Gen$, $Fes$, $archive$
represents generation number, number of fitness evaluations, and archive,
respectively. The archive is used to store optimum values. Steps 4-9 are
performed until maximum fitness evaluations are reached. Algorithm 1 is
invoked for creating multiple species out of the global population in step 5.
After the formation of multiple species, there is a need to balance the
species as each species could be representing different niches in the function
landscape. To this aim, step 5 is introduced which balances the species based
on the shape and size of niches. Now, for each species that is present in
multi-species, steps 7-9 are carried out. In step 7, modified ODE (section
3.1.3) is applied to each species and returns the local best for that species.
Local search (section 3.1.5) is a method to refine the accuracy of the best
solution obtained so far. Step 8 tries to find the better solution nearby the
$localbest$ achieved. It may so happen that the $localbest$ obtained is very
very close to one of the optima already present in the archive and both are
part of the same niche. Hence, there is a need to check the redundancy of the
optimizers obtained. To this aim, step 9 is introduced.

In the next section, we describe the two-level speciation procedure using the
application of Algorithm 1 on the global population.

####  3.1.1 Two-level speciation

This section corresponds to step-4 of the multi-species framework described in
Algorithm 2. In this section, we describe how the global population is divided
into local sub-populations using Algorithm 1. In Algorithm 1, while using
mNBCFBK-DE  for speciation, it may so happen that the root node ($e_{r}$) and
follower node ($e_{f}$) belong to the same niche. In such a case, the existing
mNBC will cut off the edge between $e_{r}$ and $e_{f}$ resulting in the
formation of two different species even though they are part of the same
species. This will lead to wastage of fitness evaluations. Hence, to improve
the efficacy of mNBC, an additional check is applied in Algorithm 1. Step 12
of Algorithm 1 represents the check that checks for the presence of valley
between $e_{r}$ and $e_{f}$ by using an extra fitness evaluation. Even though
an extra fitness evaluation is used, it will save many unnecessary fitness
evaluations. We have applied two-level speciation to capture the narrow
regions of a peak. First-level speciation is performed with higher $minsize$
and hence it may associate multiple regions of a peak in a single species. To
identify those close species, we employ a second level of speciation within
the large species to segregate the nearby species. Second-level speciation is
performed with lower $minsize$. Hence, it introduces four parameters i.e.
$\varphi_{1},\varphi_{2},minsize_{1},minsize_{2}$. The parameters
$\varphi_{1},minsize_{1}$ are associated with first-level speciation and
$\varphi_{2},minsize_{2}$ are associated with second-level speciation.

In the next section we describe the procedure to balance the species that are
obtained after two-level speciation.

####  3.1.2 Species Balance Strategy

In this section, we describe a method to balance the species by redistributing
members across the different species. Step-5 in Algorithm 2 corresponds to the
species balancing strategy. After multiple species are obtained, some species
with a narrow and small basin of attraction could have many individuals and
some species with a wide and large basin of attraction could have fewer
individuals, hence there is a need to balance such unbalanced species. Since
the niches could be of different shapes and sizes, the species balance
strategy should be dynamic. In algorithm 3, initially, a minimum threshold of
10 individuals has to be present in the species which is taken care of by
steps 3-4. It is done to be able to compute the covariance matrix of $t$ best
individuals in the species indicated in step 6.

Algorithm 3 balanceSpecies

1:Input: $multi$-$species$

2:for $species\in multi$-$species$ do

3: if $speciessize\leq dim$ or $speciessize\leq$ 10 then

4: $c=max(dim-popsize,10)$

5: Initialize $c$ individuals around the species seed

6: $t=max(speciessize/Gen,10)$

7: Compute covariance matrix of $t$ best members of species

8: Compute variance of species using: $var=\sum_{i=1}^{k}eigen_{i}$

9: Average species size$(avgsize)=\frac{NP}{\text{No of Species}}$

10: end if

11: if $speciessize\textgreater\delta*avgsize$ then

12: Remove the $(speciessize-\delta*avgsize)$ worst individuals from the
species

13: end if

14:end for

15:$sortedspecies=$Sort the species in descending order of variances.

16:for $species\in sortedspecies$ do

17: if $(species\textless\delta*avgsize)$ then

18: Generate ($\delta*avgsize$-$speciessize$) number of individuals around the
species seed using covariance matrix computed in step 7 and Gaussian
distribution with standard deviation given by equation (3.8)

19: end if

20:end for

The variance or spread of a species is calculated by summing up the
eigenvalues of the covariance matrix. The eigen values are represented using
$eigen_{i}$ (step 7) where i runs from 1 to $k$ and $k$ represents the
dimensionality of the problems. The average species size is computed in step 8
which considers the global population size(NP). Step 10 checks if the size of
the species($speciessize$) is greater than the product of $\delta$ and the
average size of the species($avgsize$) then the individuals with the least
fitness values are removed from the species indicated by step 11. $\delta$ is
a hyperparameter here that controls the number of new individuals generated in
a species. The underbalanced species with large variance is balanced first and
hence step 14. For each unbalanced species, the new individuals are generated
in the vicinity of species seed using the Gaussian distribution. To this aim,
Steps 15-17 are introduced.

In the next section we present the modified opposition DE for evolving the
species that we obtained after balancing them.

####  3.1.3 Modified Opposition DE

In this section, we introduce the modified opposition DE to evolve the sub-
population. This procedure corresponds to step 7 mentioned in Algorithm 2
(section 3.1). Algorithm 4 consists of the application of modified opposition
DE to each species.
The notations used in algorithm 4 are as below:
$F_{1}set$: Set to store successful $\vec{F_{1}}$ values, $F_{2}set$:Set to
store successful $\vec{F_{2}}$ values, $CRset$: Set to store successful
$\vec{CR}$ values, $Max_{F}ES$: Maximum Fitness Evaluations, $FES$: Current
fitness evaluations, $Gen$: Current Generation Number, $MaxGen$: Maximum
Generation Number, $\vec{FB}$: First Best Member, $\vec{SB}$: Second Best
Member, $\vec{TB}$: Third Best Member, $\vec{X_{rk}}$: Population Member,
$pr$: Probability, NP: Population size, $\vec{U}$: Trial vector,
$\vec{F}_{1}$: Mutation factor, $\vec{F}_{2}$: Mutation factor, $\vec{CR}$:
Crossover rate, $JR$: Jumping Rate, $D$: Dimensions,$MIN_{j}$: Minimum in
$j^{th}$ dimension , $MAX_{j}$: Maximum in $j^{th}$ dimension, $OP_{ij}$:
Opposite member, $P_{ij}$: Current member , $P$: Current Population, $OP$:
Opposite Population
Three sets are taken to store the successful values of $\vec{F_{1}}$,
$\vec{F_{2}}$, and $\vec{CR}$. The successful values here mean the values
which lead to the production of better offspring. $pr$ represents the
probability with which the exploitation needs to be done. In early
generations, more exploration of search space needs to be done while in the
later phases more exploitation needs to be done. Hence, to this aim steps 7-18
are introduced. Step 9 represents the perturbation of lower degree as there is
only one mutation component while step 11 represents the perturbation to be of
higher degree as two mutation components are introduced. The moderate and
strong exploitation are introduced by Step 16 and 18 respectively.

Algorithm 4 Modified Opposition DE

1:$species,speciesfitness$

2:$F_{1}set=\phi$, $F_{2}set=\phi$, $CRset=\phi$

3:while $FES\textless MaxFES$ and $Gen\textless MaxGen$ do

4: for member$\in$ species do

5: Find first($\vec{FB}$), second($\vec{SB}$) and third($\vec{TB}$) best
members in the species.

6: Randomly sample 5
members$(\vec{X_{r1}},\vec{X_{r2}},\vec{X_{r3}},\vec{X_{r4}},\vec{X_{r5}})$
from the species

7: $pr=Gen/MaxGen$

8: if $pr\leq$0.33 then

9: if random(0,1)$\leq$ 0.75 then

10: $\vec{V}=\vec{X}_{r1}+\vec{F}_{1}*(\vec{X}_{r2}-\vec{X}_{r3})$

11: else

12:
$\vec{V}=\vec{X}_{r1}+\vec{F}_{1}*(\vec{X}_{r2}-\vec{X}_{r3})+\vec{F}_{2}*(\vec{X}_{r3}-\vec{X}_{r4})$

13: end if

14: else

15: if $pr\leq$0.67 then

16: Randomly select 2 members ($\vec{X_{k1}},\vec{X_{k2}}$) out of
$\frac{NP}{2}$ best members

17: $\vec{V}=\vec{FB}+\vec{F}_{1}*(\vec{X}_{k1}-\vec{X}_{K2})$

18: else

19: $\vec{V}=\vec{FB}+\vec{F_{1}}*(\vec{SB}-\vec{TB})$

20: end if

21: end if

22: Apply binomial crossover operation to get $\vec{U}$

23: Check for the bounds of $\vec{U}$ in each dimension using equation (5)

24: Evaluate the child fitness

25: FES+=1

26: if childfitness>parentfitness then

27: $F_{1}set\cup\vec{F_{1}},F_{2}set\cup\vec{F_{2}},CRset\cup\vec{CR},$

28: end if

29: end for

30: $JR=Gen/MaxGen$

31: if $JR\leq 1$ and $JR\textgreater 0.67$ then

32: Find the minimum, maximum bounds of the current species in each dimension.

33: for (i=0;i<NP;i++) do

34: if random(0,1)$\textless$ 0.33 then

35: for (j=0;j<D;j++) do

36: $OP_{i,j}=MIN_{j}+MAX_{j}-P_{i,j}$

37: end for

38: else

39: for (j=0;j<D;j++) do

40: $OP_{i,j}$=$MIN_{j}+random(0,1)*(MAX_{j}-MIN_{j}$)

41: end for

42: end if

43: end for

44: end if

45: Evaluate the opposite population

46: FES+=NP

47: Restart the candidates randomly that are stuck for 10 successive
generations.

48: Select the NP candidates with best fitness from $\\{P\cup OP\\}$

49: Update $\vec{F_{1}},\vec{F_{2}}$, and $\vec{CR}$

50: $Gen+=1$

51: Return the best member of the species

52:end while

Step 16 uses the candidate with the highest fitness value within that species
represented by $\vec{FB}$ with a random differential vector to create the
target vector. Step 18 uses the best 3 different candidates of the species.
The bounds are checked using the following equation.

$P_{i,d}=\begin{cases}min(ub[d],2*lb[d]-P_{i,d})&\text{if}(P_{i,d}\textless
lb[d])\\\ max(lb[d],2*ub[d]-P_{i,d})&\text{if}(P_{i,d}\textgreater
ub[d])\end{cases}$ (7)

where ’d’ represents dimension, ub[d] and lb[d] represents upper and lower
bounds in dimension ’d’. $P_{i,d}$ represents the value of $d^{th}$ dimension
of $i^{th}$ member in the population. The successful parameters such as
mutation factor and crossover rate are stored in their respective sets using
step 26. The jumping rate JR determines the degree of opposition applied to
the species to create the opposite population. The opposite population is
generated during the 33% of the evolution process to enable better
exploitation. The key point to note here is that the mutation factors
($\vec{F1},\vec{F2}$) and crossover rate ($\vec{CR}$) are calculated for each
dimension for every generation in the population. We have used the vector
representations of mutation factors and crossover rate to be able to capture
the degree of perturbation in every generation individually. This helps in
creating diverse offspring for better exploration.

While creating the opposite population, we have introduced a tweak while
creating the offspring. Steps 33 - 39 indicate the modification. Instead of
creating the opposite population in the shrunken space, some amount of
randomization is introduced 67% of the time to help in the exploration of the
function space. For the rest of the 33% of time, the opposite member is
calculated using step 35. During the evolution process, the population may get
trapped in local optima. To avoid the trap, the population members are checked
for improvement in fitness values, if they do not improve over some k (let say
10) successive generations as compared to global best member, then it is an
indication of potential trap at local optima. To this end, step 46 is
introduced. Step 48 is introduced at the end of every generation to update the
algorithm parameters using the adaptive parameter strategy described in the
next section.

####  3.1.4 Adaptive Parameter Strategy

In this section, we describe the adaptive parameter strategy to update the
mutation factors and crossover rate in such a way that it adapts to the
landscapes of different shapes and sizes and helps in the efficient evolution
of the species. The evolution process i.e. modified ODE is discussed in
section 3.1.3. The adaptive parameter strategy is used by modified ODE and it
indirectly corresponds to step 7 of Algorithm 2 (section 3.1). To update the
values of $\vec{F}_{1},\vec{F}_{2},\vec{CR}$, we have used a mechanism
inspired from SHADESHADE . The complete pseudocode is given by Algorithm 5 and
Algorithm 6. We have used two mutation factors ($\vec{F}_{1}$ and
$\vec{F}_{2}$) to guide the population towards the optima based on the stage
of evolution. Weighted power mean is used in the calculation of $\vec{F1}$ and
$\vec{CR}$ whereas weighted Lehmer mean is used in the calculation of the
$\vec{F2}$. Power mean is used to interpolate between the minimum and maximum
values using arithmetic mean and harmonic mean. Lehmer mean is used to capture
the non-linearity of the moving averages of the parameter values. Weights are
used to capture the improvement in influencing parameter adaptation. Weighted
Lehmer mean and weighted power mean are computed using equations (9) and (8)
respectively. The $S_{\vec{F}}$ represents the set of mutation factors that
produces a better offspring. $w_{k}$, computed using equation (10), represents
the weight, and $\Delta{f_{k}}$ represents the fitness difference between
parent and offspring. The consideration of all the dimensions while
calculating the $\vec{F1}$, $\vec{F2}$, $\vec{CR}$ for next-generation helps
predict the degree of perturbation to be introduced in each dimension for
better exploration. $random(0,1)$ represents the random number between 0 and 1
including them. $wlm$, $wpm$ represents weighted Lehmer mean and weighted
power mean respectively. $\vec{F_{old}}$ represents the $\vec{F}$ of the
previous generation. $\vec{max}$, $\vec{min}$ represents the boundaries of the
current species (sub-population) to which DE operations are applied.
$\vec{U}$, $\vec{L}$ represents the upper and lower bounds of the objective
function respectively. So, essentially the calculation of $\vec{F}$ considers
the older values of $\vec{F}$, the shape and size of the basin of attraction,
and the information about the stage of evolution to estimate the parameter
values of the algorithm for the next generation.

$Mean_{W,Power}(\vec{F}_{success})=\left(\frac{1}{\mathinner{\\!\left\lvert
S_{\vec{F}}\right\rvert}}\sum_{k=1}^{\mathinner{\\!\left\lvert
S_{\vec{F}}\right\rvert}}w_{k}*{S_{\vec{F},k}}\right)^{\frac{1}{1.5}}$ (8)

$\operatorname{Mean}_{W,Lehmer}\left(\vec{F}_{success}\right)=\frac{\sum_{k=1}^{\mathinner{\\!\left\lvert
S_{\vec{F}}\right\rvert}}w_{k}\cdot
S_{\vec{F},k}^{2}}{\sum_{k=1}^{\mathinner{\\!\left\lvert
S_{\vec{F}}\right\rvert}}w_{k}\cdot S_{\vec{F},k}}$ (9)

$w_{k}=\frac{\Delta f_{k}}{\sum_{k=1}^{\mathinner{\\!\left\lvert
S_{CR}\right\rvert}}\Delta f_{k}}$ (10)

where $Mean_{W,Power}(\vec{F}_{success})$,
$\operatorname{Mean}_{W,Lehmer}\left(\vec{F}_{success}\right)$ represents
weighted power mean and weighted Lehmer mean. $S_{\vec{F},k}$ represents a set
of successful mutation factors while $\mathinner{\\!\left\lvert
S_{\vec{F}}\right\rvert}$ represents the size of the set which contains
successful mutation factors. $\Delta f_{k}$ represents the fitness difference
between the parent and successful offspring.

Algorithm 5 Adaptive Parameter Strategy($\vec{F}$)

1:$wf=0.8+0.2*random(0,1)$

2:Compute weighted power mean for $\vec{F1}$ and weighted Lehmer mean for
$\vec{F2}$ using equation (8) and (9) respectively.

3:$\vec{F}=0.25*\vec{F}_{old}+0.25*(\vec{max}-\vec{min})/(\vec{U}-\vec{L})+0.5*(1-\frac{FES}{MaxFes})$

4:$\vec{F}=wf*\vec{F}_{old}+(1-wf)*wpm$(or $wlm$)

Algorithm 6 Adaptive Parameter Strategy($\vec{CR}$)

1:$wf=0.9+0.1*random(0,1)$

2:Compute weighted power mean using equation (8).

3:$\vec{CR}=wf*\vec{CR_{old}}+(1-wf)*wpm$

The adaptive parameter strategy for mutation factors is described in Algorithm
5. Algorithm 6 introduces an adaptive parameter strategy for crossover rate.
In the next section, we define the local search method to improve the best
solution obtained after applying modified ODE in the hope to further improve
its fitness.

####  3.1.5 Local Search

Local search methods are generally used to refine the accuracies of the
obtained solutions. It is being used in step 8 of Algorithm 2. Our proposed
local search method is described in Algorithm 7. Essentially, we generate some
members in the vicinity of the best candidate obtained so far in the species
in the hope to find a further better candidate. We initialize variances
($vars$) using step 3 described in algorithm 7. The key idea is to keep
incrementing the variances ($vars$) if the better offspring are not found in
the nearby region of the best individual of the species. dim represents the
dimensionality of the problem. $dirvec$, computed in step 18, can be
visualized as a direction vector and it is used to perform a guided search
nearby the $localbest$ candidate. The degree of the shift from the $localbest$
is continuously updated in the hope of finding a better solution by
incrementing the variances. In Algorithm 7, Gaussian distribution has been
used to generate offsprings around the $localbest$ in the current species. For
efficient search in the nearby region of the $localbest$ candidate, we
generate offspring by differential vector perturbation and estimating the
distribution of the handful of better candidates in the species. Mean and
standard deviation are computed using (9), (10).

$\mu_{i}^{d}=\frac{1}{M}\sum_{j=1}^{M}X_{j}^{d}$ (11)
$\delta_{i}^{d}=\sqrt{\frac{1}{M-1}\sum_{j=1}^{M}\left(X_{j}^{d}-\mu_{i}^{d}\right)^{2}}$
(12)

where $\mu_{i}=\left[\mu_{i}^{1},\ldots,\mu_{i}^{d},\ldots,\mu_{i}^{D}\right]$
and $\delta_{i}=$
$\left[\delta_{i}^{1},\ldots,\delta_{i}^{d},\ldots,\delta_{i}^{D}\right](1\leq
i\leq s)$ are, respectively, the mean and standard deviation (std) vectors of
the $i^{th}$ niche,
$\boldsymbol{X}_{j}=\left[X_{j}^{1},\ldots,X_{j}^{d},\ldots,X_{j}^{D}\right]$
is the $j^{th}$ individual in the $i^{th}$ niche and $D$ is the dimension size
of the multimodal problem.

Algorithm 7 localSearch

1::$localbest$,$species$

2:$dirvec$=None

3:k=10

4:$vars=random(0.001,0.01)$

5:while (k) do

6: $mean=localbest$

7: if $dirvec$ is not None and $random(0,1)\leq$ 0.5 then

8: offspring[k]=$localbest$+$vars$*$dirvec$

9: else

10: $mbest$=max(speciessize/4,10)

11: Find the covariance matrix of the $mbest$ candidates of species using
equation (11)

12: if dim$\textgreater$ 1 then

13: offspring[k]=Apply multivariate normal distribution with mean as $mean$
and covariance matrix as computed in step 10

14: else

15: offspring[k]=Apply univariate normal distribution with mean as $mean$ and
standard deviation is as computed using equation (10)

16: end if

17: end if

18: if $\vec{O}_{fitness}$ >$\vec{X}_{bestfitness}$ then

19: $dirvec$=offspring[k]-$localbest$

20: $localbest$=offspring[k]

21: else

22: for d in dim do

23: $vars[d]$+=$random(0.001,0.01)$

24: end for

25: end if

26: k=k-1

27:end while

28:Return the $localbest$

The covariance matrix is used to capture the behavior of variances of each
dimension with respect to the other dimensions. Covariance between two
variables of dimension D is computed as follows:

$\mathrm{COV_{X,Y}}=\frac{\sum_{i=1}^{D}\left(X_{i}-\bar{x}\right)\left(Y_{i}-\bar{y}\right)}{D-1}$
(13)

where $\bar{x},\bar{y}$ represents X mean and Y mean.
$\vec{O}_{fitness},\vec{X}_{bestfitness}$ represents the offspring fitness and
member with the best fitness respectively. In the next section, the
mergeArchive procedure is discussed that aims at removing the peaks belonging
to the same niche.

####  3.1.6 Merge Archive

In this section, we describe the procedure to deal with duplicate peaks in the
$archive$. This section corresponds to step 9 of Algorithm 2. Since we are
storing the peaks in an $archive$, it is important to store the peaks only.
Otherwise, we may end up storing the entire candidates of a sub-population
(species) and subsequently the entire population. This will lead to the
consumption of more computational resources since the $archive$ will contain
huge points. Hence, given two individuals, it is essential to check whether
they belong to similar species or different species. It is critical to
algorithm performance in locating all the global optima. In algorithm 8, the
inputs are $localbest$, $archive$ and it needs to be checked whether the
$localbest$ is a new species seed or an old species seed. The nearest peak to
the given local best is found out, and the midpoint between them is
calculated. If the midpoint has lower fitness than the nearest peak and
$localbest$, it indicates the presence of a valley between them. Hence, local
best can be considered as the species seed of new species, and it is added to
the $archive$. In the next section, we describe the experiments performed and
the results obtained.

Algorithm 8 mergeArchive

1:$archive,localbest$

2:Find the nearest peak to $localbest$ present in the $archive$

3:Find the mid-point between $localbest$ and nearest peak

4:if (midpointfitness >$localbest$ and midpointfitness >nearest peak) then

5: Add $localbest$ to the $archive$

6:else

7: Replace the better fit individual between $localbest$ and nearest peak with
mid-point.

8:end if

###  3.2 Illustrative Example

In order to explain the components of EODE, an example is presented in this
section. We have considered low dimensional (2-D) composition function
$F_{11}$.

Figure 2: Initialization

Figure 3: Niching and Balancing

Figure 4: ODE and Local Search

As it can be seen from Figure 2, the population is randomly initialized in the
function space. Now we will try to map the steps given in Algorithm 2 to the
processes that happen in Fig. 2-4.

  1. 1.

Initially, the randomly distributed population needs to be divided into sub-
populations that represent regions around the peaks. The multiple species
marked as black circles in Fig 3 are shown as S1-S7. In species S4, the blue
circles represent two species that got located as a single species S4. To
avoid such merging, two-level speciation is introduced. To this aim, step 4
(Algorithm 2) is introduced.

  2. 2.

As the multiple species are formed, some species may have huge population size
and some may have small population size as can be seen from Figure 3. Hence,
the species balance strategy is introduced that aims to generate offspring
using Gaussian distribution. Also, the species balance strategy is discussed
in Algorithm 3 which maps to step 5 in Algorithm 2.

  3. 3.

Modified ODE (section 3.1.3) acts as an evolutionary process that tries to
further converge the sub-populations towards the peak they represent. The
species are represented in Figure 4 as P1-P8. Further, it is more effective to
converge towards global optima rather than local optima. It can be seen in
Figure 4 in which P7 and P8 are local optima. This process maps to step 7 in
Algorithm 2.

  4. 4.

After the evolution process completes, the population gets converged near the
optima but may not be able to reach the optima accurately as can be seen in
Figure 4. Hence, we proposed a local search method to search in the vicinity
of the best individual of the species as the actual peak would be nearby the
best value in the species. The local search method is described in Algorithm 7
and it maps to step 8 of Algorithm 2.

  5. 5.

Since the speciation process is repeated after every certain number of
generations, the best solutions returned by the local search method could be
part of the same species that has been explored in the previous generations.
Hence, there is a need to detect the best individuals that belong to the same
species and they should not be added to the archive. Therefore, Algorithm 8 is
introduced that maps to step 9 of Algorithm 2.

Table 1: Information of the benchmark problems and the population size Index | Function | NKP | Peak height | $r$ | MaxFEs | $NP$
---|---|---|---|---|---|---
1 | $F_{1}(1D)$ | 2 | $200.0$ | $0.01$ | $5.0\mathrm{E}+4$ | 250
2 | $F_{2}(1D)$ | 5 | $1.0$ | $0.01$ | $5.0\mathrm{E}+4$ | 250
3 | $F_{3}(1D)$ | 1 | $1.0$ | $0.01$ | $5.0\mathrm{E}+4$ | 250
4 | $F_{4}(2D)$ | 4 | $200.0$ | $0.01$ | $5.0\mathrm{E}+4$ | 250
5 | $F_{5}(2D)$ | 2 | $1.03163$ | $0.5$ | $5.0\mathrm{E}+4$ | 250
6 | $F_{6}(2D)$ | 18 | $186.731$ | $0.5$ | $2.0\mathrm{E}+5$ | 2000
7 | $F_{7}(2D)$ | 36 | $1.0$ | $0.2$ | $2.0\mathrm{E}+5$ | 2000
8 | $F_{6}(3D)$ | 81 | $2709.0935$ | $0.5$ | $4.0\mathrm{E}+5$ | 3000
9 | $F_{7}(3D)$ | 216 | $1.0$ | $0.2$ | $4.0\mathrm{E}+5$ | 4000
10 | $F_{8}(2D)$ | 12 | $-2.0$ | $0.01$ | $2.0\mathrm{E}+5$ | 1000
11 | $F_{9}(2D)$ | 6 | 0 | $0.01$ | $2.0\mathrm{E}+5$ | 1000
12 | $F_{10}(2D)$ | 8 | 0 | $0.01$ | $2.0\mathrm{E}+5$ | 1000
13 | $F_{11}(2D)$ | 6 | 0 | $0.01$ | $2.0\mathrm{E}+5$ | 1000
14 | $F_{11}(3D)$ | 6 | 0 | $0.01$ | $4.0\mathrm{E}+5$ | 1000
15 | $F_{12}(3D)$ | 8 | 0 | $0.01$ | $4.0\mathrm{E}+5$ | 1000
16 | $F_{11}(5D)$ | 6 | 0 | $0.01$ | $4.0\mathrm{E}+5$ | 1000
17 | $F_{12}(5D)$ | 8 | 0 | $0.01$ | $4.0\mathrm{E}+5$ | 2000
18 | $F_{11}(10D)$ | 6 | 0 | $0.01$ | $4.0\mathrm{E}+5$ | 1000
19 | $F_{12}(10D)$ | 8 | 0 | $0.01$ | $4.0\mathrm{E}+5$ | 1000
20 | $F_{12}(20D)$ | 8 | 0 | $0.01$ | $4.0\mathrm{E}+5$ | 800
Table 2: Parameters in EODE Parameters | Values
---|---
$\varphi_{1}$ | $1$
$\varphi_{2}$ | $1$
$minsize_{1}$ | $-1$
$minsize_{2}$ | $5$
$\delta$ | $1.0$
$F1$ | $(0,1)$
$F2$ | $(0,1)$
$CR$ | $(0,1)$
$MaxGen$ | $40(D<=10)$
| $60(D\textgreater 10)$

##  4 Experiments and Results

We have performed the experiments on a computer system with RAM 8GB, 1.8GHz
CPU and MacOS 11 operating system. In this section, EODE is independently run
50 times for each function. The algorithmic results are computed in five
levels of accuracy $\epsilon$ =
{$1\mathrm{e}{-1},1\mathrm{e}{-2},1\mathrm{e}{-3},1\mathrm{e}{-4},1\mathrm{e}{-5}$}.

Table 3: Results on accuracy levels 1e-1,1e-2,1e-3,1e-4,1e-5 | 1E-1 | 1E-2 | 1E-3 | 1E-4 | 1E-5
---|---|---|---|---|---
Functions | (PR,SR) | (PR,SR) | (PR,SR) | (PR,SR) | (PR,SR)
F1(1D) | (1,1) | (1,1) | (1,1) | (1,1) | (1,1)
F2(1D) | (1,1) | (1,1) | (1,1) | (1,1) | (1,1)
F3(1D) | (1,1) | (1,1) | (1,1) | (1,1) | (1,1)
F4(2D) | (1,1) | (1,1) | (1,1) | (1,1) | (1,1)
F5(2D) | (1,1) | (1,1) | (1,1) | (1,1) | (1,1)
F6(2D) | (1,1) | (0.995,0.9) | (0.995,0.9) | (0.995,0.9) | (0.9,0.822)
F7(2D) | (0.805,0) | (0.805,0) | (0.805,0) | (0.805,0) | (0.805,0)
F6(3D) | (0.886,0) | (0.852,0) | (0.852,0) | (0.845,0) | (0.832,0)
F7(3D) | (0.509,0) | (0.505,0) | (0.505,0) | (0.505,0) | (0.442,0)
F8(2D) | (1,1) | (1,1) | (1,1) | (1,1) | (1,1)
F9(2D) | (1,1) | (1,1) | (1,1) | (1,1) | (1,1)
F10(2D) | (0.975,0.8) | (0.975,0.8) | (0.975,0.8) | (0.975,0.8) | (0.975,0.8)
F11(2D) | (1,1) | (1,1) | (1,1) | (1,1) | (1,1)
F11(3D) | (0.8,0) | (0.8,0) | (0.8,0) | (0.8,0) | (0.8,0)
F12(3D) | (0.8,0) | (0.8,0) | (0.8,0) | (0.8,0) | (0.77,0)
F11(5D) | (0.733,0) | (0.730,0) | (0.730,0) | (0.730,0) | (0.728,0)
F12(5D) | (0.7,0) | (0.7,0) | (0.7,0) | (0.684,0) | (0.684,0)
F11(10D) | (0.7,0) | (0.7,0) | (0.7,0) | (0.684,0) | (0.684,0)
F12(10D) | (0.525,0) | (0.525,0) | (0.525,0) | (0.520,0) | (0.505,0)
F12(20D) | (0.25,0) | (0.25,0) | (0.25,0) | (0.25,0) | (0.25,0)

For different problems, the value of NP is shown in Table 1. To decide upon
the values of NP, the algorithm is run for 50 times with different population
sizes for each function and that NP is chosen which provides best value of
average PR over 50 runs for the given function. The results are depicted by
the Table 4. The $NP_{1},NP_{2},NP_{3},NP_{4}$ and $NP_{5}$ represents
different population sizes for each function and the
$PR_{1},PR_{2},PR_{3},PR_{4}$ and $PR_{5}$ represents the corresponding
average PR values for the given population sizes $NP_{i}$.

Table 4: PR values on different population sizes Functions | $NP_{1}$ | $PR_{1}$ | $NP_{2}$ | $PR_{2}$ | $NP_{3}$ | $PR_{3}$ | $NP_{4}$ | $PR_{4}$ | $NP_{5}$ | $PR_{5}$
---|---|---|---|---|---|---|---|---|---|---
F1(1D) | 50 | 1 | 100 | 1 | 250 | 1 | 500 | 1 | 1000 | 1
F2(1D) | 50 | 1 | 100 | 1 | 250 | 1 | 500 | 1 | 1000 | 1
F3(1D) | 50 | 1 | 100 | 1 | 250 | 1 | 500 | 1 | 1000 | 1
F4(2D) | 50 | 0.6 | 100 | 0.93 | 250 | 1 | 500 | 1 | 1000 | 1
F5(2D) | 50 | 1 | 100 | 1 | 250 | 1 | 500 | 1 | 1000 | 1
F6(2D) | 500 | 0.966 | 1000 | 0.968 | 2000 | 0.995 | 3000 | 0.985 | 5000 | 0.911
F7(2D) | 500 | 0.466 | 1000 | 0.695 | 2000 | 0.805 | 3000 | 0.797 | 5000 | 0.659
F6(3D) | 500 | 0.506 | 1000 | 793 | 2000 | 0.829 | 3000 | 0.845 | 8000 | 0.817
F7(3D) | 1000 | 0.264 | 2000 | 0.479 | 4000 | 0.505 | 8000 | 0.493 | 10000 | 0.427
F8(2D) | 250 | 0.65 | 500 | 0.895 | 1000 | 1 | 2000 | 0.969 | 3000 | 0.912
F9(2D) | 100 | 0.633 | 250 | 0.998 | 500 | 1 | 1000 | 1 | 2000 | 1
F10(2D) | 100 | 0.475 | 250 | 0.781 | 500 | 0.919 | 1000 | 0.975 | 2000 | 0.958
F11(2D) | 100 | 0.366 | 250 | 0.766 | 500 | 0.933 | 1000 | 1 | 2000 | 0.993
F11(3D) | 100 | 0.366 | 250 | 0.688 | 500 | 0.788 | 1000 | 0.8 | 2000 | 0.697
F12(3D) | 100 | 0.4 | 250 | 0.795 | 500 | 0.639 | 1000 | 0.8 | 2000 | 0.795
F11(5D) | 500 | 0.466 | 1000 | 0.730 | 2000 | 0.719 | 5000 | 0.685 | 8000 | 0.644
F12(5D) | 500 | 0.433 | 1000 | 0.684 | 2000 | 0.598 | 5000 | 0.515 | 8000 | 0.489
F11(10D) | 500 | 0.466 | 1000 | 0.684 | 2000 | 0.681 | 5000 | 0.589 | 8000 | 0.512
F12(10D) | 500 | 0.375 | 1000 | 0.520 | 2000 | 0.505 | 5000 | 0.318 | 8000 | 0.305
F12(20D) | 500 | 0.175 | 800 | 0.250 | 1000 | 0.218 | 2000 | 0.198 | 5000 | 0.195

The parameters of EODE are listed in Table 2. $\varphi_{1}$, $\varphi_{2}$ are
used for creating multiple species out of the global population by applying
mNBC twice and their values are set as 1, 1 respectively. $minsize_{1}$,
$minsize_{2}$ are parameters used in two-level application of mNBC and their
values are set as -1, 5 respectively. We define the $minsize_{1}$ as -1 to
adapt the first level species as per the dimension and generation of the
population. $minsize_{2}$ is set to be 5 to identify the peaks that have
narrow basins of attraction as the species size would be very low in such
regions. $\delta$ is used in the species balance strategy and it is set to 1.

The mutation factors ($\vec{F1},\vec{F2}$) and crossover rate ($\vec{CR}$) are
initially chosen from the closed interval [0,1] after which the algorithm
learns to adapt their values. $MaxGen$ represents the number of generations
for which the ODE needs to run for the given species. The experimental results
of EODE for all benchmark problems are listed in Table 3 at all five accuracy
levels. These results show that EODE is very stable. In addition, EODE finds
all peaks on the simple functions and most of the peaks in the low-dimensional
composition problems. EODE can find more than 80% of all peaks on the
functions containing a large number of global peaks, except for F7(3D).

###  4.1 Comparison with other algorithms

We have compared the results of various algorithms (including EODE) in this
section. We have compared the results at $\epsilon=1\mathrm{e}{-4}$, which are
commonly adopted in Wangg  and Yangg . To have an extensive comparison, we
have evaluated the performance of EODE and compared it with 15 popular
algorithms such as CDE crowding , SDEDE5 , NCDE, NSDE Qu-DE , MOMMOP MOMMOP ,
LoICDE, LoISDE ada2 , PNPCDE pnpcde , LIPS Qu , and the recently proposed
algorithms of $DE_{cl}$ DEcl , LMCEDA, LMSEDA EDA , FBK-DEFBK-DE , LBPADEZhao
, MaHDEMaHDE . Table 4, 5, 6 shows the different PRs and SRs at the accuracy
level $\epsilon=1\mathrm{e}{-4}$, where the bottom row “bprs” represents the
number of functions that achieved best PR values for these algorithms. The
algorithm(s) with the best PR value for a given function is marked in bold.
All the algorithms used their own default population sizes. The results of the
compared algorithms have been referred from their respective papers. From
Table 4, 5, 6, it is clear that EODE obtains the best PR values for most of
the functions among the compared algorithms. The detailed analysis is given
below.
1) For functions 1-5, most of the algorithms including EODE can find all
global optimal solutions.
2) For functions 6-9 which consists of a large number of global peaks, EODE
does not perform equally well as compared to all the algorithms. However, it
achieved better results as compared to many algorithms. For problem 6, the
difference is of order 0.005 which is very small. Also, for the $10^{th}$
problem, EODE can find all the optima.
3) For functions 11,12,13 and 15 that are low-dimensional composition
functions, EODE obtained the best results. Although on the $14^{th}$ function
EODE does not produce optimal results, the difference between EODE and the
best algorithm is small (around 10%). It is important to note that EODE and
FBK-DE are the two algorithms that are able to find all global optima on the
$13^{th}$ function while the other algorithms do not.
4) For the functions $16^{th}$–$20^{th}$ that are composition functions having
relatively high dimensionality, EODE achieves the best results for all
functions except the $20^{th}$ function. It is worth noting that on the 5-D
problems, the result obtained by EODE exceeds 5% of the best result obtained
by other algorithms. For function 18, EODE performed a little better than the
best algorithm while for the $19^{th}$ function it performed equivalently
well.

Table 5: Comparison of EODE with other algorithms on accuracy level 1e-4 Function | EODE | CDE | SDE | NCDE | NSDE | MOMMOP
---|---|---|---|---|---|---
Index | PR | SR | PR | SR | PR | SR | PR | SR | PR | SR | PR | SR
1 | 1 | 1 | 1 | 1 | 0.657 | 0.373 | 1 | 1 | 1 | 1 | 1 | 1
2 | 1 | 1 | 1 | 1 | 0.737 | 0.529 | 1 | 1 | 0.776 | 0.667 | 1 | 1
3 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
4 | 1 | 1 | 1 | 1 | 0.284 | 0 | 1 | 1 | 0.240 | 0 | 1 | 1
5 | 1 | 1 | 1 | 1 | 0.922 | 0.843 | 1 | 1 | 0.745 | 0.490 | 1 | 1
6 | 0.995 | 0.9 | 1 | 1 | 0.056 | 0 | 0.305 | 0 | 0.056 | 0 | 1 | 1
7 | 0.805 | 0 | 0.861 | 0 | 0.054 | 0 | 0.873 | 0 | 0.053 | 0 | 1 | 1
8 | 0.845 | 0 | 0 | 0 | 0.015 | 0 | 0.001 | 0 | 0.013 | 0 | 1 | 1
9 | 0.505 | 0 | 0.474 | 0 | 0.011 | 0 | 0.461 | 0 | 0.006 | 0 | 1 | 1
10 | 1 | 1 | 1 | 1 | 0.147 | 0 | 0.989 | 0.863 | 0.098 | 0 | 1 | 1
11 | 1 | 1 | 0.330 | 0 | 0.314 | 0 | 0.729 | 0.059 | 0.248 | 0 | 0.716 | 0.020
12 | 0.975 | 0.8 | 0.002 | 0 | 0.208 | 0 | 0.252 | 0 | 0.135 | 0 | 0.939 | 0.549
13 | 1 | 1 | 0.141 | 0 | 0.297 | 0 | 0.667 | 0 | 0.225 | 0 | 0.667 | 0
14 | 0.8 | 0 | 0.026 | 0 | 0.216 | 0 | 0.667 | 0 | 0.190 | 0 | 0.667 | 0
15 | 0.8 | 0 | 0.005 | 0 | 0.108 | 0 | 0.319 | 0 | 0.125 | 0 | 0.618 | 0
16 | 0.730 | 0 | 0 | 0 | 0.108 | 0 | 0.667 | 0 | 0.170 | 0 | 0.650 | 0
17 | 0.684 | 0 | 0 | 0 | 0.076 | 0 | 0.250 | 0 | 0.108 | 0 | 0.505 | 0
18 | 0.684 | 0 | 0.167 | 0 | 0.026 | 0 | 0.500 | 0 | 0.163 | 0 | 0.497 | 0
19 | 0.520 | 0 | 0 | 0 | 0.105 | 0 | 0.348 | 0 | 0.098 | 0 | 0.223 | 0
20 | 0.250 | 0 | 0 | 0 | 0 | 0 | 0.250 | 0 | 0.123 | 0 | 0.125 | 0
bprs | 14 | 7 | 1 | 5 | 2 | 10

Table 6: Comparison of EODE with other algorithms on accuracy level 1e-4 Function | EODE | LoICDE | LoISDE | PNPCDE | LIPS | $DE_{cl}$
---|---|---|---|---|---|---
Index | PR | SR | PR | SR | PR | SR | PR | SR | PR | SR | PR | SR
1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0.833 | 0.686 | 1 | 1
2 | 1 | 1 | 1 | 1 | 0.235 | 0.039 | 1 | 1 | 1 | 1 | 1 | 1
3 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0.961 | 0.961 | 1 | 1
4 | 1 | 1 | 0.975 | 0.902 | 0.250 | 0 | 1 | 1 | 0.990 | 0.961 | 1 | 1
5 | 1 | 1 | 1 | 1 | 0.667 | 0.333 | 1 | 1 | 1 | 1 | 1 | 1
6 | 0.995 | 0.9 | 1 | 1 | 0.056 | 0 | 0.537 | 0 | 0.246 | 0 | 0.942 | 0.340
7 | 0.805 | 0 | 0.705 | 0.02 | 0.029 | 0 | 0.874 | 0 | 0.4 | 0 | 0.986 | 0.64
8 | 0.845 | 0 | 0 | 0 | 0.012 | 0 | 0 | 0 | 0.084 | 0 | 0.999 | 0.9
9 | 0.505 | 0 | 0.187 | 0 | 0.005 | 0 | 0.472 | 0 | 0.104 | 0 | 0.726 | 0
10 | 1 | 1 | 1 | 1 | 0.083 | 0 | 1 | 1 | 0.748 | 0 | 1 | 1
11 | 1 | 1 | 0.66 | 0 | 0.167 | 0 | 0.66 | 0 | 0.974 | 0.843 | 0.667 | 0
12 | 0.975 | 0.8 | 0.495 | 0 | 0.125 | 0 | 0 | 0 | 0.574 | 0 | 0.943 | 0.58
13 | 1 | 1 | 0.51 | 0 | 0.167 | 0 | 0.461 | 0 | 0.794 | 0.176 | 0.667 | 0
14 | 0.8 | 0 | 0.657 | 0 | 0.167 | 0 | 0.592 | 0 | 0.644 | 0 | 0.667 | 0
15 | 0.8 | 0 | 0.299 | 0 | 0.125 | 0 | 0.258 | 0 | 0.336 | 0 | 0.623 | 0
16 | 0.730 | 0 | 0.559 | 0 | 0.167 | 0 | 0 | 0 | 0.304 | 0 | 0.667 | 0
17 | 0.684 | 0 | 0.223 | 0 | 0.076 | 0 | 0 | 0 | 0.162 | 0 | 0.42 | 0
18 | 0.684 | 0 | 0.219 | 0 | 0.157 | 0 | 0.147 | 0 | 0.098 | 0 | 0.667 | 0
19 | 0.520 | 0 | 0.037 | 0 | 0.027 | 0 | 0 | 0 | 0 | 0 | 0.357 | 0
20 | 0.250 | 0 | 0.123 | 0 | 0.088 | 0 | 0 | 0 | 0 | 0 | 0.212 | 0
bprs | 14 | 6 | 2 | 6 | 2 | 6

Table 7: Comparison of EODE with other algorithms on accuracy level 1e-4 Function | EODE | LMCEDA | LMSEDA | FBK-DE | LBPADE | MaHDE
---|---|---|---|---|---|---
Index | PR | SR | PR | SR | PR | SR | PR | SR | PR | SR | PR | SR
1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
2 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
3 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
4 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
5 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
6 | 0.995 | 0.9 | 0.99 | 0.843 | 0.972 | 0.588 | 0.990 | 0.820 | 1 | 1 | 1 | 1
7 | 0.805 | 0 | 0.734 | 0 | 0.673 | 0 | 0.813 | 0 | 0.889 | 0 | 0.804 | 0
8 | 0.845 | 0 | 0.367 | 0 | 0.613 | 0 | 0.824 | 0 | 0.575 | 0 | 0.983 | 0.291
9 | 0.505 | 0 | 0.284 | 0 | 0.248 | 0 | 0.425 | 0 | 0.476 | 0 | 0.351 | 0
10 | 1 | 1 | 1 | 1 | 0.998 | 0.98 | 1 | 1 | 1 | 1 | 1 | 1
11 | 1 | 1 | 0.667 | 0 | 0.892 | 0.392 | 1 | 1 | 0.674 | 0 | 0.725 | 0.078
12 | 0.975 | 0.8 | 0.75 | 0 | 0.99 | 0.922 | 0.935 | 0.480 | 0.750 | 0 | 0.650 | 0
13 | 1 | 1 | 0.667 | 0 | 0.667 | 0 | 1 | 1 | 0.667 | 0 | 0.667 | 0
14 | 0.8 | 0 | 0.667 | 0 | 0.667 | 0 | 0.907 | 0.460 | 0.667 | 0 | 0.667 | 0
15 | 0.8 | 0 | 0.696 | 0 | 0.738 | 0 | 0.730 | 0 | 0.654 | 0 | 0.648 | 0
16 | 0.730 | 0 | 0.667 | 0 | 0.667 | 0 | 0.707 | 0 | 0.667 | 0 | 0.667 | 0
17 | 0.684 | 0 | 0.456 | 0 | 0.62 | 0 | 0.630 | 0 | 0.532 | 0 | 0.352 | 0
18 | 0.684 | 0 | 0.657 | 0 | 0.66 | 0 | 0.667 | 0 | 0.667 | 0 | 0.663 | 0
19 | 0.520 | 0 | 0.451 | 0 | 0.458 | 0 | 0.520 | 0 | 0.475 | 0 | 0.455 | 0
20 | 0.250 | 0 | 0.059 | 0 | 0.248 | 0 | 0.450 | 0 | 0.275 | 0 | 0.250 | 0
bprs | 14 | 6 | 5 | 10 | 7 | 7

###  4.2 Different Components of EODE

In this section, we have discussed the different components embedded in EODE
and their effect on the algorithm. We have also discussed the mutation
operators that are used based on evolution stages. This section also contains
a brief discussion about the effect of $\varphi_{1}$ and $\varphi_{2}$ that
are key parameters in deciding the formation of the total number of species.
We briefly describe the effects of JR and the convergence of population
towards global optima.

####  4.2.1 Different Mutation Operators

We have compared the results of the algorithms using different mutation
operators one at a time. The balancing ability of the algorithm is shown by
comparing the four different algorithms denoted as EODE, EODE-r, EODE-b, EODE-
rb. DE/rand/1 and DE/rand/2 are used as mutation operators in EODE-r to evolve
the species. Besides, EODE-b uses both DE/best/1 and DE/best/2 while EODE-rb
uses DE/rand/1, DE/rand/2, DE/best/1 and DE/best/2 for mutation, which can be
referred from mutation .

Table 8: Experimental results of different mutation operators on the benchmark problems at the accuracy level $\epsilon=1\mathrm{e}{-4}$ Function | EODE | EODE-r | EODE-b | EODE-rb
---|---|---|---|---
Index | PR | SR | PR | SR | PR | SR | PR | SR
1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
2 | 1 | 1 | 1 | 1 | 0.88 | 0.4 | 1 | 1
3 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
4 | 1 | 1 | 1 | 1 | 1 | 1 | 0.75 | 0.4
5 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
6 | 0.995 | 0.9 | 1 | 1 | 0.997 | 0 | 1 | 1
7 | 0.805 | 0 | 0.734 | 0 | 0.733 | 0 | 0.788 | 0
8 | 0.845 | 0 | 0.567 | 0 | 0.306 | 0 | 0.809 | 0
9 | 0.505 | 0 | 0.314 | 0 | 0.452 | 0 | 0.425 | 0
10 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
11 | 1 | 1 | 1 | 1 | 0.933 | 0.4 | 1 | 1
12 | 0.975 | 0.8 | 0.925 | 0 | 0.8 | 0 | 1 | 1
13 | 1 | 1 | 0.735 | 0 | 0.55 | 0 | 0.667 | 0
14 | 0.8 | 0 | 0.667 | 0 | 0.505 | 0 | 0.667 | 0
15 | 0.8 | 0 | 0.580 | 0 | 0.612 | 0 | 0.785 | 0
16 | 0.730 | 0 | 0.409 | 0 | 0.441 | 0 | 0.549 | 0
17 | 0.684 | 0 | 0.186 | 0 | 0.279 | 0 | 0.667 | 0
18 | 0.684 | 0 | 0.147 | 0 | 0 | 0 | 0.558 | 0
19 | 0.520 | 0 | 0 | 0 | 0 | 0 | 0.395 | 0
20 | 0.250 | 0 | 0 | 0 | 0 | 0 | 0.125 | 0
bprs | 18 | 8 | 5 | 8

The two different mutation operators are used with a probability of 0.5 each.
The other components except mutation operators are not changed to note down
the effects of mutation operators. The results corresponding to the usage of
different mutation operators are shown in Table 8 at accuracy level $\epsilon$
= $1\mathrm{e}{-4}$. It is clear from Table 8 that EODE produces the best PR
values for all benchmark functions except for the $6^{th}$ and $12^{th}$
functions. The performances of EODE-r and EODE-b are not so good as no peaks
are found on the relatively high-dimensional functions.

Besides, Table 8 also shows that most of the compared algorithms perform well
on functions 1-5, while EODE-b and EODE-rb are unable to find all the global
peaks for functions 2 and 4 respectively. EODE-r finds all peaks in the
functions 1-5. Here, we specifically compare the results of EODE and EODE-rb
to illustrate the effect of three stage-wise mutation operators. As it can be
referred from the previous experiment, the list functions are categorized in
four different ways and the results on them are discussed below.
1) On the functions 1-5 which are considered to be simple, both EODE and EODE-
rb locates all the optima except for function 4.
2) For functions 6-10 that consist of many global optima, it is clear that the
performance of EODE is significantly better than that of EODE-rb except for
function 6, which shows that the usage of mutation operators as per the stage
of evolution process is highly effective.
3) For functions 11-15 which are low-dimensional composition functions, the PR
values for EODE are superior to those of EODE-rb except for function 12. For
functions 11-12, EODE-rb either performs better or equivalently well. For
problems 13-15, EODE-rb does not perform well as compared to EODE.
4) On the relatively high-dimensional composition problems (the
$16{th}–20^{th}$ problems), EODE performs better than EODE-rb, indicating that
the mutation operations selected as per the stage of the evolution process can
improve the overall performance of the algorithm.

####  4.2.2 Different values of $\varphi_{1},\varphi_{2}$

In this section we have compared the results for different values of
$\varphi_{1},\varphi_{2}$. In Section 3.1.1, $\varphi_{1},\varphi_{2}$ are
used to form local sub-populations from global population. Thus, we set
($\varphi_{1},\varphi_{2}$) to three different values (1,1), (0.6,0.6), and
(2,1) respectively. From Table 7 it can be seen that for
$\varphi_{1}=1,\varphi_{2}=1$, EODE achieves the best overall results. For
$\varphi_{1}=0.6,\varphi_{2}=0.6$, EODE achieves best result for function 15.
For functions 1-5, $\varphi_{1}=2,\varphi_{2}=1$ achieves better results as
compared to that of $\varphi_{1}=0.6,\varphi_{2}=0.6$. For functions 6-9
having large number of optima, $\varphi_{1}=0.6,\varphi_{2}=0.6$ achieves
better results than $\varphi_{1}=2,\varphi_{2}=1$. However, for higher
dimensional functions 16-20, EODE with $\varphi_{1}=2,\varphi_{2}=1$ seems to
perform better except for functions 16 and 20. The behaviour is attributed to
the number of species formed out of the global population. If the number of
optima is large, then lower value of $\varphi_{1},\varphi_{2}$ is effective
while if the number of optima is small, then higher value of
$\varphi_{1},\varphi_{2}$ is effective.

###  4.3 Different values of Jumping Rate (JR)

Here, different values of JR are compared and it can be referred to from Table
9. In section 3.1.3, JR is used to form the opposite population. We modified
the JR values to see if the opposite population in the exploration stage helps
or they help in the exploitation stage.

Table 9: Experimental results of different values of $\varphi_{1},\varphi_{2}$ on the benchmark problems at the accuracy level $\epsilon=1\mathrm{e}{-4}$ Function |  $\varphi_{1}$=1,$\varphi_{2}$=1 |  $\varphi_{1}$=0.6,$\varphi_{2}$=0.6 |  $\varphi_{1}$=2,$\varphi_{2}$=1
---|---|---|---
Index | PR | SR | PR | SR | PR | SR
1 | 1 | 1 | 1 | 1 | 1 | 1
2 | 1 | 1 | 0.68 | 0.2 | 0.92 | 0.6
3 | 1 | 1 | 1 | 1 | 1 | 1
4 | 1 | 1 | 0.5 | 1 | 0.9 | 0.6
5 | 1 | 1 | 1 | 1 | 1 | 1
6 | 0.995 | 0.9 | 1 | 1 | 1 | 1
7 | 0.805 | 0 | 0.744 | 0 | 0.688 | 0
8 | 0.845 | 0 | 0.723 | 0 | 0.650 | 0
9 | 0.505 | 0 | 0.476 | 0 | 0.420 | 0
10 | 1 | 1 | 1 | 1 | 1 | 1
11 | 1 | 1 | 1 | 1 | 0.95 | 0.9
12 | 0.975 | 0.8 | 0.95 | 0.4 | 1 | 1
13 | 1 | 1 | 0.733 | 0 | 0.95 | 0.9
14 | 0.8 | 0 | 0.733 | 0 | 0.667 | 0
15 | 0.8 | 0 | 1 | 1 | 0.8 | 0
16 | 0.730 | 0 | 0.709 | 0 | 0.667 | 0
17 | 0.684 | 0 | 0.515 | 0 | 0.580 | 0
18 | 0.684 | 0 | 0.342 | 0 | 0.620 | 0
19 | 0.520 | 0 | 0.125 | 0 | 0.418 | 0
20 | 0.250 | 0 | 0.175 | 0 | 0.125 | 0
bprs | 19 | 7 | 7

JR is dependent upon the current generation and the maximum generation allowed
using equation (12):

$JR=\frac{Gen}{MaxGen}$ (14)

$0\leq JR\leq 0.5$ indicates that the opposite population is generated in the
early stage of evolution i.e. exploration while $0.67\leq JR\leq 1$ indicates
that the opposite population is generated in the later stage of evolution i.e.
exploitation. $0\leq JR\leq 1$ indicates that the opposite population is
generated across the evolution stages. It can be seen that the opposite
population generated in the later phase of evolution turns out to be helpful
in locating more peaks. For the functions 1-5, EODE with $0.67\leq JR\leq 1$
and $0\leq JR\leq 1$ performs equivalently well. For functions 6-9, EODE with
$0\leq JR\leq 1$ performs better than EODE with $0\leq JR\leq 0.5$ except the
function 8. For functions 17 and 20, EODE with $0\leq JR\leq 1$ achieves best
results similar to that of EODE with $0.67\leq JR\leq 1$.

Table 10: Experimental results of different values of JR on the benchmark problems at the accuracy level $\epsilon=1\mathrm{e}{-4}$ Function | 0.67 $\leq$ JR $\leq$ 1 | 0 $\leq$ JR $\leq$ 0.5 | 0 $\leq$ JR $\leq$ 1
---|---|---|---
Index | PR | SR | PR | SR | PR | SR
1 | 1 | 1 | 1 | 1 | 1 | 1
2 | 1 | 1 | 0.84 | 0.6 | 1 | 1
3 | 1 | 1 | 1 | 1 | 1 | 1
4 | 1 | 1 | 0.95 | 0.9 | 1 | 1
5 | 1 | 1 | 1 | 1 | 1 | 1
6 | 0.995 | 0.9 | 1 | 1 | 1 | 1
7 | 0.805 | 0 | 0.722 | 0 | 0.768 | 0
8 | 0.845 | 0 | 0.701 | 0 | 0.682 | 0
9 | 0.505 | 0 | 0.433 | 0 | 0.488 | 0
10 | 1 | 1 | 1 | 1 | 0.9 | 0.4
11 | 1 | 1 | 1 | 1 | 0.933 | 0.6
12 | 0.975 | 0.8 | 0.95 | 0.4 | 0.925 | 0.4
13 | 1 | 1 | 0.866 | 0.2 | 0.9 | 0.4
14 | 0.8 | 0 | 0.667 | 0 | 0.667 | 0
15 | 0.8 | 0 | 0.667 | 0 | 0.6 | 0
16 | 0.730 | 0 | 0.667 | 0 | 0.670 | 0
17 | 0.684 | 0 | 0.480 | 0 | 0.684 | 0
18 | 0.684 | 0 | 0.486 | 0 | 0.470 | 0
19 | 0.520 | 0 | 0.225 | 0 | 0.345 | 0
20 | 0.250 | 0 | 0.150 | 0 | 0.250 | 0
bprs | 19 | 6 | 8

###  4.4 Convergence of population

It is an important property of an evolutionary algorithm to converge the
population towards optima. The algorithm needs to maintain the found optima
towards the end of evolution and still search for other optima. To assess the
performance of EODE on converging the population towards optima, the solution
distributions of some functions (i.e., F2, F4, F6, F7, F10, F11, F12, and
F13), are presented on some specific generations. Figures 5-7 show the
solution distribution of F2 with different generations (i.e., the generations
are 1, 2, and 3).

Figure 5: Gen-1 (F2)

Figure 6: Gen-2 (F2)

Figure 7: Gen-3 (F2)

Figure 8: Gen-1 (F4)

Figure 9: Gen-5 (F4)

Figure 10: Gen-9 (F4)

Figure 11: Gen-1 (F6)

Figure 12: Gen-2 (F6)

Figure 13: Gen-3 (F6)

Figure 14: Gen-1 (F7)

Figure 15: Gen-8 (F7)

Figure 16: Gen-15 (F7)

Figure 17: Gen-1 (F10)

Figure 18: Gen-5 (F10)

Figure 19: Gen-11 (F10)

Figure 20: Gen-1 (F11)

Figure 21: Gen-8 (F11)

Figure 22: Gen-14 (F11)

Figure 23: Gen-1 (F12)

Figure 24: Gen-5 (F12)

Figure 25: Gen-11 (F12)

Figure 26: Gen-1 (F13)

Figure 27: Gen-5 (F13)

Figure 28: Gen-11 (F13)

It can be clearly seen that EODE has located all of the global optima when the
generation is 3 on F2 [i.e., Figure 7]. Hence, it can be concluded that EODE
can locate the global optima quickly which indicates that EODE has a good
ability of global search. It is important to note that the generation here
refers to the number of times a multi-species framework is performed given in
Algorithm 2. It does not refer to the generation of DE operations directly.
Figures 8-10 show the solution distribution of F4 with different generations
(i.e., 1, 5, and 9). It can be inferred from Figure 10 that the population
converges completely in the $9^{th}$ generation. Similarly, for F7, F10, EODE
obtained complete convergence. However, for F6 (figs 11-13), F11( figs 20-22),
F12(figs 23-25), and F13 (figs 26-28) there is some dispersion from the peak
but still, the convergence is nearly complete. Hence, it can be concluded that
as the generation increases the solutions move towards a convergence state.
Also, it is clear that the population does not diverge from the optima with
the increasing evolution in EODE. Hence, our EODE can converge towards the
global optima by the end of evolution.

##  5 Conclusion

In this paper, we have proposed EODE to handle MMOPs, where several species
techniques have been adopted to simultaneously find multiple optima. In EODE,
two-level NBC-minsize has been used to divide the population into several
species to locate the optima. To eliminate the drawbacks of NBC-minsize an
additional check is also performed to improve its efficacy. Besides, a species
balance strategy has been proposed to balance the species having uneven sizes.
The modified opposition differential evolution algorithm acts as the
evolutionary process. The evolutionary process is further divided virtually
into three stages and mutation operations for each of the stages are embedded
individually. Besides, an adaptive parameter strategy has also been proposed
to learn the parameters of the algorithm dynamically. We have also proposed a
local search method to improve the quality of candidate solutions. Finally,
EODE has been compared with several state-of-the-art algorithms. The
experimental results show that EODE performed better than other algorithms on
most multimodal benchmark problems. However, in our experiments, we have
achieved better results for problems having a moderate number of global
optima, but it still needs to improve on highly multimodal functions. Hence,
we aim to improve EODE in our future work and test it on highly multimodal
problems.

##  6 Acknowledgement

The authors would like to thank Xin Lin and Wenjian Luo for providing source
code of NBC-Minsize  that provided some helpful insights into the design and
working of the proposed algorithm.

##  7 Declarations

###  7.1 Funding and/or Conflicts of interests/Competing interests/Data
Availability Statements

This work is part of an academic course and it has not been funded by any
organization. The author(s) declare no competing interests. The dataset can be
found at https://github.com/mikeagn/CEC2013.

## References

  *   * (1) Qing, L., Gang, W., Qiuping, W.: Restricted evolution based multimodal function optimization in holographic grating design. In: 2005 IEEE Congress on Evolutionary Computation, vol. 1, pp. 789–7941 (2005). https://doi.org/10.1109/CEC.2005.1554763
  * (2) Woo, D.-K., Choi, J.-H., Ali, M., Jung, H.-K.: A novel multimodal optimization algorithm applied to electromagnetic optimization. IEEE Transactions on Magnetics 47(6), 1667–1673 (2011). https://doi.org/%****␣sn-article.bbl␣Line␣75␣****10.1109/TMAG.2011.2106218
  * (3) Wong, K.-C., Leung, K., Wong, M.H.: Protein structure prediction on a lattice model via multimodal optimization techniques, vol. 2010, pp. 155–162 (2010). https://doi.org/10.1145/1830483.1830513
  * (4) Boughanem, M., Tamine, L.: A study on using genetic niching for query optimisation in document retrieval. In: Crestani, F., Girolami, M., van Rijsbergen, C.J. (eds.) Advances in Information Retrieval, pp. 135–149. Springer, Berlin, Heidelberg (2002). https://doi.org/10.1007/3-540-45886-7_10
  * (5) Covantes Osuna, E., Sudholt, D.: Runtime analysis of crowding mechanisms for multimodal optimization. IEEE Transactions on Evolutionary Computation 24(3), 581–592 (2020). https://doi.org/10.1109/TEVC.2019.2914606
  * (6) Li, Y., Zhan, Z.-H., Lin, S., Zhang, J., Luo, X.: Competitive and cooperative particle swarm optimization with information sharing mechanism for global optimization problems. Information Sciences 293, 370–382 (2015). https://doi.org/10.1016/j.ins.2014.09.030
  * (7) Yao, J., Kharma, N., Grogono, P.: Bmpga: a bi-objective multi-population genetic algorithm for multi-modal function optimization. In: 2005 IEEE Congress on Evolutionary Computation, vol. 1, pp. 816–8231 (2005). https://doi.org/10.1109/CEC.2005.1554767
  * (8) Yang, Q., Chen, W.-N., Yu, Z., Gu, T., Li, Y., Zhang, H., Zhang, J.: Adaptive multimodal continuous ant colony optimization. IEEE Transactions on Evolutionary Computation 21(2), 191–205 (2017). https://doi.org/10.1109/TEVC.2016.2591064
  * (9) Li, Y., Zhan, Z.-H., Lin, S., Zhang, J., Luo, X.: Competitive and cooperative particle swarm optimization with information sharing mechanism for global optimization problems. Information Sciences 293, 370–382 (2015). https://doi.org/10.1016/j.ins.2014.09.030
  * (10) Cao, Y., Zhang, H., Li, W., Zhou, M., Zhang, Y., Chaovalitwongse, W.A.: Comprehensive learning particle swarm optimization algorithm with local search for multimodal functions. IEEE Transactions on Evolutionary Computation 23(4), 718–731 (2019). https://doi.org/10.1109/TEVC.2018.2885075
  * (11) Wang, Z.-J., Zhan, Z.-H., Lin, Y., Yu, W.-J., Wang, H., Kwong, S., Zhang, J.: Automatic niching differential evolution with contour prediction approach for multimodal optimization problems. IEEE Transactions on Evolutionary Computation 24(1), 114–128 (2020). https://doi.org/10.1109/TEVC.2019.2910721
  * (12) Mahfoud, S.W.: Niching methods for genetic algorithms. PhD thesis, USA (1996). UMI Order No. GAX95-43663
  * (13) De jong, k. a. (1975). an analysis of the behavior of a class of genetic adaptive systems (doc-. (2007)
  * (14) Petrowski, A.: A clearing procedure as a niching method for genetic algorithms. In: Proceedings of IEEE International Conference on Evolutionary Computation, pp. 798–803 (1996). https://doi.org/10.1109/ICEC.1996.542703
  * (15) Holland, J.H.: Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence. The MIT Press, ??? (1992). https://doi.org/10.7551/mitpress/1090.001.0001. https://doi.org/10.7551/mitpress/1090.001.0001
  * (16) Della Cioppa, A., De Stefano, C., Marcelli, A.: Where are the niches? dynamic fitness sharing. IEEE Transactions on Evolutionary Computation 11(4), 453–465 (2007). https://doi.org/10.1109/TEVC.2006.882433
  * (17) Stoean, C., Preuss, M., Stoean, R., Dumitrescu, D.: Multimodal optimization by means of a topological species conservation algorithm. IEEE Transactions on Evolutionary Computation 14(6), 842–864 (2010). https://doi.org/10.1109/TEVC.2010.2041668
  * (18) Gao, W., Yen, G.G., Liu, S.: A cluster-based differential evolution with self-adaptive strategy for multimodal optimization. IEEE Transactions on Cybernetics 44(8), 1314–1327 (2014). https://doi.org/10.1109/TCYB.2013.2282491
  * (19) Yao, J., Kharma, N., Grogono, P.: Bi-objective multipopulation genetic algorithm for multimodal function optimization. IEEE Transactions on Evolutionary Computation 14(1), 80–102 (2010). https://doi.org/10.1109/TEVC.2009.2017517
  * (20) Biswas, S., Kundu, S., Das, S.: Inducing niching behavior in differential evolution through local information sharing. IEEE Transactions on Evolutionary Computation 19(2), 246–263 (2015). https://doi.org/10.1109/TEVC.2014.2313659
  * (21) Li, X.: Niching without niching parameters: Particle swarm optimization using a ring topology. IEEE Transactions on Evolutionary Computation 14(1), 150–169 (2010). https://doi.org/%****␣sn-article.bbl␣Line␣375␣****10.1109/TEVC.2009.2026270
  * (22) Rahnamayan, S., Tizhoosh, H.R., Salama, M.M.A.: Opposition-based differential evolution. IEEE Transactions on Evolutionary Computation 12(1), 64–79 (2008). https://doi.org/10.1109/TEVC.2007.894200
  * (23) Preuss, M.: Niching the cma-es via nearest-better clustering, pp. 1711–1718 (2010). https://doi.org/10.1145/1830761.1830793
  * (24) Lin, X., Luo, W., Xu, P.: Differential evolution for multimodal optimization with species by nearest-better clustering. IEEE Transactions on Cybernetics 51(2), 970–983 (2021). https://doi.org/10.1109/TCYB.2019.2907657
  * (25) Zaharie, D.: A multipopulation differential evolution algorithm for multimodal optimization (2004)
  * (26) Lin, X., Luo, W., Xu, P.: Differential evolution for multimodal optimization with species by nearest-better clustering. IEEE Transactions on Cybernetics 51(2), 970–983 (2021). https://doi.org/10.1109/TCYB.2019.2907657
  * (27) Tanabe, R., Fukunaga, A.: Success-history based parameter adaptation for differential evolution. In: 2013 IEEE Congress on Evolutionary Computation, pp. 71–78 (2013). https://doi.org/10.1109/CEC.2013.6557555
  * (28) Wang, Z.-J., Zhan, Z.-H., Lin, Y., Yu, W.-J., Yuan, H.-Q., Gu, T.-L., Kwong, S., Zhang, J.: Dual-strategy differential evolution with affinity propagation clustering for multimodal optimization problems. IEEE Transactions on Evolutionary Computation 22(6), 894–908 (2018). https://doi.org/10.1109/TEVC.2017.2769108
  * (29) Yang, Q., Chen, W.-N., Yu, Z., Gu, T., Li, Y., Zhang, H., Zhang, J.: Adaptive multimodal continuous ant colony optimization. IEEE Transactions on Evolutionary Computation 21(2), 191–205 (2017). https://doi.org/10.1109/TEVC.2016.2591064
  * (30) Thomsen, R.: Multimodal optimization using crowding-based differential evolution. In: Proceedings of the 2004 Congress on Evolutionary Computation (IEEE Cat. No.04TH8753), vol. 2, pp. 1382–13892 (2004). https://doi.org/10.1109/CEC.2004.1331058
  * (31) Li, X.: Efficient differential evolution using speciation for multimodal function optimization, pp. 873–880 (2005). https://doi.org/10.1145/1068009.1068156
  * (32) Qu, B.Y., Suganthan, P.N., Liang, J.J.: Differential evolution with neighborhood mutation for multimodal optimization. IEEE Transactions on Evolutionary Computation 16(5), 601–614 (2012). https://doi.org/10.1109/TEVC.2011.2161873
  * (33) Wang, Y., Li, H.-X., Yen, G.G., Song, W.: Mommop: Multiobjective optimization for locating multiple optimal solutions of multimodal optimization problems. IEEE Transactions on Cybernetics 45(4), 830–843 (2015). https://doi.org/10.1109/TCYB.2014.2337117
  * (34) Biswas, S., Kundu, S., Das, S.: An improved parent-centric mutation with normalized neighborhoods for inducing niching behavior in differential evolution. IEEE Transactions on Cybernetics 44(10), 1726–1737 (2014). https://doi.org/10.1109/TCYB.2013.2292971
  * (35) Qu, B.Y., Suganthan, P.N., Das, S.: A distance-based locally informed particle swarm model for multimodal optimization. IEEE Transactions on Evolutionary Computation 17(3), 387–402 (2013). https://doi.org/10.1109/TEVC.2012.2203138
  * (36) Bošković, B., Brest, J.: Clustering and differential evolution for multimodal optimization. In: 2017 IEEE Congress on Evolutionary Computation (CEC), pp. 698–705 (2017). https://doi.org/%****␣sn-article.bbl␣Line␣600␣****10.1109/CEC.2017.7969378
  * (37) Zhao, H., Zhan, Z.-H., Lin, Y., Chen, X., Luo, X.-N., Zhang, J., Kwong, S., Zhang, J.: Local binary pattern-based adaptive differential evolution for multimodal optimization problems. IEEE Transactions on Cybernetics 50(7), 3343–3357 (2020). https://doi.org/10.1109/TCYB.2019.2927780
  * (38) Hong, Z., Chen, Z.-G., Liu, D., Zhan, Z.-H., Zhang, J.: A multi-angle hierarchical differential evolution approach for multimodal optimization problems. IEEE Access 8, 178322–178335 (2020). https://doi.org/10.1109/ACCESS.2020.3027559
  * (39) Islam, S.M., Das, S., Ghosh, S., Roy, S., Suganthan, P.N.: An adaptive differential evolution algorithm with novel mutation and crossover strategies for global numerical optimization. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 42(2), 482–500 (2012). https://doi.org/10.1109/TSMCB.2011.2167966