An Analysis of Control Parameters of MOEA/D Under
Two Different Optimization Scenarios
Ryoji Tanabe, Hisao Ishibuchi∗
Shenzhen Key Laboratory of Computational Intelligence, Department of Computer Scienceand Engineering, Southern University of Science and Technology, Shenzhen, 518055, China

Abstract
An unbounded external archive (UEA), which stores all nondominated solutions
found during the search process, is frequently used to evaluate the performanceof multi-objective evolutionary algorithms (MOEAs) in recent studies. A recent
benchmarking study also shows that decomposition-based MOEA (MOEA/D)
is competitive with state-of-the-art MOEAs when the UEA is incorporated into
MOEA/D. However, a parameter study of MOEA/D using the UEA has not
yet been performed. Thus, it is unclear how control parameter settings influence the performance of MOEA/D with the UEA. In this paper, we present ananalysis of control parameters of MOEA/D under two performance evaluation
scenarios. One is a final population scenario where the performance assessmentof MOEAs is performed based on all nondominated solutions in the final population, and the other is a reduced UEA scenario where it is based on a pre-specified
number of selected nondominated solutions from the UEA. Control parametersof MOEA/D investigated in this paper include the population size, scalarizing
functions, and the penalty parameter of the penalty-based boundary intersection (PBI) function. Experimental results indicate that suitable settings of thethree control parameters significantly depend on the choice of an optimization
scenario. We also analyze the reason why the best parameter setting is totallydifferent for each scenario.
Keywords: multi-objective optimization, decomposition based evolutionary
algorithms, parameter study, external archive∗ Corresponding author

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

Preprint submitted to Applied Soft Computing

April 23, 2019

1. Introduction
An unconstrained (bound-constrained) multi-objective optimization problem (MOP) can be formulated as follows:
Tminimize f (x) = f1 (x), ..., fM (x)(1)

