1

A Review of Evolutionary Multi-modal
Multi-objective Optimization

Solution space

Objective space

Ryoji Tanabe, Member, IEEE,and Hisao Ishibuchi, Fel-

low, IEEE
Fig. 1: Illustration of a situation where the four solutions are identical
Abstract—Multi-modal multi-objective optimization aims tofind all Pareto optimal solutions including overlapping solutionsin the objective space. Multi-modal multi-objective optimization
has been investigated in the evolutionary computation community
since 2005. However, it is difficult to survey existing studies
in this field because they have been independently conducted
and do not explicitly use the term “multi-modal multi-objective
optimization”. To address this issue, this paper reviews existingstudies of evolutionary multi-modal multi-objective optimization,
including studies published under names that are different from“multi-modal multi-objective optimization”. Our review also
clarifies open issues in this research area.
Index Terms—Multi-modal multi-objective optimization, evolutionary algorithms, test problems, performance indicators

I. I NTRODUCTION
A multi-objective evolutionary algorithm (MOEA) is anefficient optimizer for a multi-objective optimization problem(MOP) [1]. MOEAs aim to find a non-dominated solution
set that approximates the Pareto front in the objective space.
The set of non-dominated solutions found by an MOEA is
usually used in an “a posteriori” decision-making process [2].
A decision maker selects a final solution from the solution set
according to her/his preference.
Since the quality of a solution set is usually evaluated in
the objective space, the distribution of solutions in the solution
space has not received much attention in the evolutionary
multi-objective optimization (EMO) community. However, thedecision maker may want to compare the final solution to otherdissimilar solutions that have an equivalent quality or a slightly
inferior quality [3], [4]. Fig. 1 shows a simple example. In Fig.
1, the four solutions xa , xb , xc , and xd are far from each otherin the solution space but close to each other in the objective
space. xa and xb have the same objective vector. xc and xa
are similar in the objective space. xd is dominated by thesesolutions. This kind of situation can be found in a numberof real-world problems, including functional brain imaging
problems [3], diesel engine design problems [5], distillation
plant layout problems [6], rocket engine design problems [7],and game map generation problems [8].
If multiple diverse solutions with similar objective vectors
like xa , xb , xc , and xd in Fig. 1 are obtained, the decision
maker can select the final solution according to her/his preference in the solution space. For example, if xa in Fig. 1
becomes unavailable for some reason (e.g., material shortages,
R. Tanabe and H. Ishibuchi are with Shenzhen Key Laboratory of Computational Intelligence, University Key Laboratory of Evolving Intelligent Systemsof Guangdong Province, Department of Computer Science and Engineering,Southern University of Science and Technology, Shenzhen 518055, China.
e-mail: (rt.ryoji.tanabe@gmail.com, hisao@sustc.edu.cn). (Corresponding author: Hisao Ishibuchi)

or close to each other in the objective space but are far from eachother in the solution space (a minimization problem).

mechanical failures, traffic accidents, and law revisions), thedecision maker can select a substitute from xb , xc , and xd .
A practical example is given in [4], which deals with twoobjective space mission design problems. In [4], Schütze et al.
considered two dissimilar solutions x1 = (782, 1288, 1788)Tand x2 = (1222, 1642, 2224)T for a minimization problem,
whose objective vectors are f (x1 ) = (0.462, 1001.7)T andf (x2 ) = (0.463, 1005.3)T , respectively. Although x1 dominates x2 , the difference between f (x1 ) and f (x2 ) is smallenough. The first design variable is the departure time fromthe Earth (in days). Thus, the departure times of x2 and x1
differ by 440 days (= 1222 − 782). If the decision maker
accepts x2 with a slightly inferior quality in addition to x1 ,the two launch plans can be considered. If x1 is not realizablefor some reason, x2 can be the final solution instead of x1 .
As explained here, multiple solutions with almost equivalentquality support a reliable decision-making process. If thesesolutions have a large diversity in the solution space, they canprovide insightful information for engineering design [3], [5].
A multi-modal multi-objective optimization problem(MMOP) involves finding all solutions that are equivalent
to Pareto optimal solutions [3], [9], [10]. Below, we explain
the difference between MOPs and MMOPs using the twoobjective and two-variable Two-On-One problem [11]. Figs.
2 (a) and (b) show the Pareto front F and the Pareto optimal
solution set O of Two-On-One, respectively. Two-On-One
has two equivalent Pareto optimal solution subsets O 1 andO 2 that are symmetrical with respect to the origin, where
O = O 1 ∪ O 2 . Figs. 2 (c) and (d) show O 1 and O 2 ,respectively. In Two-On-One, the three solution sets O, O 1 ,and O 2 (Figs. 2 (b), (c) and (d)) are mapped to F (Fig. 2(a)) by the objective functions. On the one hand, the goal of
MOPs is generally to find a solution set that approximates
the Pareto front F in the objective space. Since O 1 andO 2 are mapped to the same F in the objective space, it issufficient for MOPs to find either O 1 or O 2 . On the other
hand, the goal of MMOPs is to find the entire equivalent
Pareto optimal solution set O = O 1 ∪ O 2 in the solution
space. In contrast to MOPs, it is necessary to find both O 1
and O 2 in MMOPs. Since most MOEAs (e.g., NSGA-II[12] and SPEA2 [13]) do not have mechanisms to maintainthe solution space diversity, it is expected that they do notwork well for MMOPs. Thus, multi-modal multi-objective
evolutionary algorithms (MMEAs) that handle the solution
space diversity are necessary for MMOPs.
This paper presents a review of evolutionary multi-modal

5
4
3
2
1
0

x2

f2

2

8 10 12 14 16 18 20
f1

2
1
0
−1
−2
−2 −1

0
x1

(c) O 1

1

2

1

2

(b) O

x2

x2

(a) F

2
1
0
−1
−2
−2 −1

0
x1

1

2

2
1
0
−1
−2
−2 −1

0
x1

(d) O 2

Fig. 2: (a) The Pareto front F and (b) the Pareto optimal solution
set O of Two-On-One [11]. Figs. (c) and (d) show the two Pareto
optimal solution subsets O 1 and O 2 , respectively.

multi-objective optimization. This topic is not new and has
been studied for more than ten years. Early studies include [3],[5], [11], [14]–[16]. Unfortunately, most existing studies wereindependently conducted and did not use the term “MMOPs”(i.e., they are not tagged). For this reason, it is difficult tosurvey existing studies of MMOPs despite their significant
contributions. In this paper, we review related studies of
MMOPs including those published under names that weredifferent from “multi-modal multi-objective optimization”. We
also clarify open issues in this field. Multi-modal singleobjective optimization problems (MSOPs) have been wellstudied in the evolutionary computation community [10].
Thus, useful clues to address some issues in studies of MMOPs
may be found in studies of MSOPs. We discuss what can belearned from the existing studies of MSOPs.
This paper is organized as follows. Section II gives definitions of MMOPs. Section III describes MMEAs. Section
IV presents test problems for multi-modal multi-objective
optimization. Section V explains performance indicators forbenchmarking MMEAs. Section VI concludes this paper.

