A Niching Indicator-Based Multi-modal Many-objective
Optimizer
Ryoji Tanabe, Hisao Ishibuchi∗
Shenzhen Key Laboratory of Computational Intelligence, University Key Laboratory of
Evolving Intelligent Systems of Guangdong Province, Department of Computer Science andEngineering, Southern University of Science and Technology, Shenzhen 518055, China

Abstract
Multi-modal multi-objective optimization is to locate (almost) equivalent
Pareto optimal solutions as many as possible. Some evolutionary algorithms for
multi-modal multi-objective optimization have been proposed in the literature.
However, there is no efficient method for multi-modal many-objective optimization, where the number of objectives is more than three. To address this issue,
this paper proposes a niching indicator-based multi-modal multi- and manyobjective optimization algorithm. In the proposed method, the fitness calculation is performed among a child and its closest individuals in the solution space
to maintain the diversity. The performance of the proposed method is evaluatedon multi-modal multi-objective test problems with up to 15 objectives. Results
show that the proposed method can handle a large number of objectives andfind a good approximation of multiple equivalent Pareto optimal solutions. The
results also show that the proposed method performs significantly better thaneight multi-objective evolutionary algorithms.
Keywords: Multi-modal multi-objective optimization, many-objective
optimization, indicator-based evolutionary algorithms, niching methods

1. Introduction
A bound-constrained multi-objective optimization problem (MOP) can beformulated as follows:
T
Minimize f (x) = f1 (x), ..., fM (x)(1)
subject to x ∈ S ⊆ RD ,∗ Corresponding author

Email addresses: rt.ryoji.tanabe@gmail.com (Ryoji Tanabe), hisao@sustc.edu.cn(Hisao Ishibuchi)

Preprint submitted to Swarm and Evolutionary Computation

June 2, 2019

