1

A Framework to Handle Multi-modal
Multi-objective Optimization in
Decomposition-based Evolutionary Algorithms
Ryoji Tanabe, Member, IEEE,and Hisao Ishibuchi, Fellow, IEEE

Solution space

Abstract—Multi-modal multi-objective optimization is to locate(almost) equivalent Pareto optimal solutions as many as possible.
While decomposition-based evolutionary algorithms have good
performance for multi-objective optimization, they are likely
to perform poorly for multi-modal multi-objective optimization
due to the lack of mechanisms to maintain the solution spacediversity. To address this issue, this paper proposes a framework
to improve the performance of decomposition-based evolutionary
algorithms for multi-modal multi-objective optimization. Our
framework is based on three operations: assignment, deletion,
and addition operations. One or more individuals can be assigned
to the same subproblem to handle multiple equivalent solutions.
In each iteration, a child is assigned to a subproblem basedon its objective vector, i.e., its location in the objective space.
The child is compared with its neighbors in the solution spaceassigned to the same subproblem. The performance of improved
versions of six decomposition-based evolutionary algorithms by
our framework is evaluated on various test problems regarding
the number of objectives, decision variables, and equivalent
Pareto optimal solution sets. Results show that the improved
versions perform clearly better than their original algorithms.
Index Terms—Multi-modal multi-objective optimization,
decomposition-based evolutionary algorithms, reference vectorbased evolutionary algorithms, solution space diversity

I. I NTRODUCTIONM

ULTI-OBJECTIVE optimization problems (MOPs) appear in real-world applications. Since no solution x
can simultaneously minimize multiple objective functions ingeneral, the goal of MOPs is usually to find a Pareto optimal
solution preferred by a decision maker [1]. When the decision
maker’s preference information is unavailable a priori, an “aposteriori” decision making is conducted. The decision makerselects the final solution xfinal from a set of solutions thatapproximates the Pareto front in the objective space.
An evolutionary multi-objective optimization algorithm(EMOA) is frequently used for the “a posteriori” decision
making. Since EMOAs are population-based optimizers, they
are likely to find a set of solutions in a single run. A number
of EMOAs have been proposed in the literature. They can beclassified into the following three categories: dominance-based
EMOAs (e.g., NSGA-II [2]), indicator-based EMOAs (e.g.,
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)

Objective space

Fig. 1: Illustration of a situation where three solutions are almost
the same in the objective space but dissimilar in the solution space.

IBEA [3]), and decomposition-based EMOAs (e.g., MOEA/D[4] and NSGA-III [5]). In particular, decomposition-based
EMOAs have shown promising performance [6].
In the EMO community, it has been implicitly assumed
that the decision maker is interested only in the distributionof solutions in the objective space. Thus, the distributionof solutions in the solution space has not received muchattention. However, after the decision maker has selected thefinal solution xfinal based on its objective vector f (xfinal ),she/he may want to examine other dissimilar solutions withequivalent quality or slightly inferior quality [7]–[9].
Fig. 1 shows a situation where the four solutions xa , xb ,xc , and xd are far from each other in the solution space but
close to each other in the objective space. Although xd is
a dominated solution, it may be acceptable for the decision
maker. This is because the difference between f (xa ), f (xb ),f (xc ), and f (xd ) is small enough. If the decision makerhas obtained multiple dissimilar solutions with similar quality,
she/he can select xfinal according to her/his preference inthe solution space. For example, suppose that xa in Fig.
1 becomes unavailable due to some accident (e.g., material
shortages and mechanical failures) after the decision makerhas selected xa as xfinal . In such a case, she/he can substituteone of xb , xc , and xd for xa . Schütze et al. give a morepractical example on space mission design problems [9]. Sebag
et al. also demonstrate the importance of multiple equivalentsolutions for the decision maker on functional brain imaging
problems [7]. Previous studies on other real-world problems
with multiple equivalent solutions include diesel engine design
problems [10], distillation plant layout problems [11], androcket engine design problems [12].
A multi-modal multi-objective problem (MMOP) is to locate (almost) equivalent Pareto optimal solutions as many
as possible [7], [13]. On the one hand, it is sufficient for

2

MOPs to find one of xa , xb , and xc in Fig. 1 since their
quality is almost the same. On the other hand, all of xa ,xb , xc , and xd should be found for MMOPs. Since multiple
equivalent solutions play crucial roles in the reliable decision
making process as explained above, MMOPs are important inpractice. Evolutionary multi-modal multi-objective optimization algorithms (EMMAs) are specially designed optimizers
for MMOPs. Unlike EMOAs, EMMAs have mechanisms tohandle multiple equivalent solutions. Representative EMMAs
include PQ, -MOEA [9], Omni-optimizer [13], DIOP [14],Niching-CMA [15], and MO Ring PSO SCD [16].
In our earlier study [17], we proposed MOEA/D-AD, which
is a specially designed MOEA/D for MMOPs. MOEA/D-AD
can assign multiple individuals to each subproblem in orderto handle equivalent solutions. For each iteration, a child u isassigned to a subproblem whose weight vector is closest to its
objective vector f (u), in terms of the perpendicular distance.
Then, u is compared to the individuals assigned to the same
subproblem based on the weighted Tchebycheff function andthe contribution to the solution space diversity.
In this paper, we extend our earlier study [17] and propose ageneral framework (called ADA) to handle multi-modal multiobjective optimization in a variety of decomposition-based
EMOAs (including reference-based algorithms such NSGAIII [5] and θ-DEA [18]). The proposed ADA frameworkuses assignment, deletion, and addition operations. Whiledecomposition-based EMOAs work well for multi-objective
optimization, they are likely to perform poorly for multimodal multi-objective optimization. This is because most ofthem do not have mechanisms to maintain the solution spacediversity. In this paper, we show that efficient EMMAs could
be realized by incorporating a solution space diversity maintenance mechanism into decomposition-based EMOAs. The proposed ADA framework facilitates the solution space diversity
maintenance of existing EMOAs. We incorporate ADA into
the following six representative decomposition-based EMOAs:
MOEA/D-AGR [19], MOEA/D-DU [20], eMOEA/D [21],NSGA-III [5], θ-DEA [18], and RVEA [22]. We examine
the ability of the ADA versions of those EMOAs (MOEA/DAGR-ADA, MOEA/D-DU-ADA, eMOEA/D-ADA, NSGAIII-ADA, θ-DEA-ADA, and RVEA-ADA) to locate multiple
equivalent Pareto optimal solutions on various test problems.
We also analyze the behavior of those ADA versions.
The differences between this paper and our previous study[17] are as follows:
1) We propose ADA. While MOEA/D-AD [17] is a single
algorithm for MMOPs, ADA is a framework to improve
the performance of the existing decomposition-based
EMOAs for MMOPs. In contrast to MOEA/D-AD, eachconfiguration (e.g., the assignment operation) in ADA
depends on an EMOA to be combined.
2) We demonstrate that ADA can be combined with thesix decomposition-based EMOAs, including so-calledreference vector-based EMOAs. MOEA/D-AD does nothave such flexibility.
3) We introduce a practical two-phase decision making
method with a solution set found by an ADA-based algorithm. We also introduce two post-processing methods

for benchmarking an ADA-based algorithm.
4) While we used only two-objective and two-variableproblems in [17], we use problems with a various number of objectives M , design variables D, and equivalent
Pareto optimal solution subsets O to investigate thescalability of EMMAs. Such a scale-up study has notbeen performed in the literature. We also use problems
with distance-related variables [23], [24].
5) We compare the six ADA-based algorithms to threestate-of-the-art EMMAs, including TriMOEA-TA&R[24] proposed after the publication of [17].
The rest of this paper is organized as follows. Section II provides some preliminaries of this paper. Section III introduces
ADA and explains how to incorporate it into decompositionbased EMOAs. Section IV describes the experimental setup.
Section V examines the performance of the six ADA variants.
Section VI concludes this paper.
II. P RELIMINARIESFirst, Subsections II-A and II-B give definitions of MOPs
and MMOPs, respectively. Then, Subsections II-C and II-D
describe decomposition-based EMOAs and reference vectorbased EMOAs, respectively. Subsections II-C and II-D mainly
explain components of the six existing EMOAs used in ADA.
Algorithms S.1–S.6 in the supplementary file provide theoverall procedures of the six EMOAs. Section III introduces
how each component is used in ADA.
A. Definition of MOPs
A continuous MOP is to find a solution x ∈ S ⊆ RD
that minimizes a given
vector f : S →QDobjective function
max
RM , x 7→ f (x). S = j=1 [xmin
,
x
]
is
the
D-dimensionaljjmax
solution space where xmin
≤
x
≤
x
for
each
index j ∈jjj
{1, ..., D}. RM is the M -dimensional objective space.
A solution x1 is said to dominate x2 if and only if 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 iscalled 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 for the “a posteriori” decision making
is to find a non-dominated solution set that approximates thePareto front in the objective space.
B. Definition of MMOPs
Although the term “multi-modal multi-objective optimization” was firstly coined in [7], [25] in 2005, the definitionof MMOPs has not been explicitly given in most previous
studies. We consider the following two types of MMOPs: (i)Type1-MMOP is to locate all Pareto optimal solutions. (ii)
Type2-MMOP is to locate all Pareto optimal solutions and nonPareto optimal solutions which have acceptable quality for thedecision maker. Those non-Pareto optimal solutions should be
far from the Pareto optimal solutions and the other non-Pareto
optimal solutions in the solution space.
For example, xa , xb , and xc in Fig. 1 should be found for
Type1-MMOPs. In addition, the non-Pareto optimal solution

3