2) Definitions of MMOPs: The term “MMOP” was firstcoined in [3], [14] in 2005. However, “MMOP” was not usedin most studies from 2007 to 2012. Terms that represent
MMOPs were not explicitly defined in those studies. For
example, MMOPs were referred to as problems of obtaining
a diverse solution set in the solution space in [17]. It seems
that “multi-modal multi-objective optimization” has been usedagain as of 2016. Apart from these instances, MMOPs were
denoted as “Multi-objective multi-global optimization” and“Multi-modal multi-objective wicked problems” in [18] and[19], respectively.
Although MMOPs have been addressed for more than ten
years, the definition of an MMOP is still controversial. In
this paper, we define an MMOP using a relaxed equivalency
introduced by Rudolph and Preuss [17] as follows:
Definition 1. An MMOP involves finding all solutions thatare equivalent to Pareto optimal solutions.
Definition 2. Two different solutions x1 and x2 are said tobe equivalent iff kf (x1 ) − f (x2 )k ≤ δ.
where kak is an arbitrary norm of a, and δ is a non-negative
threshold value given by the decision maker. If δ = 0, theMMOP should find all equivalent Pareto optimal solutions. If
δ > 0, the MMOP should find all equivalent Pareto optimal
solutions and dominated solutions with acceptable quality. The
main advantage of our definition of an MMOP is that thedecision maker can adjust the goal of the MMOP by changing
the δ value. Most existing studies (e.g., [9], [20], [21]) assume
MMOPs with δ = 0. MMOPs with δ > 0 were discussed in[3], [4], [19], [22]. For example, xa , xb , and xc in Fig. 1
should be found for MMOPs with δ = 0. In addition, thenon-Pareto optimal solution xd should be found for MMOPs
with δ > 0 if kf (xd ) − f (xa )k ≤ δ.
Although there is room for discussion, MMOPs with δ > 0
may be more practical in real-world applications. This isbecause the set of solutions of an MMOP with δ > 0 can
provide more options for the decision maker than that ofan MMOP with δ = 0. While it is usually assumed in the
EMO community that the final solution is selected from nondominated solutions, the decision maker may also be interested
in some dominated solutions in practice [3], [4]. Below, we
use the term “MMOP” regardless of the δ value for simplicity.
III. MMEA S

II. D EFINITIONS OF MMOP S
1) Definition of MOPs: A continuous MOP involves finding a solution x ∈ S ⊆ RD that minimizes a given objectivefunction vector f : S → RM . Here, S is the D-dimensional
solution space, and RM is the M -dimensional objective space.
A solution x1 is said to dominate x2 iff fi (x1 ) ≤ fi (x2 ) for
all i ∈ {1, ..., M } and fi (x1 ) < fi (x2 ) for at least one index
i. If x∗ is not dominated by any other solutions, it is called a
Pareto optimal solution. The set of all x∗ is the Pareto optimal
solution set, and the set of all f (x∗ ) is the Pareto front. The
goal of MOPs is generally to find a non-dominated solution
set that approximates the Pareto front in the objective space.

This section describes 12 dominance-based MMEAs, 3decomposition-based MMEAs, 2 set-based MMEAs, and a
post-processing approach. MMEAs need the following threeabilities: (1) the ability to find solutions with high quality,(2) the ability to find diverse solutions in the objective space,and (3) the ability to find diverse solutions in the solution
space. MOEAs need the abilities (1) and (2) to find a solution
set that approximates the Pareto front in the objective space.
Multi-modal single-objective optimizers need the abilities (1)and (3) to find a set of global optimal solutions. In contrast,
MMEAs need all abilities (1)–(3). Here, we mainly describe
mechanisms of each type of MMEA to handle (1)–(3).

31) Pareto dominance-based MMEAs: The most representative MMEA is Omni-optimizer [9], [14], which is an NSGAII-based generic optimizer applicable to various types of problems. The differences between Omni-optimizer and NSGA-II
are fourfold: the Latin hypercube sampling-based population
initialization, the so-called restricted mating selection, the dominance-based non-dominated sorting, and the alternative
crowding distance. In the restricted mating selection, an individual xa is randomly selected from the population. Then, xa
and its nearest neighbor xb in the solution space are comparedbased on their non-domination levels and crowding distancevalues. The winner among xa and xb is selected as a parent.
The crowding distance measure in Omni-optimizer takes
into account both the objective and solution spaces. For the ith individual xi in each non-dominated front R, the crowding
distance in the objective space cobj
is calculated in a similarimanner to NSGA-II. In contrast, the crowding distance valueof xi in the solution space csol
is calculated in a differentimanner. First, for each j ∈ {1, ..., D}, a “variable-wise”
crowding distance value of xi in the j-th decision variablecsoli,j is calculated as follows:
 

x
−xi,j

if xi,j = xmin
2 xi+1,jmax −xmin
j

j
  j

xi,j −xi−1,j
sol
,
(1)
ci,j = 2 xmax −xmin
else if xi,j = xmax
j

j
j


x
−x
i+1,j
i−1,j
 max minotherwise
x
−x
j

j