where f : S → RM is an objective function vector that consists of M potentially conflicting objective functions, and RM is the objective space. Here,minmax
x = (x1 , ..., xD )T is a D-dimensional solution vector, and S = ΠD]
j=1 [xj , xj
max
is the bound-constrained solution space where xmin
≤
x
≤
x
for
each
index
j
j
j
j ∈ {1, ..., D}.
We say that x1 dominates x2 iff fi (x1 ) ≤ fi (x2 ) for all i ∈ {1, ..., M } andfi (x1 ) < fi (x2 ) for at least one index i. Here, x∗ is a Pareto-optimal solution if
there exists no x ∈ S such that x dominates x∗ . In this case, f (x∗ ) is a Paretooptimal objective vector. The set of all x∗ in S is the Pareto-optimal solution
set (PS), and the set of all f (x∗ ) is the Pareto front (PF). In general, the goalof MOPs is to find a set of non-dominated solutions that are well-distributedand close to the PF in the objective space.
A multi-objective evolutionary algorithm (MOEA) is an efficient approachfor solving MOPs [1]. However, representative Pareto dominance-based MOEAs(e.g., NSGA-II [2]) do not work well on MOPs with a large number of objectives[3], where MOPs with four or more objectives are referred to as many-objective
optimization problems (MaOPs). This is because most individuals in the population are non-dominated with each other in MaOPs, resulting a weak selection
pressure to the PF [3]. Since MaOPs appear in real-world problems [4], thepoor performance of MOEAs on MaOPs is critical. For this reason, a numberof MOEAs for MaOPs have been proposed in the literature [4].
Some recent studies (e.g., [3, 5]) show that indicator-based MOEAs perform well on MaOPs. In this framework, a so-called fitness value is assigned toeach individual in the population using quality indicators [6]. Then, the mating and/or environmental selections are performed based on the fitness values.
Representative indicator-based MOEAs include IBEA [7], SMS-EMOA [8], andHypE [9].
There are multiple Pareto optimal solutions that have (almost) the same
objective vector in some real-world problems, such as diet design problems [10],space mission design problems [11], rocket engine design problems [12], andfunctional brain imaging problems [13]. Fig. 1 explains a situation where threesolutions are identical or close to each other in the objective space but far
from each other in the solution space. On the one hand, the performance of
MOEAs is generally evaluated based only on the distribution of non-dominatedsolutions obtained in the objective space [6]. Thus, the distribution of nondominated solutions in the solution space is usually not taken into account toevaluate the performance of MOEAs [14]. On the other hand, diverse nondominated solutions in the solution space are helpful for decision-making inpractice [11, 10, 15, 16]. If multiple non-dominated solutions with (almost) thesame objective vector are obtained, users can make a final decision according
to their preference which cannot be represented by the objective functions. For
example, in Fig. 1, if xa becomes unavailable for some reasons (e.g., materials
shortages and traffic accidents), xb and xc can be considered as candidates forthe final solution instead of xa . Diverse solutions are also helpful to figure out
properties of problems [13, 17].

2

Objective space

Solution space

Figure 1: Illustration of a situation where three solutions (xa , xb , and xc ) are identicalor close to each other in the objective space (right) but far from each other in thesolution space (left). This figure was made using [16, 18] as reference.

The above-discussed issue has been addressed as a multi-modal MOP (MMOP)[13, 15, 19]. While the goal of MOPs is to achieve a good approximation of thePF, the goal of MMOPs is to obtain as many as possible equivalent Pareto optimal solutions. For example, in Fig. 1, it is sufficient to find one of xa , xb , andxc for multi-objective optimization. This is because their objective vectors arealmost the same. In contrast, it is necessary to find all of xa , xb , and xc for
multi-modal multi-objective optimization.
Since most MOEAs (e.g., NSGA-II [2] and SPEA2 [20]) do not have mechanisms to find diverse solutions in the solution space, they are unlikely to perform
well for MMOPs. Thus, multi-modal multi-objective evolutionary algorithms(MMEAs) that handle the solution space diversity are necessary for MMOPs.
MMEAs need the following three abilities: (i) the ability to find solutions withhigh quality, (ii) the ability to find diverse solutions in the objective space, and(iii) the ability to find diverse solutions in the solution space. MOEAs need theabilities (i) and (ii) to find a solution set that approximates the Pareto front
in the objective space. Multi-modal single-objective evolutionary algorithms[21, 18] need the abilities (i) and (iii) to find a set of global optimal solutions.
In contrast, an efficient MMEA should be able to handle all three abilities (i)–(iii) in order to find all equivalent Pareto optimal solutions.
Wheres multi-modal single-objective optimization is a hot research topic in
the evolutionary computation community, MMOPs have not been well studied[18]. As reviewed in [19], some MMEAs have been proposed, such as Omnioptimizer [15], Niching-CMA [16], and MO Ring PSO SCD [22]. In addition,
most existing MMEAs use the Pareto dominance relation for the mating andenvironmental selections. Since dominance-based MOEAs do not work wellon MaOPs as mentioned above, it is expected that existing MMEAs are notcapable of handling many objectives. In this paper, an MMOP with four or moreobjectives is denoted as a multi-modal many-objective optimization problem(MMaOP). While MMaOPs appear in real-world problems (e.g., five-objective
rocket engine design problems [12]), there is no method that can efficiently locate
multiple equivalent Pareto optimal solutions in problems with more than three3

objectives.
To address this issue, this paper proposes a niching indicator-based multimodal many-objective optimizer (NIMMO). Niching is a loosely defined term.
Usually, in evolutionary multi-modal single-objective optimization, the term“niching” indicates a mechanism to find multiple local optimal solutions and/or
maintain the diversity of the population [18]. An evolutionary algorithm with
the niching method is usually capable of finding multiple local optima. Also,
as mentioned above, some indicator-based MOEAs can handle a large numberof objectives. Therefore, it is expected that an efficient multi-modal manyobjective optimization method can be realized by incorporating a niching mechanism in the solution space into an indicator-based MOEA.
Our contributions in this paper are as follows:
1. We address MMaOPs. Although most previous studies address only twoand three-objective MMOPs, MMaOPs can be found in real-world applications and important problems.
2. We propose NIMMO. While most existing MMEAs were designed for twoand three-objective MMOPs, NIMMO is designed for MMaOPs. NIMMO
handles the diversity in both the objective and solution spaces in a novel
manner.
3. We investigate the performance of NIMMO on various test problems.
We demonstrate that NIMMO can find multiple equivalent solutions on
MMaOPs with up to 15 objectives. We also compare NIMMO with stateof-the-art MMEAs (TriMOEA-TA&R [23], MO Ring PSO SCD [22], andOmni-optimizer [15]).
4. We analyze the performance of NIMMO. We examine the influence of acontrol parameter and the population size on the performance of NIMMO.
In addition, we investigate how the fitness assignment scheme affects theperformance of NIMMO.
The rest of this paper is organized as follows. Section 2 describes related
work. Section 3 introduces NIMMO. Section 4 describes experimental settings.
Section 5 reports experimental results of NIMMO on MMaOPs. Section 6 concludes this paper with discussions on directions for future work.
2. Related work
Since the proposed NIMMO is an indicator-based MMEA, this section describes indicator-based MOEAs and MMEAs as related work. Subsection 2.1
explains indicator-based MOEAs. Subsection 2.2 describes MMEAs.
2.1. Indicator-based MOEAs for multi-objective optimization
Unary quality indicators map a solution set A to a real value that quantitatively represents one or more properties of A in the objective space [6, 24]. For
example, the generational distance (GD) [25] evaluates the convergence of A to
the PF in the objective space. The hypervolume [26] evaluates both convergence4and diversity of A in the objective space. In general, quality indicators play acentral role in evaluating the performance of MOEAs.
Indicator-based MOEAs directly use a quality indicator in the search [7].
For each iteration, so-called fitness values are assigned to individuals in the
population using a quality indicator. Then, the fitness values are used in themating and/or environmental selection. The first indicator-based MOEA is
IBEA [7]. IBEA can use any binary Pareto-compliant indicator. The additive 
indicator (I+ ) and the hypervolume difference indicator (IHD ) are proposed in[7] for IBEA. Subsection 3.1 explains I+ in detail. The IHD (y, x) value is thevolume of the space dominated by x but not by y with respect to a reference
point [7]. A binary version of the R2 indicator [27] (IR2 ) is also proposed in[28] for IBEA. The R2 indicator uses a utility function with a weight vector set
that maps a solution set to a scalar value. IR2 (y, x) is calculated from the R2
values of {y} and {y, x}.
One of the main issues of IBEA with I+ is that the distribution of individuals in the objective space is biased [3, 29, 30]. To address this issue, somevariants of IBEA with I+ have been proposed in the literature. Two Arch2[31] uses two archives: the convergence archive and the diversity archive. Whileindividuals in the convergence archive are updated based on the selection in
IBEA with I+ , individuals in the diversity archive are updated based on their
Lp norm-based crowding distance values. SRA [32] simultaneously uses I+ and
ISDE , where ISDE is the shift-based density estimation in SPEA2+SDE [33].
SRA uses the stochastic ranking [34] to sort individuals based on their I+ and
ISDE values. The BCE framework [35] aims to address the biased distributionof individuals in the population in IBEA with I+ and MOEA/D [36]. BCE
uses a bounded external archive that maintains well-distributed non-dominatedindividuals in the objective space. The selection is performed from individuals
in the population and the external archive. The above-mentioned Two Arch2,
SRA, and IBEA with BCE have shown good performance for many-objective
optimization.
Another representative indicator-based MOEA is SMS-EMOA [8], which isa steady-state algorithm. In the environmental selection in SMS-EMOA, theprimary and secondary criteria are based on the non-domination levels andthe hypervolume contribution, respectively. The worst individual regarding thehypervolume contribution is removed from the last front. Since the hypervolume
calculation is expensive on problems with many objectives, it is difficult to apply
SMS-EMOA to MaOPs. HypE [9] addresses this issue by using Monte Carlosimulation for the hypervolume approximation. HypE uses an approximatedhypervolume value instead of its exact value. The results presented in [9] showthat HypE has good performance on problems with up to 50 objectives. Some
MOEAs use a similar framework to SMS-EMOA, such as R2-EMOA [37], MyODEMR [38], and AR-MOEA [39]. Some other hypervolume-based MOEAs havealso been proposed, such as FV-MOEA [40] and LIBEA [41].
In addition, a number of indicator-based MOEAs have been proposed in
the literature. GD-MOEA [42] and IGD+ -EMOA [43] use the GD indicator[25] and the IGD+ indicator [44], respectively. DDE [45] uses the averaged
5

Hausdorff distance indicator (∆p ) [46] for the search, where ∆p is calculated
from GD and IGD [47]. Although MOMBI-II [29] uses the R2 indicator similar
to IBEA with IR2 , the ranking method in MOMBI-II differs from that in IBEA.
AGE [48] uses an unbounded external archive that maintains all non-dominatedsolutions found so far. In AGE, the fitness value of an individual is based on
the contribution to the additive  indicator [6] (not I+ in IBEA [7]) using theunbounded external archive. Unlike most indicator-based MOEAs, MIHPS [49]adaptively uses multiple indicators. DIVA [50] uses an indicator that integrates
the solution space diversity into the hypervolume indicator.
2.2. MMEAs
As reviewed in [19], multi-modal multi-objective optimization has been studied since around 2005 in the evolutionary computation community. Real-world
applications of MMEAs include functional brain imaging problems [13], molecular drug design problems [51], and military operational planning problems [52].
Here, we explain some MMEAs. For an exhaustive survey of MMEAs, refer to[19].
Omni-optimizer [15] is the most representative MMEA. Omni-optimizer is
an NSGA-II-based algorithm which can be applied to various problem domains,
including MOPs and MMOPs. The population in Omni-optimizer is initializedusing the Latin hypercube sampling method so that individuals are distributed
as uniformly as possible. In the environmental selection, Omni-optimizer uses
the -dominance-based non-dominated sorting and the alternative crowding distance. The crowding distance in Omni-optimizer is based on both the objective
and solution spaces. Roughly speaking, the alternative crowding distance valueof an individual is the maximum value between its crowding distance values in
the objective and solution spaces. In contrast to NSGA-II, Omni-optimizer uses
the restricted tournament mating selection. First, an individual xa is randomly
selected from the population. The nearest neighborhood individual xb of xa inthe solution space is also selected from the population. Then, xa and xb arecompared based on their non-domination levels. Ties are broken by the alternative crowding distance values of xa and xb . The winner can be a parent. This
procedure is iterated to select another parent.
Niching-CMA [16] and MO Ring PSO SCD [22] are CMA-ES and PSO algorithms for multi-modal multi-objective optimization, respectively. NichingCMA uses a niching strategy with the dynamic peak identification method[53] in the environmental selection. The number of niches and the niche radius are adaptively adjusted. For each niche, better individuals regarding nondomination levels survive to the next iteration. Niching-CMA uses the aggregated distance in both the objective and solution spaces. Although the crowding
distance metric in MO Ring PSO SCD is similar to that of Omni-optimizer,
MO Ring PSO SCD handles extreme individuals in the objective space in adifferent manner. MO Ring PSO SCD also uses the index-based ring topology [54] to maintain high diversity in the solution space. SMPSO-MM [55]is based on MO Ring PSO SCD. SMPSO-MM uses the self-organizing map to6

select neighborhood best particles. The new positions of the particles are perturbed by the mutation operation with a pre-defined probability. In addition,some dominance-based MMEAs have been proposed, such as SPEA2+ [56],4D-Miner [13], PQ, -MOEA [11], and DN-NSGA-II [57].
Some decomposition-based MMEAs have been proposed, such as a multistart decomposition-based approach [10]. MOEA/D-AD [58] can assign one ormore individuals to each decomposed subproblem. This mechanism is to handle
equivalent solutions in the framework of MOEA/D [36]. First, a child is assigned
to the subproblem whose weight vector is closest to the objective vector of thechild, with respect to the perpendicular distance. If the child is close to oneof the individuals assigned to the same subproblem in the solution space, they
are compared based on their scalarizing function values. Otherwise, the childis assigned to the subproblem with no comparison. Unlike other MMEAs, thepopulation size in MOEA/D-AD is adaptively adjusted.
TriMOEA-TA&R [23] consists of various components, including the twoarchives-based method as in Two Arch2 [31] and the decision variable-decomposition
method in MOEA/DVA [59]. At the beginning of the search, TriMOEA-TA&R
decomposes design variables of a given problem into position-related variables
and distance-related variables. While position-related variables affect only the
position of the objective vector on the PF, distance-related variables affect only
the distance between the objective vector and the PF. Then, TriMOEA-TA&R
exploits the properties of position-related and distance-related variables in an
explicit manner. In TriMOEA-TA&R, a child is assigned to the subproblem
whose weight vector is closest to the objective vector of the child, with respect to their angle as in RVEA [60]. Although the methods of assigning thechild to the subproblem in TriMOEA-TA&R and MOEA/D-AD are similar, thepopulation size in TriMOEA-TA&R is constant. While MOEA/D-AD uses therelative distance-based neighborhood criterion, TriMOEA-TA&R uses the absolute distance-based neighborhood criterion. In TriMOEA-TA&R, if the distance
between two individuals is less than a pre-defined value, they are said to be in
the same niche.
3. Proposed method
This section introduces the proposed NIMMO. NIMMO can use any indicatorbased fitness assignment scheme. In this paper, we use the fitness assignment
scheme with the additive epsilon indicator I+ in IBEA [7, 61]. As reviewed
in Subsection 2.1, a number of indicator-based MOEAs have been proposed in
the literature. Since most indicator-based MOEAs ignore the diversity of the
population in the solution space, it is likely that they cannot locate multiple
equivalent Pareto optimal solutions. For example, in Fig. 1, even if all of xa ,xb , and xc are included in the current population, two of them are very likely tobe removed from the population after the environmental selection. NIMMO addresses this issue by using a niching mechanism in the solution space. Although
DIVA [50] uses the indicator that handles the solution space diversity, the indicator calculation in DIVA is expensive. For this reason, DIVA has been applied
7

only to two-objective problems. In contrast, NIMMO handles the solution spacediversity in a computationally cheap manner.
MMEAs need to handle the diversity in both the objective and solution
spaces. While most MMEAs aggregate the objective and solution space diversity
metrics (e.g., Omni-optimizer and SMPSO-MM), NIMMO handles the objective and solution space diversity in a two-phase manner similar to MOEA/D-AD
and TriMOEA-TA&R. NIMMO also uses the relative distance-based neighborhood criterion similar to MOEA/D-AD. As described in Subsection 2.2, first,
MOEA/D-AD (TriMOEA-TA&R) assigns the child to the subproblem basedon the perpendicular distance (angle) between the objective vector of the childand the weight vector. Then, the child is compared to individuals in the same
niche in the solution space. In NIMMO, first, individuals to be compared withthe child are selected based on the distribution in the solution space. Then, thefitness values are assigned to the child and the selected individuals based on the
distribution in the objective space.
Multi-modal single-objective optimization has been well studied in the evolutionary computation community. As reviewed in [18, 21], a number of niching
techniques have been proposed for multi-modal single-objective optimization.
Traditional approaches include fitness sharing [62], crowding [63], and clustering [64]. Recent approaches include the adaptive radius-based method [65, 66],
the nearest-better clustering [67], and the taboo region-based method [68]. The
niching method used in NIMMO is similar to the deterministic crowding method[69]. While the deterministic crowding method compares the child to its closest
parent in the solution space, NIMMO compares the child to its T closest individuals in the solution space. In addition, the environmental selection in NIMMO
is based on the indicator values of the child and those T closest individuals.
3.1. The procedure of NIMMO
Algorithm 1 shows the procedure of NIMMO. Similar to SMS-EMOA [8],NIMMO is a steady state MMEA and generates a child u for each iteration t.
After the initialization of the population P (line 1), the following operations
are iteratively applied until a termination condition is satisfied.
At the beginning of each iteration t, the mating and the reproduction are
performed (lines 3–5). First, parent individuals xa and xb are randomly selectedfrom P (line 3). Then, a child u is generated by recombining xa and xb (line4). A mutation operator is also applied to u (line 5).
After the generation of u, the environmental selection is performed (lines
6–13). Since NIMMO is based on the (µ + 1)-selection, only one individual is discarded from the population P . First, for each individual xi in P ,the distance di between xi and u is calculated (lines 6–7). The function
normalizedEuclideanDistance(x, y) in line 7 returns the normalized Euclidean
distance between x and y using the upper and lower bounds for each variable.
Next, T closest individuals to u in the solution space are selected from P basedon their distances and stored into R (line 8), where T is the neighborhood size.
The T individuals stored in R are temporarily removed from P (line 9), and u
is inserted into R (line 10). Thus, R is a set of T neighbors of u in the solution
8

Algorithm 1: The proposed NIMMOt ← 1, initialize the population P = {x1 , ..., xµ };while The termination criteria are not met do
3
Randomly select the parent individuals xa and xb from P such thata 6= b;4
Generate the child u by recombining xa and xb ;5
Apply the mutation operator to u;
6
for i ∈ {1, ..., µ} do
7di ← normalizedEuclideanDistance(xi , u);
1

2

8
9
10
11
12
13
14

Select T closest individuals to u based on the distance values d1 , ..., dµand store them into R;
P ← P \R;R ← R ∪ {u};Calculate the fitness value F (x) of x ∈ R using (2);Remove the worst individual xworst from R;
P ← P ∪ R;t ← t + 1;

space and u. The environmental selection in NIMMO is performed only among
T + 1 individuals in R. The worst individual in R cannot survive to the next
iteration.
In line 11, a fitness value F (x) is assigned to each individual x in R (notP ). Although any fitness assignment method of indicator-based MOEAs canbe incorporated into NIMMO, we use that of IBEA [7], which is the most basic
one. The fitness assignment in IBEA is performed as follows:X
I(y, x)F (x) =exp −,(2)
κ I maxy ∈R\{x}where I is a Pareto dominance preserving binary indicator. If x dominates y,
I(x, y) < I(y, x). The maximum absolute indicator value I max is defined asfollows: I max = maxx,y ∈R |I(x, y)|. The scale factor κ in (2) is usually set to0.05. In (2), x and every y in R\{x} are compared with each other according
to their objective vectors. A smaller fitness value represents a better individual.
We use the additive epsilon indicator I+ [7, 61] for I in (2). As describedin Subsection 2.1, some variants of IBEA with I+ (e.g., Two Arch2 and SRA)show good performance on MaOPs. Therefore, it is expected that NIMMO
with I+ is capable of handling many objectives. The indicator value I+ (y, x)
is calculated as follows: 0I+ (y, x) = max
fi (y) − fi0 (x) ,(3)i∈{1,...,M }where fi0 (x) is the normalized objective value of x with fimin = minx∈R {fi (x)}and fimax = maxx∈R {fi (x)} as follows: fi0 (x) = (fi (x) − fimin )/(fimax − fimin ).
9In (3), I+ (y, x) represents the minimal translation such that y dominates x [61].
Although the original I+ is for comparing two sets of solutions, (3) describesI+ as an indicator for comparing two solutions for the sake of simplicity.
The worst individual xworst regarding the fitness value is discarded fromR (line 12). After the removal operation, the remaining individuals in R are
returned into P (line 13).
3.2. The effect of the niching mechanism in NIMMO on MMOPs
We explain the effect of the niching mechanism in NIMMO. Fig. 2 shows
examples of fitness assignment results in IBEA and NIMMO on a two-objective
two-variable problem.
Figs. 2(a) and (b) show the distribution of 12 individuals {A, ..., L} in thesolution and objective spaces, respectively. The fitness value of each individualin IBEA is also shown, where it can be calculated by replacing R with P in (2).
In Fig. 2, let us assume that L is a newly generated child. Our next task
is to remove a single solution from the 12 individuals. As seen from Figs. 2(a)and (b), L and H are close to each other in the objective space but far fromeach other in the solution space. Therefore, it is desirable that both L and H
remain in the population in the next iteration for multi-modal multi-objective
optimization.
Fig. 2(b) shows that L has the worst fitness value (0.97) in IBEA due
to its poor diversity in the objective space. As a result, L is discarded fromthe population after the environmental selection. Since L contributes to thediversity of the population in the solution space more than some individuals,
such a removal operation is not appropriate.
NIMMO addresses this issue by using the niching mechanism in the solution
space. If T = 3, neighborhood individuals of the child L in the solution spaceare D, E, and K. In NIMMO, the environmental selection is performed only
among the four individuals (L, D, E, and K). Fig. 2(c) shows the distributionof the four individuals in the objective space and their fitness values obtainedby (2). Since the fitness calculation in (2) is relative, the fitness value of eachindividual in NIMMO is different from that in IBEA.
Fig. 2(c) indicates that D is the worst individual and thus removed from the
population in NIMMO. As a result, the individual L, which does not survive
to the next iteration in IBEA, can remain in the population. As shown here,
NIMMO can maintain the solution space diversity in the population by using
the niching mechanism.
4. Experimental settings
4.1. Test problems
We use the following 13 two-objective MMOPs: the Two-On-One problem[70], the Omni-test problem [15], the 3 SYM-PART problems [10], and the 8
MMF problems [57]. Also, we evaluate the performance of MMEAs on two
MMaOPs: the Polygon problem [71] and its rotated version named “RPolygon
10

J

CH
I

A

6

K

G

4
F

2
0
0

2

x1

L

D
E

B

4

F (A) = 0.62 F (G) = 0.73

68(a) Solution space10

F (B) = 0.57 F (H) = 0.96
10 A
F (C) = 0.61 F (I) = 0.43
8 B
F (D) = 0.63 F (J) = 0.51
F (E) = 0.65 F (K) = 0.59
C6
D F (F) = 0.68 F (L) = 0.97
E
4
F
G
H
2
L I
J
K
0
0 2 4 6 8 10
f1(b) Objective space in IBEA

f2

x2

8

f2

10

10
8
6
4
2
0

F (D) = 0.18
F (E) = 0.03
F (K) = 0.01
F (L) = 1.48 × 10−4

D
E
L
K

0

2

4

6
f1

8

10

(c) Objective space in NIMMO

Figure 2: Illustration of the fitness assignment in IBEA and NIMMO with T = 3 on a
two-objective two-variable problem.

problem”. The number of Pareto optimal solutions with equivalent quality
N same is 2 for Two-On-One, 9 for the three SYM-PART problems, and 360 for
the Omni-test problem. For the eight MMF problems, N same is 2 for MMF1,
MMF2, MMF3, and MMF7, and N same is 4 for the other problems. The numberof objectives M can be set to an arbitrary number in the Polygon problem. Thus,
the Polygon problem is suitable for evaluating the scalability of MMEAs withrespect to M . According to [72], we set M of the Polygon problem as follows:M ∈ {3, 5, 8, 10, 15}. In total, we use 23 problem instances. In our study, N samein the Polygon problem is set to nine. The number of design variables D is fivefor the Omni-test problem and two for the other problems.
4.2. Performance indicators
We use the inverted generational distance (IGD) [6], IGDX [14], Pareto sets
proximity (PSP) [22] for performance assessment of MMEAs. All indicators
require a set of reference solutions A∗ . For all problems (except for the eight
MMF problems), we use two reference solution sets A∗sol and A∗obj of about 5 000
Pareto optimal solutions. While the solutions in A∗sol are uniformly distributedin the solution space, those in A∗obj are uniformly distributed in the objective
space. We use A∗sol for the IGDX and PSP calculations and A∗obj for the IGD
calculation. The reference point sets are available at the supplementary website(https://sites.google.com/view/nimmopt/). For the eight MMF problems,
we use reference solution sets provided by the authors of [22].
Below, A denotes a set of non-dominated solutions of the final populationof an MMEA. The IGD value [6] is the average distance from each solution inA∗ to its nearest non-dominated solution in A in the objective space as follows:


n
o
1  X
IGD(A) =min ED f (x), f (z)  ,
(4)
|A∗ |
∗ x∈A
z ∈A
where ED(x, y) is the Euclidean distance between x and y. The IGD metricevaluates the quality of A in terms of both convergence to the PF and diversity
in the objective space.
11

While IGD measures the quality of A in the objective space, IGDX evaluates
how good A approximates A∗ in the solution space. Similar to the IGD indicator
described in (4), the IGDX value of A is given as follows:


n
o
1  X
(5)
IGDX(A) =min ED x, z  .
|A∗ |
∗ x∈A
z ∈A
Similar to IGDX, the following PSP evaluates the quality of A in the solution
space:
PSP(A) =CR(A),IGDX(A)(6)

The cover rate (CR) in (6) is given as follows:
1
! 2D
D
Y
CR(A) =
δi
,

(7)

i=1

δi =

min(x∗,max, xmax) − max(x∗,min, xmin)iiiix∗,max− x∗,min
i
i

!2
,

(8)

where xi∗,min and x∗,max
in (8) are the minimum and maximum values of theii-th design variable in the Pareto solution set, respectively. If x∗,max= x∗,min
,
i
i
∗,max
∗,min
max
δi is specified as δi = 1. If xmin
≥
x
or
x
≤
x
,
δ
is
specified
as
i
i
i
i
i
δi = 0.
On the one hand, since a small IGD value indicates that A is a good approximation of the PF, the corresponding MMEA is an efficient multi-objective
optimizer. On the other hand, A with a small IGDX value and a large PSP value
is a good approximation of the Pareto optimal solutions. Thus, the MMEA finding such an A is a high performance multi-modal multi-objective optimization
method.
Fig. 3 shows a comparison of two sets of non-dominated solutions A1 and
2A in the objective (Figs. 3(a) and (c)) and solution (Figs. 3(b) and (d)) spaceson the Two-On-One problem, where A1 and A2 are evenly generated on thePS, and |A1 | = |A2 | = 20. The IGD and IGDX values of A1 and A2 are also
shown in Fig. 3. In the Two-On-One problem, there are two equivalent Pareto
optimal solution sets that are symmetrical with respect to the origin.
Figs. 3(a) and (c) show that A2 covers the PF more densely than A1 .
However, A2 does not cover the whole PS while A1 covers the whole PS as
shown in Figs. 3(b) and (d). Thus, A1 and A2 are good approximations of the
PS and the PF, respectively. While the IGD value of A2 is approximately twotimes better than that of A1 , the IGDX value of A2 is approximately 10 times
worse than that of A1 . It is likely that if a set of non-dominated solutions has
a very good IGDX value (e.g., A1 ), its IGD value is good. However, a set ofnon-dominated solutions with a very good IGD value (e.g., A2 ) can have a bad
IGDX value.
12

5
3

1

x2

f2

4

2

IGD(A1) = 0.355

2

−1

1
0

0

−2
−2

8 10 12 14 16 18 20
f1

IGD(A2) = 0.183

3

1

2

1

2

0
−1

1
0

0
x1

2

x2

f2

4

−1

(b) Solution Space(a) Objective Space

5

IGDX(A1) = 0.059−2−2

8 10 12 14 16 18 20
f1

IGDX(A2) = 0.561
−1

0
x1

1

2

(d) Solution Space(c) Objective Space

Figure 3: Comparison of two sets of uniformly generated non-dominated solutions A1(upper) and A2 (lower) on the Two-On-One problem.

4.3. MMEAs
We compare NIMMO with the following MMEAs and MOEAs: Omnioptimizer [15], MO Ring PSO SCD [22], TriMOEA-TA&R [23], IBEA [7], NSGAII [2], MOEA/D [36], NSGA-III [72], and MOEA/DD [73]. NSGA-II, IBEA, and
MOEA/D are representative dominance-, indicator-, and decomposition-based
MOEAs, respectively. NSGA-III and MOEA/DD are well-known MOEAs for
MaOPs. Omni-optimizer is a representative MMEA for MMOPs. MO Ring PSO SCD
is a recently proposed multi-modal multi-objective PSO algorithm. TriMOEATA&R is the latest MMEA.
We used available source code through the Internet for algorithm implementation. A set of weight/reference vectors of MOEA/D, MOEA/DD, and
NSGA-III were generated using the simplex-lattice design for M ≤ 5 and its
two-layered version [72] for M ≥ 8. For all algorithms (except for NSGA-II
and Omni-optimizer), the population size µ was set to 200, 210, 210, 156, 210,
230, and 135 for M = 2, 3, 5, 8, 9, 10, and 15, respectively. The number ofweight vectors cannot be specified arbitrarily due to the combinatorial natureof the simplex-lattice design. For this reason, we selected the population size µsuch that µ is close to 200 as possible. Since µ values in NSGA-II and Omnioptimizer must be multiples of two and four respectively, they were set to theclosest multiples of two and four which are larger than the above µ values respectively. The number of maximum function evaluations was set to 10 000, and31 runs were performed.

13

We set T of NIMMO to b0.1µc. To investigate the contribution of the nichingmechanism of NIMMO, we evaluated the performance of a steady-state IBEA
with I+ and the randomized mating selection, which is equivalent to NIMMO
with T = µ. Since our preliminary results show that this version of IBEA
outperforms the original IBEA regarding the IGD, IGDX, and PSP indicators,
we show only the results of the former. According to [73], the neighborhood
size of MOEA/D and MOEA/DD was set to 20. We used the Tchebycheff and
PBI functions for MOEA/D and MOEA/DD, respectively. The SBX crossoverand the polynomial mutation were used in all MMEAs and MOEAs, except for
MO Ring PSO SCD. Their control parameters were set as follows: pc = 1, ηc =20, pm = 1/D, and ηm = 20. All the control parameters of MO Ring PSO SCD
and TriMOEA-TA&R were set according to [22, 23].

14

Table 1: Results of the 9 algorithms on the 23 test problem instances. The mean IGD values for 31 runs are shown. The numbers inparenthesis indicate the ranks of the nine algorithms.

15

NIMMO

TriMOEATA&R

MO Ring
PSO SCD

Omnioptimizer

IBEA

NSGA-II

MOEA/D

NSGA-III

MOEA/DD

Two-On-One
Omni-test
SYM-PART1SYM-PART2SYM-PART3

0.0250 (3)0.0769 (7)0.0246 (4)0.0365 (6)0.0298 (4)0.0627− (9)0.1365− (8)0.0379− (7)0.0396− (7)0.0798− (9)0.0348− (7)0.2957− (9)0.0386− (9)0.0446− (9)0.0493− (8)0.0254− (5)0.0268+ (2)0.0137+ (2)0.0183+ (2)0.0413≈ (7)0.0233+ (1)0.0315+ (3)0.0169+ (3)0.0199+ (3)0.0196+ (2)0.0245≈ (2)0.0266+ (1)0.0104+ (1)0.0120+ (1)0.0118+ (1)0.0264≈ (6)0.0622+ (5)0.0341− (6)0.0361≈ (5)0.0369− (5)0.0251≈ (4)0.0525+ (4)0.0262− (5)0.0284+ (4)0.0279≈ (3)0.0429− (8)0.0700≈ (6)0.0383− (8)0.0413− (8)0.0389− (6)

MMF1
MMF2
MMF3
MMF4
MMF5
MMF6
MMF7
MMF8

0.0050 (7)0.0216 (5)0.0190 (6)0.0048 (7)0.0039 (6)0.0044 (7)0.0030 (7)0.0092 (8)0.0092≈ (8)0.0294− (8)0.0230− (8)0.0496− (9)0.0087− (8)0.0087− (8)0.0064− (9)0.0113− (9)0.0037+ (5)0.0216≈ (6)0.0154≈ (5)0.0037+ (5)0.0038+ (5)0.0035+ (5)0.0037− (8)0.0047+ (6)0.0033+ (3)0.0178+ (2)0.0134+ (3)0.0029+ (4)0.0032+ (3)0.0031+ (4)0.0029≈ (6)0.0031+ (4)0.0043+ (6)0.0286− (7)0.0215≈ (7)0.0057− (8)0.0042− (7)0.0043+ (6)0.0027+ (4)0.0079+ (7)0.0028+ (1)0.0179+ (3)0.0114+ (2)0.0026+ (3)0.0028+ (1)0.0027+ (2)0.0026+ (3)0.0026+ (2)0.0152− (9)0.1246− (9)0.0966− (9)0.0047+ (6)0.0104− (9)0.0108− (9)0.0027+ (5)0.0044+ (5)0.0034+ (4)0.0181+ (4)0.0144+ (4)0.0023+ (2)0.0034+ (4)0.0031+ (3)0.0025+ (2)0.0026+ (3)0.0030+ (2)0.0118+ (1)0.0104+ (1)0.0020+ (1)0.0029+ (2)0.0026+ (1)0.0023+ (1)0.0022+ (1)3-Polygon
5-Polygon
8-Polygon
10-Polygon
15-Polygon
3-RPolygon
5-RPolygon
8-RPolygon
10-RPolygon
15-RPolygon

0.0025 (2)0.0044 (2)0.0069 (2)0.0064 (2)0.0106 (2)0.0025 (2)0.0044 (2)0.0069 (2)0.0064 (2)0.0105 (2)0.0040− (8)0.0149− (9)0.0180− (8)0.0204− (8)0.0211− (8)0.0046− (9)0.0149− (9)0.0190− (8)0.0185− (8)0.0226− (8)0.0034− (6)0.0057− (6)0.0092− (7)0.0087− (7)0.0136− (7)0.0034− (7)0.0058− (7)0.0093− (7)0.0086− (7)0.0140− (7)0.0028− (4)0.0051− (4)0.0082− (4)0.0074− (3)0.0122− (4)0.0028− (4)0.0052− (4)0.0083− (3)0.0075− (3)0.0125− (5)0.0022+ (1)0.0036+ (1)0.0057+ (1)0.0053+ (1)0.0085+ (1)0.0022+ (1)0.0036+ (1)0.0057+ (1)0.0053+ (1)0.0085+ (1)0.0028− (3)0.0052− (5)0.0083− (5)0.0075− (4)0.0123− (5)0.0028− (3)0.0052− (5)0.0083− (4)0.0075− (4)0.0123− (4)0.0034− (7)0.0051− (3)0.0084− (6)0.0080− (6)0.0106≈ (3)0.0034− (6)0.0051− (3)0.0086− (6)0.0082− (6)0.0107− (3)0.0028− (5)0.0057− (7)0.0081− (3)0.0077− (5)0.0128− (6)0.0029− (5)0.0058− (6)0.0083− (5)0.0078− (5)0.0131− (6)0.0042− (9)0.0136− (8)0.0279− (9)0.0348− (9)0.0406− (9)0.0043− (8)0.0133− (8)0.0279− (9)0.0321− (9)0.0413− (9)Table 2: Results of the 9 algorithms on the 23 test problem instances. The mean IGDX values for 31 runs are shown.

16

NIMMO

TriMOEATA&R

MO Ring
PSO SCD

Omnioptimizer

IBEA

NSGA-II

MOEA/D

NSGA-III

MOEA/DD

Two-On-One
Omni-test
SYM-PART1SYM-PART2SYM-PART3

0.0226 (1)1.6052 (1)0.0530 (1)0.0838 (1)0.1418 (1)0.5699− (9)5.2291− (9)0.7040− (6)9.8161− (9)5.5682− (8)0.0296− (3)2.4142− (6)0.1673− (2)0.2303− (2)0.8697− (3)0.0261− (2)1.7443− (2)0.3452− (3)0.4325− (3)0.8538− (2)0.0548− (6)2.4403− (7)3.8899− (8)3.7481− (7)3.8662− (7)0.0405− (4)1.9110− (4)1.0533− (7)1.8697− (6)2.3685− (6)0.2209− (8)3.5174− (8)6.6705− (9)7.5890− (8)6.6383− (9)0.0461− (5)1.8711− (3)0.6525− (5)1.5404− (5)2.3069− (5)0.0727− (7)2.1113− (5)0.5588− (4)1.5147− (4)2.2809− (4)

MMF1
MMF2
MMF3
MMF4
MMF5
MMF6
MMF7
MMF8

0.0611 (3)0.0440 (3)0.0440 (4)0.0372 (3)0.0846 (1)0.1037 (3)0.0216 (1)0.2064 (3)0.1369− (8)0.0697− (7)0.0747− (7)0.1628≈ (8)0.1993− (8)0.1641− (8)0.0591− (8)0.4460− (5)0.0490+ (1)0.0430≈ (2)0.0282+ (1)0.0275+ (1)0.0869≈ (2)0.0740+ (1)0.0263− (2)0.0680+ (1)0.0529+ (2)0.0614− (6)0.0423≈ (3)0.0332+ (2)0.0980− (3)0.0835+ (2)0.0267− (3)0.1088+ (2)0.0930− (7)0.1219− (8)0.0841− (8)0.0820− (7)0.1618− (7)0.1471− (7)0.0435− (7)0.7538− (8)0.0665− (5)0.0531≈ (4)0.0448≈ (5)0.0523− (6)0.1187− (4)0.1119− (6)0.0383− (5)0.4556− (6)0.2606− (9)0.3290− (9)0.2359− (9)0.2741− (9)0.4480− (9)0.3753− (9)0.1231− (9)2.2762− (9)0.0765− (6)0.0580≈ (5)0.0488≈ (6)0.0487− (4)0.1263− (6)0.1113≈ (5)0.0414− (6)0.5011− (7)0.0657≈ (4)0.0373+ (1)0.0315+ (2)0.0519− (5)0.1191− (5)0.1054≈ (4)0.0365− (4)0.4089− (4)3-Polygon
5-Polygon
8-Polygon
10-Polygon
15-Polygon
3-RPolygon
5-RPolygon
8-RPolygon
10-RPolygon
15-RPolygon

0.0056 (1)0.0070 (1)0.0089 (1)0.0072 (1)0.0100 (1)0.0059 (1)0.0074 (1)0.0093 (1)0.0076 (1)0.0104 (1)0.0063− (2)0.0162− (6)0.0136− (2)0.0123− (5)0.0122− (2)0.0295− (6)0.0400− (7)0.0747− (6)0.0404− (7)0.0667− (6)0.0091− (4)0.0113− (4)0.0143− (4)0.0120− (4)0.0156− (4)0.0090− (3)0.0113− (3)0.0144− (3)0.0118− (3)0.0154− (2)0.0083− (3)0.0110− (3)0.0140− (3)0.0112− (3)0.0155− (3)0.0085− (2)0.0110− (2)0.0138− (2)0.0112− (2)0.0166− (3)0.0133− (7)0.0109− (2)0.0151− (5)0.0108− (2)0.0212− (6)0.0339− (7)0.0289− (5)0.0426− (5)0.0257− (5)0.0408− (5)0.0093− (5)0.0119− (5)0.0189− (7)0.0123− (6)0.0201− (5)0.0194− (4)0.0196− (4)0.0402− (4)0.0228− (4)0.0373− (4)0.0496− (9)0.0584− (8)0.1059− (8)0.0799− (8)0.0949− (8)0.1652− (9)0.1798− (9)0.1662− (9)0.1418− (9)0.1542− (9)0.0097− (6)0.0169− (7)0.0187− (6)0.0139− (7)0.0286− (7)0.0260− (5)0.0348− (6)0.0788− (7)0.0289− (6)0.0952− (7)0.0232− (8)0.0883− (9)0.1349− (9)0.0858− (9)0.1364− (9)0.0463− (8)0.1021− (8)0.1454− (8)0.1182− (8)0.1406− (8)Table 3: Results of the 9 algorithms on the 23 test problem instances. The mean PSP values for 31 runs are shown.

17

NIMMO

TriMOEATA&R

MO Ring
PSO SCD

Omnioptimizer

IBEA

NSGA-II

MOEA/D

NSGA-III

MOEA/DD

Two-On-One
Omni-test
SYM-PART1SYM-PART2SYM-PART3

40.09 (1)0.49 (1)18.73 (1)
11.47 (1)9.48 (1)0.13− (9)0.00− (9)10.35− (2)0.00− (9)0.02− (8)32.72− (3)0.33− (5)5.98− (4)4.81− (3)1.28− (3)36.84− (2)0.42− (2)7.63− (3)4.83− (2)2.15− (2)15.58− (6)0.12− (7)1.31− (8)0.13− (7)0.12− (7)24.19− (4)0.33− (4)4.73− (7)0.84− (6)0.23− (6)3.56− (8)0.00− (8)0.08− (9)0.01− (8)0.01− (9)20.89− (5)0.35− (3)5.32− (5)0.93− (5)0.25− (5)13.83− (7)0.22− (6)5.14− (6)1.21− (4)0.75− (4)

MMF1
MMF2
MMF3
MMF4
MMF5
MMF6
MMF7
MMF8

16.22 (3)20.61 (2)19.90 (5)27.06 (3)11.68 (1)9.09 (5)43.02 (1)4.97 (3)11.64− (7)13.73− (7)12.22− (7)22.98≈ (4)7.72− (6)9.49+ (3)18.86− (8)2.22− (4)19.89+ (1)19.36≈ (4)29.54+ (1)35.25+ (1)11.20− (2)13.16+ (1)36.57− (3)14.29+ (1)18.34+ (2)16.87≈ (6)21.72≈ (3)29.99+ (2)10.01− (3)11.58+ (2)36.68− (2)8.84+ (2)10.21− (8)7.48− (8)9.84− (8)12.00− (8)5.97− (8)6.09− (8)19.59− (7)0.73− (8)14.67− (5)20.07≈ (3)20.62≈ (4)19.30− (7)8.11− (4)8.80≈ (6)25.25− (5)1.48− (6)2.02− (9)0.85− (9)0.96− (9)2.42− (9)0.78− (9)0.77− (9)5.24− (9)0.11− (9)12.59− (6)18.17≈ (5)18.38≈ (6)20.12− (5)7.35− (7)8.67≈ (7)21.52− (6)1.16− (7)14.77− (4)25.75≈ (1)25.28+ (2)19.60− (6)7.99− (5)9.14≈ (4)25.88− (4)1.93− (5)3-Polygon
5-Polygon
8-Polygon
10-Polygon
15-Polygon
3-RPolygon
5-RPolygon
8-RPolygon
10-RPolygon
15-RPolygon

178.71 (1)141.30 (1)111.35 (1)138.23 (1)97.52 (1)
166.54 (1)132.04 (1)104.28 (1)126.33 (1)90.97 (1)158.39− (2)61.55− (7)72.69− (2)77.00− (6)79.39− (2)49.63− (6)28.83− (7)13.01− (7)25.64− (7)15.45− (6)109.85− (4)87.99− (4)69.98− (4)83.24− (4)63.86− (4)107.81− (3)86.99− (3)67.94− (3)82.54− (3)62.65− (2)120.17− (3)90.96− (2)71.09− (3)89.26− (3)63.97− (3)113.59− (2)88.22− (2)70.04− (2)86.48− (2)58.17− (3)104.94− (6)90.83− (3)68.43− (5)92.51− (2)54.63− (5)41.75− (7)45.59− (5)26.01− (5)45.27− (5)26.24− (4)108.30− (5)84.30− (5)58.88− (7)81.59− (5)52.91− (6)65.24− (4)54.22− (4)26.76− (4)49.44− (4)24.72− (5)31.29− (9)32.56− (8)10.85− (8)16.26− (8)13.80− (8)2.74− (9)2.47− (9)2.07− (9)3.66− (9)2.86− (9)102.75− (7)68.18− (6)60.18− (6)76.29− (7)42.42− (7)53.17− (5)31.70− (6)14.12− (6)41.00− (6)7.56− (7)50.72− (8)10.60− (9)5.86− (9)10.85− (9)5.72− (9)24.66− (8)6.34− (8)3.02− (8)5.06− (8)2.94− (8)

5. Experimental results
5.1. Performance of NIMMO
Tables 1, 2, and 3 show the comparison between the nine algorithms. Tables
1, 2, and 3 show the IGD, IGDX, and PSP values, respectively. In Tables 1,
2, and 3, M -(R)Polygon denotes the M -objective (R)Polygon problem. The
symbols +, −, and ≈ indicate that a given algorithm performs significantly
better (+), significantly worse (−), and not significantly different better or worse(≈) compared to NIMMO according to the Wilcoxon rank-sum test with p <0.05.
• Results on the two-objective MMOPs: Table 1 shows that NSGA-II
achieves good IGD values on the two-objective MMOPs. In particular, NSGAII performs the best on six problems regarding the IGD metric. MOEA/DD
has the best IGD value on six MMF problems.
Table 2 shows that NIMMO obtains the best IGDX values on the seventwo-objective MMOPs. According to the IGDX metric, IBEA is not capableof finding multiple equivalent Pareto optimal solutions. Therefore, the nichingmechanism of NIMMO mainly contributes to its good ability to locate multiple
Pareto optimal solutions with equivalent quality. Table 3 shows that NIMMO
performs the best regarding PSP on the seven two-objective MMOPs. The
results of PSP are consistent with those of IGDX in almost all cases.
Although MOEA/DD does not have any explicit diversity maintenance mechanism in the solution space, MOEA/DD performs the best on MMF2 regarding
IGDX and PSP. Fig. 4 shows the distribution of non-dominated solutions in
the final population of each algorithm in the solution space on MMF2. Resultsof a single run with a median IGDX value among 31 runs are shown. Fig. 5
shows the reference solutions used for the IGDX calculation for MMF2. Asseen from Fig. 4, MOEA/DD finds well-distributed non-dominated solutionsin the solution space. The distribution of solutions found by the four MMEAs
is worse than that by MOEA/DD. The reason for these results is not obvious. The good performance of MOEA/DD may be due to some unintended
properties of MMF2. For example, diverse solutions can be found on MMF2
without handling the solution space diversity. Further analysis of multi-modal
multi-objective test problems including MMF2 is needed as future research.
It seems that NIMMO cannot find better approximations of the PF (notPS) than NSGA-II. It should be noted that such an observation is not unique
to NIMMO and can be found in other MMEAs [16, 22]. This is because it isdifficult (or almost impossible) to obtain a set of non-dominated solutions thatminimizes both the IGD and IGDX metrics as shown in Fig. 3. As explained in
Subsection 4.2 using the illustrative example, a good solution set regarding IGD
is not always good regarding IGDX, and vice versa (for details, see Subsection
4.2). Thus, there is no single solution set that simultaneously minimizes IGD
and IGDX in most multi-modal multi-objective problems. In other words, theimprovements in IGD and IGDX have a trade-off relationship in most cases.
• Results on the MMaOPs with up to 15 objectives: Table 1 shows that
IBEA achieves the best IGD value on all the MMaOP instances, followed by
18

2.0

2.0

2.0

1.5

1.5

1.5

1.0

1.0

1.0

0.5

0.5

0.5

0.0
0.0

0.2

0.4

0.6

0.8

1.0

0.0
0.0

(a) NIMMO

0.2

0.4

0.6

0.8

1.0

0.0
0.0

(b) TriMOEA-TA&R

2.0

2.0

2.0

1.5

1.5

1.5

1.0

1.0

1.0

0.5

0.5

0.5

0.0
0.0

0.2

0.4

0.6

0.8

1.0

0.0
0.0

0.2

(d) Omni-optimizer

0.4

0.6

0.8

1.0

0.0
0.0
2.0

1.5

1.5

1.5

1.0

1.0

1.0

0.5

0.5

0.5

0.4

0.6

0.8

1.0

0.0
0.0

0.2

(g) MOEA/D

0.4

0.6

0.4

0.6

0.8

1.0

0.4

0.6

0.8

1.0

0.8

1.0

(f) NSGA-II

2.0

0.2

0.2

(e) IBEA

2.0

0.0
0.0

0.2

(c) MO Ring PSO SCD

0.8

1.0

(h) NSGA-III

0.0
0.0

0.2

0.4

0.6

(i) MOEA/DD

Figure 4: Distribution of non-dominated solutions in the final population of each algo-rithm in the solution space on MMF2. The horizontal and vertical axis represent x1
and x2 , respectively.
2.0
1.5
1.0
0.5
0.0
0.0

0.2

0.4

0.6

0.8

1.0

Figure 5: Distribution of reference solutions in the solution space on MMF2.

19

×10−2

0
×10

2
−2

4
f1

6

8

8

6

4

2
f2

8
6
4 f3
2
0
0
×10−2

Figure 6: Pareto front of the three-objective Polygon problem.

NIMMO. One may be surprised that NSGA-II outperforms NSGA-III on most
MMaOP instances (in terms of the IGD metric). MOEA/DD also shows thepoor IGD values on MMaOPs. However, such a non-intuitive result has alreadybeen reported in [74, 5]. Since NSGA-III was specially designed for manyobjective optimization, NSGA-III can be outperformed by NSGA-II on two- andthree- objective problems. The poor performance of a many-objective optimizer
is also reported in [75]. The poor performance of NSGA-III (and MOEA/DD)
on the Polygon problems with more than three objectives is mainly due to their
Pareto front shapes. Fig. 6 shows the Pareto front of the three-objective Polygon
problem. While most multi-objective test problems (e.g., DTLZ [76] and WFG[77]) have triangular Pareto fronts, Polygon has an inverted triangular Pareto
front, as shown in Fig. 6. The analysis presented in [78] shows that somedecomposition-based MOEAs (including NSGA-III and MOEA/DD) perform
poorly on problems with inverted triangular Pareto fronts. This is because theshape of the distribution of weight vectors is different from the shape of thePareto front. For the same reason, our results show that NSGA-III performs
poorly on the Polygon problems in terms of IGD.
Tables 2 and 3 indicate that NIMMO achieves the best IGDX and PSP values
on all the MMaOPs, respectively. Thus, NIMMO is capable of locating multiple
equivalent Pareto optimal solutions on MMaOPs with a large number of objectives. According to the IGDX and PSP indicators, Omni-optimizer performs
the second best as a multi-modal multi-objective optimization algorithm.
Fig. 7 shows the distribution of non-dominated solutions in the final population of each algorithm in the solution space on the 15-Polygon problems, whose
Pareto optimal solutions are on the inside of the nine regular 15-sided polygons.
Results of a single run with a median IGDX value among 31 runs are shown.
Fig. 7 shows all algorithms (except for MOEA/D and MOEA/DD) locate thenine Pareto optimal solution subsets. While the solutions by NIMMO are almost evenly distributed, those by the other algorithms are biased. Fig. 7 (b)shows that the solutions in the decomposition-based TriMOEA-TA&R are distributed on the edges of the polygons. As explained above, this result is due tothe undesirable property of decomposition-based algorithms on problems with20

10 ×10−1

10 ×10−1

10 ×10

8

8

8

6

6

6

4

4

4

2

2

0
0

2

4

6

8

(a) NIMMO

10 ×10

10

×10−1

−1

2

0
0

2

4

6

8

(b) TriMOEA-TA&R

10 ×10

10

×10−1

−1

0
0
10 ×10

8

8

6

6

6

4

4

4

2

2

2

0
0

0
0

4

6

8

(c) Omni-optimizer

10 ×10

10

×10−1

−1

2

4

6

8

(d) IBEA

10 ×10

10

×10−1

−1

10 ×10
8

6

6

6

4

4

4

2

2

2

0
0

0
0

4

6

8

10

×10−1

2

4

6

(g) NSGA-III

8

10

×10−1

6

8

2

4

6

8

(e) NSGA-II

8(f) MOEA/D

4

10

×10−1

−1

0
0

8

2

2

(c) MO Ring PSO SCD

8

2

−1

10

×10−1

−1

0
0

2

4

6

(h) MOEA/DD

8

10

×10−1

Figure 7: Distribution of non-dominated solutions in the final population of each algo-rithm in the solution space on the 15-Polygon problem. The horizontal and verticalaxis represent x1 and x2 , respectively.

irregular Pareto fronts.
• Summary: Table 4 summarizes the average rankings of the nine algorithms
on the 23 problem instances by the Friedman test [79, 80]. We used the CONTROLTEST package (https://sci2s.ugr.es/sicidm) to calculate the rankings. Table 4 (a) shows that NIMMO is outperformed by Omni-optimizer,
IBEA, and NSGA-II regarding IGD. Tables 4 (b) and (c) show that NIMMO
performs the best among the nine algorithms regarding IGDX and PSP. Although TriMOEA-TA&R is an MMEA, it is outperformed by some MOEAs(e.g., NSGA-II and NSGA-III regarding PSP). As explained in Subsection 2.2,
TriMOEA-TA&R explicitly exploits the properties of position-related and distancerelated variables. However, the 23 problem instances used in our benchmarking
study do not have position-related and distance-related variables. It is alsoquestionable whether a real-world problem has such a special design variable.
In summary, NIMMO is the best multi-modal multi- and many-objective optimizer among the competitors.
5.2. Influence of the neighborhood size T on NIMMO
Since it is important to find an appropriate control parameter setting of anovel evolutionary algorithm, we investigate how the performance of NIMMO
is influenced by the value of T . Table 5 shows the average ranks of NIMMO

21

Table 4: Average rankings of the nine algorithms by the Friedman test. “MO Ring”

stands for “MO Ring PSO SCD”.
(a) IGD
Algorithm
NIMMO
TriMOEA-TA&R
MO RingOmni-optimizer
IBEA
NSGA-II
MOEA/D
NSGA-III
MOEA/DD(b) IGDX
Rank
4.22
8.22
6.87
3.70
3.17
2.83
5.91
4.30
5.78

Algorithm
NIMMO
TriMOEA-TA&R
MO RingOmni-optimizer
IBEA
NSGA-II
MOEA/D
NSGA-III
MOEA/DD(c) PSP
Rank
1.52
6.43
2.96
2.61
6.17
5.00
8.65
5.70
5.96

Algorithm
NIMMO
TriMOEA-TA&R
MO RingOmni-optimizer
IBEA
NSGA-II
MOEA/D
NSGA-III
MOEA/DD

Rank
1.72
6.04
3.16
2.64
6.32
4.76
8.56
5.72
6.08

Table 5: Average rankings of NIMMO with the 11 T values by the Friedman test.
(a) IGD
Algorithm
b0.05µc
b0.1µcb0.2µcb0.3µcb0.4µcb0.5µcb0.6µcb0.7µcb0.8µcb0.9µc
b1µc

Rank
8.12
7.76
7.88
7.64
7.24
6.84
5.80
4.76
4.04
3.32
2.60

(b) IGDX
Algorithm
b0.05µc
b0.1µcb0.2µcb0.3µcb0.4µcb0.5µcb0.6µcb0.7µcb0.8µcb0.9µc
b1µc

Rank
1.76
1.52
5.08
6.80
6.28
7.20
6.88
6.72
7.24
8.00
8.52

(c) PSP
Algorithm
b0.05µc
b0.1µcb0.2µcb0.3µcb0.4µcb0.5µcb0.6µcb0.7µcb0.8µcb0.9µc
b1µc

Rank
1.80
1.52
4.48
6.72
6.68
6.92
7.00
7.04
7.52
8.04
8.28

with T ∈ {b0.05µc, b0.1µc, b0.2µc, ..., b0.9µc, b1.0µc}. Tables S.1–S.3 in the supplementary file show detailed results of IGD, IGDX, and PSP, respectively.
The supplementary file is also available at the supplementary website (https://sites.google.com/view/nimmopt/).
As seen from Table 5 (a), NIMMO with T = b1.0µc performs the best interms of IGD. A large T value is generally beneficial for NIMMO on most problems. This is because a large T value relaxes the restriction of the niching
behavior and makes NIMMO more efficient for multi-objective optimization. In
fact, NIMMO with T = b1.0µc (i.e., T = µ) does not consider the solution
space diversity. NIMMO with T = µ is identical with IBEA used in this study(see Subsection 4.3). Tables 5 (b) and (c) show that NIMMO with T = b0.1µc
performs the best in terms of IGDX and PSP, respectively. As seen from Tables
S.2 and S.3 in the supplementary file, NIMMO with T = b0.1µc performs significantly better than NIMMO with T = b0.05µc on all the MMaOPs regarding
IGDX and PSP.
Fig. 8 shows the comparison of NIMMO with 11 values of T on six representative test problems regarding IGDX. Figs. 8 (a)–(c) show that a small T value
is beneficial for NIMMO. We do not claim that NIMMO with smaller T values
always performs better than NIMMO with larger T values. For example, Fig.
8 (c) shows that NIMMO with T = b0.6µc performs worse than NIMMO with22

0.045
0.040
0.035
0.030
0.025
0.2

0.4 0.6 0.8
T values

1(a) Two-On-One

Mean IGDX values

Mean IGDX values

0.014
0.012
0.010
0.008
0.006
0.2

0.4 0.6 0.8
T values(d) 3-Polygon

0.4 0.6 0.8
T values

1.0

3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0

(b) Omni-test

0.016

0.004
0

0.2

1

0.040
0.035
0.030
0.025
0.020
0.015
0.010
0.005
0

0.2

0.4 0.6 0.8
T values(e) 8-Polygon

0.2

0.4 0.6 0.8
T values

1(c) SYM-PART1

Mean IGDX values

0.020
0

2.2
2.1
2.0
1.9
1.8
1.7
1.6
1.5
0.0

Mean IGDX values

Mean IGDX values

Mean IGDX values

0.050

1

0.045
0.040
0.035
0.030
0.025
0.020
0.015
0.010
0.005
0

0.2

0.4 0.6 0.8
T values

1(f) 15-Polygon

Figure 8: Comparison of NIMMO with 11 values of T . The mean IGDX value at the
final iteration among 31 runs is shown.

T = b0.7µc on SYM-PART1). Interestingly, Figs. 8 (d)–(e) show that NIMMO
with T = b0.3µc has the worst performance regarding IGDX on the Polygon
problems with 3, 8, and 15 objectives. These results suggest that an appropriate specification of T is problem-dependent. T = b0.1µc is the best value forthe three Polygon problems whereas T = b0.05µc is the best for Two-On-One
and Omni-test in Fig. 8. NIMMO on the Polygon problems with 3, 8, and 15
objectives. Although the best T value in NIMMO is different for each problem,
T = b0.1µc is used for all experiments in this paper.
5.3. The influence of the population size µ on the performance of NIMMO
Here, we investigate the influence of µ on the performance of NIMMO. As
discussed in [81], it is difficult to compare MOEAs with different settings of µ
in a straightforward manner. Since a large number of solutions are beneficial
for some indicators (e.g., IGD and IGDX), a large µ value may be evaluated as
being better than a small µ value, regardless of their actual effectiveness. To
address this issue, we compare NIMMO with different µ values based only on apre-defined number of solutions µselected . In this manner, we can fairly compare
algorithms under the same solution set size. For each solution set found byan algorithm, we selected µselected solutions from the µ solutions in the final
population P using the distance-based selection method presented in [58].
At the beginning of the procedure, A is an empty set. First, a solution
is randomly selected from the final population P and stored into A. Then,

23

Table 6: Average rankings of NIMMO with 10 values of µ by the Friedman test.
(a) IGD
Algorithmµ = 100µ = 200µ = 300µ = 400µ = 500µ = 600µ = 700µ = 800µ = 900µ = 1000

Rank
6.04
3.76
4.24
4.60
4.28
4.92
5.72
6.60
6.96
7.88

(b) IGDX
Algorithmµ = 100µ = 200µ = 300µ = 400µ = 500µ = 600µ = 700µ = 800µ = 900µ = 1000

Rank
7.64
4.32
4.48
4.20
4.60
4.96
5.20
5.88
6.60
7.12

(c) PSP
Algorithmµ = 100µ = 200µ = 300µ = 400µ = 500µ = 600µ = 700µ = 800µ = 900µ = 1000

Rank
8.64
4.84
4.68
4.08
4.44
4.68
4.92
5.68
6.12
6.92

the solution with the maximum D(x, A) value is repeatedly added to A until|A| = µselected . Here, D(x, A) is the distance between a solution x in P and
its nearest solution in A in the normalized solution space. Finally, we calculate
the IGD, IGDX, and PSP indicator values of A, not P .
Table 6 shows the comparison of NIMMO with µ ∈ {100, 200, ..., 900, 1 000}.
It should be noted that we can arbitrarily specify the population size µ in
NIMMO. We set µselected to 100. Tables S.4–S.6 in the supplementary file show
detailed results of IGD, IGDX, and PSP. Table 6 (a) shows that NIMMO with
µ = 200 performs the best regarding IGD. Tables 6 (b) and (c) show that
NIMMO with µ = 400 has the best performance regarding IGDX and PSP,
respectively. These results suggest that a little larger µ value is suitable for
NIMMO for multi-modal multi-objective optimization. However, Table 6 showsthat NIMMO with µ ≥ 700 performs poorly regarding all three indicators. In
summary, any µ value in 200 ≤ µ ≤ 600 is suitable for NIMMO.
5.4. The performance of NIMMO with other fitness assignment schemes
Although NIMMO uses the I+ indicator to assign the fitness values to individuals in the population, NIMMO can use any indicator-based fitness assignment scheme. Here, we investigate the performance of NIMMO with thefollowing three fitness assignment schemes: the hypervolume difference indicator (IHD ) [7], the binary R2 indicator (IR2 ) [28], and the ranking method in
SRA (SRA) [32]. For IHD , IR2 , and SRA, see Subsection 2.1.
We denote NIMMO with the three fitness assignment schemes as NIMMOIHD , NIMMO-IR2 , and NIMMO-SRA, respectively. In NIMMO-IHD and NIMMOIR2 , I+ in (3) is replaced with IHD and IR2 , respectively. In NIMMO-SRA,
the fitness assignment in (2) is replaced with the ranking method in SRA. Thus,
the ranks of individuals are their fitness values in NIMMO-SRA. As in [32], thecontrol parameter pc in SRA is randomly generated in the range [0.4, 0.6] foreach iteration.
Table 7 summarizes the comparison results of the four variants of NIMMO.
NIMMO-I+ is the original version of NIMMO described in Section 3. Tables
S.7–S.9 in the supplementary file show detailed results of IGD, IGDX, and PSP,
respectively. Note that we do not intend to generalize the results of NIMMO
24

Table 7: Average rankings of the four variants of NIMMO by the Friedman test.
(a) IGD
Algorithm
NIMMO-INIMMO-IHD
NIMMO-IR2
NIMMO-SRA(b) IGDX
Rank
1.83
3.26
3.35
1.57

Algorithm
NIMMO-INIMMO-IHD
NIMMO-IR2
NIMMO-SRA(c) PSP
Rank
1.61
3.22
3.35
1.83

Algorithm
NIMMO-INIMMO-IHD
NIMMO-IR2
NIMMO-SRA

Rank
1.65
3.26
3.35
1.74

with the four fitness assignment schemes. Even though NIMMO-I+ performs
better than NIMMO-IR2 , it does not mean that I+ performs better than IR2 .
It only means that I+ is more suitable in NIMMO than IR2 for multi-modal
many-objective optimization.
Table 7 (a) shows that NIMMO-SRA performs the best regarding IGD.
Thus, NIMMO-SRA is the best configuration for multi-objective optimization.
The good performance of NIMMO-SRA regarding IGD is consistent with theresults reported in [32]. In contrast, Tables 7 (b) and (c) show that NIMMO-I+
performs the best regarding IGDX and PSP, respectively. Although NIMMOIHD and NIMMO-SRA outperforms NIMMO-I+ on some test problems (e.g.,some MMF problems) regarding IGDX and PSP, NIMMO-I+ has the best
performance on the Polygon problems with up to 15 objectives. I+ and thesolution space diversity maintenance mechanism in NIMMO may be compatible.
In summary, our results indicate that I+ is most suitable in NIMMO for multimodal multi-objective optimization.
6. Conclusion
We proposed NIMMO for multi-modal many-objective optimization. The
performance of NIMMO was investigated on various problems. Our results
showed that NIMMO is capable of finding multiple equivalent Pareto optimal
solutions on MMaOPs with up to 15 objectives.
We uploaded source code of NIMMO at the supplementary website (https://sites.google.com/view/nimmopt/). Source code of the 23 test problem
instances used in our benchmarking study is also available at the supplementary
website. Although MMOPs have not been well studied [18, 19], we hope that
this paper and its supplement encourage further development of efficient multimodal multi- and many-objective optimizers.
It was shown in [82] that there are multiple equivalent Pareto optimal solutions in multi-objective knapsack problems. An analysis of the performance of
NIMMO on combinatorial MMOPs and MMaOPs is an avenue for future work.
Acknowledgement
This work was supported by National Natural Science Foundation of China(Grant No. 61876075), the Program for Guangdong Introducing Innovative andEnterpreneurial Teams (Grant No. 2017ZT07X386), Shenzhen Peacock Plan
25

(Grant No. KQTD2016112514355531), the Science and Technology Innovation
Committee Foundation of Shenzhen (Grant No. ZDSYS201703031748284), and
the Program for University Key Laboratory of Guangdong Province (Grant No.
2017KSYS008).
References[1] K. Deb, Multi-Objective Optimization Using Evolutionary Algorithms,
John Wiley & Sons, 2001.
[2] K. Deb, S. Agrawal, A. Pratap, T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE TEVC 6 (2) (2002) 182–197.
[3] T. Wagner, N. Beume, B. Naujoks, Pareto-, Aggregation-, and IndicatorBased Methods in Many-Objective Optimization, in: EMO, 2007, pp. 742–
756.
[4] B. Li, J. Li, K. Tang, X. Yao, Many-objective evolutionary algorithms: A
survey, ACM C. S. 48 (1) (2015) 13:1–13:35.
[5] R. Tanabe, H. Ishibuchi, A. Oyama, Benchmarking Multi- and Manyobjective Evolutionary Algorithms under Two Optimization Scenarios,
IEEE Access 5 (2017) 19597–19619.
[6] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, V. G. da Fonseca, Performance assessment of multiobjective optimizers: an analysis and review,
IEEE TEVC 7 (2) (2003) 117–132.
[7] E. Zitzler, S. Künzli, Indicator-based selection in multiobjective search, in:
PPSN, 2004, pp. 832–842.
[8] N. Beume, B. Naujoks, M. T. M. Emmerich, SMS-EMOA: multiobjective
selection based on dominated hypervolume, EJOR 181 (3) (2007) 1653–
1669.
[9] J. Bader, E. Zitzler, HypE: An Algorithm for Fast Hypervolume-Based
Many-Objective Optimization, Evol. Comput. 19 (1) (2011) 45–76.
[10] G. Rudolph, B. Naujoks, M. Preuss, Capabilities of EMOA to Detect andPreserve Equivalent Pareto Subsets, in: EMO, 2007, pp. 36–50.
[11] O. Schütze, M. Vasile, C. A. C. Coello, Computing the Set of EpsilonEfficient Solutions in Multiobjective Space Mission Design, JACIC 8 (3)(2011) 53–70.
[12] F. Kudo, T. Yoshikawa, T. Furuhashi, A study on analysis of design variables in Pareto solutions for conceptual design optimization problem ofhybrid rocket engine, in: IEEE CEC, 2011, pp. 2558–2562.

26

[13] M. Sebag, N. Tarrisson, O. Teytaud, J. Lefèvre, S. Baillet, A MultiObjective Multi-Modal Optimization Approach for Mining Stable SpatioTemporal Patterns, in: IJCAI, 2005, pp. 859–864.
[14] A. Zhou, Q. Zhang, Y. Jin, Approximating the Set of Pareto-Optimal Solutions in Both the Decision and Objective Spaces by an Estimation ofDistribution Algorithm, IEEE TEVC 13 (5) (2009) 1167–1189.
[15] K. Deb, S. Tiwari, Omni-optimizer: A generic evolutionary algorithm forsingle and multi-objective optimization, EJOR 185 (3) (2008) 1062–1087.
[16] O. M. Shir, M. Preuss, B. Naujoks, M. T. M. Emmerich, Enhancing Decision Space Diversity in Evolutionary Multiobjective Algorithms, in: EMO,
2009, pp. 95–109.
[17] T. Hiroyasu, S. Nakayama, 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.
[18] X. Li, M. G. Epitropakis, K. Deb, A. P. Engelbrecht, Seeking Multiple
Solutions: An Updated Survey on Niching Methods and Their Applications,
IEEE TEVC 21 (4) (2017) 518–538.
[19] R. Tanabe, H. Ishibuchi, A Review of Evolutionary Multi-modal Multiobjective Optimization, IEEE TEVC (2019 (in press)) .
[20] E. Zitzler, M. Laumanns, L. Thiele, SPEA2: Improving the Strength Pareto
Evolutionary Algorithm, Tech. Rep. 103, ETHZ (2001).
[21] S. Das, S. Maity, B. Qu, P. N. Suganthan, Real-parameter evolutionary
multimodal optimization - A survey of the state-of-the-art, SWEVO 1 (2)(2011) 71–88.
[22] C. Yue, B. Qu, J. Liang, A Multiobjective Particle Swarm Optimizer Using Ring Topology for Solving Multimodal Multiobjective Problems, IEEETEVC 22 (5) (2018) 805–817.
[23] Y. Liu, G. G. Yen, D. Gong, A Multi-Modal Multi-Objective Evolutionary
Algorithm Using Two-Archive and Recombination Strategies, IEEE TEVC
(2019 (in press)) .
[24] M. Li, X. Yao, Quality Evaluation of Solution Sets in Multiobjective Optimisation: A Survey, ACM CSUR (2019 (in press)) .
[25] D. A. V. Veldhuizen, G. B. Lamont, Multiobjective Evolutionary Algorithm
Research: A History and Analysis, Tech. Rep. TR-98-03, Air Force Instituteof Technology (1998).
[26] E. Zitzler, L. Thiele, Multiobjective Optimization Using Evolutionary Algorithms - A Comparative Case Study, in: PPSN, 1998, pp. 292–304.

27

[27] M. P. Hansen, A. Jaszkiewicz, Evaluating the quality of approximations of
the non-dominated set, Tech. Rep. IMM-REP-1998-7, Technical Universityof Denmark (1998).
[28] D. H. Phan, J. Suzuki, R2-IBEA: R2 indicator based evolutionary algorithm for multiobjective optimization, in: IEEE CEC, 2013, pp. 1836–1845.
doi:10.1109/CEC.2013.6557783.
[29] R. H. Gómez, C. A. C. Coello, Improved Metaheuristic Based on the R2
Indicator for Many-Objective Optimization, in: GECCO, 2015, pp. 679–686.
[30] W. Li, E. Özcan, R. John, J. H. Drake, A. Neumann, M. Wagner, A modified indicator-based evolutionary algorithm (mIBEA), in: IEEE CEC, 2017,
pp. 1047–1054.
[31] H. Wang, L. Jiao, X. Yao, Two Arch2: An Improved Two-Archive Algorithm for Many-Objective Optimization, IEEE TEVC 19 (4) (2015) 524–
541.
[32] B. Li, K. Tang, J. Li, X. Yao, Stochastic Ranking Algorithm for ManyObjective Optimization Based on Multiple Indicators, IEEE TEVC 20 (6)(2016) 924–938.
[33] M. Li, S. Yang, X. Liu, Shift-Based Density Estimation for Pareto-Based
Algorithms in Many-Objective Optimization, IEEE TEVC 18 (3) (2014)
348–365.
[34] T. P. Runarsson, X. Yao, Stochastic ranking for constrained evolutionary
optimization, IEEE TEVC 4 (3) (2000) 284–294.
[35] M. Li, S. Yang, X. Liu, Pareto or Non-Pareto: Bi-Criterion Evolution in
Multiobjective Optimization, IEEE TEVC 20 (5) (2016) 645–665.
[36] Q. Zhang, H. Li, MOEA/D: A multiobjective evolutionary algorithm basedon decomposition, IEEE TEVC 11 (6) (2007) 712–731.
[37] D. Brockhoff, T. Wagner, H. Trautmann, R2 Indicator-Based Multiobjective Search, Evol. Comput. 23 (3) (2015) 369–395.
[38] R. Denysiuk, L. A. Costa, I. A. Espı́rito-Santo, Many-objective optimization using differential evolution with variable-wise mutation restriction, in:
GECCO, 2013, pp. 591–598.
[39] Y. Tian, R. Cheng, X. Zhang, F. Cheng, Y. Jin, An Indicator-Based Multiobjective Evolutionary Algorithm With Reference Point Adaptation for
Better Versatility, IEEE TEVC 22 (4) (2018) 609–622.
[40] S. Jiang, J. Zhang, Y. Ong, A. N. Zhang, P. S. Tan, A Simple and Fast Hypervolume Indicator-Based Multiobjective Evolutionary Algorithm, IEEE
Trans. Cyber. 45 (10) (2015) 2202–2213.
28

[41] S. Z. Martı́nez, A. L. Jaimes, A. Garcı́a-Nájera, LIBEA: A Lebesgue
Indicator-Based Evolutionary Algorithm for multi-objective optimization,
SWEVO 44 (2019) 404–419.
[42] A. Menchaca-Mendez, C. A. C. Coello, GD-MOEA: A New Multi-Objective
Evolutionary Algorithm Based on the Generational Distance Indicator, in:
EMO, 2015, pp. 156–170.
[43] E. M. Lopez, C. A. C. Coello, IGD+ -EMOA: A multi-objective evolutionary algorithm based on IGD+ , in: IEEE CEC, 2016, pp. 999–1006.
[44] H. Ishibuchi, H. Masuda, Y. Tanigaki, Y. Nojima, Modified Distance Calculation in Generational Distance and Inverted Generational Distance, in:
EMO, 2015, pp. 110–125.
[45] C. A. R. Villalobos, C. A. C. Coello, A new multi-objective evolutionary
algorithm based on a performance assessment indicator, in: GECCO, 2012,
pp. 505–512.
[46] O. Schütze, X. Esquivel, A. Lara, C. A. C. Coello, Using the Averaged
Hausdorff Distance as a Performance Measure in Evolutionary Multiobjective Optimization, IEEE TEVC 16 (4) (2012) 504–522.
[47] C. A. C. Coello, M. R. Sierra, A Study of the Parallelization of a Coevolutionary Multi-objective Evolutionary Algorithm, in: MICAI, 2004, pp.
688–697.
[48] K. Bringmann, T. Friedrich, F. Neumann, M. Wagner, ApproximationGuided Evolutionary Multi-Objective Optimization, in: IJCAI, 2011, pp.
1198–1203.
[49] J. G. Falcón-Cardona, C. A. C. Coello, A multi-objective evolutionary
hyper-heuristic based on multiple indicator-based density estimators, in:
GECCO, 2018, pp. 633–640.
[50] T. Ulrich, J. Bader, E. Zitzler, Integrating decision space diversity into
hypervolume-based multiobjective search, in: GECCO, 2010, pp. 455–462.
[51] J. W. Kruisselbrink, A. Aleman, M. T. M. Emmerich, A. P. IJzerman,
A. Bender, T. Bäck, E. V. D. Horst, Enhancing search space diversity
in multi-objective evolutionary drug molecule design using niching, in:
GECCO, 2009, pp. 217–224.
[52] F. Zeng, J. Decraene, M. Y. Low, S. Zhou, W. Cai, Evolving Optimal andDiversified Military Operational Plans for Computational Red Teaming,
IEEE Sys. Jour. 6 (3) (2012) 499–509.
[53] O. M. Shir, T. Bäck, Niching with derandomized evolution strategies in
artificial and real-world landscapes, Nat. Comput. 8 (1) (2009) 171–196.

29

[54] X. Li, Niching Without Niching Parameters: Particle Swarm Optimization
Using a Ring Topology, IEEE TEVC 14 (1) (2010) 150–169.
[55] J. Liang, Q. Guo, C. Yue, B. Qu, K. Yu, A Self-organizing Multi-objective
Particle Swarm Optimization Algorithm for Multimodal Multi-objective
Problems, in: ICSI, 2018, pp. 550–560.
[56] M. Kim, T. Hiroyasu, M. Miki, S. Watanabe, SPEA2+: Improving the
Performance of the Strength Pareto Evolutionary Algorithm 2, in: PPSN,
2004, pp. 742–751.
[57] J. J. Liang, C. T. Yue, B. Y. Qu, Multimodal multi-objective optimization:
A preliminary study, in: IEEE CEC, 2016, pp. 2454–2461.
[58] R. Tanabe, H. Ishibuchi, A Decomposition-Based Evolutionary Algorithm
for Multi-modal Multi-objective Optimization, in: PPSN, 2018, pp. 249–
261.
[59] X. Ma, F. Liu, Y. Qi, X. Wang, L. Li, L. Jiao, M. Yin, M. Gong, A Multiobjective Evolutionary Algorithm Based on Decision Variable Analyses forMultiobjective Optimization Problems With Large-Scale Variables, IEEETEVC 20 (2) (2016) 275–298.
[60] R. Cheng, Y. Jin, M. Olhofer, B. Sendhoff, A Reference Vector Guided
Evolutionary Algorithm for Many-Objective Optimization, IEEE TEVC
20 (5) (2016) 773–791.
[61] M. Basseur, E. K. Burke, Indicator-based multi-objective local search, in:IEEE CEC, 2007, pp. 3100–3107.
[62] D. E. Goldberg, J. Richardson, Genetic algorithms with sharing for multimodal function optimization, in: ICGA, 1987, pp. 41–49.
[63] K. D. Jong, An Analysis of the Behavior of a Class of Genetic Adaptive
Systems, Ph.D. thesis, University of Michigan (1975).
[64] X. Yin, N. Germay, A fast genetic algorithm with sharing scheme using
cluster analysis methods in multimodal function optimization, in: Artificial
Neural Nets and Genetic Algorithms, 1993, pp. 450–457.
[65] S. Bird, X. Li, Adaptively choosing niching parameters in a PSO, in:
GECCO, 2006, pp. 3–10.
[66] M. G. Epitropakis, X. Li, E. K. Burke, A dynamic archive niching differential evolution algorithm for multimodal optimization, in: IEEE CEC, 2013,
pp. 79–86.
[67] M. Preuss, Improved Topological Niching for Real-Valued Global Optimization, in: EvoApplications, 2012, pp. 386–395.

30

[68] A. Ahrari, K. Deb, M. Preuss, Multimodal Optimization by CovarianceMatrix Self-Adaptation Evolution Strategy with Repelling Subpopulations,Evol. Comput. 25 (3) (2017) 439–471.
[69] S. W. Mahfoud, Crowding and Preselection Revisited, in: PPSN, 1992, pp.
27–36.
[70] M. Preuss, B. Naujoks, G. Rudolph, Pareto Set and EMOA Behavior for
Simple Multimodal Multiobjective Functions, in: PPSN, 2006, pp. 513–522.
[71] H. Ishibuchi, N. Akedo, Y. Nojima, A many-objective test problem forvisually examining diversity maintenance behavior in a decision space, in:
GECCO, 2011, pp. 649–656.
[72] K. Deb, H. Jain, An evolutionary many-objective optimization algorithm
using reference-point-based nondominated sorting approach, part I: solving
problems with box constraints, IEEE TEVC 18 (4) (2014) 577–601.
[73] K. Li, K. Deb, Q. Zhang, S. Kwong, An Evolutionary Many-Objective
Optimization Algorithm Based on Dominance and Decomposition, IEEETEVC 19 (5) (2015) 694–716.
[74] H. Ishibuchi, R. Imada, Y. Setoguchi, Y. Nojima, Performance comparisonof NSGA-II and NSGA-III on various many-objective test problems, in:IEEE CEC, 2016, pp. 3045–3052.
[75] M. Wagner, F. Neumann, A fast approximation-guided evolutionary multiobjective algorithm, in: GECCO, 2013, pp. 687–694.
[76] K. Deb, L. Thiele, M. Laumanns, E. Zitzler, Scalable Test Problems for
Evolutionary Multi-Objective Optimization, in: Evolutionary Multiobjective Optimization. Theoretical Advances and Applications, Springer, 2005,
pp. 105–145.
[77] S. Huband, P. Hingston, L. Barone, R. L. While, A review of multiobjective
test problems and a scalable test problem toolkit, IEEE TEVC 10 (5) (2006)
477–506.
[78] H. Ishibuchi, Y. Setoguchi, H. Masuda, Y. Nojima, Performance ofDecomposition-Based Many-Objective Algorithms Strongly Depends on
Pareto Front Shapes, IEEE TEVC 21 (2) (2017) 169–190.
[79] S. Garcı́a, A. Fernández, J. Luengo, F. Herrera, Advanced nonparametrictests for multiple comparisons in the design of experiments in computational intelligence and data mining: Experimental analysis of power, Inf.
Sci. 180 (10) (2010) 2044–2064.
[80] J. Derrac, S. Garcı́a, D. Molina, F. Herrera, A practical tutorial on the useof nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms, SWEVO 1 (1) (2011) 3–18.
31

[81] H. Ishibuchi, H. Masuda, Y. Nojima, Comparing solution sets of different
size in evolutionary many-objective optimization, in: IEEE CEC, 2015, pp.
2859–2866.
[82] Y. Nojima, K. Narukawa, S. Kaige, H. Ishibuchi, Effects of Removing Overlapping Solutions on the Performance of the NSGA-II Algorithm, in: EMO,
2005, pp. 341–354.

32