Noname manuscript No.
(will be inserted by the editor)

Review and Analysis of Three Components of Differential
Evolution Mutation Operator in MOEA/D-DE
Ryoji Tanabe · Hisao Ishibuchi

Received: date / Accepted: date

Abstract Abstract A decomposition-based multi-objective
Keywords Multi-objective optimization · Decomevolutionary algorithm with a differential evolution variposition based evolutionary algorithms · Differential
ation operator (MOEA/D-DE) shows high performanceevolution operators · Implementation of algorithms
on challenging multi-objective problems (MOPs). The
DE mutation consists of three key components: a mutation strategy, an index selection method for parent in1 Introduction
dividuals, and a bound-handling method. However, theconfiguration of the DE mutation operator that should
A multi-objective evolutionary algorithm (MOEA) is
be used for MOEA/D-DE has not been thoroughly inan efficient approach for solving multi-objective optivestigated in the literature. This configuration choicemization problems (MOPs) [10]. Since MOEAs use a set
confuses researchers and users of MOEA/D-DE. To adof individuals for the search, it is expected that well disdress this issue, we present a review of the existing contributed nondominated solutions can be found by a sinfigurations of the DE mutation operator in MOEA/Dgle run. MOEAs have been successfully applied to realDE and systematically examine the influence of eachworld problems, such as aerodynamic wing design probcomponent on the performance of MOEA/D-DE. Ourlems [30], vehicle design problems [28], oil well problems
review reveals that the configuration of the DE muta[24], and groundwater monitoring design problems [22].
tion operator differs depending on the source code of
A variation operator is used to generate a child in
MOEA/D-DE. In our analysis, a total of 30 configuMOEAs. Although genetic algorithm (GA) operators
rations (three index selection methods, two mutation
are frequently used in the evolutionary multi-objective
strategies, and five bound handling methods) are invesoptimization community, some recent studies report the
tigated on 16 MOPs with up to five objectives. Results
superiority of differential evolution (DE) operators [35].
show that each component significantly affects the perFor example, the performance of NSGA-II [11], SPEA2
formance of MOEA/D-DE. We also present the most[46], and IBEA [45] with GA and DE operators is invessuitable configuration of the DE mutation operator,
tigated in [37]. The results reported in [37] show thatwhich maximizes the effectiveness of MOEA/D-DE.
better nondominated solutions can be found by using
the DE operator. A further examination of the impact
Keywords Multi-objective optimization · Decomof both of these operators on MOEAs is presented in [5].
position based evolutionary algorithms · Differential
In [42], the performance of NSGA-III [12] with GA andevolution operators · Implementation of algorithms
DE operators is investigated on MOPs with more thanfour objectives. A comparison of SMS-EMOA [4] with
R. Tanabe and H. Ishibuchi
Shenzhen Key Laboratory of Computational Intelligence,various operators is performed in [2]. The results preUniversity Key Laboratory of Evolving Intelligent Systems ofsented in [2] show that the DE operator is more suitable
Guangdong Province, Department of Computer Science andthan the GA operator for SMS-EMOA in most cases.
Engineering, Southern University of Science and Technology,
MOEA/D-DE [25] is one of the most successful exShenzhen 518055, China
E-mail: rt.ryoji.tanabe@gmail.com, hisao@sustc.edu.cn (Coramples of the DE variation operator. MOEA/D [43,responding author: Hisao Ishibuchi)
36] is an MOEA that decomposes a given MOP into

2

Ryoji Tanabe, Hisao Ishibuchi

multiple single-objective sub-problems and solves theminfluence of configurations of the DE mutation operasimultaneously. MOEA/D-DE is a variant of MOEA/D
tor on MOEA/D-DE has not been thoroughly investiin which the GA operator is replaced by the DE operagated. Since the configuration of the DE mutation opertor1 . The results reported in [25] show that MOEA/Dator that should be used is unclear, algorithm designers
DE is capable of handling complicated Pareto solution
and users of MOEA/D-DE face confusion regarding thesets. Some improved MOEA/D-DE algorithms have been configuration choice. In fact, different implementations
proposed in the literature, and representative examare used in the uploaded source code of MOEA/D-DE.
ples include MOEA/D-DRA [44], MOEA/D-MAB [26],For example, while the current/1 strategy is incorpoMOEA/D-STM [27], MOEA/D-GRA [40], and MOEA/D- rated into MOEA/D-DE in the jMetal source code2 ,TPN [21].
the rand/1 strategy is used in the MOEA Framework
source code3 . Different comparison results could be obIn general, the term “DE operator” denotes a comtained simply depending on the choice of a library if
plex differential mutation operator and a crossover meMOEA/D-DE in an MOEA library performs signifithod [35]. In the procedure of the DE operator, a mucantly better than that in another MOEA library. Thistant vector is first generated by applying the DE musituation is undesirable for researchers. Although users
tation operator to some individuals in the population.
of MOEAs generally want to apply a well-performing
Then, a child is generated by recombining the mutant
algorithm to their MOPs, the best configuration of thevector and a parent individual. Binomial crossover is
DE mutation operator for MOEA/D-DE is unknown.
widely used in MOEA/D-DE-type algorithms. HowTo address these issues, we present a review of theever, to handle nonseparability, the procedure of theexisting configurations of the DE mutation in MOEA/Dcrossover operator is not actually performed in mostDE. We also analyze the influence of each element on
MOEA/D-DE-type algorithms by setting the crossoverthe effectiveness of MOEA/D-DE. The purpose of thisrate C ∈ [0, 1] to 1 [25].
review is to clarify the current situation. We examine
In contrast to the crossover method, there is large
the three components that can be found in the litflexibility in implementing the DE mutation operator.
erature, source code, and MOEA libraries. Then, we
The DE mutation operator consists of the following
present a large-scale experimental study to investigate
three key components:
the effect of the three components on MOEA/D-DE. A– A type of mutation strategy (e.g., the rand/1, rand/2, total of 30 configurations (three index selection methand current-to-rand/1 strategies),
ods, two mutation strategies, and five bound-handling– An index selection method for parent individuals
methods) are investigated in our study. The influence(i.e., how to select indices r1 , r2 , and r3 ), andof each component on the performance of MOEA/D– A bound-handling method to repair an infeasible
DE is carefully examined in a component-wise manner.
solution.
We also suggest the best configuration of the DE mutation operator, which maximizes the performance ofOn the one hand, previous studies [1, 29, 32, 39] reMOEA/D-DE.
veal that the performance of DE algorithms for singleThe reminder of this paper is organized as follows:
objective optimization is significantly influenced by the
Section 2 describes the background of MOPs and MOEA/Dabove-listed three components. For example, the results
DE. Section 3 reviews existing configurations of the DE
presented in [1] show that the choice of bound-handling
mutation operator in MOEA/D-DE. Section 4 describes
methods has a large impact on the performance of theexperimental settings. Section 5 shows results of MOEA
basic DE [35]. For this reason, most researchers in the/D-DE with various configurations of the DE mutation
DE community carefully select each component of theoperator. Finally, Section 6 concludes this paper.
DE mutation operator to maximize the performance ofDE.
2 Preliminaries
On the other hand, in contrast to the DE community, the details of the DE mutation operator have
2.1 Definition of MOPs
not received considerable attention in the evolutionary multi-objective optimization community. Thus, theA bound-constrained MOP, which is addressed in this1
paper, can be formulated as follows:
Strictly speaking, the differences between MOEA/D-DE[25] and the original MOEA/D [43] are as follows: (i) the parT
minimize f (x) = f1 (x), ..., fM (x) ,(1)
ent individuals are selected from the whole population withsome probability, (ii) the number of individuals replaced by achild is restricted, and (iii) the DE variation operator is usedin [25].

23

http://jmetal.sourceforge.net/
http://moeaframework.org/index.html