where we assume that all individuals in R are sorted based on
their j-th decision variable values in descending order. In (1),xmin
= minx∈R {xj } and xmax
= maxx∈R {xj }. Unlike thejjcrowding distance in the objective space, an infinitely large
value is not given to a boundary individual.
Then, an “individual-wise” crowding
distance value csolisi
PD solsol
calculated as follows: ci = ( j=1 ci,j )/D. The averagevalue csol
avg of all individual-wise crowding distance values is
P|R| solalso calculated as follows: csol
avg = (
i=1 ci )/|R|. Finally,
the crowding distance value ci of xi is obtained as follows:
(
sol
sol
sol
max{cobj
if cobj> cobjavg or ci > cavg
i , ci }
i
ci =
, (2)
obj solmin{ci , ci } otherwise
where cobj
avg is the average value of all crowding distance values
in the objective space. As shown in (2), ci in Omni-optimizer
is the combination of cobj
and csoli . Due to its alternativei
crowding distance, the results presented in [9] showed that
Omni-optimizer finds more diverse solutions than NSGA-II.
In addition to Omni-optimizer, two extensions of NSGAII for MMOPs have been proposed. DNEA [23] is similar toOmni-optimizer but uses two sharing functions in the objective
and solution spaces. DNEA requires fine-tuning of two sharing
niche parameters for the objective and solution spaces. The
secondary criterion of DN-NSGA-II [24] is based on thecrowding distance only in the solution space. DN-NSGA-II
uses a solution distance-based mating selection.
The following are other dominance-based MMEAs. An
MMEA proposed in [25] utilizes DBSCAN [26] and the rake
selection [27]. DBSCAN, which is a clustering method, isused for grouping individuals based on the distribution ofindividuals in the solution space. The rake selection, which isa reference vector-based selection method similar to NSGA-III[28], is applied to individuals belonging to each niche for theenvironmental selection. SPEA2+ [5], [15] uses two archives
Aobj and Asol to maintain diverse non-dominated individuals
in the objective and solution spaces, respectively. While theenvironmental selection in Aobj is based on the density of
individuals in the objective space similar to SPEA2 [13], thatin Asol is based on the density of individuals in the solution
space. For the mating selection in SPEA2+, neighborhoodindividuals in the objective space are selected only from Aobj .
PQ, -MOEA [4], 4D-Miner [3], [29], and MNCA [19] arecapable of handling dominated solutions for MMOPs with
δ > 0. PQ, -MOEA uses the -dominance relation [30] sothat an unbounded archive can maintain individuals with acceptable quality according to the decision maker. Unlike other
MMEAs, PQ, -MOEA does not have an explicit mechanism tomaintain the solution space diversity. 4D-Miner was speciallydesigned for functional brain imaging problems [3]. The
population is initialized by a problem-specific method. 4DMiner maintains dissimilar individuals in an external archive,
whose size is ten times larger than the population size. The
environmental selection in 4D-Miner is based on a problemspecific metric. Similar to DIOP [22] (explained later), MNCA
simultaneously evolves multiple subpopulations P 1 , ..., P S ,where S is the number of subpopulations. In MNCA, theprimary subpopulation P 1 aims to find an approximationof the Pareto front that provides a target front for othersubpopulations P 2 , ..., P S . While the update of P 1 is basedon the same selection mechanism as in NSGA-II, the updateof P 2 , ..., P S is performed with a complicated method thattakes into account both the objective and solution spaces.
Although the above-mentioned MMEAs use genetic variation operators (e.g., the SBX crossover and the polynomial
mutation [12]), the following MMEAs are based on otherapproaches. Niching-CMA [20] is an extension of CMAES [31] for MMOPs by introducing a niching mechanism.
The number of niches and the niche radius are adaptively
adjusted in Niching-CMA. An aggregate distance metric in
the objective and solution spaces is used to group individuals into multiple niches. For each niche, individuals with
better non-domination levels survive to the next iteration.
MO Ring PSO SCD [21], a PSO algorithm for MMOPs,
uses a diversity measure similar to Omni-optimizer. However,
MO Ring PSO SCD handles the boundary individuals in theobjective space in an alternative manner. In addition, an indexbased ring topology is used to create niches.
Two extensions of artificial immune systems [32] havebeen proposed for MMOPs: omni-aiNet [18] and cob-aiNet[33]. These two methods use a modified version of thepolynomial mutation [12]. The primary and secondary criteriaof omni-aiNet are based on -nondomination levels [30] anda grid operation, respectively. In addition, omni-aiNet uses
suppression and insertion operations. While the suppression
operation deletes an inferior individual, the insertion operation
adds new individuals to the population. The population size
is not constant due to these two operations. The primaryand secondary criteria of cob-aiNet are based on the fitness

4

TABLE I: Properties of 18 MMEAs. µ and nmax denote theSet

Decomp.

Dominance

population size and the maximum number of evaluations used ineach paper, respectively. “δ > 0” indicates whether each method canhandle MMOPs with δ > 0. “U” means whether each method has an
unbounded population/archive. Initial µ values are reported for omniaiNet, cob-aiNet, PQ, -MOEA, and MOEA/D-AD. µ and nmax usedin the post-processing step are shown for a method in [17].

P.

assignment method in SPEA2 [13] and a diversity measurewith a sharing function in the solution space, respectively. The
maximum population size is introduced in cob-aiNet.
2) Decomposition-based MMEAs: A three-phase multistart method is proposed in [16]. First, (1, λ)-ES is carriedout on each M objective functions K times to obtain M × K
best-so-far solutions. Then, an unsupervised clustering method
is applied to the M × K solutions to detect the number of
equivalent Pareto optimal solution subsets s. Finally, s runsof (1, λ)-ES are performed on each N single-objective subproblem decomposed by the Tchebycheff function. The initialindividual of each run is determined in a chained manner.
The best solution found in the j-th subproblem becomes aninitial individual of (1, λ)-ES for the j + 1-th subproblem(j ∈ {1, ..., N − 1}). It is expected that s equivalent solutions
are found for each N decomposed subproblems.
Two variants of MOEA/D [34] for MMOPs are proposedin [35], [36]. MOEA/D decomposes an M -objective problem
into N single-objective subproblems using a set of weight vectors, assigning a single individual to each subproblem. Then,
MOEA/D simultaneously evolves the N individuals. Unlike
MOEA/D, the following two methods assign one or moreindividuals to each subproblem to handle the equivalency.
The MOEA/D algorithm presented in [35] assigns K individuals to each subproblem. The selection is conducted basedon a fitness value combining the PBI function value [34]and two distance values in the solution space. K dissimilar
individuals are likely to be assigned to each subproblem.
The main drawback of the above methods [16], [35] is thedifficulty in setting a proper value for K, because it is problem
dependent. MOEA/D-AD [36] does not need such a parameterbut requires a relative neighborhood size L. For each iteration,a child u is assigned to the j-th subproblem whose weightvector is closest to f (u), with respect to the perpendiculardistance. Let X be a set of individuals already assigned to thejth-subproblem. If x in X is within the L nearest individuals
from the child u in the solution space, x and u are comparedbased on their scalarizing function values g(x) and g(u). If
g(u) ≤ g(x), x is deleted from the population and u enters
the population. u also enters the population when no x in X
is in the L neighborhood of u in the solution space.
3) Set-based MMEAs: DIOP [22] is a set-based MMEA
that can maintain dominated solutions in the population. In
the set-based optimization framework [37], a single solution
in the upper level represents a set of solutions in the lower
level (i.e., a problem). DIOP simultaneously evolves an archive
A and a target population T . While A approximates only the
Pareto front and is not shown to the decision maker, T obtains
diverse solutions with acceptable quality by maximizing thefollowing G indicator: G(T ) = wobj Dobj (T ) + wsol Dsol (T ).
Here, wobj + wsol = 1. Dobj is a performance indicator in
the objective space, and Dsol is a diversity measure in thesolution space. In [22], Dobj and Dsol were specified by thehypervolume indicator [38] and the Solow-Polasky diversity
measure [39], respectively. Meta-individuals in T that are dominated by any meta-individuals in A are excluded for the
calculation of the G metric. At the end of the search, T
is likely to contain meta-individuals (i.e., solution sets of a

MMEAs
SPEA2+ [5], [15]

Year
2004

µ
100

nmax
50 000Omni-optimizer [9], [14]

2005

1 000

500 000

4D-Miner [3], [29]

2005

200

8 000

omni-aiNet [18]

2006

400

40 000Niching-CMA [20]

2009

50

50 000A method in [25]2010

Not clearly reported

PQ, -MOEA [4]

2011

200

5 000

cob-aiNet [33]

2011

100

40 000

MNCA [19]

2013

100

100 000

DN-NSGA-II [24]

2016

800

80 000MO Ring PSO SCD [21]

2017

800

80 000DNEA [23]

2018

210

63 000A method in [16]

2007

10

20 000A method in [35]

2018

1 120

89 600MOEA/D-AD [36]

2018

100

30 000

DIOP [22]

2010

50

100 000

A method in [40]

2012

200

400 000

A method in [17]

2009

20

2 000

δ>0

U

X
X

X

X

X

X
X

problem) -nondominated by meta-individuals in A.
Another set-based MMEA is presented in [40]. Unlike
DIOP, the proposed method evolves only a single population.
Whereas DIOP maximizes the weighted sum of values of Dobj
and Dsol , the proposed method treats Dobj and Dsol as meta
two-objective functions. NSGA-II is used to simultaneously
maximize Dobj and Dsol in [40].
4) A post-processing approach: As pointed out in [17], it
is not always necessary to locate all Pareto optimal solutions.
Suppose that a set of non-dominated solutions A has alreadybeen obtained by an MOEA (e.g., NSGA-II) but not an MMEA
(e.g., Omni-optimizer). After the decision maker has selected
the final solution xfinal from A according to her/his preferencein the objective space, it is sufficient to search solutions whose
objective vectors are equivalent to f (xfinal ).
A post-processing approach is proposed in [17] to handle this problem. First, the proposed approach formulates ameta constrained two-objective minimization problem where
f1meta = kf (x) − f (xfinal )k2 , f2meta = −kx − xfinal k2 , andg meta (x) = f1meta (x) − θ < 0. The meta objective functions
f1meta and f2meta represent the distance between x and xfinal in
the objective and solution spaces. Thus, smaller f1meta (x) andf2meta (x) indicate that x is similar to xfinal in the objective
space and far from xfinal in the solution space, respectively.
The constraint g meta with θ > 0 prevents f2meta (x) frombecoming an infinitely small value in unbounded problems.
NSGA-II is used as a meta-optimizer in [17].
5) Open issues: Table I summarizes the properties of the
18 MMEAs reviewed in this section.
While some MMEAs require an extra parameter (e.g.,
L in MOEA/D-AD), Omni-optimizer does not require such