xd should be found for Type2-MMOPs if its quality f (xd ) is
acceptable to the decision maker. While most existing studies(e.g., [13], [16]) assume Type1-MMOPs, only a few studies(e.g., [7], [9], [14]) address Type2-MMOPs. Although there
is room for discussion, Type2-MMOPs may be more practicalthan Type1-MMOPs. Diverse solutions with similar quality to
the final solution xfinal are beneficial to the decision makereven if their quality is slightly worse than xfinal [7], [9].
Nevertheless, we focus only on Type1-MMOPs in thispaper. The main reason is due to the difficulty in benchmarking
EMMAs on Type2-MMOPs. In contrast to Type1-MMOPs,
Type2-MMOPs are loosely defined. This is because the terms“acceptable quality” and “far from” in Type2-MMOPs significantly depend on the decision maker. It is difficult to define thetwo key factors in Type2-MMOPs for benchmarking purposesin a fair manner. How to evaluate the performance of EMMAs
for Type2-MMOPs itself is another research topic.
C. Decomposition-based EMOAs (MOEA/D-type algorithms)
Although some frameworks of decomposition-based
EMOAs have been proposed in the literature, MOEA/D-type
algorithms show the promising performance [6]. Here, we
describe MOEA/D-type algorithms. First, we explain the mostbasic MOEA/D [4] and its modified version called MOEA/DDE [26]. The framework of MOEA/D-DE is used in most
MOEA/D-type algorithms. Then, we describe componentsof three MOEA/D-type algorithms: MOEA/D-AGR [19],MOEA/D-DU [20], and eMOEA/D [21].
1) MOEA/D: The most basic MOEA/D [4] decomposes an
M -objective MOP into N single-objective subproblems usinga scalarizing function g : RM → R and a set of uniformly
distributed weight vectors W = {w1 , ..., w
PNM}. For eachi ∈ {1, ..., N }, wi = (wi,1 , ..., wi,M )T , and k=1 wi,k = 1.
The j-th individual xj in the population P is assigned to the jth subproblem with wj . Thus, the population size µ is always
equal to the number of weight vectors N . The j-th subproblem
also has its neighborhood index list B j = {bj,1 , ..., bj,S },which consists of indices of the S closest weight vectors to
wj in the weight vector space.
After the initialization of P and W , the following steps arerepeatedly performed until a termination condition is satisfied.
In each iteration t, for the j-th subproblem (j ∈ {1, ..., N }),an index list T is set to B j . Two indices a and b of theparent individuals xa and xb are randomly selected from T .
u is generated by applying variation operators to xa and xb .
The SBX crossover and the polynomial mutation [27] areused in the original MOEA/D. After u has been generated,
the replacement selection is applied to each neighborhood
subproblem k ∈ T . The current solution xk of the k-th
subproblem is replaced with u if g(u|wk ) ≤ g(xk |wk ).
2) MOEA/D-DE: The differences between MOEA/D and
MOEA/D-DE are threefold. First, the differential evolution(DE) operator [28] is used instead of the SBX crossover [27].
Second, two types of index lists (B j and {1, ..., N }) are used.
The index list T is set to B j with a probability of δ ∈ [0, 1]
and {1, ..., N } with a probability of 1 − δ. Thus, T can be
set to all individual indices. Third, the maximum number ofindividuals replaced by the child u is restricted to nrep .

3) MOEA/D-AGR: The difference between MOEA/D-AGR
and MOEA/D-DE is caused by a replacement method. In
MOEA/D-AGR, indices of subproblems to be updated by a
newly generated solution u are set to K neighborhood indicesof the j-th subproblem in the weight vector space where
j is the subproblem index whose scalarizing function valueg(u|wj ) is the best among all N subproblems:
j = argmin g(u|wk ) .
(1)
k∈{1,...,N }The weighted Tchebycheff function g tch is used in (1):g tch (x|w) = max
wi |fi (x) − zi∗ | ,(2)
i∈{1,...,M }∗∗ T= (z1∗ , ..., zM)

where zis the ideal point. Since finding
the true ideal point is difficult in general, its approximation isused in (2). The i-th element of the approximation of z ∗ isthe minimum objective value found during the search process.
The replacement neighborhood size K plays a crucial role inbalancing exploration and exploitation in MOEA/D-AGR in asimilar manner to S in MOEA/D. A large K value encourages
exploitation. In MOEA/D-AGR, the K value deterministicallyincreases with the number of iterations. The results presented
in [19] show that a scheduling method based on the sigmoid
function is suitable for MOEA/D-AGR.
4) MOEA/D-DU: In contrast to MOEA/D-AGR, the selection of subproblems that need to be updated is based on
the distance in the objective and weight vector spaces in
MOEA/D-DU. In the normalized objective space, the perpendicular distance between the normalized objective vector
f 0 (u) and wj is calculated for each j ∈ {1, ..., N }. The replacement is performed for K subproblems with the minimum
perpendicular distance. Similar to MOEA/D-AGR, K controls
the balance between exploration and exploitation.
MOEA/D-DU uses the division version of the weighted
Tchebycheff function g dtch [29]:|fi (x) − zi∗ |
dtch
,
(3)
g
(x|w) = maxwii∈{1,...,M }if wi = 0, it is set to 10−6 to avoid division by zero. It isreported in [29] that the distribution of the search directionsof MOEA/D with g dtch is more uniform than that with g tch .
5) eMOEA/D: In eMOEA/D, the replacement method isapplied to K subproblems with the minimum scalarizing
function values. MOEA/D-AGR and eMOEA/D are the same
regarding the use of scalarizing function values to selectsubproblems. The following multiplicative scalarizing function(MSF) [21] is used in eMOEA/D:1+αmaxi∈{1,...,M } w1i |fi (x − zi∗ )|α , (4)g msf (x|w) = mini∈{1,...,M } w1i |fi (x − zi∗ )|where α controls the size of the so-called improvement region.
α plays a similar role in the penalty value θ of the PBI functionin (7), which will be explained later. Even if f (x) is far
from wj regarding the the perpendicular distance, x would
be evaluated as being superior using a sufficiently large α
value. In (4), g msf with α = 0 is identical to g dtch .

4According to a general rule of thumb “emphasize diversity
and convergence at the early and later stages, respectively”, αdecreases linearly with the number of iterations t:

 

t
α = β 1 − maxM
min {wk } ,(5)tk∈{1,...,M }where tmax is the maximum number of iterations. The recommended value of β is 1.
D. Reference vector-based EMOAs
Representative reference vector-based EMOAs include
NSGA-III [5], θ-DEA [18], RVEA [22], VaEA [30], and
SPEA/R [31]. While w ∈ W is called the weight vector indecomposition-based EMOAs, it is referred to the reference
vector in reference vector-based EMOAs. Below, we explain
NSGA-III, θ-DEA, and RVEA.
1) NSGA-III: NSGA-III is an improved version of NSGAII [2] for many-objective optimization. For each iteration t,
children Q are generated by applying variation operators
to randomly selected pairs of individuals from P . In theenvironmental selection, µ individuals for the next iterationt + 1 are selected from the union of P and Q. The primaryand secondary criteria are based on the non-domination levels
and the reference vector-based niching method, respectively.
Individuals in the union P ∪ Q are grouped as F 1 , F 2 , ...
according to their non-domination levels. First, P and the front
index i are initialized as P = ∅ and i = 1, respectively.
Then, individuals in F i are added to P and i is incremented
until |P | + |F i | ≥ µ. After this operation, P = F 1 ∪ ... ∪F l−1 , where l is the index of the last front. If |P | < µ, otherµ − |P | individuals are selected from the last front F l using
the niching method described below. At the beginning of theniching procedure, an individual x in P and F l is assigned tothe j-th subproblem with the minimum perpendicular distance
between its normalized objective vector f 0 (x) and wj :no(6)
j = argmin PD f 0 (x), wk ,k∈{1,...,N }where the function PD returns the perpendicular distance
between two input vectors. After the assignments of all individuals, other individuals in the next iteration are selectedfrom F l based on the number of individuals assigned to each
subproblem and their perpendicular distance.
2) θ-DEA: The environmental selection in θ-DEA is similar
to that of NSGA-III. However, the secondary criterion in
θ-DEA is based on the so-called θ-dominance. First, eachindividual in the union of P and F l is assigned to the jth subproblem using (6). Then, individuals assigned to each
subproblem are ranked based on their θ-dominance levels.
Let us assume that two individuals x and y are assigned
to the same j-th subproblem. x is said to θ-dominate y if
g pbi (x|wj ) < g pbi (y|wj ). The PBI function g pbi is given as:g pbi (x|w) = d1 + θ d2 ,

(7)
∗ T

k (f (x) − z ) wk
,
kwk


w
∗
d2 = f (x) − z + d1
,
kwk

d1 =

(8)
(9)

where kak indicates the Euclidean norm of a. The distance d1
represents how close the objective vector f (x) is to the Pareto
front, and d2 is the perpendicular distance between f (x) andw. The penalty parameter θ balances the convergence (d1 )and the diversity (d2 ). The recommended θ value is 106 forM vectors w with the axis directions (e.g., (1, 0, ..., 0)T ) and
5 for all the other vectors w.
3) RVEA: RVEA uses a set of unit reference vectors
V = {v 1 , ..., v N }, instead of a set of reference vectors W ={w1 , ..., wN }, where v i = wi /kwi k for each i ∈ {1, ..., N }.
RVEA adaptively adjusts V based on the current P .
After children Q have been generated, the environmental
selection is applied to the union of P and Q. First, for eachindividual in P ∪ Q, f (x) is transformed as f 0 (x) = f (x) −z ∗ . Then, x in P ∪ Q is assigned to the j-th subproblem withthe minimum angle between f 0 (x) and v j :
n
o
(10)
j = argmin angle f 0 (x), v k ,k∈{1,...,N }where the function angle(a, b) in (10) returns the angle
between the two input vectors a and b.
Then, individuals assigned to the same subproblem are compared based on their angle-penalized distance (APD) values.
For each j ∈ {1, ..., N }, the best individual with the minimum
APD value can survive to the next iteration. The APD valueof x with the unit reference vector v is given as:APD(x) = 1 + P (x, v) kf 0 (x)k,
(11)
!α
!
0
angle f (x), vt,(12)
P (x, v) = M
max
t
γ(v)

γ(v) = minangle(v, s) ,(13)
s∈V \{v }where P (x, v) is a penalty value for x. The larger the angle
between f 0 (x) and v is, the larger the penalty value is given tox. γ(v) is used to normalize the angle. t is the current numberof iterations, and tmax is the maximum number of iterations.
The influence of the angle-based penalty scheme increases as
the search progresses. The recommended setting of α is 2.
III. P ROPOSED ADA FRAMEWORKThis section explains the proposed ADA framework and thesix ADA versions. Unlike EMOAs for MOPs, EMMAs for
MMOPs need to maintain the diversity of the population inboth the objective and solution spaces. For example, Omnioptimizer [13] uses an aggregate crowding distance metricin the objective and solution spaces. While most EMMAs(e.g., [14]–[16]) aggregate the objective and solution spacediversity metrics similar to Omni-optimizer, ADA uses themin a two-phase manner similar to TriMOEA-TA&R [24]. Morespecifically, ADA handles the objective space diversity by anassignment method in an original EMOA and the solution
space diversity by a simple niching criterion.
On the one hand, a single individual is assigned to each
subproblem in most decomposition-based EMOAs. Thus, thepopulation size µ is equal to the number of weight/reference
vectors N (i.e., µ = N ). On the other hand, one or moreindividuals can be assigned to each subproblem in the ADA