Review and Analysis of Three Components of Differential Evolution Mutation Operator in MOEA/D-DE

where S is the solution space. f : S → RM is an objective function vector that consists of M potentiallyconflicting objective functions, and RM is the objective
space. In (1), x = (x1 , ..., xD )T is a solution vector,
QD
max
and S = j=1 [xmin] is the bound-constrained soj , xj
minlution space where xj ≤ xj ≤ xmax
for each indexjj ∈ {1, ..., D}.
We say that x1 dominates x2 if and only if fi (x1 ) ≤fi (x2 ) for all i ∈ {1, ..., M } and fi (x1 ) < fi (x2 ) for atleast one index i. Here, x∗ is a Pareto-optimal solution
if no x ∈ S exists such that x dominates x∗ . In thiscase, f (x∗ ) is a Pareto-optimal objective vector. The
set of all x∗ in S is the Pareto-optimal solution set,
and the set of all f (x∗ ) in RM is the Pareto front. In
general, no solution can simultaneously minimize all objective functions f1 , ..., fM in an MOP. Thus, the goalof multi-objective optimization is to find a set of nondominated solutions that are well distributed and close
to the Pareto front in the objective space.

3

Algorithm 1: The procedure of MOEA/D-DE
1
2
3
4
5
6
7
8
9
10

11

12
13

14
15
16
17
18

19
20

21

t ← 1, initialize the population P = {x1 , ..., xµ };for i ∈ {1, ..., µ} do
Set the neighborhood index list B i = {i1 , ..., iT };while The termination criteria are not met dofor i ∈ {1, ..., µ} do
if rand[0, 1] ≤ δ thenR ← Bi ;else
R ← {1, ..., µ};Select parent indices from R with an index
selection method (Subsection 3.2);
Generate the mutant vector v i using a mutation
strategy (Subsection 3.1);
if v i ∈/ S thenRepair v i using a bound-handling method(Subsection 3.3);
Generate the child ui by crossing xi and v i ;Apply a GA mutation operator to ui ;c ← 1;
while c ≤ nrep and R 6= ∅ do
Randomly select an index j from R, andR ← R\{j};if g(ui |wj , z ∗ ) ≤ g(xj |wj , z ∗ ) thenxj ← ui , c ← c + 1;t ← t + 1;

2.2 MOEA/D-DE
MOEA/D-type algorithms, including MOEA/D-DE [25],decompose an M -objective MOP defined in (1) into
µ single-objective sub-problems g1 (x|w1 ), ..., gµ (x|wµ )using a set of uniformly distributed weight vectors W ={w1 , ..., wµ } and a scalarizing function g : RM → R,
i T
where µ is the population size, wi = (w1i , ..., wM)PMifor each i ∈ {1, ..., µ}, and
j=1 wj = 1. For eachi ∈ {1, ..., µ}, an individual xi is assigned to the ith sub-problem. MOEA/D-type algorithms attempt tofind the optimal solutions of all sub-problems simultaneously.
Algorithm 1 shows the procedure of MOEA/D-DE.
After initialization (lines 1–3), the following steps areiteratively performed. For each individual, an index list
is selected for use in the mating and replacement selections (lines 5–9). Then, the mating and reproduction
operations are performed (lines 10–15). Finally, the replacement selection is applied to the child and the individuals in the population (lines 17–20). We explain
each step of MOEA/D-DE in detail below.
At the beginning of the search, all individuals in the
population P = {x1 , ..., xµ } are randomly generated
in the solution space (line 1). For each sub-problem
index i ∈ {1, ..., µ}, an index list B i = {i1 , ..., iT } isinitialized (lines 2-3): B i consists of indices of the T
closest weight vectors to wi in the weight vector space,
where T is the neighborhood size. After initialization,
the following steps (lines 5–20) are repeatedly applied
to each sub-problem until a termination condition issatisfied.