5

a parameter. This parameter-less property is an advantageof Omni-optimizer. However, Omni-optimizer is a Pareto
dominance-based MMEA. Since dominance-based MOEAs
perform poorly on most MOPs with more than three objectives[28], Omni-optimizer is unlikely to handle many objectives.
In addition to MMEAs, some MOEAs handling the solution
space diversity have been proposed, such as GDEA [41],DEMO [42], DIVA [43], “MMEA” [44], DCMMMOEA [45],and MOEA/D-EVSD [46]. Note that solution space diversity
management in these MOEAs aims to efficiently approximate
the Pareto front for MOPs. Since these methods were notdesigned for MMOPs, they are likely to perform poorly for
MMOPs. For example, “MMEA”, which stands for a modelbased multi-objective evolutionary algorithm, cannot find multiple equivalent Pareto optimal solutions [44]. Nevertheless,helpful clues for designing an efficient MMEA can be foundin these MOEAs.
The performance of MMEAs has not been well analyzed.
The post-processing method may perform better than MMEAs
when the objective functions of a real-world problem arecomputationally expensive. However, an in-depth investigationis necessary to determine which approach is more practical.
Whereas the population size µ and the maximum number of
evaluations nmax were set to large values in some studies,
they were set to small values in other studies. For example,
Table I shows that µ = 1 000 and nmax = 500 000 for
Omni-optimizer, while µ = 50 and nmax = 50 000 for
Niching-CMA. It is unclear whether an MMEA designed withlarge µ and nmax values works well with small µ and nmax
values. While MMOPs with four or more objectives appear
in real-world applications (e.g., five-objective rocket engine
design problems [7]), most MMEAs have been applied to only
two-objective MMOPs. A large-scale benchmarking study isnecessary to address the above-mentioned issues.
The decision maker may want to examine diverse dominatedsolutions. As explained in Section I, dominated solutions
found by PQ, -MOEA support the decision making in space
mission design problems [4]. The results presented in [29]showed that diverse solutions found by 4D-Miner help neuroscientists analyze brain imaging data. Although most MMEAs
assume MMOPs with δ = 0 as shown in Table I, MMEAs
that can handle MMOPs with δ > 0 may be more practical.
Since most MMEAs (e.g., Omni-optimizer) remove dominatedindividuals from the population, they are unlikely to finddiverse dominated solutions. Some specific mechanisms arenecessary to handle MMOPs with δ > 0 (e.g., the multiple
subpopulation scheme in DIOP and MNCA).
As explained at the beginning of this section, MMEAs
need the three abilities (1)–(3). While the abilities (1) and(2) are needed to approximate the Pareto front, the ability(3) is needed to find equivalent Pareto optimal solutions.
Most existing studies (e.g., [9], [20], [21], [36]) report thatthe abilities (1) and (2) of MMEAs are worse than those of
MOEAs. For example, the results presented in [36] showedthat Omni-optimizer, MO Ring PSO SCD, and MOEA/DAD perform worse than NSGA-II in terms of IGD [47](explained in Section V). If the decision maker is not interested
in the distribution of solutions in the solution space, it would