5framework (i.e., µ ≥ N ). This mechanism is to maintainmultiple equivalent individuals in each subproblem. The µvalue is adaptively adjusted in ADA.
Algorithm 1 shows the ADA framework. Algorithms S.7–
S.12 in the supplementary file show the six ADA-based algorithms. Algorithms S.7–S.12 are almost the same to Algorithm
1. Lines 3, 11, and 16 are different between Algorithms S.7–
S.12 and Algorithm 1. At the beginning of the search, µ is
set to N (line 1). The population P and the weight/reference
vector set W are also initialized. The i-th individual xi in P is
assigned to the i-th subproblem (lines 2–3). This operation is
unnecessary for those EMOAs with a reassignment procedureof individuals to subproblems such as RVEA. After the child u
has been generated by the reproduction operations (line 6–8),
the population P is updated using the assignment (line 10),deletion (lines 16–17), and addition (lines 18–19) operations.
The assignment and deletion operations differ depending on
an EMOA to be combined with the ADA framework. Afterthe normalization of the µ + 1 objective vectors (line 9), u isassigned to the j-th subproblem (line 10). X denotes a set ofindividuals that have been assigned to the j-th subproblem andare in the neighborhood of u in the solution space (line 11).
Thus, ADA requires a neighborhood criterion in the solution
space. X is explained later in detail using Fig. 2. Two Boolean
variables bwinner and bexplorer (line 12) are used in the addition
operation (lines 18–19). u enters P if it satisfies either of thefollowing two conditions. One is that there is no individual in
X (lines 13–14). The other is that u is better than at least one
individual in X in the deletion operation (lines 15–17).
The following Subsections (from III-A to III-F) explain
each step of ADA in detail. Subsection III-G presents aneffective decision making based on a solution set found by
the proposed approach. Apart from the decision making,Subsection III-H introduces two post-processing methods forbenchmarking. Subsection III-I discusses the applicability of
ADA. Subsection III-J discusses the originality of ADA.
A. Reproduction operation
ADA uses the basic GA operators (i.e., the SBX crossoverand the polynomial mutation [27]) and the simplest methodof selecting parents. First, two parents xa and xb are randomly selected from P such that a 6= b, regardless of thesubproblem to which each individual has been assigned. Then,u is reproduced by applying SBX to xa and xb . Finally,
the polynomial mutation is applied to u. All ADA-basedalgorithms use the same mating selection scheme, regardlessof the mating selection schemes in their original EMOAs.
In principle, any variation operators can be incorporated into
ADA, including DE operators [26], PSO operators [16], andmodel-based methods [15], [32]. Any parent selection methods can also be used in ADA, including the distance-based
selection methods in the solution space [33]. The performanceof ADA can be improved by using these more sophisticatedmethods. However, if such effective methods are used in ADA,
it is unclear which algorithmic component mainly contributes
to the overall performance of ADA-based algorithms. Since wewant to investigate the effectiveness of the ADA framework

Algorithm 1: The ADA frameworkt ← 1, µ ← N , initialize the population P = {x1 , ..., xµ }and the weight/reference vector set W = {w1 , ..., wN };2 for i ∈ {1, ..., N } do
3
Assign xi to the i-th subproblem;
1

while The termination criteria are not met doµ ← |P |;6
Randomly select a and b from {1, ..., µ} such that a 6= b;
7Generate the child u by applying the crossover operation
to xa and xb ;8
Apply the mutation operation to u;
9Normalize the objective vectors f (x1 ), ..., f (xµ ), f (u);10
Assign u to the j-th subproblem;11
X ← {x ∈ P |x has been assigned to the j-th
subproblem and is in the neighborhood of u in thesolution space};12
bexplorer ← FALSE and bwinner ← FALSE;
13
if X = ∅ then14
bexplorer ← TRUE;
4
5

15
16
17

for x ∈ X do
if x is worse than u thenP ← P \{x} and bwinner ← TRUE

19

if bwinner = TRUE or bexplorer = TRUE thenP ← P ∪ {u};20t ← t + 1;

18

21

return P for the decision making (Subsection III-G) or
benchmarking (Subsection III-H);

in an isolated manner, we use the simplest variation operatorsand parent selection method in this paper.
B. Normalization
Since the objective functions of most real-world problems
are differently scaled, the normalization method is mandatory.
In ADA, the normalization method depends on an EMOA to
be combined. For example, NSGA-III-ADA uses the interceptbased normalization method in NSGA-III.
If the original EMOA does not have a normalization method(e.g., MOEA/D-AGR), we use the following simple normalization method. The objective vector f (x) is normalized using∗ T
the approximated ideal point z ∗ = (z1∗ , ..., zM) and the worst
worst T
worst
worst
, ..., zM) . The i-th
point in the population z= (z1element of the normalized objective vector f 0 (x) is given as:
fi0 (x) = (fi (x) − zi∗ )/(ziworst − zi∗ ).
C. Assignment operation
The child u is assigned to the j-th subproblem based on
its normalized objective vector (line 10 in Algorithm 1). The
assignment operation plays a crucial role in ADA to maintainthe diversity of the population in the objective space.
Each EMOA has a different method to select the index j of
the subproblem to which u is assigned. MOEA/D-AGR-ADA
and eMOEA/D-ADA select the j-th subproblem with the minimum scalarizing function value g(u|wj ) as in (1). MOEA/DDU-ADA, NSGA-III-ADA, and θ-DEA-ADA assign u to thej-th subproblem with the minimum perpendicular distance
between f 0 (u) and wj as in (6). RVEA selects the subproblem

6with the minimum angle between f 0 (u) and the unit reference
vector v j as in (10). ADA requires only a single index j for u
whereas MOEA/D-AGR, eMOEA/D, and MOEA/D-DU selectK subproblem indices from {1, ..., N }.

Solution space

Objective space

D. Neighborhood criterion in the solution space
ADA requires the neighborhood criterion in the solution
space (line 11 in Algorithm 1). While the objective space diversity is maintained by the assignment operation, the solution
space diversity is controlled by the neighborhood criterion.
We use a simple relative distance-based neighborhood criterion presented in [17]. First, the normalized Euclidean distance
between each individual x in P and u is calculated in thesolution space. The upper and lower bounds for each decisionvariable of a problem are used for the normalization. Then, allµ individuals in P are sorted based on their distance values
in descending order. If x is within the L nearest individuals
from u, x is said to be a neighbor of u in the solution space.
L is a control parameter in this neighborhood criterion.
In addition to the relative distance-based neighborhood
criterion, any neighborhood criterion can be incorporated into
ADA. A number of niching methods have been proposed in
the multi-modal single-objective optimization community [34].
Modern niching methods include the adaptive radius-based
method [35] and the nearest-better clustering [36]. However,
our preliminary results show that such sophisticated niching
methods do not work well in ADA. The main reason for the
failure is that an appropriate parameter specification for sucha modern niching method is difficult due to the difficulty inunderstanding the fitness landscape of an MOP [37]. For thisreason, we use the simple neighborhood criterion.
E. Deletion operation
Let X be a set of individuals that are in the neighborhoodof u in the solution space among the individuals assigned
to the same j-th subproblem as u (line 11 in Algorithm 1).
Fig. 2 shows an example of X. In Fig. 2, µ = 10, N = 4,
L = 3, and u has been assigned to the third-subproblem (w3 ).
A set of neighborhood individuals of u in the solution space
is Y = {xb , xc , xd }. A set of individuals assigned to thethird subproblem in the objective space is Z = {xa , xb , xc }.
In this case, X = Y ∩ Z = {xb , xc }. Whereas xd is theneighborhood individual of u, xd has been assigned to thesecond subproblem (w2 ). xa has been assigned to the third
subproblem with u, but xa is not in the neighborhood of u.
Since xa is dominated by u with respect to the two objectives
in Fig. 2, xa is deleted in most EMOAs. In contrast, xa cansurvive in P in ADA. This is because xa is not a neighbor ofu (and any other individual assigned to the third subproblem).
Thus, xa is not compared to u. Although the quality of xa ispoor, xa contributes to the solution space diversity of P .
A paired comparison between u and each x in X is
performed (lines 15–17 in Algorithm 1). If x is evaluated as
being worse than u (by the evaluation criterion explained in
the next paragraph), x is deleted from P . The deletion ofsuch x is reasonable since it is based on both the quality inthe objective space and the diversity in the solution space.

Individuals assigned tothe same subproblem

Neighborhood individuals
in the solution space

Fig. 2: Example of X. In this example, X = {xb , xc }. Since xa
and xd are not compared to the child u in ADA, xa and xd survive
to the next iteration independent of their quality.

The comparison criterion depends on the environmental
selection in each original EMOA. In MOEA/D-type algorithms, the comparison is based on the scalarizing function g.
MOEA/D-AGR-ADA, MOEA/D-DU-ADA, and eMOEA/DADA use g tch in (2), g dtch in (3), and g msf in (4), respectively. If g(x|wj ) ≥ g(u|wj ) on the j-th subproblem, x isremoved from P . In NSGA-III-ADA and θ-DEA-ADA, theprimary comparison between x and u is based on the Paretodominance relation. In NSGA-III-ADA, ties are broken by
the perpendicular distance between the normalized objective
vector and wj . In θ-DEA-ADA, the tie-breaker is the θdominance. RVEA-ADA uses the APD value given in (11).
F. Addition operation
The child u is added to the population P if either of thefollowing two conditions is met (lines 18–19 in Algorithm 1).
One is that no individual exists in X. An empty X meansthat there is no neighborhood individual of u in P in thesolution space (or the objective space). If the first condition
is met, u enters P without any comparison. Although u withinferior quality is likely to enter P , it helps P to maintainthe solution space diversity (or the objective space diversity).
While a dominated individual is unlikely to survive to the next
iteration in most EMOAs, it can remain in P due to the first
criterion in the addition operation.
The other is that u performs better than at least one
individual in X in the deletion operation (lines 16–17 in
Algorithm 1). Since u has good quality in its neighborhoodin the solution space, it should be added to P .
G. Decision making
We introduce an effective decision making with a solution set found by an ADA-based algorithm. Fig. 3 shows
an example of the decision making. Fig. 3 (a) exhibits the
distribution of all objective vectors in the final population P
of MOEA/D-AGR-ADA on the two-objective and two-variable
SYM-PART1 problem [38]. As shown in Fig. 3 (d), nine
equivalent Pareto optimal solution sets are on the nine lines in
SYM-PART1. The experimental setting is described in Section
IV later. The decision making in ADA is the following twophase method based on Aprimary and Asecondary .

7

4

4

3

3

2
1
0
0

1

2

3

(a) P in the obj. space

4

4
3

1.70

2

2

1.66

1

1

The objective vector of
the 77th subproblem

001
primary

(b) A

2

3

4

1.68
1.64
1.62
0.51 0.53 0.55

001secondaryin the obj. space (c) A

2

3

15
10
5
0
−5
−10
−15
−15−10 −5 0
4secondaryin the obj. space (d) A

5 10 15in the solution space

Fig. 3: Example of the decision making based on multiple equivalent solutions found by the proposed approach. Solutions obtained by

MOEA/D-AGR-ADA on SYM-PART1 are shown. These figures show results of a single run with a median IGD+ value among 31 runs. In
Fig. 3 (a), (b), and (c), the x and y axes represent f1 and f2 , respectively. In Fig. 3 (d), the x and y axes represent x1 and x2 , respectively.
The nine gray lines in Fig. 3 (d) are the equivalent Pareto optimal solution sets.

1) Aprimary -based decision making: Fig. 3 (a) shows that
P contains a large number of solutions. This is undesirable forthe decision maker, because she/he usually wants to examine
a small number of well-distributed solutions in the objective
space [4]. To address this issue, the best solution for the j-th
subproblem is selected from P based on the same criterion
in the deletion operation (j ∈ {1, ..., N }). Then, only nondominated solutions are selected from the N best solutions.
Let Aprimary be a set of the non-dominated solutions obtainedby this procedure. Fig. 3 (b) shows non-dominated solutions in
Aprimary in the objective space. The decision maker examines
objective vectors in Aprimary , rather than P .
When the decision maker wants to examine N or less nondominated solutions (e.g., only 10 solutions), other solution selection methods can be used to select Aprimary from P . If thenumber of objectives M is less than 4, efficient hypervolumebased selection approaches (e.g., [39]) are available. If M ≥ 4,
computationally cheap distance-based selection methods (e.g.,[40]) can be used.
2) Asecondary -based decision making: After the decision
maker has examined Aprimary in the objective space, she/heselects the final solution xfinal from Aprimary . Suppose thatthe decision maker has selected a solution on the 77-th
subproblem as xfinal in Fig. 3 (b). In ADA, she/he can examine
other dissimilar solutions with similar quality to xfinal . Let
Asecondary be a set of all solutions assigned to the same
77-th subproblem with xfinal . Figs. 3 (c) and (d) show the
distribution of all nine solutions in Asecondary in the objective
and solution spaces, respectively. Although the nine solutionsin Asecondary have almost the same quality in the objective
space, they are dissimilar in the solution space. In addition
to xfinal , the decision maker can examine other candidates in
Asecondary based on her/his preference in the solution space.