subject to x ∈ S ⊆ RD ,where f : S → RM is an objective function vector that consists of M potentiallyconflicting objective functions, and RM is the objective function space. Here,minmax
x = (x1 , ..., xD )T is a D-dimensional solution vector, and S = ΠD]
j=1 [xj , xj
minmax
is the bound-constrained search space where xj ≤ xj ≤ xj
for each indexj ∈ {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 at least one index i. Here, x∗ is a Paretooptimal solution if there exists no x ∈ S such that x dominates x∗ . In this case,f (x∗ ) is a Pareto-optimal objective function vector. The set of all x∗ in S is thePareto-optimal solution set (PS), and the set of all f (x∗ ) is the Pareto front(PF). Usually, no solution can simultaneously minimize all objective functions
f1 , ..., fM in MOPs. Thus, the goal of MOPs is to find a set of nondominatedsolutions that are well-distributed and close to the PF in the objective function
space. MOPs frequently appear in engineering problems such as aerodynamic
wing design problems [1], financial and economic problems [2], oil well problems[3], and unit commitment problems [4]. In general, it is difficult to find a set ofgood nondominated solutions on MOPs with large values of M and/or D. Thisis because the objective function and solution spaces exponentially grow with
M and D, respectively.
A multi-objective evolutionary algorithm (MOEA) is an efficient approachfor solving MOPs [5]. Since MOEAs use a set of individuals (solutions of a given
MOP) for the search, it is expected that well-distributed nondominated solutions
can be found in a single run. A number of MOEAs have been proposed in theevolutionary computation community [6, 7]. MOEAs can be roughly classified
into the following three categories: dominance-based MOEAs, indicator-based
MOEAs, and decomposition-based MOEAs. Dominance-based MOEAs (e.g.,
NSGA-II [8], SPEA2 [9], and -MOEA [10]) mainly use the Pareto dominanceor relaxed dominance relations for the mating and environmental selections.
An indicator-based MOEA assigns a so-called fitness value to each individualin the population using a quality indicator [11]. Representative indicator-based
MOEAs include IBEA [11], SMS-EMOA [12], and HypE [13]. A decompositionbased MOEA decomposes a given MOP with M objectives into multiple singleobjective sub-problems1 using a scalarizing function g : RM → R and tries tofind good solutions for all the subproblems. Well-known decomposition based
1 Some decomposition based MOEAs (e.g., MOEA/D-M2M [14]) decompose a given MOP
into a set of simple MOPs.

2

MOEAs are MOGLS [15], C-MOGA [16], MSOPS [17], and MOEA based on decomposition (MOEA/D) [18]. In particular, recent studies report the promising
performance of MOEA/D-type algorithms [19].
MOEA/D [18] decomposes an M -objective MOP defined in equation (1) into
µ single-objective sub-problems g1 (x|w1 ), ..., gµ (x|wµ ) using a set of weight veci T
tors W = {w1 , ..., wµ } and a scalarizing function g, where wi = (w1i , ..., wM)PMii(i ∈ {1, ..., µ}) and j=1 wj = 1. MOEA/D assigns each individual x (i ∈{1, ..., µ}) to each sub-problem and tries to find the optimal solution of all subproblems simultaneously. Unlike other MOEAs, in MOEA/D, the mating andenvironmental selections are performed only in a set of neighborhood individualsof each weight vector wi . Recent studies show that improved MOEA/D-type
algorithms are capable of finding good nondominated solutions on MOPs withcomplex PFs [20, 21, 22, 23]. Although MOEA/D was originally designed for
MOPs with up to four objectives, its variants can efficiently handle a large
number of objectives [24, 25]. Also, MOEA/D-type algorithms are successfully
applied to real-world problems [19].
In general, the performance of evolutionary algorithms significantly depends
on control parameter settings [26, 27]. MOEAs including MOEA/D are not an
exception. Therefore, it is important to understand how each control parameter
influences the performance of MOEA/D. General rules of thumb are helpful tousers for tuning control parameters of MOEA/D (e.g., the population size µshould be set to approximately 28 on three-objective multimodal MOPs). For
these reasons, some parameter studies have been performed for MOEA/D as
briefly reviewed below:
Population size µ: The population size µ is an important parameter for all
MOEAs. However, there exist only a few studies that investigate the impact of µ on the performance of MOEA/D. Experimental results in theoriginal MOEA/D paper [18] show that MOEA/D with a small µ value(µ = 20) can successfully find well-approximated nondominated solutions
close to the PF. The performance of NSGA-II and MOEA/D with variousµ values is investigated in [28]. Results on multi-objective knapsack problems show that MOEA/D with a very large µ value works well on most of
the problems, while such µ values make the convergence speed of NSGAII slow. In [29], the performance of two reference-vector based MOEAs(NSGA-III [30] and θ-DEA [31]) and two MOEA/D variants (MOEA/Dand MOEA/DD [24]) with three µ values are evaluated on three- and fiveobjective MOPs. Experimental results in [29] show that the best µ value
is different for each MOEA.
Scalarizing function g: In MOEA/D, a fitness value of an individual on each
sub-problem j ∈ {1, ..., µ} is given by a pre-defined scalarizing function
g. Therefore, the performance of MOEA/D is affected by the scalarizing function g used for the search [18, 32, 33]. Typical scalarizing functions for MOEA/D include the weighted sum, Chebyshev, and Penaltybased Boundary Intersection (PBI) functions [18]. It is pointed out that
MOEA/D with the weighted sum function does not have an ability to han3dle MOPs with non-convex PFs [18]. In [33], the behavior of MOEA/Dwith the three representative scalarizing functions is investigated on multiobjective knapsack problems with up to 10 objectives. Some studies (e.g.,[18, 30]) report that MOEA/D with the PBI function can find evenly
distributed nondominated solutions compared to the Chebyshev function.
Since an appropriate scalarizing function is different problem dependent,
adaptive selection strategies have also been proposed for these [32, 34].
Penalty parameter θ: The efficiency of the PBI function is significantly influenced by the setting of the penalty parameter θ. The impact of θ values
on the performance of MOEA/D is examined in [35, 36]. A deterministic
control method for θ is also proposed in [22].
It should be noted that the above-described parameter studies (except for[29]) are based on only the final population scenario [37], where all nondominatedsolutions in the final population are used for the performance assessment. The
final population scenario is the most widely used optimization scenario andadopted in almost all previous studies of MOEA/D (e.g., [21, 22, 23, 24, 30,31]). While the final population scenario is commonly used in the evolutionary
computation community, it is not always a practically desirable optimization
scenario [38]. In the final population scenario, if the number of nondominatedsolutions found in the search exceeds the predefined population size µ, they
are removed from the population to keep the population size constant. This
operation is undesirable if a user wants to know the entire PF using a large
number of nondominated solutions. A good potential solution found in the
search is also possibly to be discarded from the population [38]. Such problems
are easily addressed by using an unbounded external archive (UEA) [39, 40]
that stores all nondominated solutions found during the search process. For suchreasons, the UEA is frequently used in recent work (e.g., [37, 38, 41, 42, 43, 44]).
The recent COCO platform2 with the BBOB-biobj functions [45] also adopts the
UEA for the performance assessment of multi-objective optimization methods.
One may think that decision makers usually want to know only a small number of representative, well-distributed nondominated solutions, and thus the useof a large number of nondominated solutions in the UEA for performance assessment does not make sense. Fortunately, such an issue can be easily addressedby applying a selection method of a small number of nondominated solutions(e.g., [37, 38, 46, 47, 48, 49, 50]) to the UEA. Therefore, there is no particularreason not to incorporate the UEA into MOEAs especially if the computational
cost for archive maintenance is sufficiently small in comparison with the objective function evaluation of each solution as in many real-world application
problems3 . According to [37], an optimization scenario that uses a pre-specified
number of selected nondominated solutions from the UEA for the performanceassessment is called the reduced UEA scenario in this paper. A benchmarking
2 http://coco.gforge.inria.fr/doku.php?id=algorithms-biobj
3 A simulation run that takes a long time is required for some real-world problems toevaluate a single solution [51, 52].

4

study in [37] shows that MOEA/D and its variants (MOEA/D-DE [20] andMOEA/D-DRA [49]) are competitive with state-of-the-art methods under the
reduced UEA scenario.
This paper presents an analysis of control parameters of MOEA/D on MOPs
with up to five objectives under the final population and reduced UEA scenarios. Since MOEA/D performs well under the reduced UEA scenario as reported
in [37], such an analysis is worth performing and helpful to both its users andalgorithm designers who try to develop more efficient MOEA/D. We investigate the following three control parameters of MOEA/D: the population size µ,scalarizing functions g, and the penalty parameter θ of the PBI function. As far
as we know, parameter studies of MOEA/D for the reduced UEA scenario havenot been well performed. As mentioned above, most of the previous analyticalstudies are based only on results for the traditional final population scenario.
Although the reduced UEA scenario is frequently used in recent studies, it is unclear whether MOEA/D with a suitable parameter setting for the final population scenario works similarly under the reduced UEA scenario. Some MOEA/D
algorithms with an external archive have also been proposed (e.g., [18, 53, 54]),
but parameter studies are not performed in the literature. Whereas parametertuning studies (not analytical studies) of MOEAs for the final population scenario (e.g., [55, 56]) and the UEA scenario4 (e.g., [43, 57]) have been presented
individually, they do not investigate how different the tuned parameter values
are between the two different optimization scenarios. The main contributionsof this paper can be summarized as follows:• We carefully examine proper settings of the three control parameters (µ, g,and θ) which realize the best performance of MOEA/D on the four DTLZ
[58] and nine WFG [59] test problems in a component-wise manner.
• In addition to the one-by-one analysis, we investigate dependencies between two control parameters of MOEA/D.
• By comparing experimental results between the final population and reduced UEA scenarios, we also reveal how different appropriate settings ofthe three control parameters of MOEA/D are for each scenario.
• Furthermore, the reason why the best parameter setting on some MOPs
is totally different for each scenario is analyzed in this paper.
This paper is organized as follows: Section 2 describes the basic procedureof MOEA/D. Experimental settings are introduced in Section 3. Section 4shows the analysis of the three control parameters (µ, g, and θ) of MOEA/D,
separately. Section 5 presents further analysis of the three control parameters.
Finally, Section 6 concludes this paper and discusses the future work.
4 The UEA scenario was named in [37]. All nondominated solutions in the UEA are used
for the performance assessment under the UEA scenario.

5

Algorithm 1: The procedure of MOEA/Dt ← 1, initialize the population P = {x1 , ..., xµ };2 for i ∈ {1, ..., µ} do
3
Set the neighborhood index list B i = {i1 , ..., iT };1

while The termination criteria are not met dofor i ∈ {1, ..., µ} do
6Randomly select two indices k and l from B i ;7
Generate a child ui by crossing xk and xl ;8
Apply a mutation operator to ui ;9
for j ∈ B i do
10
if g(ui |wj , z ∗ ) ≤ g(xj |wj , z ∗ ) then11
xj ← ui ;4512t ← t + 1;

2. MOEA/D
Here, MOEA/D [18] is briefly explained. Algorithm 1 shows the overall
procedure of MOEA/D. MOEA/D decomposes an MOP with M objectives
into µ single-objective sub-problems using a set of uniformly distributed weightvectors W = {w1 , ..., wµ }. In our study, Das and Dennis’s systematic approach[60] was used to generate the weight vectors W .
At the beginning of the search, all individuals in the population are randomly
generated in the search space S (line 1). For each subproblem index i ∈ {1, ..., µ},an index list B i = {i1 , ..., iT }, which is used for the mating and replacement
selections, is initialized: B i consists of indices of the T closest weight vectors to
wi in the weight vector space (lines 2-3) where T is the neighborhood size.
After the initialization, the following steps are repeatedly applied for each
subproblem i ∈ {1, .., µ} until the search termination criteria are met. For eachi, the parent indices k and l are randomly selected from B i (line 6). Then,a child ui is generated by crossing xk and xl (line 7). A mutation operator
is applied to the child ui if necessary (line 8). After ui has been generated,
the replacement selection is performed using a predefined scalarizing function g(lines 9–11). For each j ∈ B i , the individual xj is compared with the child ui
based on g. If ui is better than xj according to their scalarizing function values
based on the weight vector wj , xj is replaced by ui (lines 10–11).
Since the replacement of individuals is based on their scalarizing function
values, g plays a crucial role in MOEA/D. Although there are a number of scalarizing functions as reviewed in [61], in this paper we investigated the following
three scalarizing functions: the two Chebyshev functions (a multiplication version g chm [18] and a division version g chd [62, 63, 25]) and the PBI function g pbi[18]. Since the three scalarizing functions (g chm , g chd , and g pbi ) are most widelyused ones in the literature, they are worth investigating. The three scalarizing

6

functions are defined as follows:g chm (x|w, z ∗ ) ={wi |fi (x) − zi∗ |},|fi (x) − zi∗ |chd∗g (x|w, z ) = max,wii∈{1,...,M }maxi∈{1,...,M }g pbi (x|w, z ∗ ) = d1 + θ d2 ,

(2)
(3)
(4)

∗ T

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

(5)
(6)

where all the three scalarizing functions in equations (2), (3), and (4) should∗ T) is the ideal point. Since it is difficult tobe minimized. The z ∗ = (z1∗ , ..., zM∗
obtain the actual ideal point z of a given MOP, its approximated point thatconsists of the minimum function value for each objective fi (i ∈ {1, ..., M })found during the search process is usually used for the calculation of equations(2) – (6). For g chd in equation (3), if wi = 0, it was set to 10−6 for its implementation to avoid division by zero. While the Chebyshev function g chm definedin equation (2) is one of the most frequently used scalarizing functions, somestudies (e.g., [62, 25]) use its alternative version g chd . This is because the search
directions of MOEA/D can be more evenly distributed by using g chd [63].
In equations (5) and (6), d1 denotes how close the objective function vector
f (x) is to the PF, and d2 is the perpendicular distance between f (x) and w.
The two distance measures d1 and d2 evaluate the convergence and diversity ofthe solution x in the objective function space, respectively. The PBI function
value calculated by equation (4) is the sum of d1 and θ d2 . The penalty parameter θ > 0 controls the balance between the convergence (d1 ) and diversity (d2 )in the population. While a small θ encourages convergence toward the PF, alarge θ value emphasizes the importance of diversity in the population [33].
3. Experimental settings
This section introduces our experimental settings. Experimental results arereported in Section 4 and 5. Subsection 3.1 describes test problems and aperformance indicator used in our study. In our study, we used the average
performance score (APS) [13] to aggregate the performance of MOEA/D withvarious configurations on 13 MOPs. The calculation method of the APS value
is described in Subsection 3.2. Subsection 3.3 introduces the parameter settingsof MOEA/D.
3.1. Test problems and performance indicator
The four DTLZ [58] and nine WFG [59] test problems with M ∈ {2, 3, 4, 5}
were used in our analysis study. Table 1 summarizes their properties. The

7

Table 1: Properties of the DTLZ and WFG test problems.
Problem

Shape of PF

DTLZ1LinearDTLZ2Nonconvex

DTLZ3

Nonconvex

DTLZ4

Nonconvex

WFG1

Mixed

WFG2

Discontinuous

WFG3

Partially
Degenerate

WFG4

Nonconvex

WFG5

Nonconvex

WFG6

Nonconvex

WFG7

Nonconvex

WFG8

Nonconvex

WFG9

Nonconvex

Multimodality

Nonseparability

OthersX

X
BiasedBiasedX

X
X

X
Deceptive
X
BiasedX

X

BiasedX

Deceptive, Biased

shapes of the PFs of the WFG1, WFG2, and WFG3 test problems are complicated, discontinuous, and partially degenerate [64], respectively. The DTLZ1
problem has a linear PF. The shapes of the PFs of the other problems arenonconvex PFs. According to [58], for the DTLZ problems, the number of theposition variables k was set to k = 5 for the DTLZ1 problem and k = 10 for theother DTLZ problems, where the number of variables is D = M + k − 1. Also,
as suggested in [59], for the WFG test problems, the number of the position
variables k was set to k = 2 (M − 1), and the number of the distance variables
l was set to l = 20, where D = k + l.
The hypervolume (HV) indicator [65] was used for evaluating the quality ofa set of obtained nondominated solutions A. Before calculating the HV value,
the objective function vector f (x) of each x ∈ A was normalized using the ideal
point and the nadir point of the target MOP. As suggested in [66], the reference
point for the HV calculation was set to (1.1, ..., 1.1)T . In this setting, the HV
range for all of the WFG test problems is [0, 1.1M ]. We further normalized HVvalues ∈ [0, 1.1M ] to the range [0, 1] by dividing by 1.1M . The HV value of A
was calculated for every 2 000 function evaluations.
As mentioned in Section 1, all nondominated solutions in the population P
are used for the HV calculation under the final population scenario. For the
reduced UEA scenario, a selection method of a small number of nondominatedsolutions from the UEA is necessary. Although there are some computationallycheap selection methods (e.g., [38, 49, 50]), we used a distance-based selection
method described in [37]. In this selection method, a pre-defined number ofsolutions b is selected from the UEA. First, M extreme solutions having theminimum objective function values for fi (i ∈ {1, ..., M }) are selected for B.
8

0.4
f 0.6

0.8

0.0 0.2

0.4
f 0.6
1

0.8

1.0

1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8 f 21.0

0.0 0.2

0.4
f 0.6
1

0.8

0.0 0.2

0.4
f 0.6
1

0.8

(i) Final Population scenario

f3

1

1.0

1.0

1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8 f 2
1.0

1.0

1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8 f 2
1.0

f3

0.0 0.2

1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8 f 21.0

0.0 0.2

0.4
f 0.6
1

0.8

f3

1.0

1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8 f 2
1.0

1.0

1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8 f 2
1.0

f3

1.0

1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8 f 2
1.0

FEvals = 5 × 104
f3

FEvals = 1 × 104
f3

FEvals = 2 × 103

0.0 0.2

0.4
f 0.6
1

0.8

0.0 0.2

0.4
f 0.6
1

0.8

1.0

1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8 f 2
1.0

f3

f3

1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8 f 21.0

0.0 0.2

0.4
f 0.6
1

0.8

1.0

f3

(ii) UEA scenario

0.0 0.2

0.4
f 0.6
1

0.8

(iii) Reduced UEA scenario
Figure 1: Distribution of nondominated solutions found by MOEA/D with µ = 210 and g chm
on the three-objective WFG4 problem under (i) the final population scenario, (ii) the UEA
scenario, and (iii) the reduced UEA scenario. Results at 2 × 103 (left), 1 × 104 (center), and5×104 (right) function evaluations (FEvals) are shown. In the UEA scenario, all nondominatedsolutions in the UEA are used for the performance assessment [37]. The number of selected
nondominated solutions b for the reduced UEA scenario was set to 210 for M = 3. Results ofa single run are shown.

After that, a nondominated solution which is farthest from a set of alreadyselected ones in the objective function space is repeatedly added to B until
the size of B is equal to b. It is expected that a set of uniformly distributed
nondominated solutions in the objective function space are obtained by the
distance-based selection method. In our study, the number of selected nondominated solutions b was set to 200, 210, 220, and 210 for M = 2, 3, 4, and 5,
respectively.
For the sake of explanation, we show the distribution of nondominated solutions found by MOEA/D with µ = 210 and g chm on the three-objective WFG4
problem under (i) the final population scenario, (ii) the UEA scenario, and (iii)the reduced UEA scenario in Figure 1. Results at 2 × 103 , 1 × 104 , and 5 × 104
function evaluations are shown. As shown in Figure 1, the distribution of nondominated solutions in the population and the UEA is significantly different
after 1 × 104 function evaluations. In the final population scenario, nondominated solutions to be preserved in the population in the next iteration aredetermined according to the environmental selection of MOEA/D (lines 9–11 in
9Table 2: Settings for the three control parameters of MOEA/D (the population size µ, thescalarizing function g, and the penalty value θ of the PBI function g pbi ). Unless explicitly
noted, the default settings were used in our experimental study. Examined settings in thetable denote parameter values which are investigated in our analytical study, and each result
can be found in its corresponding section.

Controlparameters

Default
settings

Examinedsettings

Corresponding
sectionsµ for M = 2{
µ for M = 3{
µ for M = 4{
µ for M = 5{

200{
210{
220{
210{

{25, 50, 100, 200, 300, 400}
{28, 55, 105, 210, 300, 406}
{35, 56, 120, 220, 286, 455}
{126, 210, 330, 495}

Section 4.1

g

g chm

{g chm , g chd , g pbi }Section 4.2

θ5{0.1, 0.5, 1, 2, 3, 5, 8, 10}Section 4.3

Algorithm 1). In contrast, in the UEA scenario, if a newly generated solution is
nondominated with respect to existing solutions in the UEA, it enters the UEA.
Then, solutions dominated by the newly inserted solution are discarded fromthe UEA. As seen from Figure 1(ii), the number of nondominated solutions in
the UEA gradually increases as the search progresses. At 5×104 function evaluations, nondominated solutions in the UEA cover the entire PF. In the reduced
UEA scenario, a limited number of nondominated solutions are selected fromthe UEA using the above-mentioned distance-based selection method. Sparsely
distributed nondominated solutions can be found in Figure 1(iii). In summary,the selection of nondominated solutions in the (reduced) UEA scenario is performed independently from the environmental selection of MOEA/D.
3.2. Average Performance Score (APS)
Here, the APS calculation method [13] is introduced. Suppose that n algorithms A1 , ..., An are compared for a given problem instance based on the HVvalues obtained in multiple runs. For each i ∈ {1, ..., n} and j ∈ {1, ..., n}\{i},let δi,j = 1, if Aj significantly outperforms Ai using the Wilcoxon rank-sum
test with p < 0.05, otherwisePnδi,j = 0. Then, the performance score P (Ai ) isdefined as follows: P (Ai ) = j∈{1,...,n}\{i} δi,j . The score P (Ai ) represents thenumber of algorithms outperforming 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 of the target algorithm is better than othercompared algorithms.
3.3. Parameter settings for MOEA/D
Table 2 shows the settings of the three control parameters (µ, g, and θ) of
MOEA/D. Unless explicitly noted, the default settings in Table 2 were used in10our experimental study. Then, the influence from each control parameter on theperformance of MOEA/D is investigated in a component-wise manner. That is,
we analyzed the effect of µ, g, and θ on MOEA/D individually.
We implemented MOEA/D by modifying the source code of MOEA/D-DE[20] which was downloaded from the jMetal website5 . According to [18], theneighborhood size T of MOEA/D was set to 20. The SBX crossover and thepolynomial mutation were used in MOEA/D. The control parameters of thevariation operators were set as follows: pc = 1, ηc = 20, pm = 1/D, andηm = 20. A simple normalization strategy described in [18] was introduced into
MOEA/D to handle differently scaled objective function values.
4. Experimental results and discussion
Here, we report and discuss the experimental results of MOEA/D with various control parameter settings. See Table 2 for the organization of this section.
Subsections 4.1, 4.2, and 4.3 are also organized as follows: First, we summarize the overall performance of MOEA/D with various configurations on the
13 MOPs based on the APS (Subsections 4.1.1, 4.2.1, and 4.3.1). Then, experimental results on each problem are described (Subsections 4.1.2, 4.2.2, and4.3.2). Finally, we discuss and analyze our results of MOEA/D with each control
parameter setting (Subsections 4.1.3, 4.2.3, and 4.3.3).
Since some real-world problems require the execution of a simulation thattakes a long time to evaluate the solution, the maximum number of functions
evaluations is dependent on the user’s available time [51, 52]. For this reason,
a well-performing MOEA should be able to return a set of good nondominatedsolutions to a user at any time [43]. Also, as pointed out in [44], the comparison
based on the end-of-the-run results does not provide sufficient information about
the performance of MOEAs. Therefore, we mainly discuss our experimental
results based on the anytime performance of MOEA/D, rather than the end-ofthe-run results.
On the one hand, when the number of nondominated solutions obtained by
MOEA/D exceeds a predefined µ value, some of them are discarded from the
population in order to keep µ constant. Therefore, the quality of solutions,regarding the HV metric, directly depends on µ and the environmental selection under the final population scenario. On the other hand, when using the
UEA, all nondominated solutions generated by MOEA/D are stored into the
UEA independently from the environmental selection. Thus, the performanceof MOEA/D is indirectly influenced by µ and the environmental selection under
the reduced UEA scenario. Also, a large number of nondominated solutions canbe examined and selected for the HV calculation for the reduced UEA scenario.
Due to this reason, in most cases, the HV value under the reduced UEA scenario is better than that under the final population scenario in our experimental
results. For example, see Figure 3.
5 http://jmetal.sourceforge.net/

11

4

M =2

2
0

0

10000 20000 30000 40000 50000

4

M =3

APS values

20010000 20000 30000 40000 50000

4
2
0

M =4
0

µ = 25 4µ = 50µ = 400µ = 100 2µ = 300µ = 200µ = 28µ = 55 3µ = 105µ = 210 2µ = 300 1µ = 406µ = 35µ = 56 3µ = 120µ = 220 2µ = 286 1µ = 455

10000 20000 30000 40000 50000

M =2

0

M =3

0

0

µ = 2101µ = 330 1

0

M =5
0

10000 20000 30000 40000 50000

M =4µ = 126 22

10000 20000 30000 40000 50000µ = 406µ = 300µ = 210µ = 28µ = 55µ = 105µ = 455µ = 35µ = 286µ = 56µ = 220µ = 120

10000 20000 30000 40000 50000

M =5µ = 495µ = 126µ = 210µ = 495µ = 330
010000 20000 30000 40000 50000µ = 400µ = 300µ = 25µ = 200µ = 100µ = 50010000 20000 30000 40000 50000

Number of function evaluations

Number of function evaluations(a) Final Population scenario(b) Reduced UEA scenario

Figure 2: Overall performance of MOEA/D with various µ settings on the 13 problems (thefour DTLZ and nine WFG problems) with M ∈ {2, 3, 4, 5} (lower is better). The horizontal
and vertical axes represent the number of function evaluations and the APS values, respectively. A small APS value indicates that a corresponding setting performs relatively better
than the other settings. For the APS value [13], see Subsection 3.2.

4.1. Influence of the population size µ on the performance of MOEA/D under
two scenarios
4.1.1. Overall performance
Figure 2 shows the overall performance of MOEA/D with various µ settings
on the 13 problems with M ∈ {2, 3, 4, 5} under (a) the final population scenarioand (b) the reduced UEA scenario. It should be noted that results for (a) the
final population scenario and (b) the reduced UEA scenario shown in this paper
are obtained from the identical MOEA/D runs.
First, we describe results of MOEA/D for the final population scenario shownin Figure 2(a). For M = 2, while MOEA/D with small µ values (µ ∈ {25, 50})
perform better than that with larger µ values (µ ∈ {100, 200, 300, 400}) justafter the beginning of the search, their APS values gradually deteriorate as thesearch progresses. When M is increasing, such a tendency is more noticeable.
For M ≥ 3, after several function evaluations, the larger µ value is used for
MOEA/D, the better APS value is achieved. The best performance of MOEA/D
is obtained with the largest µ value (µ = 406, 455, and 495 for M = 3, 4, and5) for M ≥ 3, and MOEA/D with the smallest µ value (µ = 25, 28, 35, and 126
12

0.35

Median hypervolume values across all 31 runs

0.30
0.25
0.6

WFG4 (M = 2)
10000 20000 30000 40000 50000

0.5
0.4
0.8
0.7
0.6
0.5
0.4
0.3
0.8
0.7
0.6
0.5
0.4

WFG4 (M = 3)
10000 20000 30000 40000 50000

WFG4 (M = 4)
10000 20000 30000 40000 50000

WFG4 (M = 5)

µ = 200 0.35
µ = 300µ = 100µ = 400 0.30
µ = 50µ = 250.25µ = 406 0.6µ = 300µ = 210µ = 105 0.5µ = 55µ = 28 0.4µ = 455 0.8µ = 286 0.7µ = 220 0.6µ = 120 0.5µ = 560.4µ = 350.3µ = 495 0.8

WFG4 (M = 2)
10000 20000 30000 40000 50000

WFG4 (M = 3)
10000 20000 30000 40000 50000

WFG4 (M = 4)
10000 20000 30000 40000 50000µ = 330 0.7
0.6µ = 2100.5µ = 1260.4µ = 25µ = 50µ = 100µ = 200µ = 300µ = 400µ = 28µ = 55µ = 105µ = 210µ = 300µ = 406µ = 35µ = 56µ = 120µ = 220µ = 286µ = 455µ = 210µ = 330µ = 126

WFG4 (M = 5)10000 20000 30000 40000 50000

10000 20000 30000 40000 50000

Number of function evaluations

Number of function evaluations(a) Final Population scenario(b) Reduced UEA scenarioµ = 495

Figure 3: Performance of MOEA/D with various µ settings on the WFG4 problem withM ∈ {2, 3, 4, 5}. The horizontal and vertical axes represent the number of function evaluations
and the HV values, respectively. The shaded area indicates 25-75 percentiles.

for M = 2, 3, 4, and 5) performs worst among all the configurations for all M .
Next, the results of MOEA/D for the reduced UEA scenario in Figure 2(b)are described below. In contrast to the results for the final population scenario,
MOEA/D with the largest µ value always shows the worst APS value. For
M ∈ {2, 3}, small µ values (µ ∈ {50, 100} and µ ∈ {55, 105}, respectively) aresuitable for MOEA/D at any time. For M ∈ {4, 5}, MOEA/D with a smallµ value also shows a good performance within a short period. After several
function evaluations, MOEA/D with µ = 120 and µ ∈ {210, 330} outperform
other configurations on the four- and five-objective MOPs, respectively.
4.1.2. Results on each problem
Figure 3 shows the anytime performance of MOEA/D with various µ settings
on the WFG4 problem with M ∈ {2, 3, 4, 5} under two scenarios. Although we
do not show results on all the 13 problems here due to the space constraints,
they can be found in Figures S.1 – S.13 in the supplementary file6 of this paper.
6 The supplementary file is also available online (https://sites.google.com/site/moeadanalysis/home/ti-moead-asoc17-supp.pdf)

13

As shown in Figure 3(a), under the final population scenario, while MOEA/Dwith µ ∈ {100, 200} show good performance for M = 2, MOEA/D with thelargest µ value outperforms other configurations after several function evaluations for M ≥ 3. A larger µ value is beneficial to MOEA/D for the final
population scenario. Similar results can be found on other problems shown in
Figures S.1(a) – S.13(a). Of course, there are some exceptions. For example,
results on the WFG1 problem (Figure S.5) show that MOEA/D with a small µvalue exhibits a good performance under the final population scenario.
While good HV values are obtained with a large µ value in most cases forthe final population scenario, MOEA/D with a small µ value performs the best
under the reduced UEA scenario as seen from Figure 3(b). It is interesting to
notice that the performance rank of MOEA/D with the smallest and largestµ values is totally different between the two scenarios on the three- and fourobjective WFG4 problems. For example, for M = 3, while MOEA/D with
µ = 406 performs significantly better than MOEA/D with µ = 28 under the
final population scenario, their performance rank is inverted even on the sameproblem instance in the case of the reduced UEA scenario. Similar conclusions
can be found on most of other problems shown in Figures S.1(b) – S.13(b),except for the results on the DTLZ4, WFG2, WFG3, and WFG8 problems.
A large µ value is suitable for MOEA/D on the DTLZ4, WFG2, WFG3, andWFG8 test problems even when we use the reduced UEA scenario.
According to the percentiles in Figure 3(a) and (b), there is no large variationin the HV values obtained by using all the µ values for M = 2. However, thevariation of the HV values becomes large with increasing M . In particular,a larger variation is observed in results of MOEA/D with small µ values for
M = 5. Thus, it is likely that the performance of MOEA/D with small µ values
on MOPs with a large M is significantly different for each run.
4.1.3. Discussion
We here discuss why the best population size µ is dependent on the choice ofan optimization scenario, the number of objectives M , and the characteristics ofa given problem. As described in Subsections 4.1.1 and 4.1.2, while MOEA/Dwith the largest µ value outperforms other configurations in most cases under
the final population scenario, a small µ value is suitable for MOEA/D on someproblems for the reduced UEA scenario. In particular, the performance rankof MOEA/D with various µ values on the WFG4 problem with M ∈ {3, 4} is
totally different depending on the considered optimization scenario.
This seems to be simply because MOEA/D with a large µ value is capableof keeping a large number of nondominated solutions in the population. As
discussed in [67], such a large number of nondominated solutions are always
beneficial for the HV calculation. For this reason, the results in Subsections
4.1.1 and 4.1.2 show that MOEA/D with a large µ value performs well on mostproblems. However, nondominated solutions in the population do not directly
influence the performance of MOEA/D under the reduced UEA scenario, because solutions selected from the UEA are used for the HV calculation. Also,
MOEA/D can store nondominated solutions whose size is over µ in the UEA.
14

µ = 28µ = 406

0.0

0.2

0.2

0.6

0.6

0.4
0.2
0.0
0.0

0.4
0.2
0.0
0.0

0.0

0.4

0.4

f1 0.6

f3

1.0
0.8

f3

1.0
0.8

0.8

0.8
1.0

0.2

0.6 f 20.2
0.4

f1 0.6

1.0

0.8

0.8
1.0

0.4
0.6 f 2

1.0(a) Final Population scenarioµ = 28µ = 406

0.0

0.2

0.2
0.4

f1 0.6

0.8

0.8
1.0

f3

1.0
0.8

f3

1.0
0.8
0.6

0.6

0.4
0.2
0.0
0.0

0.4
0.2
0.0
0.0

0.0

0.4
0.6 f 2

0.2

0.2
0.4

f1 0.6

1.0

0.8

0.8
1.0

0.4
0.6 f 2

1.0(b) Reduced UEA scenario
Figure 4: Distribution of nondominated solutions in the objective function space on the threeobjective WFG4 problem for the final and reduced UEA scenarios. Nondominated solutions
shown here are found by MOEA/D with µ = 28 and µ = 406 at 5 × 104 function evaluations.
The number of selected nondominated solutions b for the reduced UEA scenario was set to
210 for M = 3. Results of a single run with the median HV value are shown.

Due to this reason, MOEA/D with a small µ value works well under the reduced
UEA scenario.
Figure 4 shows the distribution of nondominated solutions in the objectivefunction space for each scenario found by MOEA/D with µ = 28 and µ =406 on the three-objective WFG4 problem. For the final population scenario,
while MOEA/D with µ = 406 achieves a set of diverse nondominated solutions,
MOEA/D with µ = 28 can provide only a small number of nondominatedsolutions. Unlike the results for the final population scenario, MOEA/D withboth µ settings successfully maintain the densely distributed solutions covering
the entire PF in the UEA for the reduced UEA scenario. It is likely that
MOEA/D with a small µ value cannot keep good nondominated solutions in
the population but has an ability to generate well-distributed solutions.
We do not intend to argue that a small µ value is always beneficial for
MOEA/D under the reduced UEA scenario. Our results in Subsection 4.1.2
show that a (relatively) large µ value is necessary for MOEA/D to find good15

nondominated solutions on some problems, especially MOPs with a large M
or an irregular shape of the PF. As seen from the results on the five-objective
WFG4 problem (Figure 3), while MOEA/D with the smallest µ value (µ =126) achieves the highest HV value until 14 000 function evaluations, it clearly
performs worst after that. Since the objective function space is exponentially
increasing with M , a large population size is necessary for the search if M islarge.
Results on the DTLZ4, WFG2, WFG3, and WFG8 problems (Figure S.4,
S.6, S.7, and S.12, respectively) show that a large µ value is beneficial for
MOEA/D for all M independently of the choice of an optimization scenario.
The shape of the PF of the WFG2 and WFG3 problems are discontinuous andpartially degenerate [64], respectively. Since the Pareto optimal solutions do
not exist in most of the decomposed subproblems of MOPs having discontinuous and degenerate PFs [68], a large µ value, which makes a large numberof subproblems, is necessary for MOEA/D to find well-approximated solutions.
For this reason, MOEA/D with a small µ value performs poorly on the WFG2
and WFG3 problems. Also, the mapping from the solution space to the objective function space f : S → RM is biased in the DTLZ4 and WFG8 problems.
The existence of the bias makes the problem more difficult, and thus MOEA/Dwith a small µ value cannot find good solutions on the DTLZ4 and WFG8
problems.
4.2. Impact of the scalarizing function g on the performance of MOEA/D under
two scenarios
4.2.1. Overall performance
Figure 5 shows the overall performance of MOEA/D with the three scalarizing functions (g chm , g chd , and g pbi with θ = 5) over all the 13 test problems.
For details of the three functions, see Section 2. The influence of θ on theperformance of MOEA/D using g pbi is investigated in Subsection 4.3.
The results on the two- and three-objective 13 problems for the final population scenario show that the performance of MOEA/D with the two Chebyshev
functions (g chm and g chd ) is better than that with g pbi at anytime (Figure5(a)). For M = 4, MOEA/D with g pbi outperforms the other configurations after 38 000 function evaluations. For M = 5, MOEA/D with g pbi always exhibits
the best performance.
The results for the reduced UEA scenario (Figure 5(b)) are significantly different from the results for the final population scenario described above. As
seen from Figure 5(b), MOEA/D with the two Chebyshev functions are competitive with each other. The g chm and g chd are the best scalarizing functionsfor M ∈ {2, 4, 5} and M = 3, respectively. Interestingly, while MOEA/D withg pbi exhibits a good performance on MOPs with more than four objectives after
several function evaluations under the final population (Figure 5(a)), its performance is worst among the three configurations for all M under the reduced
UEA scenario (Figure 5(b)).

16

2.0

g pbi 2.0

1.5

1.5

M =2

1.0

g chd 0.5

0.0

0.0

0

1.5

M =3

1.0

APS values

10000 20000 30000 40000 50000

g pbi 2.0

g chd 0.5
0.0

010000 20000 30000 40000 50000

1.5

M =4
0

10000 20000 30000 40000 50000

1.5

0

g pbi

M =4

g chm 2.0

M =5

0

0

10000 20000 30000 40000 50000

g pbi

M =5

g pbi 0.5
0.0

10000 20000 30000 40000 50000

g chd
g chm

g chd 1.0

0.5

10000 20000 30000 40000 50000

g pbi 0.5
0.0

g chm
g chd

1.5

1.0
0.0

M =3

g chd 1.0

0.5
2.0

g chm 2.0

10000 20000 30000 40000 50000

g pbi

1.5

1.0
0.0

0

g chm 1.0

0.0

g chd
g chm

1.5

0.5
2.0

M =2

g chm 1.0

0.5
2.0

g pbi

g chd
g chm

0

10000 20000 30000 40000 50000

Number of function evaluations

Number of function evaluations(a) Final Population scenario(b) Reduced UEA scenario

Figure 5: Overall performance of MOEA/D with the three scalarizing functions (g chm , g chd ,and g pbi with θ = 5) on the 13 problems (the four DTLZ and nine WFG problems) withM ∈ {2, 3, 4, 5} (lower is better). The horizontal and vertical axes represent the number offunction evaluations and the APS values, respectively. For the APS value [13], see Subsection
3.2.

4.2.2. Results on each problem
Figure 6 shows the anytime performance of MOEA/D with the three scalarizing functions on the WFG4 problem with M ∈ {2, 3, 4, 5}. Results on all the
13 problems can be found in Figures S.14 – S.26 in the supplementary file.
As shown in Figure 6, the results on the WFG4 problem are similar to
the aggregated results on all of the 13 problems shown in Figure 5. For bothscenarios, MOEA/D with the two Chebyshev functions perform significantly
better than that with g pbi on the two- and three-objective WFG4 problems. For
the final population scenario, with increasing M , MOEA/D with g pbi achieves
the highest HV values among the three configurations in the later search stage.
However, for the reduced UEA scenario, the worst HV values are obtainedby using g pbi , and MOEA/D with the two Chebyshev functions show better
performance at any time. As seen from the percentiles in Figure 6(a) and (b),the variation of the HV values achieved by using the two Chebyshev functions
is large for M = 5. In contrast, the performance of MOEA/D with g pbi is notmuch different for each run. The quality of nondominated solutions found by

17

0.35

0.35

g chd

Median hypervolume values across all 31 runs

0.30

g chm

0.25
0.20
0.6

WFG4 (M = 2)010000 20000 30000 40000 50000

0.5

g chm

WFG4 (M = 3)
0

g pbi

0.6

g chm

0.5

0.8
0.7
0.6
0.5
0.4
0.3

g pbi

10000 20000 30000 40000 50000

0.7

0.4

0.6

WFG4 (M = 2)
0

g pbi

10000 20000 30000 40000 50000

g chd

0.5

0.4
0.3

g chm

0.25

g pbi 0.20
g chd

g chd

0.30

WFG4 (M = 4)
0

g chd

10000 20000 30000 40000 50000

g pbi
g chd

WFG4 (M = 5)
0

g chm

10000 20000 30000 40000 50000

g chm

0.4
0.3

WFG4 (M = 3)
0

0.7

g chm

0.6

g chd

0.5
0.4
0.8
0.7
0.6
0.5
0.4
0.3

g pbi

10000 20000 30000 40000 50000

WFG4 (M = 4)
0

g pbi

10000 20000 30000 40000 50000

g chm
g chd

WFG4 (M = 5)
0

g pbi

10000 20000 30000 40000 50000

Number of function evaluations

Number of function evaluations(a) Final Population scenario(b) Reduced UEA scenario

Figure 6: Performance of MOEA/D with the three scalarizing functions (g chm , g chd , and g pbi
with θ = 5) on the WFG4 problem with M ∈ {2, 3, 4, 5}. The horizontal and vertical axes
represent the number of function evaluations and the HV values, respectively. The shaded
area indicates 25-75 percentiles.

using g pbi is likely to be stable even on MOPs with a large M .
Results on the DTLZ2, WFG5, WFG6, WFG7, WFG8, and WFG9 problems are similar to results on the WFG4 problem. On the DTLZ1 problem,
MOEA/D with g chd and g pbi show good performance (Figure S.14). On the
DTLZ3 problem, g chd is the best scalarizing function for all M (Figure S.16).
On the DTLZ4 problem with M ≥ 3, MOEA/D with g pbi performs the best
(Figure S.17). On the WFG2 and WFG3 problems with all M , the performanceof MOEA/D with the two Chebyshev functions is better than that with g pbi
even for the final population scenario (Figure S.19 and S.20).
4.2.3. Discussion
A number of previous work report that MOEA/D with g chm is not capableof obtaining well-distributed nondominated solutions in the objective function
space (e.g., [18, 30, 54, 69]). Recent studies also report the promising performance of MOEA/D with g pbi on MOPs with more than four objectives [24, 30].
Our results for the final population described in Subsections 4.2.1 and 4.2.2 areconsistent with these previous studies. However, the performance of MOEA/D

18

Normalized objective function values

1.0
0.8
0.6
0.4
0.2
0.0

1.0
0.8
0.6
0.4
0.2
0.0

g chm

1

2

3

4

Objective number

5

2

3

4

Objective number

g chd

1

2

3

4

5

1.0
0.8
0.6
0.4
0.2
0.0

Objective number(a) Final Population scenario

g chm

1

1.0
0.8
0.6
0.4
0.2
0.0

5

1.0
0.8
0.6
0.4
0.2
0.0

g chd

1

2

3

4

5

Objective number(b) Reduced UEA scenario

1.0
0.8
0.6
0.4
0.2
0.0

g pbi

1

2

3

4

5

4

5

Objective number

g pbi

1

2

3

Objective number

Figure 7: Parallel coordinates of the objective function values of nondominated solutions on
the three-objective WFG4 problem for the final and reduced UEA scenarios. Nondominatedsolutions shown here are found by MOEA/D with the three scalarizing functions (g chm , g chd ,and g pbi with θ = 5) at 5 × 104 function evaluations. The number of selected nondominatedsolutions b for the reduced UEA scenario was set to 210 for M = 5. The horizontal andvertical axes represent the objective number and the normalized objective function values,
respectively. Results of a single run with the median HV value are shown.

with the two Chebyshev functions (g chm and g chd ) is significantly better thanthat with g pbi on most problem instances under the reduced UEA scenario.
In summary, the best scalarizing function for MOEA/D is significantly dependent on the choice of an optimization scenario. Below, we discuss this result(mainly with respect to the comparison between the PBI function and the two
Chebyshev functions).
Figure 7 shows the parallel coordinates of the objective function values of
nondominated solutions found by MOEA/D with the three scalarizing functions
on the five-objective WFG4 problem under (a) the final population scenario and(b) the reduced UEA scenario. It should be recalled that g pbi and g chm are
the best scalarizing functions on the five-objective WFG4 problem for the final
population and reduced UEA scenarios, respectively (see Figure 6).
The number of nondominated solutions in Figure 7(a) is 146 for g chm , 150
for g chd , and 210 for g pbi . That is, MOEA/D with g pbi obtains the largest
number of nondominated solutions among the three configurations for the final
population scenario. The nondominated solutions found by MOEA/D withg pbi are also uniformly distributed in the objective function space. For these
two reasons, it seems that MOEA/D with g pbi achieves the highest HV value
at the end of the search. In contrast to the results for the final population
scenario, as shown in Figure 7(b), MOEA/D with g chm obtains well-distributed
nondominated solutions in the objective function space under the reduced UEA

19

g chm

g chd

g pbi

0.0

0.2

0.4

f10.6

0.8

0.8
1.0

0.4
0.6f 2

0.2

0.0

0.2

0.4

f10.6

1.0

0.8

0.8
1.0

0.4
0.6f 2

0.2

f3

1.0
0.8
0.6
0.4
0.2
0.0
0.0

f3

1.0
0.8
0.6
0.4
0.2
0.0
0.0

f3

1.0
0.8
0.6
0.4
0.2
0.0
0.0

0.0

0.2

0.4

f10.6

1.0

0.8

0.8
1.0

0.4
0.6f 20.2

1.0(a) Final Population scenario
g chd

g pbi

1.0
0.8
0.6
0.4
0.2
0.0
0.0

0.0

0.2

0.4

f10.6

0.8

0.8
1.0

1.0

0.4
0.6f 2

0.2

0.0

0.2

0.4

f10.6

0.8

0.8
1.0

0.4
0.6f 2

0.2

1.0

f3

1.0
0.8
0.6
0.4
0.2
0.0
0.0

f3

1.0
0.8
0.6
0.4
0.2
0.0
0.0

f3

g

chm

0.0

0.2

0.4

f10.6

0.8

0.8
1.0

0.4
0.6f 20.2

1.0(b) Reduced UEA scenario

Figure 8: Distribution of nondominated solutions in the objective function space on the threeobjective WFG4 problem for the final and reduced UEA scenarios. Nondominated solutions
shown here are found by MOEA/D with g chm , g chd , and g pbi at 5 × 104 function evaluations.
The number of selected nondominated solutions b for the reduced UEA scenario was set to
210 for M = 3. Results of a single run with the median HV value are shown.

scenario. The uniformity of the distribution of nondominated solutions found by
MOEA/D with g chm is slightly better than those obtained by the other methods.
For further discussion, we show the distribution of nondominated solutionsin the objective function space for both the final and reduced UEA scenarios on
the three-objective WFG4 problem in Figure 8. The distribution on the threeobjective WFG4 problem in Figure 8 is similar to the distribution on the fiveobjective WFG4 problem in Figure 7. In the final population scenario (Figure8(a)), while the uniformly distributed nondominated solutions are obtained byusing g pbi and g chd , the distribution of solutions found by MOEA/D with g chm isbiased. However, it seems that MOEA/D with all the three scalarizing functions
find evenly distributed nondominated solution under the reduced UEA scenario(Figure 8(b)).
As pointed out in [37, 54], it seems that MOEA/D with g chm is capable ofgenerating well-approximated solutions but cannot maintain them in its population. Such a problem of g chm is addressed by incorporating the UEA, which
automatically stores good nondominated solutions independently from the population, into MOEA/D. Due to this reason, MOEA/D with g chm performs thebest on the five-objective WFG4 problem under the reduced UEA scenario, re-20

garding the HV metric. This observation can be found in some of the other test
problems.
4.3. Effect of the penalty value θ for the PBI function on the performance of
MOEA/D4.3.1. Overall performance
Figure 9 shows the overall performance of MOEA/D using the PBI function
g pbi with θ ∈ {0.1, 0.5, 1, 2, 3, 5, 8, 10} on the 13 problems with M ∈ {2, 3, 4, 5}
under (a) the final population scenario and (b) the reduced UEA scenario. Unlike the results of µ and g presented in Subsections 4.1 and 4.2, the performance
rank of each configuration is not drastically changed during the search process
for both scenarios. That is, a configuration that performs well at the beginningof the search exhibits good performance also at the end of the search.
As shown in Figure 9(a), good APS values are obtained by θ ∈ {2, 3} forM ∈ {2, 3} for the final population scenario. The best penalty value is θ ∈ {3, 5}
and θ ∈ {5, 8} for M = 4 and M = 5, respectively. It seems that MOEA/D
using g pbi with θ = 5 works well for all M . In contrast, θ ∈ {0.1, 0.5, 1} areclearly inappropriate penalty values for all M .
Figure 9(b) shows results for the reduced UEA scenario. It seems that theperformance rank under the reduced UEA scenario is not significantly different
from that under the final population scenario. Although the smallest APS value
is obtained by θ = 0.1 until 6 000 function evaluations for M ≥ 3, MOEA/D
using g pbi with θ ∈ {3, 5} perform well among all the configurations for all M .
While the performance rank of θ = 0.1 is improved M ≥ 3, MOEA/D usingg pbi with θ ∈ {0.5, 1} performs poorly even for the reduced UEA scenario.
4.3.2. Results on each problem
Figure 10 shows the anytime performance of MOEA/D using g pbi with various θ values on the WFG4 problem with M ∈ {2, 3, 4, 5}. Figures S.27 – S.39
in the supplementary file show results on all the 13 problems.
As seen from Figure 10(a), for the final population scenario, θ = 1, 2, 2, and 3
are suitable values on the WFG4 problem with M = 2, 3, 4, and 5, respectively.
That is, the best penalty value θ increases with respect to the number objectives
M on the WFG4 problem. Interestingly, as shown in Figure 10(b), in contrast
to the results for the final population scenario, MOEA/D using g pbi with θ = 0.1
achieves the good HV value under the reduced UEA scenario for M ∈ {2, 3}.
According to the percentiles in Figure 10(a) and (b), for M = 3, the variationof the HV values obtained by using θ ∈ {0.1, 0.5} is larger than that by using
other θ values. Interestingly, in Figure 10, the 25th and 75th percentile values
are not evenly distributed. For example, on the results for M = 3 in Figure10(a), while the median and the 25th percentile values are almost the same,the difference between the median and 75th percentile values is large. Only for
M = 5, a large variation is observed in results of MOEA/D with θ = 3. While
MOEA/D with θ = 3 obtains the best HV value under both scenarios regarding
the median HV value, it performs similarly to MOEA/D with θ ∈ {0.1, 0.5, 1, 2}
for a particular run as seen from its 75th percentile values.
21

7

M =2

5
3
1
0

0

7

10000 20000 30000 40000 50000

M =3

5

APS values

3
1
0

0

7

10000 20000 30000 40000 50000

M =4

5
3
1
0

0

7

10000 20000 30000 40000 50000

M =5

5
3
1
0

0

10000 20000 30000 40000 50000

θ = 0.5
θ = 0.1
θ = 10
θ=8
θ=5
θ=1
θ=2
θ=3

7

θ = 0.5
θ=1
θ = 0.1
θ = 10
θ=8
θ=5
θ=3
θ=2

7

θ = 0.5
θ=1
θ = 0.1
θ = 10
θ=8
θ=2
θ=5
θ=3

7

θ = 0.5
θ=1
θ=2
θ = 0.1
θ=3
θ = 10
θ=8
θ=5

7

M =2

5
3
1
0

0

10000 20000 30000 40000 50000

M =3

5
3
1
0

0

10000 20000 30000 40000 50000

M =4

5
3
1
0

0

10000 20000 30000 40000 50000

M =5

5
3
1
0

0

10000 20000 30000 40000 50000

Number of function evaluations

Number of function evaluations(a) Final Population scenario(b) Reduced UEA scenario

θ = 0.1
θ = 10θ = 0.5
θ=8
θ=5
θ=1
θ=2
θ=3
θ=1
θ = 0.5
θ = 10
θ=8
θ = 0.1
θ=5
θ=3
θ=2
θ=1
θ = 0.5
θ = 10θ = 0.1
θ=8
θ=2
θ=5
θ=3
θ=1
θ = 0.5
θ=2
θ = 0.1
θ=3
θ = 10θ=8
θ=5Figure 9: Overall performance of MOEA/D using the PBI function g pbi with various θ values
on the 13 problems with M ∈ {2, 3, 4, 5}. The horizontal and vertical axes represent the
number of function evaluations and the APS values, respectively. For the APS value [13], see
Subsection 3.2.

Due to space constraints, results on other test problems are not described
here. Similar to the results on the WFG4 problem, the best θ value significantly
depends on the considered scenario on most test problems.
4.3.3. Discussion
Our results described in Subsections 4.3.1 and 4.3.2 show that the best setting of θ is significantly dependent on the characteristics of a given problem.
For example, while MOEA/D using g pbi with θ = 2 performs well on the WFG4
problem with M = 3 (Figure 10(a)), θ = 10 is the best setting for solving thethree-objective DTLZ1 problem (Figure S.27). Such conclusions are consistent
with previous studies [33, 35, 36, 70]. For example, experimental results in [33]show that MOEA/D using g pbi with θ = 0.1 performs better than that withother θ values on multi- and many-objective knapsack problems with the convex PFs. However, a small penalty value like θ = 0.1 is not appropriate for
MOEA/D to solve MOPs with nonconvex PFs [36].
In addition to the characteristics of a given problem, our results in Subsections 4.3.1 and 4.3.2 newly reveal that a suitable θ value also depends on atarget optimization scenario. For example, as shown in Figure 10, the best θ
22

0.35
0.30

Median hypervolume values across all 31 runs

0.25
0.20
0.15
0.6
0.4

WFG4 (M = 2)
10000 20000 30000 40000 50000

WFG4 (M = 3)

0.2
0.0
0.7
0.5
0.3
0.1
0.0
0.8

10000 20000 30000 40000 50000

WFG4 (M = 4)

10000 20000 30000 40000 50000

0.6
0.4

WFG4 (M = 5)

0.2
0.0

10000 20000 30000 40000 50000

θ = 1 0.35
θ=2
θ = 3 0.30
θ = 5 0.25
θ=8
θ = 10 0.20
θ = 0.5
θ = 0.1 0.15
θ = 2 0.6
θ=3
θ = 5 0.4
θ=8
θ = 10θ = 0.1 0.2
θ = 0.5
θ = 1 0.0
θ = 2 0.7
θ=3
θ = 5 0.5
θ=8
θ = 10 0.3
θ = 0.1
θ = 0.5 0.1
θ = 1 0.0
θ = 3 0.8
θ=5
θ = 8 0.6
θ = 10 0.4
θ = 0.1
θ = 1 0.2
θ = 0.5
θ = 2 0.0

WFG4 (M = 2)
10000 20000 30000 40000 50000

WFG4 (M = 3)
10000 20000 30000 40000 50000

WFG4 (M = 4)
10000 20000 30000 40000 50000

WFG4 (M = 5)10000 20000 30000 40000 50000

Number of function evaluations

Number of function evaluations(a) Final Population scenario(b) Reduced UEA scenario

θ = 0.5
θ = 0.1
θ=1
θ=2
θ=3
θ=5
θ=8
θ = 10θ = 0.1
θ=2
θ=3
θ=5
θ=8
θ = 10θ = 0.5
θ=1
θ=2
θ=3
θ=5
θ=8
θ = 10θ = 0.1
θ = 0.5
θ=1
θ=3
θ=5
θ=8
θ = 10θ = 0.1
θ = 0.5
θ=2
θ=1

Figure 10: Performance of MOEA/D using the PBI function g pbi with various θ values on
the WFG4 problem with M ∈ {2, 3, 4, 5}. The horizontal and vertical axes represent the
number of function evaluations and the HV values, respectively. The shaded area indicates25-75 percentiles.

value is totally different on the WFG4 problem (nonconvex) for each scenario.
Below, we discuss the reason why the θ = 0.1 is the best penalty value on thethree-objective WFG4 problem for the reduced UEA scenario.
Figure 11 shows the distribution of nondominated solutions found by MOEA/D
using g pbi with θ = 0.1 and θ = 2 on the three-objective WFG4 problem forthe final and reduced UEA scenarios. It should be noted that θ = 2 is thebest parameter setting on the WFG4 problem with M = 3 for the final population scenario (see Figure 10). As seen from Figure 11(a), while MOEA/D with
θ = 0.1 can find only the three extreme solutions under the final population scenario, well-distributed nondominated solutions are obtained by θ = 2. However,
Figure 11(b) shows that in addition to the extreme solutions, MOEA/D with
θ = 0.1 achieves some nondominated solutions distributed on the rim and insideregions of the PF. Note that the HV value of nondominated solutions obtainedby using θ = 0.1 is higher than that by θ = 2 in Figure 11(b).
It is pointed out that the HV indicator does not always assess the uniformity and diversity of the distribution of nondominated solutions [71]. For
further investigation, we used the inverted generational distance (IGD) [65],which evaluates both convergence and diversity performance of MOEAs. The

23

θ = 0.1

θ=2

0.0

0.2

0.2

0.6

0.6

0.4
0.2
0.0
0.0

0.4
0.2
0.0
0.0

0.0

0.4

0.4

f1 0.6

f3

1.0
0.8

f3

1.0
0.8

0.8

0.8
1.0

0.2

0.6 f 20.2
0.4

f1 0.6

1.0

0.8

0.8
1.0

0.4
0.6 f 2

1.0(a) Final Population scenario
θ = 0.1

θ=2

0.0

0.2

0.2
0.4

f1 0.6

0.8

0.8
1.0

f3

1.0
0.8

f3

1.0
0.8
0.6

0.6

0.4
0.2
0.0
0.0

0.4
0.2
0.0
0.0

0.0

0.4
0.6 f 2

0.2

0.2
0.4

f1 0.6

1.0

0.8

0.8
1.0

0.4
0.6 f 2

1.0(b) Reduced UEA scenario
Figure 11: Distribution of nondominated solutions in the objective function space on the threeobjective WFG4 problem for the final and reduced UEA scenarios. Nondominated solutions
shown here are found by MOEA/D using the PBI function g pbi with θ = 0.1 (left) and with
θ = 2 (right) at 5 × 104 function evaluations. Results of a single run with the median HVvalue are shown.

generational distance (GD) and maximum spread (MS) [65] indicator values of
nondominated solutions are also reported to discuss convergence and diversity
performance of MOEA/D, individually.
Below, let A be a set of nondominated solutions selected from the populationor the UEA. A set of reference vectors are denoted as A∗ , where each elementof A∗ is a Pareto-optimal solution. The GD value is the average distance fromeach solution in A to its nearest reference vector in A∗ in the objective function
space:


n
o
X

1 
GD(A) =
min ED f (x), f (z)  ,
(7)
|A|
z ∈A∗
x∈A
where ED(x, y) is the Euclidean distance between x and y. In contrast to the
GD metric in equation (7), the IGD value is the average distance from eachreference vector in A∗ to its nearest solution in A in the objective function24space as follows:



o
n

1 
IGD(A) =min ED f (x), f (z)  .
|A∗ |
∗ x∈A
z ∈A
X

(8)

The MS value of A is the average of the maximum and minimum objectivefunction values of solutions in A for each objective as follows:
v
uM
2
uX
MS(A) = t
max {fi (x)} − min {fi (x)} .
(9)x∈Ai=1 x∈A
Figure 12 shows the MS, GD, and IGD indicator values of nondominatedsolutions found by MOEA/D using g pbi with various θ values on the threeobjective WFG4 problem under (a) the final population scenario and (b) the
reduced UEA scenario. Before calculating each indicator value, all objectivefunction values were normalized as described in Subsection 3.1. One may wonder that GD values for the final population scenario are significantly better thanthose for the reduced UEA scenario. This is because the number of nondominated solutions in the population is usually smaller than that in the UEA (e.g.,
see Figure 11), and a small number of solutions is possibly beneficial for GD as
discussed in [67].
As seen from Figure 12, MOEA/D with θ ∈ {0.1, 2, 3, 5, 8, 10} achieve highMS values for both scenarios. The best GD values are also obtained by using θ =0.1 for the reduced UEA scenario. Thus, nondominated solutions in the UEA
obtained by using θ = 0.1 are close to the PF. This is because the contour linesof g pbi with θ = 0.1 and the weighted sum function are similar [33]. Therefore,
MOEA/D using g pbi with θ = 0.1 shows good HV values on the WFG4 problem
with M = 3.
It is interesting to notice that the performance rank based on the IGD metricis slightly different from that based on the HV indicator. Although the best HVvalues are obtained by using θ = 0.1 under the reduced UEA scenario (Figure10), MOEA/D with θ = 2 achieves the lowest IGD values (Figure 12). Since
the IGD metric assesses the uniformity of the distribution, the IGD value of
nondominated solutions found by MOEA/D with θ = 0.1, whose distribution isbiased to specific regions as shown in Figure 11, is worse than that by MOEA/Dwith θ = 2.
5. Further analysis
On the one hand, in Section 4, we examined the influence of each control parameter on MOEA/D by keeping other parameters default (see Table 2). Thisis because we wanted to examine the effect of changing each parameter in acomponent-wise manner. If two or more parameters are perturbed simultaneously, which parameter improves or degrades the performance of MOEA/D may
be unclear.
25θ = 10 1.75
θ=8
θ=5
θ=2
θ = 3 1.70
θ = 0.1
θ = 0.5
θ=1

2.0

MS

1.5
1.0
0.5
0.0

GD

0.10
0.05

IGD

0.00

10000 20000 30000 40000 50000

θ=1
θ = 10
θ=8
θ = 0.5
θ=5
θ=3
θ=2
θ = 0.1
10000 20000 30000 40000 50000

θ = 1 0.14
θ = 0.5
θ = 0.1
θ = 10θ = 8 0.07
θ=5
θ=3
θ=2

0.6
0.4
0.2
0.0

10000 20000 30000 40000 50000

θ = 10 0.15
θ=8
θ=5
θ = 3 0.10
θ=2
θ = 0.5 0.05
θ = 0.1
θ=1

0.15

0.00

1.65

10000 20000 30000 40000 50000

θ=8
θ = 10
θ=5
θ=2
θ = 0.1
θ=3
θ = 0.5
θ=1

0.00

10000 20000 30000 40000 50000

θ = 0.5
θ=1
θ = 0.1
θ = 10
θ=8
θ=5
θ=3
θ=2
10000 20000 30000 40000 50000

Number of function evaluations

Number of function evaluations(a) Final Population scenario(b) Reduced UEA scenario

Figure 12: The MS, GD, and IGD indicator values of nondominated solutions found by
MOEA/D using g pbi with various θ values on the three-objective WFG4 problem under (a)the final population scenario and (b) the reduced UEA scenario. A low GD value and ahigh MS value indicate that the corresponding method has good convergence and diversity
performance, respectively. A method showing a low IGD value performs well in terms of bothconvergence and diversity. The shaded area indicates 25-75 percentiles.

On the other hand, analysis of dependencies between multiple control parameters is important and interesting. It is expected that some control parameters(e.g., µ and g) are correlated with each other. Also, we used the default parameters (see Table 2), which are commonly used in the literature, for our analysis.
However, the effect of the three control parameters (µ, g, and θ) on MOEA/D
may be different when well-tuned parameters are used instead of the default
parameters. Thus, whether our analysis results in Section 4 are representative
or not is not guaranteed.
This section presents further analysis of MOEA/D in order to address theabove-mentioned issues. Dependencies between two control parameters are examined in Subsection 5.1. Analysis of each control parameter with the best
configuration is presented in Subsection 5.2.
5.1. Analysis of dependencies between two control parameters
While we discussed the results based on the anytime performance of MOEA/D
in Section 4, we report only the end-of-the-run results in this subsection. Although the end-of-the-run results do not provide sufficient information as de26

Hypervolume

0.35

WFG4 (M = 2)

0.34
0.33
0.32

WFG4 (M = 3)0.56
g chm 0.54
g chd 0.52
0.50
g pbi
0.48
0.46

25 100

200

300

400

g chd 0.70
g chm
g pbi

28 105210 300µ

406

WFG4 (M = 5)

0.65

g pbi
g chm

0.60
0.55

g chd

0.50
126 210

µ

330

495

µ

(a) Final Population scenario

Hypervolume

0.35

WFG4 (M = 2)

g chm

0.58

WFG4 (M = 3)0.56
g chd 0.54

0.34

g

0.33

pbi

0.52
0.50

0.32

25 100

200

300

400

0.48

WFG4 (M = 5)g chd 0.775
0.750
g chm 0.725
0.700
0.675
pbi
g 0.6500.625
28 105210 300µ

406µ

126 210

g chm
g chd
g pbi

330

495

µ

(b) Reduced UEA scenario

Figure 13: Influence of µ on the performance of MOEA/D with the three scalarizing functions(g chm , g chd , and g pbi ) on the WFG4 problem with M ∈ {2, 3, 5}. The median HV value at
50 000 evaluations among 31 runs is shown.

scribed in Section 4, analyzing dependencies between two parameters based on
the anytime performance is not a trivial task. Such an analysis is our future
research.
5.1.1. Dependencies between µ and g
Figure 13 shows the influence of µ on the performance of MOEA/D with thethree scalarizing functions (g chm , g chd , and g pbi ) on the WFG4 problem withM ∈ {2, 3, 5}. Results on other test problems can be found in Figures S.40 –S.52 in the supplementary file.
As seen from Figure 13, the best µ value is different depending on scalarizing
functions. For example, as shown in Figure 13(a), for the final population
scenario, MOEA/D with the two Chebyshev functions need a relatively large µ
value for all M . In contrast, small µ values are suitable for g pbi . In addition, as
discussed in Subsection 4.1, suitable µ values also depend on the choice of an
optimization scenario. Figure 13(b) indicates that small µ values are appropriate
for all the three scalarizing functions on the WFG4 problem under the reduced
UEA scenario. For M ∈ {2, 3}, in most cases, the increase of µ deteriorates theperformance of MOEA/D under the reduced UEA scenario. In summary, µ andg are correlated with each other, and the relation between them is dependent
on the choice of a scenario.

27

Hypervolume

WFG4 (M = 2)
0.350
0.325
0.300
0.275
0.250
0.225
0.200
0.175
25 100200

300

θ = 1 0.6
θ = 2 0.5
θ=3
θ = 5 0.4
θ = 8 0.3
θ = 10θ = 0.5 0.2
θ = 0.1 0.1

400

WFG4 (M = 3)

28 105µ

210 300

θ = 1 0.8
θ = 2 0.7
θ = 3 0.6
θ = 5 0.5
θ = 8 0.4
θ = 10 0.3
θ = 0.1 0.2
θ = 0.5 0.1

406

WFG4 (M = 5)

126 210

µ

330

θ=3
θ=5
θ=8
θ = 10θ = 0.1
θ=2
θ=1
θ = 0.5495µ(a) Final Population scenario
WFG4 (M = 2)Hypervolume

0.35
0.34
0.33
0.32

θ = 0.5 0.6
θ = 0.1
θ = 1 0.5
θ=2
θ = 3 0.4
θ=5
θ = 8 0.3
θ = 1025 100

200

300

µ

400

0.2

WFG4 (M = 3)

28 105210 300406µ

θ = 0.1
θ = 1 0.7
θ=2
θ = 3 0.6
θ = 5 0.5
θ=8
θ = 10 0.4
θ = 0.5
0.3

WFG4 (M = 5)

126 210

330

θ=3
θ=5
θ=8
θ = 10θ = 0.1
θ=2
θ = 0.5
θ=1

495

µ

(b) Reduced UEA scenario

Figure 14: Influence of µ on the performance of MOEA/D using g pbi with various θ values on
the WFG4 problem with M ∈ {2, 3, 5}. The median HV value at 50 000 evaluations among
31 runs is shown.

5.1.2. Dependencies between µ and θ
Figure 14 shows the impact of µ on the performance of MOEA/D using g pbi
with various θ values on the WFG4 problem with M ∈ {2, 3, 5}. Results onother test problems are shown in Figures S.53 – S.65 in the supplementary file.
Unlike the analysis results of µ and g, the effect of µ values on the performanceof MOEA/D with various θ values is relatively moderate, except for some results (e.g., results for θ = 1 on the three-objective WFG4 problem). Thus,
the performance rank of MOEA/D with different θ values is not significantly
influenced by µ.
5.1.3. Comparison of the two Chebyshev functions and the PBI function withvarious θ values
Figure 15 provides a comparison of the two Chebyshev functions (g chm andg chd ) and g pbi with θ ∈ {0.1, ..., 10} on the WFG4 problem with M ∈ {2, 3, 5}.
For g chm and g chd , the same HV values are shown for each θ value in Figure
15. Results on other test problems can be found Figures S.79 – S.91 in the
supplementary file.
As seen from Figure 15, the performance of MOEA/D with g pbi significantly depends on the θ value for both optimization scenarios. Although such
results under the final population scenario have already been reported in previous work (e.g., [35, 70]), similar results are observed under the reduced UEA
scenario. MOEA/D with a particular θ value performs better than that withthe two Chebyshev functions. Whereas the poor results are obtained by the28Hypervolume

WFG4 (M = 2)
0.350
0.325
0.300
0.275
0.250
0.225
0.200
0.175

g chd 0.6
0.5

WFG4 (M = 3)g chm 0.4
0.3
g pbi 0.2
0.1
0.1 2.0

5.0
θ

8.0 10.0

0.1 2.0

5.0
θg chd 0.8
0.7
0.6
chm
g 0.5
0.4
g pbi 0.3
0.2
0.1

8.0 10.0

WFG4 (M = 5)

g pbi
g chd
g chm

0.1 2.0

5.0
θ

8.0 10.0(a) Final Population scenario

Hypervolume

0.35

WFG4 (M = 2)

g chd 0.6

WFG4 (M = 3)

g chd

0.34

g chm 0.5

g chm

0.33

g pbi

0.4

g pbi

0.32

0.1 2.0

5.0
θ

8.0 10.0

0.3

0.1 2.0

5.0
θ

8.0 10.0

WFG4 (M = 5)
0.8

g chm

0.7

g chd

0.6
0.5

g pbi

0.4
0.3

0.1 2.0

5.0
θ

8.0 10.0(b) Reduced UEA scenario

Figure 15: Comparison of the two Chebyshev functions (g chm and g chd ) and g pbi with various
θ values on the WFG4 problem with M ∈ {2, 3, 5}. The median HV value at 50 000 evaluations
among 31 runs is shown.

g pbi with θ = 5 under the reduced UEA scenario in Subsection 4.2, it is likelythat MOEA/D using g pbi with a suitable θ value performs well. Interestingly,
the performance of MOEA/D regarding the HV metric is drastically improved
or degraded at around θ ∈ {1, 2, 3}. The reason is unclear, and thus its further
analysis is our future research topic.
5.2. Analysis using well-tuned parameters
In Section 4, except for a single parameter to be examined, other parameters
were set to the default values as shown in Table 2. In contrast, here, we examine
the effect of each parameter on MOEA/D by using the best parameters.
We conducted a grid search to find the best configuration for each control
parameter. For the three control parameters (µ, g, and θ), all the values in Table
2 were examined. In addition to the three control parameters (µ, g, and θ), theperformance of MOEA/D with four neighborhood sizes T ∈ {10, 20, 40, 80} wasinvestigated. We evaluated the performance of MOEA/D with all possible combinations of the four parameters (240 and 160 configurations for M ∈ {2, 3, 4}
and M = 5 respectively) on each problem instance according to the same experimental procedure described in Section 3. The best configuration for eachparameter was determined according to the best median HV value at the 50 000
evaluations on each test problem. Due to space constraints, the best parameters
are not shown, but they are different depending on each control parameter andeach problem instance.

29

Figures 16, 17, and 18 show the influence of µ, g, and θ on the overall
performance of MOEA/D with the best configurations on the 13 test problems,
respectively. For each problem instance, for each optimization scenario, thebest configuration for each control parameter to be analyzed was used. The
performance rank of MOEA/D with the best parameters (Figures 16, 17, and18) is slightly different from that with the default parameters (Figures 2, 5, and
9). For example, unlike the results of MOEA/D with the default parameters(Figure 2), MOEA/D with µ = 400 performs the best on the two-objective
problems under the final population scenario (Figure 16). The interval which
MOEA/D with g pbi performs the best for M = 5 under the final population
scenario becomes shorter (see Figures 5(a) and 17(a)). Although θ = 3 is thesecond best θ value for M = 3 under the reduced UEA scenario (Figure 9),θ = 3 is the most suitable for MOEA/D (Figure 18).
However, a qualitative difference between the results of MOEA/D with thedefault and best parameters is not significant. Thus, the effect of the threecontrol parameters on MOEA/D is almost the same even if the best parameters
are used instead of the default parameters. This result guarantees that ouranalysis results in Section 4 can be generalized to a certain degree.
6. Conclusion
In this paper, we analyzed the influence of the three control parameters (thepopulation size µ, the scalarizing function g, and the penalty value θ of the
PBI function g pbi ) on the performance of MOEA/D under the two different
optimization scenarios (the final population and reduced UEA scenarios). Our
results on the four DTLZ and nine WFG problems up to five objectives showthat a suitable setting of the three control parameters (µ, g, and θ) of MOEA/D
is totally different depending on the problem and the optimization scenario. Our
observations are summarized as follows:• While MOEA/D with a small µ value exhibits poor performance measured
by the HV indicator on most of the 13 MOPs under the final population
scenario, it performs significantly better than MOEA/D with a large µ
value under the reduced UEA scenario.
• For the final population scenario, the PBI function g pbi with θ = 5 is anappropriate scalarizing function for some MOPs, and such a conclusion isconsistent with previous studies. However, for the reduced UEA scenario,
MOEA/D with the two Chebyshev functions (g chm and g chd ) perform
significantly better than that with g pbi with θ = 5.
• While well-distributed nondominated solutions are obtained by only usinga large θ value under the final population scenario, MOEA/D with a small
θ value (e.g., θ = 0.1) performs well on some MOPs under the reduced
UEA scenario.

30

M =2

4
2
0

0

10000 20000 30000 40000 50000

M =3

4

APS values

20010000 20000 30000 40000 50000

4

M =4

2
0

µ = 25µ = 50 2µ = 100µ = 200µ = 300 1µ = 400µ = 28µ = 55 3µ = 105 2µ = 210µ = 300 1µ = 406µ = 35 3µ = 56µ = 120 2µ = 220µ = 286 1µ = 455

10000 20000 30000 40000 50000

M =2

0

0

10000 20000 30000 40000 50000

M =40µ = 1262

10000 20000 30000 40000 50000

M =3

M =5
0

µ = 330µ = 455µ = 35µ = 56µ = 286µ = 220µ = 120µ = 126µ = 2101µ = 495µ = 495µ = 330
010000 20000 30000 40000 50000µ = 406µ = 300µ = 210µ = 105µ = 55µ = 28

10000 20000 30000 40000 50000

M =5µ = 2101µ = 400µ = 25µ = 300µ = 200µ = 50µ = 100

0

10000 20000 30000 40000 50000

Number of function evaluations

Number of function evaluations(a) Final Population scenario(b) Reduced UEA scenario

Figure 16: Overall performance of MOEA/D with various µ settings on the 13 problems withM ∈ {2, 3, 4, 5} (lower is better). The horizontal and vertical axes represent the number offunction evaluations and the APS values, respectively. The best configurations (g, θ, andT ) on each problem instance are used for each µ value (i.e., the control parameters used are
different for different problems).

We also analyzed the reason why an appropriate control parameter setting
is different for each scenario. Our analytical results can be briefly summarized
as follows:• MOEA/D with a particular parameter setting (e.g., a small µ value, g chm ,and g pbi with a small θ value) is capable of generating good nondominatedsolutions but cannot keep them in the population.
• Such an issue can be addressed by incorporating the UEA into MOEA/D,
which maintains well-distributed nondominated solutions independently
from the population.
• Therefore, an MOEA/D configuration that performs poorly for the final
population scenario is probable to show a good performance for the reduced UEA scenario.
Although in this paper we presented the analysis of the three control parameters (µ, g, and θ), latest MOEA/D-type algorithms (e.g., MOEA/D-DE [20])have other control parameters, such as the neighborhood size T , the number
31

2.0

M =2

1.5

M =2

1.5

g pbi

1.0

g chm 1.0

g chd

0.5

g chd 0.5

g chm

0.0

0.0

2.0

0

10000 20000 30000 40000 50000

M =3

1.5

APS values

g pbi 2.0

g pbi 2.0010000 20000 30000 40000 50000

M =3

1.5

g pbi

1.0

g chm 1.0

g chd

0.5

g chd 0.5

g chm

0.0

0.0

2.0

0

10000 20000 30000 40000 50000

M =4

1.5

g chm 2.0010000 20000 30000 40000 50000

M =4

1.5

g pbi

1.0

g pbi 1.0

g chd

0.5

g chd 0.5

g chm

0.0

0.0

2.0

0

10000 20000 30000 40000 50000

M =5

1.5

g chm 2.0010000 20000 30000 40000 50000

M =5

1.5

g pbi

1.0

g chd 1.0

g chd

0.5

g pbi 0.5

g chm

0.0

0.0

0

10000 20000 30000 40000 50000

010000 20000 30000 40000 50000

Number of function evaluations

Number of function evaluations(a) Final Population scenario(b) Reduced UEA scenario

Figure 17: Overall performance of MOEA/D with the three scalarizing functions (g chm , g chd ,and g pbi ) on the 13 problems with M ∈ {2, 3, 4, 5} (lower is better). The horizontal andvertical axes represent the number of function evaluations and the APS values, respectively.
The best configurations (µ, θ, and T ) on each problem instance are used for each scalarizing
function.

of replacement individuals nrep , and the probability of selecting types of neighborhood δ. An analysis of remaining control parameters of MOEA/D-type
algorithms is an avenue for future work.
Acknowledgement
This work was supported by the Science and Technology Innovation Committee Foundation of Shenzhen (Grant No. ZDSYS201703031748284).
References[1] Y. S. Ong, P. B. Nair, A. J. Keane, Evolutionary optimization of computationally expensive problems via surrogate modeling, AIAA journal 41 (4)(2003) 687–696.
[2] A. Ponsich, A. L. Jaimes, C. A. C. Coello, A Survey on Multiobjective
Evolutionary Algorithms for the Solution of the Portfolio Optimization
32

7

M =2

5
3
1
0

0

7

10000 20000 30000 40000 50000

M =3

5

APS values

3
1
0

0

7

10000 20000 30000 40000 50000

M =4

5
3
1
0

0

7

10000 20000 30000 40000 50000

M =5

5
3
1
0

0

10000 20000 30000 40000 50000

θ = 0.1
θ = 0.5
θ = 10
θ=8
θ=1
θ=5
θ=3
θ=2

7

θ = 0.5
θ = 0.1
θ=1
θ = 10
θ=8
θ=5
θ=2
θ=3

7

θ=1
θ = 0.5
θ = 0.1
θ = 10
θ=2
θ=8
θ=5
θ=3

7

θ = 0.5
θ=1
θ = 0.1
θ=2
θ=3
θ = 10
θ=8
θ=5

7

M =2

5
3
1
0

0

10000 20000 30000 40000 50000

M =3

5
3
1
0

0

10000 20000 30000 40000 50000

M =4

5
3
1
0

0

10000 20000 30000 40000 50000

M =5

5
3
1
0

0

10000 20000 30000 40000 50000

Number of function evaluations

Number of function evaluations(a) Final Population scenario(b) Reduced UEA scenario

θ = 0.5
θ = 0.1
θ = 10
θ=8
θ=1
θ=5
θ=2
θ=3
θ=1
θ = 0.5
θ = 0.1
θ=2
θ = 10
θ=8
θ=5
θ=3
θ=1
θ = 0.5
θ=2
θ = 0.1
θ=3
θ = 10
θ=5
θ=8
θ=1
θ = 0.5
θ=2
θ = 0.1
θ=3
θ = 10θ=8
θ=5Figure 18: Overall performance of MOEA/D using the PBI function g pbi with various θvalues on the 13 problems with M ∈ {2, 3, 4, 5}. The horizontal and vertical axes represent
the number of function evaluations and the APS values, respectively. The best configurations(µ and T ) on each problem instance are used for each θ value.

Problem and Other Finance and Economics Applications, IEEE TEVC
17 (3) (2013) 321–344.
[3] H. K. Singh, T. Ray, R. A. Sarker, Optimum Oil Production Planning
Using Infeasibility Driven Evolutionary Algorithm, Evol. Comput. 21 (1)(2013) 65–82.
[4] A. Trivedi, D. Srinivasan, K. Pal, C. Saha, T. Reindl, Enhanced Multiobjective Evolutionary Algorithm Based on Decomposition for Solving the Unit
Commitment Problem, IEEE Trans. Industrial Informatics 11 (6) (2015)
1346–1357.
[5] K. Deb, Multi-Objective Optimization Using Evolutionary Algorithms,
John Wiley & Sons, 2001.
[6] A. Zhou, B. Qu, H. Li, S. Zhao, P. N. Suganthan, Q. Zhang, Multiobjective
evolutionary algorithms: A survey of the state of the art, Swarm and Evol.
Comput. 1 (1) (2011) 32–49.

33

[7] B. Li, J. Li, K. Tang, X. Yao, Many-objective evolutionary algorithms: A
survey, ACM C. S. 48 (1) (2015) 13:1–13:35.
[8] K. Deb, S. Agrawal, A. Pratap, T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE TEVC 6 (2) (2002) 182–197.
[9] E. Zitzler, M. Laumanns, L. Thiele, SPEA2: Improving the Strength Pareto
Evolutionary Algorithm, Tech. Rep. 103, ETHZ (2001).
[10] K. Deb, M. Mohan, S. Mishra, Evaluating the epsilon-Domination BasedMulti-Objective Evolutionary Algorithm for a Quick Computation of
Pareto-Optimal Solutions, Evol. Comput. 13 (4) (2005) 501–525.
[11] E. Zitzler, S. Künzli, Indicator-based selection in multiobjective search, in:
PPSN, 2004, pp. 832–842.
[12] N. Beume, B. Naujoks, M. T. M. Emmerich, SMS-EMOA: multiobjective
selection based on dominated hypervolume, EJOR 181 (3) (2007) 1653–
1669.
[13] J. Bader, E. Zitzler, HypE: An Algorithm for Fast Hypervolume-Based
Many-Objective Optimization, Evol. Comput. 19 (1) (2011) 45–76.
[14] H. Liu, F. Gu, Q. Zhang, Decomposition of a Multiobjective Optimization Problem Into a Number of Simple Multiobjective Subproblems, IEEETEVC 18 (3) (2014) 450–455.
[15] H. Ishibuchi, T. Murata, A multi-objective genetic local search algorithm
and its application to flowshop scheduling, IEEE Trans. SMC, Part C 28 (3)(1998) 392–403.
[16] T. Murata, H. Ishibuchi, M. Gen, Specification of Genetic Search Directions
in Cellular Multi-objective Genetic Algorithms, in: EMO, 2001, pp. 82–95.
[17] E. J. Hughes, Evolutionary many-objective optimisation: many once or one
many?, in: IEEE CEC, 2005, pp. 222–227.
[18] Q. Zhang, H. Li, MOEA/D: A multiobjective evolutionary algorithm basedon decomposition, IEEE TEVC 11 (6) (2007) 712–731.
[19] A. Trivedi, D. Srinivasan, K. Sanyal, A. Ghosh, A Survey of Multiobjective
Evolutionary Algorithms Based on Decomposition, IEEE TEVC 21 (3)(2017) 440–462.
[20] H. Li, Q. Zhang, Multiobjective Optimization Problems With ComplicatedPareto Sets, MOEA/D and NSGA-II, IEEE TEVC 13 (2) (2009) 284–302.
[21] S. Jiang, S. Yang, An Improved Multiobjective Optimization Evolutionary Algorithm Based on Decomposition for Complex Pareto Fronts, IEEE
Trans. Cybernetics 46 (2) (2016) 421–437.

34

[22] S. Yang, S. Jiang, Y. Jiang, Improving the multiobjective evolutionary
algorithm based on decomposition with new penalty schemes, Soft Comput.
21 (16) (2017) 4677–4691.
[23] H. Liu, L. Chen, Q. Zhang, K. Deb, Adaptively Allocating Search Effort in
Challenging Many-Objective Optimization Problems, IEEE TEVC (2017)(in press).
[24] K. Li, K. Deb, Q. Zhang, S. Kwong, An Evolutionary Many-Objective
Optimization Algorithm Based on Dominance and Decomposition, IEEETEVC 19 (5) (2015) 694–716.
[25] Y. Yuan, H. Xu, B. Wang, B. Zhang, X. Yao, Balancing Convergenceand Diversity in Decomposition-Based Many-Objective Optimizers, IEEETEVC 20 (2) (2016) 180–198.
[26] A. E. Eiben, R. Hinterding, Z. Michalewicz, Parameter control in evolutionary algorithms, IEEE TEVC 3 (2) (1999) 124–141.
[27] F. G. Lobo, C. F. Lima, Z. Michalewicz (Eds.), Parameter setting in evolutionary algorithms, Studies in Computational Intelligence, Springer, 2007.
[28] H. Ishibuchi, Y. Sakane, N. Tsukamoto, Y. Nojima, Evolutionary ManyObjective Optimization by NSGA-II and MOEA/D with Large Populations, in: IEEE SMC, 2009, pp. 1758–1763.
[29] H. Ishibuchi, Y. Setoguchi, H. Masuda, Y. Nojima, How to compare manyobjective algorithms under different settings of population and archive sizes,
in: IEEE CEC, 2016, pp. 1149–1156.
[30] K. Deb, H. Jain, An evolutionary many-objective optimization algorithm
using reference-point-based nondominated sorting approach, part I: solving
problems with box constraints, IEEE TEVC 18 (4) (2014) 577–601.
[31] Y. Yuan, H. Xu, B. Wang, X. Yao, A New Dominance Relation-Based
Evolutionary Algorithm for Many-Objective Optimization, IEEE TEVC
20 (1) (2016) 16–37.
[32] H. Ishibuchi, Y. Sakane, N. Tsukamoto, Y. Nojima, Simultaneous use ofdifferent scalarizing functions in MOEA/D, in: GECCO, 2010, pp. 519–526.
[33] H. Ishibuchi, N. Akedo, Y. Nojima, Behavior of Multiobjective Evolutionary Algorithms on Many-Objective Knapsack Problems, IEEE TEVC
19 (2) (2015) 264–283.
[34] R. H. Gómez, C. A. C. Coello, A hyper-heuristic of scalarizing functions,in: GECCO, 2017, pp. 577–584.
[35] A. Mohammadi, M. N. Omidvar, X. Li, K. Deb, Sensitivity analysis
of Penalty-based Boundary Intersection on aggregation-based EMO algorithms, in: IEEE CEC, 2015, pp. 2891–2898.
35[36] H. Ishibuchi, K. Doi, Y. Nojima, Characteristics of many-objective test
problems and penalty parameter specification in MOEA/D, in: IEEE CEC,
2016, pp. 1115–1122.
[37] R. Tanabe, H. Ishibuchi, A. Oyama, Benchmarking Multi- and Manyobjective Evolutionary Algorithms under Two Optimization Scenarios,
IEEE Access 5 (2017) 19597–19619.
[38] K. Bringmann, T. Friedrich, P. Klitzke, Generic Postprocessing via Subset
Selection for Hypervolume and Epsilon-Indicator, in: PPSN, 2014, pp. 518–
527.
[39] T. Hanne, On the convergence of multiobjective evolutionary algorithms,
EJOR 117 (3) (1999) 553–564.
[40] J. E. Fieldsend, R. M. Everson, S. Singh, Using unconstrained elite archives
for multiobjective optimization, IEEE TEVC 7 (3) (2003) 305–323.
[41] R. Wang, Z. Zhou, H. Ishibuchi, T. Liao, T. Zhang, Localized Weighted
Sum Method for Many-Objective Optimization, IEEE TEVC 22 (1) (2018)
3–18.
[42] M. López-Ibáñez, J. D. Knowles, M. Laumanns, On Sequential Online
Archiving of Objective Vectors, in: EMO, 2011, pp. 46–60.
[43] A. Radulescu, M. López-Ibáñez, T. Stützle, Automatically Improving theAnytime Behaviour of Multiobjective Evolutionary Algorithms, in: EMO,
2013, pp. 825–840.
[44] D. Brockhoff, T. Tran, N. Hansen, Benchmarking Numerical Multiobjective
Optimizers Revisited, in: GECCO, 2015, pp. 639–646.
[45] T. Tusar, D. Brockhoff, N. Hansen, A. Auger, COCO: The Bi-objective
Black Box Optimization Benchmarking (bbob-biobj) Test Suite, CoRR
abs/1604.00359.
[46] H. Ishibuchi, Y. Sakane, N. Tsukamoto, Y. Nojima, Selecting a small number of representative non-dominated solutions by a hypervolume-based solution selection approach, in: FUZZ-IEEE, 2009, pp. 1609–1614.
[47] M. Basseur, B. Derbel, A. Goëffon, A. Liefooghe, Experiments on Greedy
and Local Search Heuristics for d dimensional Hypervolume Subset Selection, in: GECCO, 2016, pp. 541–548.
[48] A. P. Guerreiro, C. M. Fonseca, L. Paquete, Greedy Hypervolume Subset
Selection in Low Dimensions, Evol. Comput. 24 (3) (2016) 521–544.
[49] Q. Zhang, W. Liu, H. Li, The performance of a new version of MOEA/Don CEC09 unconstrained MOP test instances, in: IEEE CEC, 2009, pp.
203–208.

36

[50] R. Wang, J. Xiong, H. Ishibuchi, G. Wu, T. Zhang, On the effect of reference point in MOEA/D for multi-objective optimization, Appl. Soft Comput. 58 (2017) 25–34.
[51] A. J. Nebro, J. J. Durillo, C. A. C. Coello, F. Luna, E. Alba, A Study ofConvergence Speed in Multi-objective Metaheuristics, in: PPSN, 2008, pp.
763–772.
[52] Y. Jin, Surrogate-assisted evolutionary computation: Recent advances andfuture challenges, Swarm and Evol. Comput. 1 (2) (2011) 61–70.
[53] H. Li, M. Ding, J. Deng, Q. Zhang, On the use of random weights in
MOEA/D, in: IEEE CEC, 2015, pp. 978–985.
[54] M. Li, S. Yang, X. Liu, Pareto or Non-Pareto: Bi-Criterion Evolution in
Multiobjective Optimization, IEEE TEVC 20 (5) (2016) 645–665.
[55] S. Wessing, N. Beume, G. Rudolph, B. Naujoks, Parameter Tuning Boosts
Performance of Variation Operators in Multiobjective Optimization, in:
PPSN, 2010, pp. 728–737.
[56] L. C. T. Bezerra, M. López-Ibáñez, T. Stützle, A Large-Scale Experimental Evaluation of High-Performing Multi-and Many-Objective Evolutionary
Algorithms, Tech. Rep. 2017-005, IRIDIA (2017).
[57] M. Andersson, S. Bandaru, A. H. C. Ng, A. Syberfeldt, Parameter Tuningof MOEAs Using a Bilevel Optimization Approach, in: EMO, 2015, pp.
233–247.
[58] K. Deb, L. Thiele, M. Laumanns, E. Zitzler, Scalable Test Problems for
Evolutionary Multi-Objective Optimization, in: Evolutionary Multiobjective Optimization. Theoretical Advances and Applications, Springer, 2005,
pp. 105–145.
[59] S. Huband, P. Hingston, L. Barone, R. L. While, A review of multiobjective
test problems and a scalable test problem toolkit, IEEE TEVC 10 (5) (2006)
477–506.
[60] I. Das, J. E. Dennis, Normal-Boundary Intersection: A New Method forGenerating the Pareto Surface in Nonlinear Multicriteria Optimization
Problems, SIAM J. Optimiz. 8 (3) (1998) 631–657.
[61] M. Pescador-Rojas, R. H. Gómez, E. Montero, N. Rojas-Morales, M. C.
Riff, C. A. C. Coello, An Overview of Weighted and Unconstrained Scalarizing Functions, in: EMO, 2017, pp. 499–513.
[62] K. Li, Q. Zhang, S. Kwong, M. Li, R. Wang, Stable Matching-Based Selection in Evolutionary Multiobjective Optimization, IEEE TEVC 18 (6)(2014) 909–923.

37

[63] Y. Qi, X. Ma, F. Liu, L. Jiao, J. Sun, J. Wu, MOEA/D with Adaptive
Weight Adjustment, Evol. Comput. 22 (2) (2014) 231–264.
[64] H. Ishibuchi, H. Masuda, Y. Nojima, Pareto Fronts of Many-Objective
Degenerate Test Problems, IEEE TEVC 20 (5) (2016) 807–813.
[65] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, V. G. da Fonseca, Performance assessment of multiobjective optimizers: an analysis and review,
IEEE TEVC 7 (2) (2003) 117–132.
[66] H. Ishibuchi, Y. Setoguchi, H. Masuda, Y. Nojima, Performance ofDecomposition-Based Many-Objective Algorithms Strongly Depends on
Pareto Front Shapes, IEEE TEVC 21 (2) (2017) 169–190.
[67] H. Ishibuchi, H. Masuda, Y. Nojima, Comparing solution sets of different
size in evolutionary many-objective optimization, in: IEEE CEC, 2015, pp.
2859–2866.
[68] M. Asafuddoula, H. K. Singh, T. Ray, An Enhanced Decomposition-Based
Evolutionary Algorithm With Adaptive Reference Vectors, IEEE Trans.
Cybernetics (2017) (in press).
[69] R. H. Gómez, C. A. C. Coello, Improved Metaheuristic Based on the R2
Indicator for Many-Objective Optimization, in: GECCO, 2015, pp. 679–686.
[70] H. Sato, Analysis of inverted PBI and comparison with other scalarizing
functions in decomposition based MOEAs, J. Heuristics 21 (6) (2015) 819–
849.
[71] H. Ishibuchi, R. Imada, Y. Setoguchi, Y. Nojima, Reference point specification in hypervolume calculation for fair comparison and efficient search,
in: GECCO, 2017, pp. 585–592.

38