be better to use MOEAs rather than MMEAs. The poor performance of MMEAs for multi-objective optimization is mainly
due to the ability (3), which prevents MMEAs from directly
approximating the Pareto front. This undesirable performance
regarding the abilities (1) and (2) is an issue in MMEAs.
• What to learn from MSOPs: An online data repository(https://github.com/mikeagn/CEC2013) that provides results of
optimizers on the CEC2013 problem suite [48] is available for
MSOPs. This repository makes the comparison of optimizers
easy, facilitating constructive algorithm development. A similar data repository is needed for studies of MMOPs.
The number of maintainable individuals in the population/archive strongly depends on the population/archive size.
However, it is usually impossible to know the number of
equivalent Pareto optimal solutions of an MMOP a priori. The
same issue can be found in MSOPs. To address this issue, thelatest optimizers (e.g., dADE [49] and RS-CMSA [50]) havean unbounded archive that maintains solutions found during
the search process. Unlike modern optimizers for MSOPs,
Table I shows that only three MMEAs have such a mechanism.
The adaptive population sizing mechanisms in omni-aiNet,
PQ, -MOEA, and MOEA/D-AD are advantageous. A general
strategy of using an unbounded (external) archive could improve the performance of MMEAs.
IV. M ULTI - MODAL MULTI - OBJECTIVE TEST PROBLEMSThis section describes test problems for benchmarking
MMEAs. Unlike multi-objective test problems (e.g., the DTLZ
[51] test suite), multi-modal multi-objective test problems
were explicitly designed such that they have multiple equivalent Pareto optimal solution subsets. The two-objective andtwo-variable SYM-PART1 [16] is one of the most representative test problems for benchmarking MMEAs: f1 (y) =(y1 +a)2 +y22 and f2 (y) = (y1 −a)2 +y22 . Here, y1 and y2 aretranslated values of x1 and x2 as follows: y1 = x1 −t1 (c+2a)and y2 = x2 − t2 b. In SYM-PART1, a controls the region of
Pareto optimal solutions, and b and c specify the positionsof the Pareto optimal solution subsets. The so-called tileidentifiers t1 and t2 are randomly selected from {−1, 0, 1}.
Fig. 3(a) shows the shape of the Pareto optimal solutions of
SYM-PART1 with a = 1, b = 10, and c = 8. As shown in
Fig. 3(a), the equivalent Pareto optimal solution subsets areon nine lines in SYM-PART1.
Other test problems include the Two-On-One [11] problem,
the Omni-test problem [9], the SYM-PART2 and SYM-PART3
problems [16], the Superspheres problem [52], the EBN
problem [53], the two SSUF problems [24], and the Polygon
problems [54]. Fig. 3 also shows the distribution of their Pareto
optimal solutions. Since there are an infinite number of Pareto
optimal solutions in the EBN problem, we do not show them.
Source codes of the ten problems can be downloaded from the
supplementary website (https://sites.google.com/view/emmo/).
In Omni-test, equivalent Pareto optimal solution subsets are
regularly located. SYM-PART2 is a rotated version of SYMPART1. SYM-PART3 is a transformed version of SYMPART2 using a distortion operation. The Superspheres problem with D = 2 has six equivalent Pareto optimal solution

6

0

0

0

0
x1

−15
−15

15

where M , D, and P denote the number of objectives, design
variables, and equivalent Pareto optimal solution subsets, respectively.
If a problem has irregularity, the shapes of its multiple equivalent
Pareto optimal solution subsets differ from each other.

−15
−8 −4

15

0
x1

4

(b) SYM-PART2(c) SYM-PART3

2
1
0
−1
−2
−2 −1

6
5
4
3
2
1
0
0 1 2 3 4 5 6
x1

5
4
3
2
1
0
0

0
x1

1

2

(e) Omni-test

x2

0

2
x1

(g) SSUF1

3

15
10
5
0
0

2

4

x1

6

(h) SSUF3

8 10
×10

π/4
x1

8

π/2

(f) Superspheres−1
20 ×10

1

−1
1

x2

(a) SYM-PART1(d) Two-On-One

x2

0
x1

x2

x2

−15
−15

x2

15

x2

15

x2

x2

TABLE II: Properties of multi-modal multi-objective test problems,
15

−1

10
8
6
4
2
0
0

2

4

x1

6

8 10

(i) Polygon

Fig. 3: Distribution of the Pareto optimal solutions for the eightproblems. Only x1 and x2 are shown on Omni-test.

subsets. However, the number of its P is unknown for D > 2.
EBN can be considered as a real-coded version of the so-calledbinary one-zero max problem. All solutions in the solution
space are Pareto optimal solutions. SSUF1 and SSUF3 areextensions of the UF problems [55] to MMOPs. There are
two symmetrical Pareto optimal solution subsets in SSUF1 and
SSUF3. Polygon is an extension of the distance minimization
problems [56] to MMOPs, where P equivalent Pareto optimal
solution subsets are inside of P regular M -sided polygons.
In addition, the eight MMF problems are presented in [21].
Similar to SSUF1 and SSUF3, the MMF problems are derived
from the idea of designing a problem that has multiple equivalent Pareto optimal solution subsets by mirroring the originalone. A bottom-up framework for generating scalable test
problems with any D is proposed in [57]. P equivalent Pareto
optimal solution subsets are in P hyper-rectangular located in
the solution space similar to the SYM-PART problems. While
the first k variables play the role of “position” parametersin the solution space, the other D − k variables represent“distance” parameters. The six HPS problem instances wereconstructed using this framework in [57].
If a given problem has the multi-modal fitness landscape, it
may have multiple non-Pareto fronts whose shapes are similar
to the true Pareto front. Such a problem (e.g., ZDT4 [58]) is
referred to as a multi-frontal test problem [59]. If the δ value(defined in Subsection II-2) is sufficiently large, a multi-frontal
test problem can be regarded as a multi-modal multi-objective
test problem. In fact, ZDT4 was used in [19] as a test problem.
The Kursawe problem [60] is a multi-modal and nonseparable
test problem with a disconnected Pareto front. The Kursawe
problem has two fronts in the objective space similar to multi-

Test problems

M

D

P

Irregularity

SYM-PART problems [16]

2

2

9

X

Two-On-One problem [11]

2

2

2

Omni-test problem [9]

2

Any

3D

Superspheres problem [52]2AnyUnknown

EBN problem [53]

2

Any

∞

Polygon problems [54]

Any

2

Any

SSUF problems [24]

2

2

2

MMF suite [21]222 or 4

HPS suite [57]

2

Any

Any

frontal problems. Thus, the Kursawe problem can be used asa multi-modal multi-objective test problem.
1) Open issues: Table II summarizes the properties ofmulti-modal multi-objective test problems reviewed here. In
Table II, P of Omni-test adheres to [22].
Table II indicates that scalable test problems do not exist,
in terms of M , D, and P . Although the SYM-PART problems
have some desirable properties (e.g., their adjustable andstraightforward Pareto optimal solution shapes), M , D, andP are constant in these problems. Only Polygon is scalable in
M . While most test problems have only two design variables,
Omni-test and HPS are scalable in D. Unfortunately, P
increases exponentially with increased D in Omni-test due
to the combinatorial nature of variables. Although the ideaof designing scalable SYM-PART and Polygon problems to
D is presented in [61], [62], they have similar issues toOmni-test. Although the HPS problems do not have suchan issue, it is questionable whether there exists a real-world
problem with design variables affecting only the distance
between the objective vectors and the Pareto front. Only SYMPART3 has irregularity. Since the shapes of the Pareto optimal
solution subsets may be different from each other in real-world
problems, we believe that test problems with the irregularity
are necessary to evaluate the performance of MMEAs. The
performance of an MMEA with an absolutely defined nichingradius (e.g., DNEA) is likely to be overestimated in test
problems without irregularity.
In addition, the relation between synthetic test problems
and real-world problems has not been discussed. The idea ofdesigning a Polygon problem based on a real-world map ispresented in [63]. However, this does not mean that such a
Polygon problem is an actual real-world problem.
• What to learn from MSOPs: Some construction methods
for multi-modal single-objective test problems are available,
such as the software framework proposed in [64], the construction method for various problems [65], and Ahrari andDeb’s method [66]. Borrowing ideas from such sophisticatedconstruction methods is a promising way to address theabove-mentioned issues of multi-modal multi-objective test

7

(4)

where A is a set of solutions obtained by an MMEA and A∗ isa set of reference solutions in the Pareto optimal solution set.
ED(x1 , x2 ) denotes the Euclidean distance between x1 andx2 . While A with a small IGD value is a good approximationof the Pareto front, A with a small IGDX approximates Pareto
optimal solutions well. Other indicators in the first category
include GDX [4], the Hausdorff distance indicator [67] in thesolution space [4], CR [21], and PSP [21]. GDX is a GD
indicator in the solution space similar to IGDX. CR is analternative version of the maximum spread [38] to measure
the spread of A. PSP is a combination of IGDX and CR.
Performance indicators in the second category include themean of the pairwise distance between two solutions [20], CS[16], SPS [16], the Solow-Polasky diversity measure [39] usedin [22], [40], and PSV [57]. CS is the number of Pareto optimal
solution subsets covered by at least one individual. SPS is thestandard deviation of the number of solutions close to each
Pareto optimal solution subset. PSV is the percentage of thevolume of A in the volume of A∗ in the solution space.
1) Open issues: Table III shows the properties of performance indicators for MMEAs reviewed in this section, where
the properties are assessed based on the description of eachindicator. While the properties of the performance indicators
for MOEAs have been examined (e.g., [38], [67]), those for
MMEAs have not been well analyzed.
Performance indicators for MMEAs should be able toevaluate the three abilities (1)–(3) explained in Section III.
Although IGDX is frequently used, it should be noted that
IGDX does not evaluate the distribution of solutions in theobjective space. Fig. 4 shows the distribution of two solution sets A1 and A2 for SYM-PART1 in the solution and