EMMAs and other optimizers. This is because the current
population P in ADA can contain an unbounded number ofsolutions. Most performance indicators cannot assess solution
sets with different sizes in a fair manner [41]. IGD+ [42] and
IGDX [32] described in Section IV later are not exceptions.
Thus, ADA requires a method of selecting a constant numberof solutions from P for benchmarking studies.
We introduce two post-processing methods that are used toevaluate the performance of ADA-based algorithms for MOPs
and MMOPs, respectively. In general, uniformly distributed
solutions are unlikely to be uniformly distributed objective
vectors due to a non-uniform mapping from the solution space
to the objective space. To address this issue, we use thetwo post-processing methods for the objective and solution
spaces, respectively. On the one hand, the same method of
selecting Aprimary from P described in Subsection III-G isused for performance indicators of MOPs (e.g., IGD [43]).
Since indicators of MOPs assess the distribution of a solution
set in the objective space, the choice of Aprimary is reasonable.
On the other hand, the solution distance-based selection
method presented in [17] is used for performance indicators for MMOPs (e.g., IGDX [32]). Let us consider thetask of selecting N sparsely distributed solutions from all
non-dominated solutions in P . Below, D(x, A) denotes the
distance between a solution x and its nearest solution in asolution set A in the normalized solution space. First, Atertiary
is set to empty. A solution is randomly selected from thenon-dominated solution set and stored into Atertiary . Then, asolution with the maximum D(x, Atertiary ) value is repeatedly
added to Atertiary until |Atertiary | = N . Unlike Aprimary andAsecondary , Atertiary is used only for benchmarking of ADAbased algorithms.
I. Applicability of ADA

H. Two Post-processing methods for benchmarking
Apart from the practical two-phase decision making method
described in Subsection III-G, we here consider benchmarkingof ADA-based EMMAs. The performance of EMOAs and
EMMAs is evaluated using performance indicators. However,a fair performance comparison is difficult between ADA-based

ADA is a framework to improve the performance of
decomposition-based EMOAs for MMOPs. We do not claimthat ADA can be combined into any decomposition-based
EMOAs. Since ADA requires a method of assigning a childu to a subproblem, ADA is inapplicable to EMOAs withno assignment mechanism. Such EMOAs include MOEA/D

8[4] and MOEA/D-DRA [44]. Also, the deletion operation
performs the pairwise comparison independently from otherindividuals. Thus, ADA is not applicable to EMOAs whoseenvironmental selection is performed for all individuals in P ,such as MOEA/D-STM [45] and VaEA [30].
In summary, ADA can be combined into EMOAs with anassignment mechanism and a pairwise comparison-based environmental selection. The six EMOAs explained in Subsections
II-C and II-D satisfy these conditions. In addition, ADA isapplicable to I-DBEA [46] and SPEA/R [31].
Fortunately, even if a decomposition-based EMOA does
not satisfy the above-mentioned two conditions, ADA can beapplied to the EMOA after some modifications. For example,
since MOEA/D does not have the assignment method, we cannot directly combine ADA in MOEA/D. However, MOEA/D
can be easily modified by using any of the three assignment
methods described in Subsection III-C. In fact, MOEA/DAD proposed in [17] is an ADA-based MOEA/D with theassignment operation in (6). Thus, the applicability of ADA
is not limited to only a few decomposition-based EMOAs.
J. Originality of ADA
In addition to MOEA/D-AD [17], a variant of MOEA/D for
MMOPs is proposed in [47]. The MOEA/D variant assigns
K individuals to each subproblem. The fitness value is basedon the PBI function value and two distance values in thesolution space. The main disadvantage of the MOEA/D variantin [47] is the difficulty in finding a proper K value. Since
the number of equivalent Pareto optimal solution subsets
is unknown a priori, fine-tuning of K is necessary for agiven problem. Although a multi-start decomposition-based
approach is proposed in [38], it has a similar disadvantage.
In contrast, ADA adaptively adjusts the number of individuals
assigned to each subproblem. Thus, ADA does not require a
problem-dependent parameter such as K.
TriMOEA-TA&R [24] consists of various advanced components, including the convergence and diversity archivesbased strategy as in Two Arch2 [48], the decision variabledecomposition method in MOEA/DVA [49], and the anglebased individual assignment in RVEA [22]. ADA and
TriMOEA-TA&R are similar in that they handle diversity inthe objective and solution spaces in a two-phase manner. While
TriMOEA-TA&R uses an absolute distance-based neighborhood criterion with σniche , ADA uses the relative distancebased neighborhood criterion with L. ADA uses a much
simpler mating scheme, and it does not use any decomposition
method of decision variables. ADA also uses the adaptive
population sizing strategy to handle equivalent solutions, asmentioned above. Whereas TriMOEA-TA&R is an algorithm
for MMOPs with distance-related variables, ADA is a generalframework for various MMOPs.
IV. E XPERIMENTAL SETTINGS
A. Test problems
We use the following eight multi-modal multi-objective test
problems: the Two-On-One problem [50], the three SYMPART problems [38], the two SSUF problems [51], the Polygon problem [52], and the Omni-test problem [13]. Table I

TABLE I: Properties of multi-modal multi-objective test problems,
where M , D, and O denote the number of objectives, design
variables, and equivalent Pareto optimal solution subsets, respectively.
Test problems

M

D

O

Two-On-One [50] and SSUF1,3 [51]

2

2

2

SYM-PART1–3 [38]

2

2

9

Polygon [52]

Any

2

Any

Omni-test [13]

2

Any

3D

shows their properties, including the number of objectives
M , the number of decision variables D, and the number of
equivalent Pareto optimal solution sets O.
Two-On-One has two equivalent Pareto optimal solution setsthat are symmetrical with respect to the origin. Equivalent
Pareto optimal solution sets are on the nine lines in SYMPART1, as shown in Fig. 3. The nine lines are rotatedin SYM-PART2. In addition, the nine lines are distorted
in SYM-PART3. SSUF1 and SSUF3 have two symmetrical
Pareto optimal solution sets. We evaluate the scalability of
EMMAs to M using Polygon. We set M of Polygon as
follows: M ∈ {3, 5, 8, 10}. Although O can be any number in
Polygon, it was set to be nine. We investigate the scalability of
EMMAs to D and O using Omni-test. We set D as follows:D ∈ {2, 3, 5, 8, 10}. O increases exponentially with increased
D in Omni-test. In summary, we use 15 test problem instances.
HPS [23] and MMMOP [24] have been recently proposed.
However, HPS and MMMOP have so-called “distance-related”
variables that affect only the distance between the objective
vector and the Pareto front. For this reason, we do not mainlyuse HPS and MMMOP for our benchmarking study. We use
HPS and MMMOP only in Subsections V-E and V-F.
B. Performance indicators
Below, A is a set of solutions obtained by an EMMA.
A∗ is also a set of reference solutions in the Pareto optimal
solution set. The size of A∗ was set to 5 000. For A∗ ofeach problem, 5 000 solutions were selected from randomly
generated 10 000 Pareto-optimal solutions using the distancebased solution selection method [17] (Subsection III-H).
We use IGD+ [42] to evaluate A in terms of both convergence to the Pareto front and diversity in the objective space:


n
o
X

1

min d f (x), f (z)  ,IGD+ (A) =
|A∗ |
∗ x∈A
z ∈A

(14)

q
2
PM
+
where d(a, b) =is a
i=1 max{ai − bi , 0} . IGD
modified version of IGD [43]. While the original IGD is Pareto
non-compliant, IGD+ is weakly Pareto compliant.
We evaluate how well A approximates the Pareto-optimal
solution set in the solution space using IGDX [32]:


n
o
1  X
IGDX(A) =min ED x, z  ,
|A∗ |
∗ x∈A
z ∈A

(15)

where ED(a, b) is the Euclidean distance between a and b.

9

EMOAs that can find A with small IGD+ and IGDX values
are efficient multi-objective optimizers and multi-modal multiobjective optimizers, respectively. In the ADA-based algorithms, Aprimary and Atertiary explained in Subsection III-H
are used for the IGD+ and IGDX calculations, respectively.
C. Average performance score
We use the average performance score (APS) [53] in orderto aggregate results on various problems. Suppose that nalgorithms A1 , ..., An are compared for a problem instancebased on the indicator values obtained in multiple runs. For
each i ∈ {1, ..., n} and j ∈ {1, ..., n} \{i}, if Aj significantly outperforms Ai using the Wilcoxon rank-sum test withp < 0.001, then δi,j = 1; otherwise, δP
i,j = 0. The scoren
P (Ai ) is defined as follows: P (Ai ) =j∈{1,...,n}\{i} δi,j .
The score P (Ai ) represents the number of algorithms that
outperform Ai . The APS value of Ai is the average of the
P (Ai ) values for all problem instances. A small APS valueof Ai indicates that Ai performs well among n algorithms.
D. EMMAs and EMOAs
We examine the performance of the six ADA-based EMMAs: MOEA/D-AGR-ADA, MOEA/D-DU-ADA, eMOEA/DADA, NSGA-III-ADA, θ-DEA-ADA, and RVEA-ADA. We
implemented the ADA-based algorithms using jMetal [54].
Their source codes can be downloaded from the supplementary
website (https://sites.google.com/view/mmoada). We compare
the ADA-based EMMAs to their original EMOAs: MOEA/DAGR [19], MOEA/D-DU [20], eMOEA/D [21], NSGA-III [5],θ-DEA [18], and RVEA [22]. Our implementations of the six
original EMOAs were based on their corresponding articles,
except for MOEA/D-AGR. We replaced the DE operator with
SBX in MOEA/D-AGR to remove the effect of variation
operators. For details, refer to the corresponding articles.
The number of maximum function evaluations was 30 000.
31 runs were performed for each test problem. A set ofweight/reference vectors were generated using the simplexlattice design method for M < 6 and its two-layered version[5] for M ≥ 6. The number of the weight/reference vectors
N was 100, 105, 210, 156, and 230, for M = 2, 3, 5, 8,
and 10, respectively. The SBX crossover and the polynomial
mutation were used in all methods, including MOEA/D-AGR.
Their control parameters were set as follows: pc = 1, ηc = 20,
pm = 1/D, and ηm = 20. According to the analysis presented
in [17], L of the neighborhood criterion in ADA was set to
L = b0.1µc. For example, L = 201 when µ = 2 019. Other
parameters were set according to the corresponding references.
V. E XPERIMENTAL RESULTS
This section shows performance analysis of the six ADAbased EMMAs. Subsection V-A describes the effect of ADA
on the six EMOAs. Subsection V-B compares the six ADAbased EMMAs to state-of-the-art EMMAs. Subsection V-C
discusses the adaptive population sizing in ADA. Subsection
V-D investigates the influence of the L value on the performance of the six ADA-based EMMAs. Subsection V-E

TABLE II: Results of the six EMOAs and their ADA versions on
the 15 test problem instances. Tables (a) and (b) show the APS valuesof the algorithms for IGD+ and IGDX, respectively. AGR and DU
stand for MOEA/D-AGR and MOEA/D-DU, respectively.
(a) IGD+
AGR-ADADU-ADAeMOEA/D-ADA
NSGA-III-ADA
θ-DEA-ADA
RVEA-ADA

Orig.
0.0 (1)0.0 (1)0.2 (1)0.0 (1)0.1 (1)0.5 (2)(b) IGDX
ADA
1.0 (2)0.8 (2)0.7 (2)0.8 (2)0.6 (2)0.3 (1)

AGR-ADADU-ADAeMOEA/D-ADA
NSGA-III-ADA
θ-DEA-ADA
RVEA-ADA

Orig.
1.0 (2)0.9 (2)1.0 (2)1.0 (2)1.0 (2)1.0 (2)ADA
0.0 (1)0.0 (1)0.0 (1)0.0 (1)0.0 (1)0.0 (1)

examines the performance of the six ADA-based EMMAs on
test problems with distance-related variables. Subsection V-Fpresents a comparison with the state-of-the-art EMMAs using
an unbounded external archive. Subsection V-G discusses theruntime of the ADA-based algorithms.
A. Effect of ADA
Tables II (a) and (b) show the paired comparison of each
EMOA and its ADA version on the 15 test problem instancesin terms of IGD+ and IGDX, respectively. Table II shows
only the APS value of each algorithm at the final iteration.
Table S.1 in the supplementary file shows detailed results. As
described in Subsection III-H, only N solutions in Aprimaryand Atertiary were used for the IGD+ and IGDX calculations.
Thus, all algorithms are compared under the same number ofsolutions. We set the p value to 0.001 for the Wilcoxon ranksum test. M -Polygon is the M -objective Polygon problem, andD-Omni-test is the D-variable Omni-test problem. Below, we
describe the results of IGD+ and IGDX.
1) IGD+ : Table II (a) shows that the original EMOAs
outperform their ADA-based EMMAs regarding IGD+ (except
for RVEA-ADA). In general, EMMAs perform slightly worse
than EMOAs for multi-objective optimization [15], [16]. While
EMOAs aim to find a good approximation of the Pareto front,
EMMAs attempt to locate all Pareto optimal solutions. For
this reason, the performance of the ADA versions regarding
IGD+ is worse than that of their original EMOAs.
However, as shown in Table S.1 in the supplementary file,the IGD+ value of the ADA versions is only 3.57 times worse
than that of their original EMOAs even in the worst case.
Also, the ADA versions perform better than their original
EMOAs on some problems. For example, eMOEA/D-ADA
obtains better IGD+ values than eMOEA/D on SYM-PART2,
SSUF3, and 10-Polygon. As discussed in [55], a mechanism
for maintaining the solution space diversity can help the ADA
versions to find high-quality solutions.
2) IGDX: Table II (b) shows that the ADA-based EMMAs
significantly outperform their original EMOAs. As shown in
Table S.1, the IGDX value of the ADA versions is 141.58
times better than that of their original EMOAs in the best case.
These results indicate that ADA improves the performance ofthe six decomposition-based EMOAs for MMOPs.
Table S.1 shows that the ADA-based EMMAs work wellon Polygon with M ∈ {3, 5, 8, 10} and Omni-test withD ∈ {2, 3, 5, 8, 10}. Since the number of equivalent Pareto

10

TABLE III: Results of the nine EMMAs on the 15 test problem

0

0.5

0.5

11

1

1

1

0.5

0.5

0.5

00 0

00 0

00

0.5

(a) eMOEA/D

0

0.5

0.5

11

0.5

11

0.5

(b) NSGA-III

0.5

11

(c) RVEA

1

1

1

0.5

0.5

0.5

00 0

(d) eMOEA/D-ADA

0.5

0.5

11

00 0

0.5

(e) NSGA-III-ADA

0.5

11

Fig. 4: Distribution of non-dominated solutions found by eachmethod in the objective space on the three-objective Polygon problem. The x, y, and z axes represent f1 , f2 , and f3 , respectively.
10

10

8

8

8

6

6

6

4

4

4

2

2

2

0
0

0
10 0

2

4

6

8

(a) eMOEA/D

2

4

6

8

10

0
0

(b) NSGA-III
10

10

8

8

8

6

6

6

4

4

4

2

2

2

0
0

0
10 0

4

6

8

(d) eMOEA/D-ADA

2

4

6

8

4

6

8

10

8

10

(c) RVEA

10

2

2

10

(e) NSGA-III-ADA

0
0

2

4

6

(a) IGD+
MOEA/D-AGR-ADA
MOEA/D-DU-ADA
eMOEA/D-ADA
NSGA-III-ADA
θ-DEA-ADA
RVEA-ADATriMOEA-TA&R
MO Ring PSO SCD
Omni-optimizer(b) IGDX
3.7 (7)2.8 (5)2.9 (6)
1.7 (2)2.7 (4)5.9 (9)2.4 (3)4.2 (8)1.6 (1)

MOEA/D-AGR-ADA
MOEA/D-DU-ADA
eMOEA/D-ADA
NSGA-III-ADA
θ-DEA-ADA
RVEA-ADATriMOEA-TA&R
MO Ring PSO SCD
Omni-optimizer

2.1 (3)2.1 (4)2.3 (5)0.8 (1)1.0 (2)4.1 (6)7.7 (9)4.3 (7)5.1 (8)00(f) RVEA-ADA

10

instances. Tables (a) and (b) show the APS values of the algorithms
for the IGD+ and IGDX indicators, respectively. The numbers inparentheses are the ranks of the algorithms based on their APS values.

verted triangular Pareto front, it is difficult for decompositionbased EMOAs to find good solutions.
For the above-mentioned reasons, no method in Fig. 4
can approximate the Pareto front of Polygon well. Nevertheless, some ADA-based algorithms (e.g., eMOEA/D-ADA
and RVEA-ADA) find better distributed objective vectors thantheir original EMOAs. This unintended effect of ADA may becaused by the diversity of the population in the solution space[55]. Although we do not claim that ADA can improve theperformance of EMOAs for MOPs in addition to MMOPs, thesolution space diversity maintenance might help the original
EMOAs to handle problems with irregular Pareto fronts.
In the Polygon problem with M = 3 and O = 9, equivalent
Pareto solution sets are inside of the nine regular triangles
in the solution space. Figs. 5 (a)–(c) show that the original
EMOAs cannot locate all equivalent Pareto solution sets well.
In contrast, Figs. 5 (d)–(f) demonstrate that their ADA versions
can locate all nine equivalent Pareto solution sets. Results onother test problems are similar to Fig. 5.