For each i, a set of individual indices R is set to
B i with a probability of δ ∈ [0, 1] or {1, ..., µ} with aprobability of 1 − δ (lines 6–9). The function rand[0, 1]
in line 6 is a randomly chosen value in the range [0, 1].
After R has been determined, MOEA/D-DE generates
a mutant vector v i (lines 10–13). First, individual indices {r1 , r2 , ...} are randomly selected from R using an
index selection method (line 10). Next, v i is generated
by applying a mutation strategy to the selected individuals {xr1 , xr2 , ...} (line 11). If an element of v i violatesa corresponding bound constraint (i.e., v i ∈/ S), then abound-handling method is applied to it such that v i ∈ S(line 13). The details of each procedure (the index selection, the mutation strategy, and the bound-handling
method) are described in Subsections 3.2, 3.1, and 3.3,
respectively.
After v i is repaired by a bound-handling method (ifneeded), a child ui is generated by recombining xi andv i (line 14). Binomial crossover [35], which is the mostbasic crossover method in DE, is defined as follows:(
vji if rand[0, 1] ≤ C or j = j rand
i
uj =
,
(2)
xij otherwise
where the crossover rate C ∈ [0, 1] in (2) controls thenumber of inherited variables from xi to ui . The decision variable index j rand in (2) is randomly selectedfrom {1, ..., D}. Since binomial crossover only exchanges
elements between the parent individual xi and the mutant vector v i , the child ui always satisfies the boundconstraints as long as xi , v i ∈ S.

4

Ryoji Tanabe, Hisao Ishibuchi

why many configurations exist in the source code of
MOEA/D-DE. Table 1 shows 15 configurations of theDE mutation operator in the source code of MOEA/DDE and its variants. We downloaded the source codefrom each website in Table 1 and carefully checked how
the differential mutation was implemented. The components of the DE mutation described here include twomutation strategies, three index selection methods forparent individuals, and five bound-handling methods.
max These components are explained in Subsections 3.1, 3.2,
where pmut ∈ [0, 1] is the mutation rate. If uij ∈/ [xmin],j , xjand 3.3, respectively.
then it is replaced by the closest value (xminor xmax).
jj
The amount of perturbation σj in (3) is defined as follows:(
3.1 Two DE mutation strategies in MOEA/D-DE
mut(2 rand[0, 1])1/(η +1) − 1
if rand[0, 1] ≤ 0.5
,
σj =
mut
1 − (2 − 2 rand[0, 1])1/(η +1) otherwise
For each target individual xi (i ∈ {1, ..., µ}) in the population P , MOEA/D-DE performs the differential mu(4)tation to generate a mutant vector v i (line 11 in Alwhere η mut > 0 in (4) is the distribution index.
gorithm 1). A number of mutation strategies exist in
After ui is generated, the replacement procedure
DE [8, 9]. Among these strategies, the following rand/1
is performed using a predefined scalarizing function gstrategy is the most basic mutation strategy:(lines 17–20). First, an index j is randomly selectedfrom R, and j is removed from R (line 18). Then, thev i = xr1 + F (xr2 − xr3 ),(6)ji
individual x is compared with the child u based on g
and the weight vector wj (line 19). If ui is better thanwhere r1 , r2 , and r3 are randomly selected from a set
xj according to their scalarizing function values, thenof individual indices R (see lines 6–9 in Algorithm 1)
xj is replaced with ui (line 20). Unlike the original
according to the index selection method, which is exMOEA/D proposed in 2007 [43], the number of indiplained in Subsection 3.2 later. The rand/1 mutation
viduals replaced by the child is limited to nrep > 0 in
strategy is widely used in DE algorithms for singleMOEA/D-DE (line 17). c is used to count the number
objective optimization and multi-objective optimization.
of individuals replaced by ui . When c reaches the maxFor example, the rand/1 operator is incorporated into
imum number of replacements nrep , the replacement
representative multi-objective DE algorithms, such asprocedure terminates (line 17).
GDE3 [23] and DEMO [33].
Since the replacement criterion is based on the scalarAlthough the rand/1 operator is commonly used inizing function g (line 19), the performance of MOEA/Dsingle- and multi-objective DE, the current/1 mutation
type algorithms significantly depends on g [18]. Although strategy is incorporated into MOEA/D-DE. Note that
there are a number of scalarizing functions, as reviewed
the term “current/1” is introduced in [14] and is notin [31], the following Tchebycheff function (g tch ) [43] isused in the original MOEA/D-DE paper [25]. The curused in MOEA/D-DE:
rent/1 mutation strategy is defined as follows:
In MOEA/D-DE, a GA mutation operator is applied to the child ui (line 15). The following polynomial
mutation, which is commonly used in the evolutionary
multi-objective optimization community, is performed
for each j ∈ {1, ..., D}:(
uij + σj (xmax− xminif rand[0, 1] ≤ pmut
j
j )
i
uj =
,
uij
otherwise(3)

g tch (x|w, z ∗ ) =
tch

max

{wi |fi (x) − zi∗ |},(5)i∈{1,...,M }∗= (z1∗ , ...,

where g in (5) should be minimized. The z∗ TzM) is the ideal point. Since obtaining the true ideal
point z ∗ for a given MOP is difficult, its approximatedpoint, which consists of the minimum function value
for each objective fi (i ∈ {1, ..., M }) found during thesearch process, is typically used.
3 A review of existing configurations of the DE
mutation operator in MOEA/D-DE
Here, we review the existing configurations of the DE
mutation operator in MOEA/D-DE. We also explainv i = xi + F (xr1 − xr2 ),(7)

where r1 and r2 are randomly selected from R. Although the base vector is a randomly selected individual from R in the rand/1 strategy, it is always identicalto the target individual xi in the current/1 strategy.
Compared to the rand/1 strategy, the mutant vector v igenerated by the current/1 strategy is generally close
to xi [14].
As mentioned above, MOEA/D-DE uses the current/1 strategy for the mutant vector generation. However, in the MOEA/D-DE paper [25], it is not explicitly
stated that the current/1 mutation strategy in (7) isincorporated into MOEA/D-DE. Instead, in [25], it isReview and Analysis of Three Components of Differential Evolution Mutation Operator in MOEA/D-DE

5

Table 1: Configurations of the DE mutation operator in the source code of MOEA/D-DE and its variants. Eightconfiguration numbers from #A to #H are used in Subsection 5.5.

MOEAs

Languages

Mutationstrategies

Index
selection
methods

Boundhandling
methods

current/1

WR

reinitializationacurrent/1

WPR

replacement#Ba

Web sites

MOEA/D-DE
described in [25]MOEA/D-DEC++

MOEA/D homepage

Configuration
number ##A

MOEA/D-DE

MATLAB

MOEA/D homepage

rand/1WOR

replacement#C

MOEA/D-DE

Java

MOEA/D homepagea rand/1

WOR

replacement#C

MOEA/D-DE
MOEA/D-DRA

Java

jMetal 4.5bcurrent/1

WPR

replacement#B

MOEA/D-DRA

Java

MOEA Frameworkc

rand/1

WPR

replacement#Dcurrent/1

WOR

r-reflection#E

current/1

WOR

no method-

current/1

WPR

replacement#Bcurrent/1

WOR

replacement#F

MOEA/D-DE
MOEA/D-DE
MOEA/D-STM
MOEA/D-STM

C++
MATLABJava
MATLABPaGMO

dPlatEMOe

Li’s website

fLi’s websitefg

MOEA/D-DE

C

COCO website

MOEA/D-DRAC++

MOEA/D homepagea current/1

ENS-MOEA/D
MOEA/D-TPN

MATLABC

Suganthan’s website

rand/1hiYang’s websiteWORr-reflection

#G

WPR

r-reflection#H

current/1

WPR

r-reflection#H

current/1

WPR

r-reflection#Ha

http://dces.essex.ac.uk/staff/zhang/webofmoead.htm
http://jmetal.sourceforge.net/c
http://moeaframework.org/index.html
dhttps://esa.github.io/pagmo2/index.html
ehttp://bimk.ahu.edu.cn/index.php?s=/Index/Software/index.html
fhttp://www.cs.bham.ac.uk/~likw/publications.html
ghttp://coco.gforge.inria.fr/doku.php?id=mo-gecco2015h
http://www3.ntu.edu.sg/home/EPNSugan/index_files/
i
http://www.tech.dmu.ac.uk/~syang/publications.htmlbreported that the rand/1 strategy is used in MOEA/DDE (equation (6) on page 9 in [25]), but the randomly
chosen index r1 is replaced by the target individual index i (Step 2.2 on page 9 in [25]). The modified rand/1
strategy is identical to the current/1 strategy.
The complicated description in [25] is somewhat difficult to understand. At first glance, one may think thatthe well-known rand/1 strategy is used in MOEA/DDE. Unless readers have carefully read all the descriptions of the DE variation operator in MOEA/D-DE,
then they will not notice that the current/1 strategy isactually used in MOEA/D-DE. In fact, the descriptionsin [25] have confused some researchers. Consequently,
the rand/1 strategy is incorporated into MOEA/D-DE
in some source code. For example, as shown in Table
1, MOEA/D-DE in MOEA Framework and the sourcecode of [6] use the rand/1 strategy, rather than thecurrent/1 strategy. Additionally, some authors describethat the rand/1 mutation strategy was used in MOEA/D-

DE-type algorithms in their articles, but the current/1
strategy was actually used in their source code. MOEA/DTPN [21] in Table 1 is such an example. Since therand/1 and current/1 strategies are essentially different, it is expected that the performance of MOEA/DDE with the two strategies is different.

3.2 Index selection methods for the parent individuals
Here, we explain three index selection methods (WOR,WR, and WPR methods) for the parent individuals in
the DE mutation operator. Algorithms 2, 3, and 4 show
the procedures of the WOR, WR, and WPR methods
for selecting two individual indices r1 and r2 for thecurrent/1 mutation strategy. In Algorithm 2 (WOR),
r1 and r2 are randomly selected such that r1 6= r2 ,r1 6= i, and r2 6= i. In contrast, i, r1 , and r2 can be thesame index (i.e., r1 = r2 = i) in Algorithm 3 (WR). In
Algorithm 4 (WPR), r1 and r2 differ from each other,

6
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

Ryoji Tanabe, Hisao Ishibuchi

void CMOEAD :: matings el e ct io n ( vector < int > & list , int cid , int size , int type ){// list : the set of the indexes of s e l e c t e d mating parents// cid : the id of current s u b p r o b l e m// size : the number of s e l e c t e d mating parents// type : 1 - n e i g h b o r h o o d ; o t h e r w i s e - whole p o p u l a t i o n
int ss= population [ cid ]. table . size () , r , p ;while ( list . size () < size ){
if ( type ==1){
r = int ( ss * rnd_uni (& rnd_uni_init ));p = population [ cid ]. table [ r ];}else
p = int ( population . size ()* rnd_uni (& rnd_uni_init ));bool flag = true ;for ( int i =0; i < list . size (); i ++){if ( list [ i ]== p ) // p is in the list{
flag = false ;break ;}}
if ( flag ) list . push_back ( p );}}Fig. 1: The “matingselection” function in the original C++ code of MOEA/D-DE, which is available on theauthors’ website. To implement the WOR selection method, line 19 should have been “if(list[i]==p && cid==p)”.

Algorithm 2: The WOR selection method.
1

doRandomly select r1 from R;
3 while r1 = i;
4 do
5Randomly select r2 from R;
6 while r2 = i and r2 = r1 ;2

Algorithm 3: The WR selection method.
1
2

Randomly select r1 from R;
Randomly select r2 from R;

Algorithm 4: The WPR selection method.
Randomly select r1 from R;
do3
Randomly select r2 from R;
4 while r2 = r1 ;1

2

but either of them can be equal to i. In the following,
the three index selection methods are described in detail.
In the original DE paper [35], the parent individual
indices r1 , r2 , and r3 for the rand/1 strategy are randomly selected from {1, ..., µ}\{i} such that they differ
from each other (i.e., r1 6= r2 , r1 6= r3 , r2 6= r3 , andrj 6= i for j ∈ {1, 2, 3}). That is, individual indices are
selected without replacement. In this paper, the traditional index selection method in DE is called a WOR
method (Algorithm 2). Although the WOR selection
method is generally used in DE, such a restriction is
not described in the original MOEA/D-DE paper [25].
In this case, the selected indices are possibly equal toeach other and also equal to i. The method that selects
individual indices with replacement is called a WR method (Algorithm 3).
However, in the original C++ code of MOEA/DDE, r1 and r2 for the current/1 strategy are selectedfrom a set of individual indices R (see Subsection 2.2)
such that they differ from each other, but there is no re-striction about i. Figure 1 shows the matingselection
function in the original C++ source code of MOEA/DDE, which is available on the authors’ website. The
matingselection function is for selecting indices forthe differential mutation. We do not explain the procedure of the matingselection function in detail, butthere is a bug in line 19 in Figure 1. The one-dimensional
vector list is for storing indices {r1 , r2 , ...} that havebeen selected in the matingselection function. In line
19, p is a randomly selected candidate that is able toenter list. Since it is checked whether p exists in list,
p may equal cid, where cid is the index of the target
sub-problem and identical to i in Algorithm 1. Thus,
while it is ensured that r1 differs from r2 , either of thetwo indices may equal i. The selection method withpartial replacement used in the original MOEA/D-DE
code is called a WPR method (Algorithm 4) in this pa-

Review and Analysis of Three Components of Differential Evolution Mutation Operator in MOEA/D-DE

per. The WPR selection method has never been foundin the DE literature.
We carefully checked the index selection method ineach source code in Table 1. As shown in Table 1, it
appears that the WR method has hardly been used in
MOEA/D-DE. The WPR selection method is used for
MOEA/D-DE in some MOEA packages (e.g., jMetal
and MOEA Framework), which may be because they
were based on the original MOEA/D-DE source code.
The WOR method is used in other source code (e.g.,
PaGMO and PlatEMO), which is because the WOR
is the most standard index selection method in the
DE community. If researchers attempt to implement
MOEA/D-DE without referring to its original C++
source code, most of them are likely to incorporate the
WOR method into MOEA/D-DE for this reason.
Wang et al. [39] investigate the impact of the index selection method on the performance of DE algorithms for single-objective optimization. Their results
show that the performance of some DE algorithms canbe significantly improved or degraded by replacing the
traditional WOR method with the WR method. Since
removing the restriction r1 6= r2 increases the probability that a child is generated near the base vector, the useof the WR method makes the search more greedy [39].
Thus, it is expected that the performance of MOEA/DDE is affected by the choice of index selection method.

3.3 Bound constraint-handling methods
It is not always ensured that the mutant vector v generated by the differential mutation is in the solution
D
max
space S = Πj=1
[xmin]. When an element of thej , xjmutant vector vj is out of the bound (j ∈ {1, ..., D}),a bound-handling method must be applied to v suchthat v ∈ S (line 13 in Algorithm 1). In our study, we
consider the following five bound-handling methods (resampling, replacement, reinitialization, reflection, andr-reflection methods), which are commonly used in theDE literature:• Resampling: The mutant vector v is repeatedly
generated using the differential mutation with different parent individuals until all elements of v satisfy thebound constraints.
• Replacement: In the case where vj is out of thebound, it is replaced with the corresponding minimum/maximum value (xminor xmax) as (8):
j
j

min

x j
vj = xmaxjvj

if vj < xminj.
if vj > xmaxjotherwise

(8)

7

max
], it is reini• Reinitialization: When vj ∈/ [xminj , xjmax
tialized in the range [xmin
,
x
]
as
(9):
j
j

min

if vj < xmin
− xmin(xmaxjj ) rand[0, 1] + xj
jor
vj > xmax
.
vj =
j

v
otherwisej(9)• Reflection: If vj violates the j-th bound constraint,
it is reflected as much as it exceeded the minimum/maximum value as (11):

min
min
min

xj + (xj − vj ) if vj < xj.
(10)
vj = xmax− vj ) if vj > xmax
+ (xmax
j
j
j


vj
otherwise• Randomized reflection (r-reflection): This method is an alternative version of the reflection method.
Unlike the reflection method, the amount of reflection
is randomly determined as (11):

min
min
min

xj + rand[0, 1] (xj − vj ) if vj < xj.
vj = xmax+ rand[0, 1] (xmax− vj ) if vj > xmax
j
j
j


vj
otherwise(11)
Table 1 shows that the replacement method is usedin the original C++ source code and most source codefor MOEA/D-DE. In addition to the replacement method, the r-reflection method is used in some implementations, as shown in Table 1. Table 1 indicates thatthe reinitialization, reflection, and resampling methods have not been incorporated into MOEA/D-DE. Although it is described in [25] that the reinitialization
method is used as the repair method (Step 2.3 on page
9 in [25]), the replacement strategy is actually incorporated into the original C++ source code. While the reflection method is one of the most popular repair methods in DE algorithms for single-objective optimization,it has not been used in MOEA/D-DE. It was shown in[1] that the resampling method is most suitable for DE
for single-objective optimization. However, the resampling method has not received much attention in theevolutionary multi-objective optimization community.
The use of a repair operator in evolutionary algorithms, including DE algorithms, influences their performance [1, 15, 16, 41]. The impact of the bound-handling
method on the basic DE algorithm is investigated in [1].
The results reported in [1] indicate that the choice ofrepair method significantly influences the performanceof the basic DE [35] for single-objective continuous optimization. Thus, it is expected that the performanceof MOEA/D-DE also depends on the choice of boundhandling method.

8

Ryoji Tanabe, Hisao Ishibuchi

Table 2: Properties of the DTLZ and WFG test problems. “Mult.” and “Nonsep.” represent the multimodality and the nonseparability, respectively.
ProblemPareto front

DTLZ1

LinearDTLZ2Nonconvex

DTLZ3

Nonconvex

DTLZ4

Nonconvex

DTLZ5Partiallydegenerate

DTLZ6Partiallydegenerate

DTLZ7

Disconnected

WFG1

Mixed

WFG2

Disconnected

WFG3

Partiallydegenerate

WFG4

Nonconvex

WFG5

Nonconvex

WFG6

Nonconvex

WFG7

Nonconvex

WFG8

Nonconvex

WFG9

Mult.

Nonconvex

Nonsep.

OthersX

Table 3: Example of the APS calculation.
Problem
I1
I2
I3
I4

A1
2
0
1
2

A2
1
0
0
0

A3
0
2
1
1

APS

1.25 (5/4)

0.25 (1/4)

1 (4/4)

X
BiasedA. Before calculating the hypervolume value, the objective vector f (x) of each x ∈ A was normalized using
ideal T
the ideal point z ideal = (z1ideal , ..., zM) and the nadir
nadir
nadir
nadir T
point z= (z1, ..., zM ) of each problem:

X
BiasedX

X
X

X
Deceptive
X
BiasedX

X

BiasedX

Deceptive
Biased4 Experimental settings
This section describes our experimental settings. Results are reported in Section 5.

4.1 Test problems
The seven DTLZ [13] and nine WFG [17] test problems with M ∈ {2, 3, 4, 5} were used in our study. Table 2 summarizes their properties. Note that the Pareto
fronts of the DTLZ5, DTLZ6, and WFG3 problems
are partially degenerate [19]. According to [13], for the
DTLZ problems, the number of position variables k was
set to k = 5 for the DTLZ1 problem, k = 20 for the
DTLZ7 problem, and k = 10 for the other DTLZ problems, where the number of variables is D = M + k − 1.
Additionally, as suggested in [17], for the WFG test
problems, the number of position variables k was set tok = 2 (M − 1), and the number of distance variables lwas set to l = 20, where D = k + l.

4.2 Performance indicators
The hypervolume indicator [47] was used for evaluating
the quality of a set of obtained nondominated solutions

finormalized (x) =fi (x) − ziideal,
zinadir − ziideal(12)

where finormalized (x) is the i-th normalized objective
value (i ∈ {1, ..., M }). The i-th element ziideal of z ideal is
the minimum value of the i-th objective function of thetrue Pareto front. Conversely, the i-th element zinadirof z nadir is the maximum value of the i-th objectivefunction of the true Pareto front. That is, the objective
space of each test problem is normalized such that theideal point and the nadir point in the normalized objective space are (0, 0, ..., 0)T and (1, 1, ..., 1)T , respectively. According to [20], the reference point for the hypervolume calculation was set to (1.1, ..., 1.1)T . In thissetting, the hypervolume value is in the range [0, 1.1M ].
The average performance score (APS) [3] was used
to aggregate the results on various problems. Supposethat n algorithms A1 , ..., An are compared for a given
problem instance based on the hypervolume values obtained in multiple runs. For each i ∈ {1, ..., n} andj ∈ {1, ..., n} \{i}, if Aj significantly outperforms Ai
using the Wilcoxon rank-sum test with p < 0.05, thenδi,j = 1; otherwise, δi,j = 0. The P
performance 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 averageof the P (Ai ) values for all the considered problem instances. In other words, the APS value of Ai represents
how good (relatively) the performance of Ai is among
the n algorithms on average over all problem instances.
A small APS value indicates that the performance ofthe target algorithm is better than other compared algorithms.
Table 3 shows a simple example of the APS calculation. Let us consider the APS values of three algorithms A1 , A2 , and A3 on four problems I1 , I2 , I3 , andI4 . Each element of Table 3 shows a performance scorevalue of the corresponding algorithm. The performance
score value represents the number of algorithms that
outperform the corresponding algorithm. For example,

Review and Analysis of Three Components of Differential Evolution Mutation Operator in MOEA/D-DE

9

Table 4: Overall performance of MOEA/D-DE with the three

Table 5: Overall performance of MOEA/D-DE with the two

index selection methods on the 16 problems. The APS value
at the final iteration is shown in the tables (lower is better).
The numbers in parentheses indicate the ranks of the three
configurations based on their APS values.

mutation strategies on the 16 problems. The APS value atthe final iteration is shown in the tables (lower is better).
The numbers in parentheses indicate the ranks of the three
configurations based on their APS values.

(a) current/1
M

WOR

WR

WPR

2
3
4
5

0.88 (3)0.56 (3)0.25 (1)0.38 (3)0.12 (2)0.06 (1)0.25 (1)0.00 (1)0.06 (1)0.31 (2)0.31 (3)0.12 (2)current/1

rand/12
3
4
5

0.12 (1)0.00 (1)0.06 (1)0.06 (1)0.62 (2)0.50 (2)0.38 (2)0.19 (2)5 Experimental results(b) rand/1
M

WOR

WR

WPR

2
3
4
5

0.56 (2)0.25 (2)0.06 (1)0.06 (1)0.25 (1)0.06 (1)0.62 (3)0.44 (3)0.06 (1)0.31 (3)0.31 (3)0.12 (2)M

A1 is outperformed by two algorithms on I1 and I4 . A3
is not outperformed by any other algorithms on I1 . The
APS value for each algorithm is calculated by dividing
the total performance score value (e.g., 5 for A1 ) by the
number of problems (i.e., 4 for all algorithms). In thisexample, the APS values of A1 , A2 , and A3 are 1.25,
0.25, and 1, respectively. It can be concluded that A2
is the best performer according to the calculated APS
values.

This section describes results of MOEA/D-DE withvarious configurations of the DE mutation operator.
First, Subsection 5.1 analyzes the effect of the index selection method in MOEA/D-DE. Next, Subsection 5.2
investigates the influence of the mutation strategy on
MOEA/D-DE. Subsection 5.3 presents a comparison of
the five bound-handling methods. Subsection 5.4 examines which configuration of the DE mutation operator isthe most suitable for MOEA/D-DE. We examine the effectiveness of a total of 30 configurations of three index
selection methods (WOR, WR, and WPR), two mutation strategies (rand/1 and current/1), and five boundhandling methods (resampling, replacement, reinitialization, reflection, and r-reflection). Finally, Subsection
5.5 compares the eight existing configurations (#A, ...,#H) shown in Table 1.

5.1 The impact of the index selection method on theperformance of MOEA/D-DE
4.3 Control parameters of MOEA/D-DE
We used the source code of MOEA/D-DE downloaded
from the jMetal website for our experiments. All thecontrol parameters (except for the population size µ) for
MOEA/D-DE were set according to the original study[25]. For the DE operators, the scale factor F and thecrossover rate C were set to 0.5 and 1, respectively.
For the polynomial mutation, pmut and η mut were set
to 1/D and 20, respectively. The other parameters of
MOEA/D-DE were set as follows: T = 20, nrep = 2, andδ = 0.9. The population size µ was set to 200, 210, 220,
and 210 for M = 2, 3, 4, and 5, respectively. The simple
normalization strategy described in [43] was introduced
into MOEA/D-DE to handle differently scaled objective function values. The maximum number of function
evaluations was set to 100 000 for all test problems, and
51 independent runs were performed.

Table 4 shows the overall performance of MOEA/DDE with the three index selection methods (the WOR,WR, and WPR methods) on the 16 problems with M ∈{2, 3, 4, 5}. Tables 4(a) and (b) show the results of MOEA/DDE with the current/1 and rand/1 strategies, respectively. We do not present the comparison results of eachproblem here due to space constraints, but they can befound in Tables S.1 and S.2 in the supplementary fileof this paper.
Table 4(a) shows that the best performance of MOEA
/D-DE with the current/1 strategy is obtained by using the WPR selection method for M = 2. Although
MOEA/D-DE with the WOR and WR methods achieve
the same APS value for M = 4, the WR method is thebest selection method for M ∈ {3, 4, 5}.
Table 4(b) indicates that the WR selection method
is the most suitable for MOEA/D-DE using the rand/1
strategy for M ∈ {2, 3}. In contrast to the results forM ∈ {2, 3}, MOEA/D-DE with the WR method per-

10

Ryoji Tanabe, Hisao Ishibuchi

Table 6: Overall performance of MOEA/D-DE with the five bound-handling methods on the 16 problems. The APS value
at the final iteration is shown in the tables (lower is better). The numbers in parentheses indicate the ranks of the fiveconfigurations.
M

replacement

reinitializationreflection

r-reflection

resampling

2
3
4
5

0.69 (1)0.62 (1)0.62 (1)0.81 (1)2.12 (5)1.56 (3)1.69 (3)1.62 (3)1.38 (3)2.31 (5)2.12 (5)2.31 (5)1.19 (2)1.88 (4)1.75 (4)1.94 (4)1.69 (4)0.81 (2)1.06 (2)1.19 (2)forms poorly for M ∈ {4, 5}. Although MOEA/D-DE
with the WOR and WPR methods show the same APS
value for M = 4, the WOR method is the best selection
method for M = 5.
In summary, our results indicate that the choice ofindex selection method significantly affects the performance of MOEA/D-DE. As described above, the bestindex selection method depends on the type of mutation
strategies used in MOEA/D-DE and on the number ofobjectives M . In addition, the most suitable index selection method also depends on the type of problems(Tables S.1 and S.2). However, on most problems, theWR and WOR methods are likely to be suitable for thecurrent/1 and rand/1 mutation strategies, respectively.

5.2 The rand/1 and current/1 mutation strategies,
which is better for MOEA/D-DE?
Table 5 shows the overall performance of MOEA/DDE with the current/1 and rand/1 mutation strategies on the 16 problems with M ∈ {2, 3, 4, 5}. Basedon the results in Subsection 5.1, the WR and WOR
methods were used for the current/1 and rand/1 mutation strategies, respectively. Results on each problem
are presented in Table S.3 in the supplementary file.
Table 5 shows that the current/1 strategy is moresuitable for MOEA/D-DE than the rand/1 strategy
for all M . As discussed in [14], the mutant vector v igenerated by the current/1 strategy is generally close
to the target individual xi , and the neighborhood region of xi is efficiently exploited, which is the reason
why MOEA/D-DE with the current/1 strategy performs well. Although the rand/1 mutation strategy isincorporated into some packages, as reviewed in Subsection 3.1, Table 5 suggests that the current/1 strategy
is more appropriate for MOEA/D-DE.

problems. Results on each MOP can be found in Table
S.4 in the supplementary file.
For all M , the best results are obtained by using
the replacement method, followed by the resampling
method. For M ≥ 3, MOEA/D-DE with the reflection method achieves the worst APS values. Although
the r-reflection method is used in some implementations of MOEA/D-DE, as reviewed in Subsection 3.3,
its APS value is poor for M ≥ 3. In summary, our results show that the replacement method is the mostsuitable bound-handling method of the DE mutation
operator in MOEA/D-DE.
The results reported in [1] show that the resampling
method is the most appropriate for the basic DE [35] forsingle-objective continuous optimization. In contrast to
the results presented in [1], MOEA/D-DE with the resampling method performs the second best on MOPs.
Thus, the most suitable bound-handling method for theDE mutation operator depends on the problem domain.
5.4 Comparison of all 30 configurations
We separately examined the effect of the three components of the DE mutation operator in Subsections
5.1, 5.2, and 5.3. The results in the three subsections
show that the WR method, the current/1 strategy, andthe replacement method are suitable in MOEA/D-DE.
However, it is still unclear whether MOEA/D-DE witha combination of the three well-performing elements
works well. To determine whether it performs well, thissection presents comparisons of MOEA/D-DE with a
total of 30 configurations (the three index selection methods, the two mutation strategies, and the five boundhandling methods) of the DE mutation operator.
First, we show the results on the 16 test problems
in Subsection 5.4.1. Then, we report the performanceof the 30 configurations on each problem type in Subsection 5.4.2.

5.3 The influence of the bound-handling method on
the performance of MOEA/D-DE

5.4.1 Results on the 16 test problems

Table 6 shows the overall performance of MOEA/DDE with the five bound-handling methods on the 16

Table 7 shows the overall performance of MOEA/D-DE
with all 30 configurations on the 16 problems. Below,

Review and Analysis of Three Components of Differential Evolution Mutation Operator in MOEA/D-DE

11

Table 7: Overall performance of MOEA/D-DE with all 30 configurations on the 16 problems with M ∈ {2, 3, 4, 5}. The rankof each configuration based on its APS values at the final iteration is shown for each M . The average rank for all M is also
reported.
Boundhandling
methods

Mutationstrategies

Index
selection
methods
WOR

current/1

WR
WPR
WOR

replacement

Configuration
number #

M =2

M =3

#F

5
1

#B
#C
#D

M =4

M =5

Avg. rank

3

1

4

3.25

1

3

1

1.5

2
11

2
6

2
5

2
2

2.0
6.0

4

4

6

4

4.5

12
20
16

7
20
8

4
22
10

6
16
11

7.25
19.5
11.25

rand/1WR

current/1

WPR
WOR
WR

rand/1

WPR
WOR
WR

15
29
22

16
26
14

18
27
16

14
19
15

15.75
25.25
16.75

current/1

WPR
WOR
WR

30
16
8

27
21
12

28
22
13

22
23
20

26.75
20.5
13.25

rand/1

WPR
WOR
WR

7
28
19

18
29
22

21
29
24

21
30
24

16.75
29.0
22.25

25
18

30
17

30
14

29
28

28.5
19.25

3

11

11

17

10.5

6
23

15
22

19
26

26
25

16.5
24.0

reinitializationreflection

WPR
WOR
current/1

#E

WR
WPR
WOR

r-reflection

#A

#H
#G

rand/1

WR

13

19

20

18

17.5

current/1

WPR
WOR
WR

24
14
10

24
10
5

25
8
7

27
8
7

25.0
10.0
7.25

rand/1

WPR
WOR
WR

9
26
21

9
28
13

9
14
12

10
12
9

9.25
20.0
13.75

WPR

27

25

17

13

20.5

resampling(b, m, s)-MOEA/D-DE denotes an MOEA/D-DE using
b, m, and s, where b, m, and s represent a boundhandling method, a mutation strategy, and an index
selection method, respectively.
Table 7 indicates that the performance rank of MOEA
/D-DE significantly depends on the configuration of theDE mutation operator. For example, (replacement, current/1, WPR)-MOEA/D-DE works well for any number of objectives, and its average rank is 2.0. In contrast,
the worst configuration is (reflection, rand /1, WOR)MOEA/D-DE, whose average rank is 29.0. Although
this configuration is frequently incorporated into DE
for single-objective continuous optimization (e.g., [34,38]), our results show that it is unsuitable for MOEA
/D-DE. As shown in Table 7, the best configuration de-pends on the number of objectives M . According to theaverage rank for all M , the best performance is obtainedby (replacement, current/1, WR)-MOEA/D-DE. Thisconfiguration consists of the three well-performing components, as shown in Subsections 5.1 – 5.3.

5.4.2 Results on each problem type
Next, we discuss the performance of the 30 configurations on each problem type. According to Table 2, weclassified the 16 test problems into the following fourgroups:Unimodal problems DTLZ2, DTLZ4, DTLZ5, DTLZ6,
WFG1, WFG3, WFG5, WFG6, WFG7, and WFG8

12

Ryoji Tanabe, Hisao Ishibuchi

Table 8: Performance of MOEA/D-DE with all 30 configurations on the four problem sets (unimodal, multimodal, separable,
and nonseparable problems).The average rank values for all M ∈ {2, 3, 4, 5} are shown.
Boundhandling
methods

Mutationstrategies

Index
selection
methods
WOR

current/1

WR
WPR
WOR

replacement

Configuration
number #Unimodal

#F

3.25

7.0

3.5

9.5

1.25

4.25

1.0

12.5

#B
#C

2.25
5.5

6.5
10.0

2.25
5.75

9.5
13.25

3.5

8.0

3.75

17.0

#D

6.25
22.5
13.0

10.5
10.25
5.75

5.75
24.0
13.5

13.25
8.25
2.5

Multimodal

Separable

Nonseparable

rand/1WR

current/1

WPR
WOR
WR

rand/1

WPR
WOR
WR

16.5
27.25
18.75

8.0
16.5
9.25

20.25
27.5
18.75

3.0
18.5
11.75

current/1

WPR
WOR
WR

28.75
20.25
12.5

16.75
21.0
20.5

29.0
20.75
11.0

16.75
16.5
20.75

rand/1

WPR
WOR
WR

15.5
27.25
20.25

17.75
26.0
24.5

16.25
28.0
18.5

16.5
24.25
25.25

#E

28.0
18.0

26.75
21.75

27.5
18.5

25.25
17.5

10.0

16.5

8.5

15.75

#H
#G

15.5
22.75

17.5
23.0

14.75
23.75

17.25
21.75

reinitializationreflection

WPR
WOR
current/1

WR
WPR
WOR

r-reflection#A

rand/1

WR

15.0

23.0

14.0

20.5

current/1

WPR
WOR
WR

23.5
9.75
8.0

25.25
12.5
6.5

24.25
10.0
7.5

22.5
8.25
7.5

rand/1

WPR
WOR
WR

9.5
19.75
13.75

10.25
19.25
17.0

10.0
19.75
13.25

8.5
17.5
17.0

WPR

21.5

18.0

21.0

17.25

resampling

Multimodal problems DTLZ1, DTLZ3, DTLZ7, WFG2,
WFG4, and WFG9
Separable problems DTLZ1, DTLZ2, DTLZ3, DTLZ4,
DTLZ5, DTLZ6, DTLZ7, WFG1, WFG4, WFG5,and WFG7
Nonseparable problems WFG2, WFG3, WFG6, WFG8,and WFG9Table 8 shows results of the 30 configurations on thefour problem types. Each value in Table 8 is the averagerank value of each configuration over all problems (M ∈{2, 3, 4, 5}) in each problem type. On the one hand, Table 8 indicates that the performance rank of each configuration for the unimodal, multimodal, and separable
problem sets is almost consistent with that for all 16
test problems in Table 7. The best average rank value isobtained by (replacement, current/1, WR)-MOEA/D-

DE. In addition, (replacement, current/1, WPR)-MOEA/DDE performs well. On the other hand, the results forthe nonseparable problem set show that the reinitialization repair method is suitable for MOEA/D-DE. (reinitialization, current/1, WR)-MOEA/D-DE and (reinitialization, current/1, WPR)-MOEA/D-DE have good
performance according to the average rank values. This
result may be because MOEA/D-DE needs to keep thediversity of the population to find good solutions on
nonseparable problems. Since the reinitialization method randomly generates an element of a child that isout of the bound, it can introduce diversity in the population. Although the best configuration also depends
on the problem type, (replacement, current/1, WR)MOEA/D-DE works well on a wide variety of problems.

Review and Analysis of Three Components of Differential Evolution Mutation Operator in MOEA/D-DE

5.5 Comparison between the eight existing
configurations
As reviewed in Section 3, the configuration of the DE
mutation operator differs depending on the source code.
According to Table 1, eight configurations (#A, ..., #H)
exist. Our results in Subsections 5.1–5.4 also show thatthe performance of MOEA/D-DE significantly depends
on components of the mutation. Here, we discuss theperformance rank of MOEA/D-DE with the eight existing configurations.
Table 7 shows that configuration #B has the best
average rank value among the eight configurations. Configuration #F is the second-best performer. Although
the average rank value of configuration #A is not goodfor all 16 problems, it is best for the nonseparable problem set, as discussed in Subsection 5.4.2 (see Table 8).
In contrast, configurations #E, #G, and #H show poor
performance. In particular, configuration #G performs
the worst among the eight configurations (#A, ..., #H).
The above results indicate that an incorrect conclusion could be obtained depending on the source codeof MOEA/D-DE used. For example, let us consider theexperiment to reproduce the results of MOEA/D-DE in
the original paper [25]. As shown in Table 1, configuration #B is used in the original C++ code of MOEA/DDE. However, configuration #G is incorporated into
the PaGMO library. If one uses the PaGMO implementation of MOEA/D-DE (#G) rather than its originalC++ source code (#B), she/he is likely to obtain results that are worse than those reported in [25]. Thatis, the performance of MOEA/D-DE is underestimatedin this case. To avoid such an undesirable comparison,
the configuration of the DE mutation operator should
receive careful attention.

6 Conclusion
We have reviewed the existing configurations of the DE
mutation in MOEA/D-DE (Section 3) and examined
the influence of each component on the performance of
MOEA/D-DE (Section 5). The main contributions of
this paper can be summarized as follows:– The DE mutation consists of three components (themutation strategy, the index selection method, andthe bound-handling method). Our review reveals
that the configuration of the DE mutation operator in MOEA/D-DE differs depending on the sourcecode. We also explained why different configurationsof the DE mutation operator are implemented in thesource code of MOEA/D-DE.

13– Although implementations of the DE mutation operator in MOEA/D-DE have not received much attention, our results shows that the performance of
MOEA/D-DE significantly depends on the configuration of the DE mutation operator. Thus, it is necessary to carefully select the configuration of the DE
mutation operator to maximize the performance of
MOEA/D-DE. Our results show that the combination of the current/1 strategy, the WR method (orthe WPR method), and the replacement method ismost suitable for MOEA/D-DE.
As we discussed in Section 5, some of our results on
MOPs are inconsistent with previous studies on singleobjective optimization problems. For example, while a
combination of the WOR method, the rand/1 strategy,
and the reflection method is frequently used in DE algorithms for single-objective optimization, our results
indicates that it is inappropriate for MOEA/D-DE. As
far as we know, such observations have not been reported in the literature. One future research topic is toinvestigate why the best configuration of the DE mutation operator depends on the problem domain.
Another future research topic is to analyze the impact of configurations of the DE mutation operator
on the performance of other DE-based MOEAs (e.g.,
DEMO [33] and GDE3[23]). We used the DTLZ and
WFG problems to examine the performance of MOEA/DDE. Comparison of various configurations of the DE
mutation in MOEA/D-DE on other test problems (e.g.,
the CEC2017 problems [7]) remains as future work.
Acknowledgement
This work was supported by the Program for Guangdong Introducing Innovative and Entrepreneurial Teams(Grant No. 2017ZT07X386), Shenzhen Peacock Plan
(Grant No. KQTD2016112514355531), the Science andTechnology Innovation Committee Foundation of Shenzhen (Grant No. ZDSYS201703031748284), the Program for University Key Laboratory of Guangdong Province(Grant No. 2017KSYS008), and National Natural Science Foundation of China (Grant No. 61876075).
This is a pre-print of an article published in Soft
Computing. The final authenticated version is available
online at: https://doi.org/10.1007/s00500-019-03842-6.
Compliance with ethical standards
Conflict of interest
Ryoji Tanabe and Hisao Ishibuchi declare that they
have no conflict of interest.

14

Ethical approval
This article does not contain any studies with humanparticipants performed by any of the authors.
References
1. J. Arabas, A. Szczepankiewicz, and T. Wroniak. Experimental comparison of methods to handle boundary constraints in differential evolution. In PPSN, pages 411–
420, 2010.
2. A. Auger, D. Brockhoff, N. Hansen, D. Tusar, T. Tusar,
and T. Wagner. The Impact of Variation Operatorson the Performance of SMS-EMOA on the Bi-objective
BBOB-2016 Test Suite. In GECCO, pages 1225–1232,
2016.
3. J. Bader and E. Zitzler. HypE: An Algorithm for
Fast Hypervolume-Based Many-Objective Optimization.
Evol. Comput., 19(1):45–76, 2011.
4. N. Beume, B. Naujoks, and M. T. M. Emmerich. SMSEMOA: multiobjective selection based on dominated hypervolume. EJOR, 181(3):1653–1669, 2007.
5. L. C. T. Bezerra, M. López-Ibáñez, and T. Stützle. To
DE or Not to DE? Multi-objective Differential Evolution
Revisited from a Component-Wise Perspective. In EMO,
pages 48–63, 2015.
6. D. Brockhoff, T. Tran, and N. Hansen. Benchmarking Numerical Multiobjective Optimizers Revisited. In
GECCO, pages 639–646, 2015.
7. R. Cheng, M. Li, Y. Tian, X. Zhang, S. Yang, Y. Jin,and X. Yao. Benchmark Functions for the CEC’2017
Competition on Many-Objective Optimization. Technicalreport, Univ. of Birmingham, 2017.
8. S. Das, S. S. Mullick, and P. N. Suganthan. Recent
advances in differential evolution - An updated survey.
Swarm and Evol. Comput., 27:1–30, 2016.
9. S. Das and P. N. Suganthan. Differential Evolution: A
Survey of the State-of-the-Art. IEEE TEVC, 15(1):4–31,
2011.
10. K. Deb. Multi-Objective Optimization Using Evolutionary Algorithms. John Wiley & Sons, 2001.
11. K. Deb, S. Agrawal, A. Pratap, and T. Meyarivan. A
fast and elitist multiobjective genetic algorithm: NSGAII. IEEE TEVC, 6(2):182–197, 2002.
12. 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, 18(4):577–601, 2014.
13. K. Deb, L. Thiele, M. Laumanns, and E. Zitzler. Scalable Test Problems for Evolutionary Multi-Objective Optimization. In Evolutionary Multiobjective Optimization.
Theoretical Advances and Applications, pages 105–145.
Springer, 2005.
14. W. Gong, Y. Wang, Z. Cai, and S. Yang. A Weighted
Biobjective Transformation Technique for Locating Multiple Optimal Solutions of Nonlinear Equation Systems.
IEEE TEVC, 21(5):697–713, 2017.
15. N. Hansen, A. S. P. Niederberger, L. Guzzella, andP. Koumoutsakos. A Method for Handling Uncertainty in
Evolutionary Optimization With an Application to Feedback Control of Combustion. IEEE TEVC, 13(1):180–
197, 2009.
16. S. Helwig, J. Branke, and S. Mostaghim. Experimental
analysis of bound handling techniques in particle swarm
optimization. IEEE TEVC, 17(2):259–271, 2013.

Ryoji Tanabe, Hisao Ishibuchi
17. S. Huband, P. Hingston, L. Barone, and R. L. While. A
review of multiobjective test problems and a scalable test
problem toolkit. IEEE TEVC, 10(5):477–506, 2006.
18. H. Ishibuchi, N. Akedo, and Y. Nojima. Behavior of Multiobjective Evolutionary Algorithms on Many-Objective
Knapsack Problems. IEEE TEVC, 19(2):264–283, 2015.
19. H. Ishibuchi, H. Masuda, and Y. Nojima. Pareto Frontsof Many-Objective Degenerate Test Problems. IEEE
TEVC, 20(5):807–813, 2016.
20. H. Ishibuchi, Y. Setoguchi, H. Masuda, and Y. Nojima.
Performance of Decomposition-Based Many-Objective
Algorithms Strongly Depends on Pareto Front Shapes.
IEEE TEVC, 21(2):169–190, 2017.
21. S. Jiang and S. Yang. An Improved Multiobjective Optimization Evolutionary Algorithm Based on Decomposition for Complex Pareto Fronts. IEEE Trans. Cyber.,46(2):421–437, 2016.
22. J. B. Kollat, P. M. Reed, and R. M. Maxwell. Manyobjective groundwater monitoring network design using
bias-aware ensemble Kalman filtering, evolutionary optimization, and visual analytics. Water Resources Research, 47(2):1–18, 2011.
23. S. Kukkonen and J. Lampinen. GDE3: the third evolution step of generalized differential evolution. In IEEE
CEC, pages 443–450, 2005.
24. M. N. Le, Y. Ong, S. Menzel, Y. Jin, and B. Sendhoff. Evolution by adapting surrogates. Evol. Comput.,21(2):313–340, 2013.
25. H. Li and Q. Zhang. Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and
NSGA-II. IEEE TEVC, 13(2):284–302, 2009.
26. K. Li, Á. Fialho, S. Kwong, and Q. Zhang. Adaptive
Operator Selection With Bandits for a Multiobjective
Evolutionary Algorithm Based on Decomposition. IEEE
TEVC, 18(1):114–130, 2014.
27. K. Li, Q. Zhang, S. Kwong, M. Li, and R. Wang. Stable Matching-Based Selection in Evolutionary Multiobjective Optimization. IEEE TEVC, 18(6):909–923, 2014.
28. X. Liao, Q. Li, X. Yang, W. Zhang, and W. Li. Multiobjective optimization for crash safety design of vehicles using stepwise regression model. Struct. Multidiscipl.
Optim., 35(6):561–569, 2008.
29. E. Mezura-Montes, J. Velázquez-Reyes, and C. A.
Coello Coello. A comparative study of differential evolution variants for global optimization. In GECCO, pages
485–492, 2006.
30. Y. S. Ong, P. B. Nair, and A. J. Keane. Evolutionary
optimization of computationally expensive problems via
surrogate modeling. AIAA journal, 41(4):687–696, 2003.
31. M. Pescador-Rojas, R. Hernández Gómez, El. Montero,N. Rojas-Morales, M. Cristina Riff, and C. A. CoelloCoello. An Overview of Weighted and Unconstrained
Scalarizing Functions. In EMO, pages 499–513, 2017.
32. A. K. Qin, V. L. Huang, and P. N. Suganthan. Differential evolution algorithm with strategy adaptation for
global numerical optimization. IEEE TEVC, 13(2):398–
417, 2009.
33. T. Robič and B. Filipič. DEMO: differential evolution
for multiobjective optimization. In EMO, pages 520–533,
2005.
34. J. Rönkkönen, S. Kukkonen, and K. V. Price. RealParameter optimization with Differential Evolution. In
IEEE CEC, pages 506–513, 2005.
35. R. Storn and K. Price. Differential Evolution - A Simple and Efficient Heuristic for Global Optimization overContinuous Spaces. J. Global Optimiz., 11(4):341–359,
1997.

Review and Analysis of Three Components of Differential Evolution Mutation Operator in MOEA/D-DE
36. A. Trivedi, D. Srinivasan, K. Sanyal, and A. Ghosh. A
Survey of Multiobjective Evolutionary Algorithms Basedon Decomposition. IEEE TEVC, 21(3):440–462, 2017.
37. T. Tusar and B. Filipic. Differential Evolution versus
Genetic Algorithms in Multiobjective Optimization. In
EMO, pages 257–271, 2007.
38. Y. Wang, Z. Cai, and Q. Zhang. Differential Evolution
With Composite Trial Vector Generation Strategies andControl Parameters. IEEE TEVC, 15(1):55–66, 2011.
39. Y. Wang, Z. Liu, J. Li, H. Li, and J. Wang. On the selection of solutions for mutation in differential evolution.
Frontiers of Computer Science, pages 1–19, 2016.
40. Z. Wang, Q. Zhang, A. Zhou, M. Gong, and L. Jiao.
Adaptive Replacement Strategies for MOEA/D. IEEE
Trans. Cyber., 46(2):474–486, 2016.
41. S. Wessing. Repair methods for box constraints revisited.
In EvoApplications, pages 469–478, 2013.
42. Y. Yuan, H. Xu, and B. Wang. An Experimental Investigation of Variation Operators in Reference-Point Based
Many-Objective Optimization. In GECCO, pages 775–
782, 2015.
43. Q. Zhang and H. Li. MOEA/D: A multiobjective evolutionary algorithm based on decomposition. IEEE TEVC,
11(6):712–731, 2007.
44. Q. Zhang, W. Liu, and H. Li. The performance of a newversion of MOEA/D on CEC09 unconstrained MOP test
instances. In IEEE CEC, pages 203–208, 2009.
45. E. Zitzler and S. Künzli. Indicator-based selection in
multiobjective search. In Proc. of PPSN, pages 832–842,
2004.
46. E. Zitzler, M. Laumanns, and L. Thiele. SPEA2: Improving the Strength Pareto Evolutionary Algorithm. Technical report, ETHZ, 2001.
47. E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, andV. Grunert da Fonseca. Performance assessment of multiobjective optimizers: an analysis and review. IEEE
TEVC, 7(2):117–132, 2003.

15