Div.

Unif.

IGDX [4], [44]

X

X

X

X

X

Hausdorff distance [4]

X

X

X

X

X

X

X

X

X

X

X

CR [21]X

PSP [21]

X

X

X

X

SPS [16]

X

X

X

X

X

Solow-Polasky [39]X

PSV [57]

1500(a) A1 in the solution space

X

X

X

X

X

15

15

Dif.

X

X

0
x1

Ref.
X

X

CS [16]

−15
−15

Spr.

X

Pairwise distance [20]

x2

(3)

Conv.
X

−15
−15

X

X

0
x1

15

(b) A2 in the solution space

4

4

3

3

2

2

f2



n
o
1  X
min ED f (x), f (z)  ,IGD(A) =|A∗ |∗ x∈ A
z ∈A


n
o
1  X
IGDX(A) =min ED x, z  ,|A∗ |∗ x∈ Az ∈A

Indicators
GDX [4]x2

V. P ERFORMANCE INDICATORS FOR MMEA S
Performance indicators play an important role in quantitatively evaluating the performance of MOEAs as well as
MMEAs. Since performance indicators for MOEAs consider
only the distribution of objective vectors (e.g., the hypervolume, GD, and IGD indicators [38], [47]), they cannot be used
to assess the ability of MMEAs to find multiple equivalent
Pareto optimal solutions. For this reason, some indicators havebeen specially designed for MMEAs. Performance indicators
for MMEAs can be classified into two categories: simple
extensions of existing performance indicators for MOEAs andspecific indicators based on the distributions of solutions.
IGDX [4], [44] is a representative example of the first
approach. The IGD and IGDX indicators are given as follows:

TABLE III: Properties of performance indicators for MMEAs(convergence to Pareto optimal solution subsets, diversity,
uniformity, spread, the use of reference solution sets, andpossibility to compare solution sets with different sizes).

f2

problems. In [64], Rönkkönen et al. present eight desirable
properties for multi-modal single-objective problem generatorssuch as scalability in D, control of the number of global andlocal optima, and regular and irregular distributions of optima.
These eight properties can be a useful guideline for designing
multi-modal multi-objective problem generators.

1

1

0

0
0

1

2

f1

3

4

(c) A1 in the objective space

0

1

2

f1

3

4

(d) A2 in the objective space

Fig. 4: Comparison of solution sets A1 and A2 for SYM-PART1.

objective spaces, where |A1 | and |A2 | are 27. While thesolutions in A1 are evenly distributed on one of the nine
Pareto optimal solution subsets, the solutions in A2 are evenly
distributed on all of them. Although A1 has 27 objective
vectors that cover the Pareto front, A2 has only 3 equivalent
objective vectors. The IGDX and IGD values of A1 and A2
are as follows: IGDX(A1 ) = 15.92, IGDX(A2 ) = 0.25,
IGD(A1 ) = 0.06, and IGD(A2 ) = 0.81. We used 5 000
Pareto optimal solutions for A∗ . Although A2 has a worse
distribution in the objective space than A1 , IGDX(A2 ) issignificantly better than IGDX(A1 ). As demonstrated here,
IGDX can evaluate the abilities (1) and (3) but cannot evaluate
the ability (2) to find diverse solutions in the objective space.
Since the other indicators in Table III do not take into accountthe distribution of objective vectors similar to IGDX, they
are likely to have the same undesirable property. For a fair
performance comparison, it is desirable to use the indicators

8

for MOEAs (e.g., hypervolume and IGD) in addition to theindicators for MMEAs in Table III.
• What to learn from MSOPs: It is desirable that the indicatorsfor multi-modal single-objective optimizers evaluate a solution
set without the knowledge of the fitness landscape such as the
positions of the optima and the objective values of the optima[68]. The same is true for indicators for MMEAs. Table III
shows that most indicators (e.g., IGDX) require A∗ . Since A∗
is usually unavailable in real-world problems, it is desirable
that indicators for MMEAs evaluate A without A∗ .
Since the archive size in modern multi-modal singleobjective optimizers is unbounded in order to store a numberof local optima [10], most indicators in this field can handle
solution sets with different sizes (e.g., the peak ratio and thesuccess rate [48]). For the same reason, it is desirable thatindicators for MMEAs evaluate solution sets with different
sizes in a fair manner. However, it is difficult to directly
use indicators for multi-modal single-objective optimizers toevaluate MMEAs.
VI. C ONCLUSION
The contributions of this paper are threefold. The firstcontribution is that we reviewed studies in this field in termsof definitions of MMOPs, MMEAs, test problems, and performance indicators. It was difficult to survey the existing studies
of MMOPs for the reasons described in Section I. Our review
helps to elucidate the current progress on evolutionary multimodal multi-objective optimization. The second contribution
is that we clarified open issues in this field. In contrast to
multi-modal single-objective optimization, multi-modal multiobjective optimization has not received much attention despite
its practical importance. Thus, some critical issues remain.
The third contribution is that we pointed out an issue associated with performance indicators for MMEAs. Reliable
performance indicators are necessary for the advancement of
MMEAs. We hope that this paper will encourage researchers
to work in this research area, which is not well explored.
ACKNOWLEDGMENT
This work was supported by the Program for Guangdong Introducing Innovative and Enterpreneurial Teams(Grant No. 2017ZT07X386), Shenzhen Peacock Plan (Grant
No. KQTD2016112514355531), the Science and Technology Innovation Committee Foundation of Shenzhen (Grant
No. ZDSYS201703031748284), the Program for University Key Laboratory of Guangdong Province (Grant No.
2017KSYS008), and National Natural Science Foundation ofChina (Grant No. 61876075).
R EFERENCES[1] K. Deb, Multi-Objective Optimization Using Evolutionary Algorithms.
John Wiley & Sons, 2001.
[2] K. Miettinen, Nonlinear Multiobjective Optimization. Springer, 1998.
[3] M. Sebag, N. Tarrisson, O. Teytaud, J. Lefèvre, and S. Baillet, “A
Multi-Objective Multi-Modal Optimization Approach for Mining Stable
Spatio-Temporal Patterns,” in IJCAI, 2005, pp. 859–864.
[4] O. Schütze, M. Vasile, and C. A. C. Coello, “Computing the Set of
Epsilon-Efficient Solutions in Multiobjective Space Mission Design,”
JACIC, vol. 8, no. 3, pp. 53–70, 2011.