(f) RVEA-ADA

Fig. 5: Distribution of non-dominated solutions found by eachmethod in the solution space on the three-objective Polygon problem.
The x and y axes represent x1 and x2 , respectively.

optimal solution sets O increases exponentially with increased
D in Omni-test, the results show that ADA can also handle
problems with a large O value. In summary, ADA has asufficient scalability to M , D, and O.
3) Distribution of solutions: Figs. 4 and 5 show the distribution of non-dominated solutions found by eMOEA/D-ADA,
NSGA-III-ADA, RVEA-ADA, and their original versions on
Polygon with M = 3 in the objective and solution spaces,
respectively. These figures show results of a single run with a
median IGD+ and IGDX values, respectively. For the ADAbased algorithms, non-dominated solutions in Aprimary andAtertiary are shown in Figs. 4 and 5, respectively. Figs. S.1
and S.2 in the supplementary file show the results of the other
ADA-based algorithms, but the results are similar to Figs. 4
and 5. It was shown in [5] that some decomposition-based
EMOAs did not work well on problems with convex Pareto
fronts. It was also presented in [56] that some decompositionbased EMOAs performed poorly on problems with invertedtriangular Pareto fronts. Since Polygon has a convex and an in-

B. Comparison with other EMMAs
In Subsection V-A, we demonstrated that ADA can improve
the performance of the six EMOAs for MMOPs. Here, we
compare the six ADA-based algorithms to the following
three EMMAs: Omni-optimizer [13], MO Ring PSO SCD[16], and TriMOEA-TA&R [24]. Omni-optimizer is the mostrepresentative EMMA. MO Ring PSO SCD and TriMOEATA&R are recently proposed methods. Default parametersettings were used for the three methods.
Table III shows the APS values of the nine EMMAs on
the 15 test problem instances. Table S.2 in the supplementary
file shows detailed results. TriMOEA-TA&R incorrectly recognizes that some problems have distance-related variables (e.g.,
SYM-PART1). In such a case, TriMOEA-TA&R generates a
set of additional solutions Aadd by recombining solutions in
the diversity archive and the distance-related variables at theend of the search. Since |Aadd | > N , methods of selecting
N solutions from Aadd are needed. We use the two postprocessing methods in ADA. In the same manner as in ADA,
Aprimary and Atertiary are used for the IGD+ and IGDX
calculations, respectively. We use the selection method in
NSGA-III-ADA to obtain Aprimary from Aadd .
Table III (a) shows that Omni-optimizer performs the best
regarding IGD+ . The good performance of Omni-optimizer

Population size µ

11

1800
1600
1400
1200
1000
800
600
400
200
0
0

2-Omni-test
3-Omni-test
5-Omni-test
8-Omni-test
10-Omni-test

10000
20000
30000
Number of evaluations

Fig. 6: Evolution of µ of NSGA-III-ADA on Omni-test with D ∈{2, 3, 5, 8, 10}. The mean µ values over 31 runs are reported.
Cumulative number

104
103
102
101
100
0

Addition (explorer)
Addition (winner)
Deletion

10000
20000
30000
Number of evaluations

Fig. 7: Cumulative number of the activations of the addition operation based on the two criteria (bexplorer and bwinner ) and the deletion
operation in NSGA-III-ADA on 3-Polygon (M = 3). Results of asingle run with a median IGDX value among 31 runs are shown.

for MOPs is consistent with the results presented in [16], [17].
NSGA-III-ADA shows the best performance regarding IGD+
among the six ADA-based algorithms. Table III (b) shows thatthe six ADA-based algorithms perform better than the threeother algorithms regarding IGDX. TriMOEA-TA&R performs
the worst regarding IGDX. This is because the 15 problem
instances have no distance-related variables. In summary, theresults indicate that the six ADA-based algorithms have better
performance than the state-of-the-art algorithms for MMOPs.
C. Adaptive population sizing in ADA
In ADA, the number of individuals assigned to each subproblem is adaptively adjusted by the deletion and addition
operations. Thus, the population size µ is not constant. Here,
we discuss the adaptive population sizing in ADA.
Fig. 6 shows the change of µ of NSGA-III-ADA on Omnitest with D ∈ {2, 3, 5, 8, 10}. Results of the other ADA
versions are similar to Fig. 6. The µ value for D = 3 isalways larger than that for D = 2. In contrast, the µ value
decreases as the D value increases from D = 3. This may be
because problems with a large D value are generally difficult
for EMMAs. Since the search does not proceed well forD ∈ {5, 8, 10}, the µ value does not significantly increase.
The difficulty of finding multiple equivalent Pareto sets has
been reported for the case of a large D value [57].
In addition, the trajectory of µ is problem-dependent. Fig.
S.3 in the supplementary file shows the change of µ on otherproblems. While the µ value is approximately 640 on Two-OnOne and 3-Polygon at the end of the search, it is approximately
960–1 090 on other problems. Fig. S.3 shows that the µ value

sharply increases in early iterations. After some evaluations,the µ value is stable.
Fig. 7 shows the cumulative numbers of the activationsof the addition operation based on the two criteria (bexplorer
and bwinner in Algorithm 1) and the deletion operation in
NSGA-III-ADA on 3-Polygon. On the one hand, the numberof activations of the addition operation based on bwinner isalmost the same as that of the deletion operation throughout
the evolution in Fig. 7. This competitive behavior of the two
operations leads to the stable µ value in Fig. S.3. On the other
hand, the addition operation based on bexplorer is frequently
performed only in an early stage of evolution in Fig. 7. Thus,
the sharp increase of the µ value in Fig. S.3 is due to the
bexplorer -based addition operation. Since only one individualis assigned to each subproblem at the beginning of the search,
the bexplorer -based addition operation is often activated.
In summary, the µ value is adaptively adjusted by thecooperative behavior of the deletion and addition operationsin ADA. In the addition operation, the roles of the two
criteria bexplorer and bwinner differ from each other. Thus, it isimportant for ADA to use both bexplorer and bwinner .
D. Impact of L on the six ADA-based algorithms
We investigate the influence of the L value on the performance of the ADA-based algorithms. Table IV shows the
APS values of the six ADA-based algorithms with six L
values on the 15 test problem instances. Tables S.3–S.8 inthe supplementary file show detailed results.
Table IV (a) shows that the six ADA-based algorithms with
L = b0.05µc have the worst performance regarding IGD+ .
Too small L values degrade the performance of the ADAbased algorithms regarding IGD+ . Table IV (b) shows thatL = b0.1µc is most suitable for all ADA-based algorithms(except for RVEA-ADA) in terms of IGDX. Too large L
values degrade the performance of the ADA-based algorithms
regarding IGDX. In summary, L = b0.1µc are suitable formost ADA-based algorithms for MMOPs.
We also investigate the influence of L on the evolution of µ.
Figs. S.4 and S.5 in the supplementary file show the change ofµ of NSGA-III-ADA with various L values on Omni-test withvarious D values and other problems, respectively. Results ofother ADA versions are similar to the results of NSGA-IIIADA. Figs. S.4 and S.5 indicate that the µ value decreases
as the L value increases. The larger the L value is, the moreindividuals can be neighbors of the child u. As a result, theniching mechanism in ADA deteriorates. These observations
are consistent with the above-mentioned results of IGDX.
E. Results on test problems with distance-related variables
This subsection shows results on the HPS and MMMOP test
problems. Although the six HPS problem instances (HPS1,..., HPS6) are proposed in [23], we use only HPS2. This isbecause all the other five problem instances are variants of
HPS2 and also because the details of only HPS2 are provided
in [23]. M , D, and O in HPS2 are as follows: M = 2, D = 7,
and O = 4. The 20 MMMOP problem instances (MMMOP1A,..., MMMOP6D) are proposed in [24]. We also use those 20

12

TABLE IV: Results of the six ADA-based algorithms with L ∈{b0.05µc, b0.1µc, b0.2µc, b0.3µc, b0.4µc, b0.5µc} on the 15 test
problem instances. Tables (a) and (b) show the APS values of thealgorithms for the IGD+ and IGDX indicators, respectively. AGR
and DU stand for MOEA/D-AGR and MOEA/D-DU, respectively.
(a) IGD+
AGR-ADADU-ADAeMOEA/D-ADA
NSGA-III-ADA
θ-DEA-ADA
RVEA-ADA

b0.05µc
1.6 (6)
1.6 (6)
1.7 (6)2.6 (6)
2.2 (6)
1.9 (6)b0.1µc0.9 (3)0.7 (2)0.9 (3)
1.5 (5)
1.3 (5)0.5 (1)

AGR-ADADU-ADAeMOEA/D-ADA
NSGA-III-ADA
θ-DEA-ADA
RVEA-ADA

b0.05µc
0.8 (2)0.8 (2)0.9 (2)1.0 (3)
1.1 (3)0.5 (1)b0.1µc0.1 (1)0.3 (1)0.3 (1)0.1 (1)0.2 (1)0.6 (2)b0.2µc0.6 (1)0.5 (1)0.5 (1)0.7 (4)0.4 (4)
1.1 (3)b0.3µc0.7 (2)0.7 (3)0.5 (2)0.1 (1)0.2 (3)0.9 (2)b0.4µc1.0 (4)1.0 (4)0.9 (4)0.1 (1)0.0 (1)
1.3 (4)b0.5µc1.3 (5)
1.3 (5)
1.1 (5)0.2 (3)0.0 (1)
1.9 (5)b0.3µc2.6 (4)2.5 (4)2.2 (4)2.2 (4)2.1 (4)2.7 (4)b0.4µc2.9 (5)2.9 (5)2.5 (5)2.7 (5)2.7 (5)2.9 (5)b0.5µc3.2 (6)
3.5 (6)
3.2 (6)
3.5 (6)
3.2 (6)
3.1 (6)(b) IGDX
b0.2µc0.8 (3)1.0 (3)
1.1 (3)0.5 (2)0.9 (2)
1.6 (3)

TABLE V: Results of the nine EMMAs on the 21 test problem instances with distance-related variables: APS value of each algorithm
for IGD+ in (a) and for IGDX in (b).
(a) IGD+
MOEA/D-AGR-ADA
MOEA/D-DU-ADA
eMOEA/D-ADA
NSGA-III-ADA
θ-DEA-ADA
RVEA-ADATriMOEA-TA&R
MO Ring PSO SCD
Omni-optimizer(b) IGDX
3.0 (5)2.2 (3)
1.5 (2)4.2 (7)3.0 (6)
4.4 (8)1.0 (1)
6.8 (9)2.3 (4)