[5] T. Hiroyasu, S. Nakayama, and M. Miki, “Comparison study of SPEA2+,
SPEA2, and NSGA-II in diesel engine emissions and fuel economy
problem,” in IEEE CEC, 2005, pp. 236–242.
[6] M. Preuss, C. Kausch, C. Bouvy, and F. Henrich, “Decision Space
Diversity Can Be Essential for Solving Multiobjective Real-World
Problems,” in MCDM, 2008, pp. 367–377.
[7] F. Kudo, T. Yoshikawa, and T. Furuhashi, “A study on analysis of design
variables in Pareto solutions for conceptual design optimization problemof hybrid rocket engine,” in IEEE CEC, 2011, pp. 2558–2562.
[8] J. Togelius, M. Preuss, and G. N. Yannakakis, “Towards multiobjective
procedural map generation,” in PCGames, 2010.
[9] K. Deb and S. Tiwari, “Omni-optimizer: A generic evolutionary algorithm for single and multi-objective optimization,” EJOR, vol. 185, no. 3,
pp. 1062–1087, 2008.
[10] X. Li, M. G. Epitropakis, K. Deb, and A. P. Engelbrecht, “Seeking
Multiple Solutions: An Updated Survey on Niching Methods and Their
Applications,” IEEE TEVC, vol. 21, no. 4, pp. 518–538, 2017.
[11] M. Preuss, B. Naujoks, and G. Rudolph, “Pareto Set and EMOA
Behavior for Simple Multimodal Multiobjective Functions,” in PPSN,
2006, pp. 513–522.
[12] K. Deb, S. Agrawal, A. Pratap, and T. Meyarivan, “A fast and elitist
multiobjective genetic algorithm: NSGA-II,” IEEE TEVC, vol. 6, no. 2,
pp. 182–197, 2002.
[13] E. Zitzler, M. Laumanns, and L. Thiele, “SPEA2: Improving the Strength
Pareto Evolutionary Algorithm,” ETHZ, Tech. Rep., 2001.
[14] K. Deb and S. Tiwari, “Omni-optimizer: A Procedure for Single andMulti-objective Optimization,” in EMO, 2005, pp. 47–61.
[15] M. Kim, T. Hiroyasu, M. Miki, and S. Watanabe, “SPEA2+: Improving
the Performance of the Strength Pareto Evolutionary Algorithm 2,” in
PPSN, 2004, pp. 742–751.
[16] G. Rudolph, B. Naujoks, and M. Preuss, “Capabilities of EMOA to
Detect and Preserve Equivalent Pareto Subsets,” in EMO, 2007, pp. 36–
50.
[17] G. Rudolph and M. Preuss, “A multiobjective approach for finding equivalent inverse images of pareto-optimal objective vectors,” in MCDM,
2009, pp. 74–79.
[18] G. P. Coelho and F. J. V. Zuben, “omni-aiNet: An Immune-Inspired
Approach for Omni Optimization,” in ICARIS, 2006, pp. 294–308.
[19] E. M. Zechman, M. H. G., and M. E. Shafiee, “An evolutionary
algorithm approach to generate distinct sets of non-dominated solutions
for wicked problems,” Eng. Appl. of AI, vol. 26, no. 5-6, pp. 1442–1457,
2013.
[20] O. M. Shir, M. Preuss, B. Naujoks, and M. T. M. Emmerich, “Enhancing
Decision Space Diversity in Evolutionary Multiobjective Algorithms,” in
EMO, 2009, pp. 95–109.
[21] C. Yue, B. Qu, and J. Liang, “A Multi-objective Particle Swarm
Optimizer Using Ring Topology for Solving Multimodal Multi-objective
Problems,” IEEE TEVC, 2018 (in press).
[22] T. Ulrich, J. Bader, and L. Thiele, “Defining and Optimizing IndicatorBased Diversity Measures in Multiobjective Search,” in PPSN, 2010,
pp. 707–717.
[23] Y. Liu, H. Ishibuchi, Y. Nojima, N. Masuyama, and K. Shang, “A
Double-Niched Evolutionary Algorithm and Its Behavior on PolygonBased Problems,” in PPSN, 2018, pp. 262–273.
[24] J. J. Liang, C. T. Yue, and B. Y. Qu, “Multimodal multi-objective
optimization: A preliminary study,” in IEEE CEC, 2016, pp. 2454–2461.
[25] O. Kramer and H. Danielsiek, “DBSCAN-based multi-objective niching
to approximate equivalent pareto-subsets,” in GECCO, 2010, pp. 503–
510.
[26] M. Ester, H. Kriegel, J. Sander, and X. Xu, “A Density-Based Algorithm
for Discovering Clusters in Large Spatial Databases with Noise,” in
KDD, 1996, pp. 226–231.
[27] O. Kramer and P. Koch, “Rake Selection: A Novel Evolutionary MultiObjective Optimization Algorithm,” in KI, 2009, pp. 177–184.
[28] K. Deb and H. Jain, “An evolutionary many-objective optimization
algorithm using reference-point-based nondominated sorting approach,part I: solving problems with box constraints,” IEEE TEVC, vol. 18,
no. 4, pp. 577–601, 2014.
[29] V. Krmicek and M. Sebag, “Functional Brain Imaging with Multiobjective Multi-modal Evolutionary Optimization,” in PPSN, 2006, pp.
382–391.
[30] M. Laumanns, L. Thiele, K. Deb, and E. Zitzler, “Combining Convergence and Diversity in Evolutionary Multiobjective Optimization,” Evol.
Comput., vol. 10, no. 3, pp. 263–282, 2002.
[31] N. Hansen and A. Ostermeier, “Completely derandomized selfadaptation in evolution strategies,” Evol. Comput., vol. 9, no. 2, pp.
159–195, 2001.