MOEA/D-AGR-ADA
MOEA/D-DU-ADA
eMOEA/D-ADA
NSGA-III-ADA
θ-DEA-ADA
RVEA-ADATriMOEA-TA&R
MO Ring PSO SCD
Omni-optimizer

2.0 (3)
1.5 (2)0.9 (1)3.0 (6)
2.2 (4)3.7 (7)2.7 (5)6.9 (9)4.1 (8)

F. Comparison using an unbounded external archive
While
the
three
EMMAs
(TriMOEA-TA&R,
MO Ring PSO SCD, and Omni-optimizer) keep the µvalue constant, the ADA-based algorithms can adaptively
adjust the µ value. One may think that the comparisons in
Subsections V-B and V-E are unfair for this reason. Here,we compare the ADA-based algorithms to the three EMMAs
using an unbounded external archive (UEA) [58], which
stores all non-dominated solutions found during the search
process. The UEA can be incorporated into any algorithms
with no changes in their algorithmic behavior. The three
EMMAs with the UEA can maintain all non-dominatedsolutions found so far. Note that only the three EMMAs use
the UEA in this section. Since an algorithm with the UEA
significantly outperforms its original version (see [58]), sucha comparison is unfair for the ADA-based algorithms.
Tables S.11 and S.12 in the supplementary file show
the comparisons with the three EMMAs using the UEA
on the 15 test problem instances with no distance-related
variables and the 21 test problem instances with distancerelated variables, respectively. The performance of the three
EMMAs is improved by using the UEA. However, Table S.11
indicates that all ADA-based algorithms (except for RVEAADA) outperform the three EMMAs with the UEA on the 15
test problems in terms of IGDX. Table S.12 also shows thateMOEA/D-ADA performs the best on the 21 test problems interms of IGDX. Thus, the results show that some ADA-basedalgorithms perform significantly better than the three EMMAs
even with the UEA. Note that the performance of the ADAbased algorithms can be further improved by using the UEA
as in the three EMMAs demonstrated here.
G. On the runtime of the ADA-based algorithms

instances. M , D, and O in MMMOP are as follows: M ∈{2, 3}, D ∈ {2, ..., 7}, and O ∈ {2, ..., 9}.
Table S.9 in the supplementary file shows results of the six
EMOAs and their ADA versions on the 21 problem instances.
The results in Table S.9 are similar to the results in Table
II. Since the solution space diversity maintenance mechanism
in ADA is not useful to handle distance-related variables, theADA-based algorithms are likely to perform poorly on the
problems with distance-related variables. However, the results
show that the ADA-based algorithms perform better than their
original EMOAs (except for NSGA-III-ADA).
Table V shows the APS values of the nine EMMAs on
the 21 problem instances. Table S.10 in the supplementary
file shows detailed results. Since TriMOEA-TA&R can explicitly exploit distance-related variables for the search similar
to MOEA/DVA [49], its performance is improved on the
21 problem instances. Table V (a) shows that TriMOEATA&R performs the best regarding IGD+ . Table V (b)shows that TriMOEA-TA&R is the fifth-ranked algorithm
regarding IGDX. MOEA/D-AGR-ADA, MOEA/D-DU-ADA,
eMOEA/D-ADA, and θ-DEA-ADA perform better than thethree competitors regarding IGDX. In summary, the ADAbased algorithms have high performance for MMOPs even on
the problems with distance-related variables.

The space complexity of ADA itself is O(µD) if D > M .
Otherwise, it is O(µM ). Although the µ value equals to the
maximum number of function evaluations nmax in ADA in
the worst case, the µ value is naturally bounded as shown in
Subsection V-C. Thus, the empirical space complexity is muchsmaller than O(nmax D) and O(nmax M ).
Below, we explain the worst-case time complexity of the
ADA-based algorithms. ADA itself in one iteration requires
O(µD), which is due to the procedure of selecting theneighbors of the child u in the solution space (line 11 in Algorithm 1). The time complexity of an ADA-based algorithm
depends on its original version. Since the time complexity
of NSGA-III is the larger value between O(µ2 logM −2 µ) andO(µ2 M ) [5], that of NSGA-III-ADA is the largest value
among O(µ2 logM −2 µ), O(µ2 M ), and O(µD). Similarly, thetime complexity of θ-DEA-ADA and RVEA-ADA is thelarger value between O(µ2 M ) [18], [22] and O(µD). The
assignment operations in MOEA/D-AGR, MOEA/D-DU, andeMOEA/D in one function evaluation require O(µM ) [19]–[21]. Thus, the time complexity of MOEA/D-AGR-ADA,
MOEA/D-DU-ADA, and eMOEA/D-ADA in one iteration is
the larger value between O(µM ) and O(µD).
Fig. 8 shows the average CPU time of each EMMA overthe five Omni-test problems with D = 2, 3, 5, 8, 10. We

13

Average CPU time

103
10

2

101
100
10−1

2 3

5

8

10

Number of decision variables D

θ-DEA-ADA
RVEA-ADA
MO Ring PSO SCD
NSGA-III-ADA
MOEA/D-AGR-ADA
MOEA/D-DU-ADA
eMOEA/D-ADA
TriMOEA-TA&R
θ-DEA
Omni-optimizer

Fig. 8: Average CPU time of the EMMAs over the five Omni-test
problems with D = 2, 3, 5, 8, 10.

obtained all results using a workstation with Xeon E5-2620
v4/2.10GHz and 8GB RAM. Fig. 8 also shows results of
θ-DEA. The six ADA-based algorithms and θ-DEA are implemented in Java, Omni-optimizer is implemented in C, andMO Ring PSO SCD and TriMOEA-TA&R are implemented
in Matlab. These implementations are based on the available
codes of these algorithms. The use of C and Matlab is dueto the original implementations of these algorithms [13], [16],[24]. Since different programming languages require different
CPU time, Fig. 8 provides rough comparison.
As shown in Fig. 8, Omni-optimizer is fastest, followed by
θ-DEA and TriMOEA-TA&R. The six ADA-based algorithms
are much slower than these three algorithms. This is because
the µ value of the ADA-based algorithms is always larger
than that of the other algorithms (see Subsection V-C). Sincesome operations require O(µ2 M ) and O(µD) as discussed
above, a larger µ value makes an algorithm slower. Although
MO Ring PSO SCD and RVEA-ADA perform similarly in
Fig. 8, we believe that MO Ring PSO SCD with an optimized source code can perform much faster than the ADAbased algorithms. Clearly, θ-DEA-ADA is slowest. This isbecause θ-DEA-ADA needs to normalize objective vectors forevery iteration in a computationally expensive manner. Since
the original θ-DEA is a generational EMOA, it can reduce thenumber of normalizations. In fact, θ-DEA is much faster than
θ-DEA-ADA. Note that θ-DEA-ADA can be speeded up by
using a simple normalization method (e.g., Subsection III-B).
As demonstrated here, the slow speed of the ADAalgorithms is their disadvantage. Speeding up the ADAalgorithms is an avenue for future work. However, this does
not mean that ADA is impractical. Some real-world problems
require a long computation time to evaluate a solution, e.g.,by the execution of an expensive computer simulation [59].
In this case, the time of function evaluations dominates thatof the other parts of optimization algorithms. For this reason,in the EMO community, the comparison is generally based on
the number of function evaluations, rather than the CPU time.
VI. C ONCLUSION
We proposed the ADA framework to improve the performance of decomposition-based EMOAs for MMOPs. ADA
adaptively assigns one or more individuals to each subproblem to locate multiple equivalent Pareto optimal solutions.
The population size µ is automatically adjusted during thesearch process. As presented in Subsection III-G, the effective