9[32] D. Dasgupta, S. Yu, and F. Niño, “Recent Advances in Artificial Immune
Systems: Models and Applications,” Appl. Soft Comput., vol. 11, no. 2,
pp. 1574–1587, 2011.
[33] G. P. Coelho and F. J. V. Zuben, “A Concentration-Based Artificial
Immune Network for Multi-objective Optimization,” in EMO, 2011, pp.
343–357.
[34] Q. Zhang and H. Li, “MOEA/D: A multiobjective evolutionary algorithm
based on decomposition,” IEEE TEVC, vol. 11, no. 6, pp. 712–731, 2007.
[35] C. Hu and H. Ishibuchi, “Incorporation of a decision space diversity
maintenance mechanism into MOEA/D for multi-modal multi-objective
optimization,” in GECCO (Companion), 2018, pp. 1898–1901.
[36] R. Tanabe and H. Ishibuchi, “A Decomposition-Based Evolutionary
Algorithm for Multi-modal Multi-objective Optimization,” in PPSN,
2018, pp. 249–261.
[37] E. Zitzler, L. Thiele, and J. Bader, “On Set-Based Multiobjective
Optimization,” IEEE TEVC, vol. 14, no. 1, pp. 58–79, 2010.
[38] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, and V. G. da Fonseca, “Performance assessment of multiobjective optimizers: an analysis
and review,” IEEE TEVC, vol. 7, no. 2, pp. 117–132, 2003.
[39] A. R. Solow and S. Polasky, “Measuring biological diversity,” Environ.
Ecol. Stat., vol. 1, no. 2, pp. 95–103, 1994.
[40] H. Ishibuchi, M. Yamane, N. Akedo, and Y. Nojima, “Two-objective
solution set optimization to maximize hypervolume and decision spacediversity in multiobjective optimization,” in SCIS, 2012, pp. 1871–1876.
[41] A. Toffolo and E. Benini, “Genetic Diversity as an Objective in MultiObjective Evolutionary Algorithms,” Evol. Comput., vol. 11, no. 2, pp.
151–167, 2003.
[42] T. Robič and B. Filipič, “DEMO: differential evolution for multiobjective
optimization,” in EMO, 2005, pp. 520–533.
[43] T. Ulrich, J. Bader, and E. Zitzler, “Integrating decision space diversity
into hypervolume-based multiobjective search,” in GECCO, 2010, pp.
455–462.
[44] A. Zhou, Q. Zhang, and Y. Jin, “Approximating the Set of ParetoOptimal Solutions in Both the Decision and Objective Spaces by an
Estimation of Distribution Algorithm,” IEEE TEVC, vol. 13, no. 5, pp.
1167–1189, 2009.
[45] H. Xia, J. Zhuang, and D. Yu, “Combining Crowding Estimation in
Objective and Decision Space With Multiple Selection and Search
Strategies for Multi-Objective Evolutionary Optimization,” IEEE Trans.
Cyber., vol. 44, no. 3, pp. 378–393, 2014.
[46] J. C. Castillo, C. Segura, A. H. Aguirre, G. Miranda, and C. León,“A multi-objective decomposition-based evolutionary algorithm withenhanced variable space diversity control,” in GECCO (Companion),
2017, pp. 1565–1571.
[47] C. A. C. Coello and M. R. Sierra, “A Study of the Parallelization ofa Coevolutionary Multi-objective Evolutionary Algorithm,” in MICAI,
2004, pp. 688–697.
[48] X. Li, A. Engelbrecht, and M. G. Epitropakis, “Benchmark Functions
for CEC’2013 Special Session and Competition on Niching Methods forMultimodal Function Optimization,” RMIT Univ., Tech. Rep., 2013.
[49] M. G. Epitropakis, X. Li, and E. K. Burke, “A dynamic archive nichingdifferential evolution algorithm for multimodal optimization,” in IEEE
CEC, 2013, pp. 79–86.
[50] A. Ahrari, K. Deb, and M. Preuss, “Multimodal Optimization by
Covariance Matrix Self-Adaptation Evolution Strategy with Repelling
Subpopulations,” Evol. Comput., vol. 25, no. 3, pp. 439–471, 2017.
[51] K. Deb, L. Thiele, M. Laumanns, and E. Zitzler, “Scalable Test Problems for Evolutionary Multi-Objective Optimization,” in Evolutionary
Multiobjective Optimization. Theoretical Advances and Applications.
Springer, 2005, pp. 105–145.
[52] M. T. M. Emmerich and A. H. Deutz, “Test problems based on lamé
superspheres,” in EMO, 2006, pp. 922–936.
[53] N. Beume, B. Naujoks, and M. T. M. Emmerich, “SMS-EMOA:
multiobjective selection based on dominated hypervolume,” EJOR, vol.
181, no. 3, pp. 1653–1669, 2007.
[54] H. Ishibuchi, Y. Hitotsuyanagi, N. Tsukamoto, and Y. Nojima, “ManyObjective Test Problems to Visually Examine the Behavior of Multiobjective Evolution in a Decision Space,” in PPSN, 2010, pp. 91–100.
[55] Q. Zhang, A. Zhou, S. Zhao, P. N. Suganthan, W. Liu, and S. Tiwari,“Multiobjective optimization Test Instances for the CEC 2009 Special
Session and Competition,” Univ. of Essex, Tech. Rep., 2008.
[56] M. Köppen and K. Yoshida, “Substitute Distance Assignments in NSGAII for Handling Many-objective Optimization Problems,” in EMO, 2007,
pp. 727–741.
[57] B. Zhang, K. Shafi, and H. A. Abbass, “On Benchmark Problems andMetrics for Decision Space Performance Analysis in Multi-Objective
Optimization,” IJCIA, vol. 16, no. 1, pp. 1–18, 2017.

[58] E. Zitzler, K. Deb, and L. Thiele, “Comparison of Multiobjective
Evolutionary Algorithms: Empirical Results,” Evol. Comput., vol. 8,
no. 2, pp. 173–195, 2000. [Online]. Available: http://dx.doi.org/10.
1162/106365600568202[59] S. Huband, P. Hingston, L. Barone, and R. L. While, “A review ofmultiobjective test problems and a scalable test problem toolkit,” IEEE
TEVC, vol. 10, no. 5, pp. 477–506, 2006.
[60] F. Kursawe, “A Variant of Evolution Strategies for Vector Optimization,”in PPSN, 1990, pp. 193–197.
[61] V. L. Huang, A. K. Qin, K. Deb, E. Zitzler, P. N. Suganthan, J. J.
Liang, M. Preuss, and S. Huband, “Problem Definitions for Performance
Assessment on Multi-objective Optimization Algorithms,” NTU, Tech.
Rep., 2007.
[62] H. Ishibuchi, M. Yamane, N. Akedo, and Y. Nojima, “Many-objective
and many-variable test problems for visual examination of multiobjective
search,” in IEEE CEC, 2013, pp. 1491–1498.
[63] H. Ishibuchi, N. Akedo, and Y. Nojima, “A many-objective test problem
for visually examining diversity maintenance behavior in a decision
space,” in GECCO, 2011, pp. 649–656.
[64] J. Rönkkönen, X. Li, V. Kyrki, and J. Lampinen, “A framework forgenerating tunable test functions for multimodal optimization,” Soft
Comput., vol. 15, no. 9, pp. 1689–1706, 2011.
[65] B. Y. Qu, J. J. Liang, Z. Y. Wang, Q. Chen, and P. N. Suganthan,“Novel benchmark functions for continuous multimodal optimization
with comparative results,” SWEVO, vol. 26, pp. 23–34, 2016.
[66] A. Ahrari and K. Deb, “A Novel Class of Test Problems for Performance
Evaluation of Niching Methods,” IEEE TEVC, vol. 22, no. 6, pp. 909–
919, 2018.
[67] O. Schütze, X. Esquivel, A. Lara, and C. A. C. Coello, “Using the
Averaged Hausdorff Distance as a Performance Measure in Evolutionary
Multiobjective Optimization,” IEEE TEVC, vol. 16, no. 4, pp. 504–522,
2012.
[68] J. Mwaura, A. P. Engelbrecht, and F. V. Nepocumeno, “Performancemeasures for niching algorithms,” in IEEE CEC, 2016, pp. 4775–4784.