decision making can be performed using multiple equivalentsolutions found by ADA. We incorporated ADA into the six
decomposition-based EMOAs. We examined the performanceof those six ADA-based algorithms on the 15 test problem
instances and the 21 test problem instances with distancerelated variables. The results show that ADA can improve theperformance of the original EMOAs for MMOPs. We also
analyzed the adaptive behavior of ADA.
Computational overhead in the ADA-based algorithms istheir disadvantage (see Subsection V-G). Thus, the unboundpopulation size of ADA is a “double-edged sword”. However,as discussed in Subsection V-G, the computational overheadin ADA is practically acceptable in real-world applications
where solution evaluations need dominant computation timein the whole optimization process.
For a fair comparison, we compared the ADA-based algorithms to the three EMMAs using the UEA in Subsection
V-F. However, the UEA is inherently a post-processing method
since an algorithm cannot take advantage of non-dominatedsolutions in the UEA during the search process as in the ADAbased algorithms. For this reason, the comparison with thethree EMMAs using the UEA is not totally fair. A generaladaptive population sizing framework which can be combinedwith any EMMA is needed for a totally fair comparison.
Although we focused on Type1-MMOPs, we believe that
ADA is applicable to Type2-MMOPs using a relaxed dominance relation or a relaxed equivalent relation. An analysis
of ADA on Type2-MMOPs is an avenue for future work.
Other post-processing methods that handle the diversity in boththe objective and solution spaces may be better than the two
post-processing methods in ADA. Further analysis of postprocessing methods is a future research topic.
ACKNOWLEDGEMENT
This work was supported by National Natural Science
Foundation of China (Grant No. 61876075), the Programfor 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), and the Program for University Key Laboratory of Guangdong Province (Grant No.
2017KSYS008).
R EFERENCES[1] K. Miettinen, Nonlinear Multiobjective Optimization. Springer, 1998.
[2] 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.
[3] E. Zitzler and S. Künzli, “Indicator-based selection in multiobjective
search,” in PPSN, 2004, pp. 832–842.
[4] Q. Zhang and H. Li, “MOEA/D: A multiobjective evolutionary algorithm
based on decomposition,” IEEE TEVC, vol. 11, no. 6, pp. 712–731, 2007.
[5] 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.
[6] A. Trivedi, D. Srinivasan, K. Sanyal, and A. Ghosh, “A Survey of Multiobjective Evolutionary Algorithms Based on Decomposition,” IEEE
TEVC, vol. 21, no. 3, pp. 440–462, 2017.
[7] M. Sebag, N. Tarrisson, O. Teytaud, J. Lefèvre, and S. Baillet, “A14Multi-Objective Multi-Modal Optimization Approach for Mining Stable
Spatio-Temporal Patterns,” in IJCAI, 2005, pp. 859–864.
[8] G. Rudolph and M. Preuss, “A multiobjective approach for finding equivalent inverse images of pareto-optimal objective vectors,” in MCDM,
2009, pp. 74–79.
[9] 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.
[10] 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.
[11] 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.
[12] 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.
[13] 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.
[14] T. Ulrich, J. Bader, and L. Thiele, “Defining and Optimizing IndicatorBased Diversity Measures in Multiobjective Search,” in PPSN, 2010,
pp. 707–717.
[15] 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.
[16] C. Yue, B. Qu, and J. Liang, “A Multiobjective Particle Swarm Optimizer Using Ring Topology for Solving Multimodal Multiobjective
Problems,” IEEE TEVC, vol. 22, no. 5, pp. 805–817, 2018.
[17] R. Tanabe and H. Ishibuchi, “A Decomposition-Based Evolutionary
Algorithm for Multi-modal Multi-objective Optimization,” in PPSN,
2018, pp. 249–261.
[18] Y. Yuan, H. Xu, B. Wang, and X. Yao, “A New Dominance RelationBased Evolutionary Algorithm for Many-Objective Optimization,” IEEE
TEVC, vol. 20, no. 1, pp. 16–37, 2016.
[19] Z. Wang, Q. Zhang, A. Zhou, M. Gong, and L. Jiao, “Adaptive
Replacement Strategies for MOEA/D,” IEEE Trans. Cyber., vol. 46,
no. 2, pp. 474–486, 2016.
[20] Y. Yuan, H. Xu, B. Wang, B. Zhang, and X. Yao, “Balancing Convergence and Diversity in Decomposition-Based Many-Objective Optimizers,” IEEE TEVC, vol. 20, no. 2, pp. 180–198, 2016.
[21] S. Jiang, S. Yang, Y. Wang, and X. Liu, “Scalarizing Functions in
Decomposition-Based Multiobjective Evolutionary Algorithms,” IEEE
TEVC, vol. 22, no. 2, pp. 296–313, 2018.
[22] R. Cheng, Y. Jin, M. Olhofer, and B. Sendhoff, “A Reference Vector Guided Evolutionary Algorithm for Many-Objective Optimization,”
IEEE TEVC, vol. 20, no. 5, pp. 773–791, 2016.
[23] 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.
[24] Y. Liu, G. G. Yen, and D. Gong, “A Multi-Modal Multi-Objective Evolutionary Algorithm Using Two-Archive and Recombination Strategies,”
IEEE TEVC, 2018 (in press).
[25] K. Deb and S. Tiwari, “Omni-optimizer: A Procedure for Single andMulti-objective Optimization,” in EMO, 2005, pp. 47–61.
[26] H. Li and Q. Zhang, “Multiobjective Optimization Problems With
Complicated Pareto Sets, MOEA/D and NSGA-II,” IEEE TEVC, vol. 13,
no. 2, pp. 284–302, 2009.
[27] K. Deb and R. B. Agrawal, “Simulated Binary Crossover for Continuous
Search Space,” Complex Systems, vol. 9, no. 2, 1995.
[28] R. Storn and K. Price, “Differential Evolution - A Simple and Efficient
Heuristic for Global Optimization over Continuous Spaces,” J. Global
Optimiz., vol. 11, no. 4, pp. 341–359, 1997.
[29] Y. Qi, X. Ma, F. Liu, L. Jiao, J. Sun, and J. Wu, “MOEA/D with Adaptive
Weight Adjustment,” Evol. Comput., vol. 22, no. 2, pp. 231–264, 2014.
[30] Y. Xiang, Y. Zhou, M. Li, and Z. Chen, “A Vector Angle-Based Evolutionary Algorithm for Unconstrained Many-Objective Optimization,”
IEEE TEVC, vol. 21, no. 1, pp. 131–152, 2017.
[31] S. Jiang and S. Yang, “A Strength Pareto Evolutionary Algorithm
Based on Reference Direction for Multiobjective and Many-Objective
Optimization,” IEEE TEVC, vol. 21, no. 3, pp. 329–346, 2017.
[32] 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.
[33] 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.
[34] 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.
[35] S. Bird and X. Li, “Adaptively choosing niching parameters in a PSO,”
in GECCO, 2006, pp. 3–10.
[36] M. Preuss, “Improved Topological Niching for Real-Valued Global
Optimization,” in EvoApplications, 2012, pp. 386–395.
[37] P. Kerschke, H. Wang, M. Preuss, C. Grimme, A. H. Deutz, H. Trautmann, and M. Emmerich, “Towards analyzing multimodality of continuous multiobjective landscapes,” in PPSN, 2016, pp. 962–972.
[38] G. Rudolph, B. Naujoks, and M. Preuss, “Capabilities of EMOA to
Detect and Preserve Equivalent Pareto Subsets,” in EMO, 2007, pp. 36–
50.
[39] K. Bringmann, T. Friedrich, and P. Klitzke, “Two-dimensional subset
selection for hypervolume and epsilon-indicator,” in GECCO, 2014, pp.
589–596.
[40] H. K. Singh, K. S. Bhattacharjee, and T. Ray, “Distance based subset
selection for benchmarking in evolutionary multi/many-objective optimization,” IEEE TEVC, 2018 (in press).
[41] H. Ishibuchi, Y. Setoguchi, H. Masuda, and Y. Nojima, “How to compare
many-objective algorithms under different settings of population andarchive sizes,” in IEEE CEC, 2016, pp. 1149–1156.
[42] H. Ishibuchi, H. Masuda, Y. Tanigaki, and Y. Nojima, “Modified
Distance Calculation in Generational Distance and Inverted Generational
Distance,” in EMO, 2015, pp. 110–125.
[43] 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.
[44] Q. Zhang, W. Liu, and H. Li, “The performance of a new version of
MOEA/D on CEC09 unconstrained MOP test instances,” in IEEE CEC,
2009, pp. 203–208.
[45] K. Li, Q. Zhang, S. Kwong, M. Li, and R. Wang, “Stable MatchingBased Selection in Evolutionary Multiobjective Optimization,” IEEE
TEVC, vol. 18, no. 6, pp. 909–923, 2014.
[46] M. Asafuddoula, T. Ray, and R. A. Sarker, “A Decomposition-Based
Evolutionary Algorithm for Many Objective Optimization,” IEEE TEVC,
vol. 19, no. 3, pp. 445–460, 2015.
[47] 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.
[48] H. Wang, L. Jiao, and X. Yao, “Two Arch2: An Improved Two-Archive
Algorithm for Many-Objective Optimization,” IEEE TEVC, vol. 19,
no. 4, pp. 524–541, 2015.
[49] X. Ma, F. Liu, Y. Qi, X. Wang, L. Li, L. Jiao, M. Yin, and M. Gong,“A Multiobjective Evolutionary Algorithm Based on Decision Variable
Analyses for Multiobjective Optimization Problems With Large-Scale
Variables,” IEEE TEVC, vol. 20, no. 2, pp. 275–298, 2016.
[50] M. Preuss, B. Naujoks, and G. Rudolph, “Pareto Set and EMOA
Behavior for Simple Multimodal Multiobjective Functions,” in PPSN,
2006, pp. 513–522.
[51] J. J. Liang, C. T. Yue, and B. Y. Qu, “Multimodal multi-objective
optimization: A preliminary study,” in IEEE CEC, 2016, pp. 2454–2461.
[52] 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.
[53] J. Bader and E. Zitzler, “HypE: An Algorithm for Fast HypervolumeBased Many-Objective Optimization,” Evol. Comput., vol. 19, no. 1, pp.
45–76, 2011.
[54] J. J. Durillo and A. J. Nebro, “jMetal: A Java framework for multiobjective optimization,” AES, vol. 42, no. 10, pp. 760–771, 2011.
[55] T. Ulrich, J. Bader, and E. Zitzler, “Integrating decision space diversity
into hypervolume-based multiobjective search,” in GECCO, 2010, pp.
455–462.
[56] H. Ishibuchi, Y. Setoguchi, H. Masuda, and Y. Nojima, “Performanceof Decomposition-Based Many-Objective Algorithms Strongly Depends
on Pareto Front Shapes,” IEEE TEVC, vol. 21, no. 2, pp. 169–190, 2017.
[57] H. Ishibuchi, Y. Peng, and K. Shang, “A Scalable Multimodal Multiobjective Test Problem,” in IEEE CEC, 2019, pp. 302–309.
[58] R. Tanabe, H. Ishibuchi, and A. Oyama, “Benchmarking Multi- andMany-objective Evolutionary Algorithms under Two Optimization Scenarios,” IEEE Access, vol. 5, pp. 19 597–19 619, 2017.
[59] T. Chugh, K. Sindhya, J. Hakanen, and K. Miettinen, “A survey on
handling computationally expensive multiobjective optimization problems with evolutionary algorithms,” Soft Comput., vol. 23, no. 9, pp.
3137–3166, 2019.

15

Ryoji Tanabe is a Research Assistant Professor
with Department of Computer Science and Engineering, Southern University of Science and Technology, China. He was a Post-Doctoral Researcher
with ISAS/JAXA, Japan, from 2016 to 2017. Hereceived his Ph.D. in Science from The University ofTokyo, Japan, in 2016. His research interests include
stochastic single- and multi-objective optimization
algorithms, parameter control in evolutionary algorithms, and automatic algorithm configuration.

Hisao Ishibuchi (M’93-SM’10-F’14) received the
B.S. and M.S. degrees in precision mechanics from
Kyoto University, Kyoto, Japan, in 1985 and 1987,
respectively, and the Ph.D. degree in computer science from Osaka Prefecture University, Sakai, Osaka, Japan, in 1992. He was with Osaka Prefecture University in 1987-2017. Since 2017, he isa Chair Professor at Southern University of Science and Technology, China. His research interests include fuzzy rule-based classifiers, evolutionary multi-objective and many-objective optimization,
memetic algorithms, and evolutionary games.
Dr. Ishibuchi was the IEEE Computational Intelligence Society (CIS) VicePresident for Technical Activities (2010-2013), an AdCom member of theIEEE CIS (2014-2019), and the Editor-in-Chief of the IEEE COMPUTATIONAL INTELLIGENCE MAGAZINE (2014-2019).