# The Cone $\epsilon$-Dominance: An Approach for Evolutionary Multiobjective
Optimization

Lucas de Souza Batista<EMAIL_ADDRESS>
Universidade Federal de Minas Gerais, Departamento de Engenharia Elétrica, Av.
Antônio Carlos 6627, 31270-010, Belo Horizonte, MG, Brasil Felipe Campelo
<EMAIL_ADDRESS>
Universidade Federal de Minas Gerais, Departamento de Engenharia Elétrica, Av.
Antônio Carlos 6627, 31270-010, Belo Horizonte, MG, Brasil Frederico Gadelha
Guimarães<EMAIL_ADDRESS>
Universidade Federal de Minas Gerais, Departamento de Engenharia Elétrica, Av.
Antônio Carlos 6627, 31270-010, Belo Horizonte, MG, Brasil Jaime Arturo
Ramírez<EMAIL_ADDRESS>
Universidade Federal de Minas Gerais, Departamento de Engenharia Elétrica, Av.
Antônio Carlos 6627, 31270-010, Belo Horizonte, MG, Brasil

###### Abstract

We propose the cone$\epsilon$-dominance approach to improve convergence and
diversity in multiobjective evolutionary algorithms (MOEAs). A
cone$\epsilon$-MOEA is presented and compared with MOEAs based on the standard
Pareto relation (NSGA-II, NSGA-II*, SPEA2, and a clustered NSGA-II) and on the
$\epsilon$-dominance ($\epsilon$-MOEA). The comparison is performed both in
terms of computational complexity and on four performance indicators selected
to quantify the quality of the final results obtained by each algorithm: the
convergence, diversity, hypervolume, and coverage of many sets metrics.
Sixteen well-known benchmark problems are considered in the experimental
section, including the ZDT and the DTLZ families. To evaluate the possible
differences amongst the algorithms, a carefully designed experiment is
performed for the four performance metrics. The results obtained suggest that
the cone$\epsilon$-MOEA is capable of presenting an efficient and balanced
performance over all the performance metrics considered. These results
strongly support the conclusion that the cone$\epsilon$-MOEA is a competitive
approach for obtaining an efficient balance between convergence and diversity
to the Pareto front, and as such represents a useful tool for the solution of
multiobjective optimization problems.

Keywords

> Multiobjective optimization, evolutionary algorithms, genetic algorithms,
> epsilon-dominance, cone epsilon-dominance.

##  1 Introduction

Applied design problems are usually expressed as optimization problems that
involve multiple and often conflicting objectives. The goal of solving these
multiobjective optimization problems (MOPs) is to find a set of trade-off
solutions, known as Pareto-optimal or nondominated solutions, in which each
solution represents different compromises for the objectives. In this context,
multiobjective evolutionary algorithms (MOEAs) have been widely recognized as
suitable approaches to solve MOPs, mainly due to their robustness and
generality. In general, researchers have found that MOEAs are more practical
compared to the preference-based classical approaches, since the user gets an
opportunity to analyze a range of other trade-off solutions before choosing
the preferable one. Besides, the MOEA search procedure is also algorithmically
efficient as the population evolves simultaneously towards many different
regions of the Pareto-optimal front (Zitzler et al.,, 2000; Coello Coello et
al.,, 2002; Deb et al., 2005a, ).

In the design of a MOEA, the algorithm is often required to find a well-
distributed set of solutions near to the global Pareto front, thus enabling
the decision-maker to choose the most suitable design by considering some
higher-level problem information. Essentially, a high-quality approximation
set should approach the true Pareto front as close as possible, and be well-
spread along its extension (Zitzler et al.,, 2001). This fact implies that the
evaluation of the quality of Pareto approximated sets achieved by evolutionary
multiobjective algorithms represents itself a multi-criteria problem, in which
convergence and diversity have to be considered simultaneously. To consider
this compromise, most MOEAs incorporate an external archive in which the
nondominated solutions estimated during the search process are stored. The
external or memory archives impose an elitist mechanism to MOEAs, and only the
solutions that are globally nondominated with respect to all the solutions
generated so far by the MOEA are maintained. As the external archive size is
usually finite, some truncation techniques have been developed in order to
enforce a good distribution of solutions into the archive, in order to obtain
a good characterization of the global Pareto front. In general, a Pareto-based
fitness assignment method is usually designed in order to guide the search
toward the global Pareto-optimal front, whereas the population diversity is
commonly promoted by employing density estimation methods such as clustering
(Zitzler and Thiele,, 1999), crowding (Deb et al.,, 2002), adaptive grids
(Knowles and Corne,, 2000), entropy (Farhang-Mehr and Azarm,, 2002), or
relaxed forms of Pareto dominance (Ikeda et al.,, 2001; Laumanns et al.,,
2002; Hernández-Díaz et al.,, 2007; Sato et al.,, 2007; Batista et al., 2011b,
). To help contextualize the scope of this work, some of these techniques are
reviewed throughout the text.

In the current paper, we propose a relaxed form of dominance called
cone$\epsilon$-dominance that is incorporated into a steady-state MOEA. The
cone$\epsilon$-MOEA is contrasted to five well known approaches, namely the
NSGA-II (two versions of this method are considered) (Deb et al.,, 2002), the
$\epsilon$-MOEA (Deb et al.,, 2003; Deb et al., 2005a, ), the SPEA2 (Zitzler
et al.,, 2001) and the clustered-NSGA-II (Deb et al.,, 2003), both in terms of
computational complexity and on a number of performance metrics selected to
quantify the convergence and diversity preservation abilities of the
algorithms. After analyzing the results obtained, both in terms of statistical
significance and magnitude of effects, it becomes clear that the proposed
cone$\epsilon$-MOEA is capable of achieving an efficient balance between
convergence and diversity, and as such represents an interesting and
competitive approach to the solution of MOPs.

The paper is organized as follows: Section 2 reviews the basic definitions of
the multiobjective optimization problem and the usual dominance criterion.
Section 3 discusses some of the most popular relaxed dominance criteria, and
Section 4 reports some practical diversity preservation methods employed in
MOEAs. Section 5 contains a detailed description and mathematical definition
of the cone$\epsilon$-dominance approach. Section 6 describes the MOEAs, the
test problems, the performance metrics, and the statistical design employed in
the evaluation of the proposed strategy. Section 7 presents a comprehensive
analysis of the performance of the proposed approach and the other algorithms
compared. Finally, the conclusions are outlined in Section 8.

##  2 Mathematical Background

###  2.1 Multiobjective Optimization Problem

Multiobjective optimization problems can be stated, without loss of
generality, as:

$\begin{split}\min\limits_{\boldsymbol{x}}&\ \
\boldsymbol{f}\left(\boldsymbol{x}\right)=\left[f_{1}(\boldsymbol{x}),\ldots,f_{m}(\boldsymbol{x})\right]^{T}\\\
\makebox{subject to: }&\begin{cases}&g_{i}\left(\boldsymbol{x}\right)\leq 0,\
i=1,\ldots,n_{g}\\\ &h_{j}\left(\boldsymbol{x}\right)=0,\ j=1,\ldots,n_{h}\\\
&\boldsymbol{x}\in\mathbb{X}\end{cases}\end{split}$ (1)

where $\boldsymbol{x}\in\mathbb{X}$ is the vector of optimization variables;
$\mathbb{X}\subset\mathbb{R}^{n}$ is the optimization domain, defined by the
Cartesian product of the domains of each optimization variable;
$\boldsymbol{f}\left(\cdot\right):\mathbb{X}\mapsto\mathbb{R}^{m}$ are the
objective functions of the problem;
$\boldsymbol{g}\left(\cdot\right):\mathbb{X}\mapsto\mathbb{R}^{n_{g}}$ and
$\boldsymbol{h}\left(\cdot\right):\mathbb{X}\mapsto\mathbb{R}^{n_{h}}$
represent the inequality and equality constraints of the problem,
respectively; and the set of feasible solutions is represented by
$\Omega\subseteq\mathbb{X}$. In this context, the goal of multiobjective
evolutionary algorithms is to obtain a diverse set of estimates of the Pareto
optimal set, which contains the nondominated solutions of the multiobjective
problem.

###  2.2 The Standard Dominance Relation

Pareto dominance (Deb,, 2001) has been the most commonly adopted criterion
used to discriminate among solutions in the multiobjective context, and
therefore it has been the basis to develop most of the MOEAs proposed so far,
e.g., SPEA (Zitzler and Thiele,, 1999) and NSGA-II (Deb et al.,, 2002).

By definition, a feasible solution $\boldsymbol{x}\in\Omega$ Pareto dominates
another point $\boldsymbol{x}^{\prime}\in\Omega$ if the following relation
holds:

$\boldsymbol{f}(\boldsymbol{x})\leq\boldsymbol{f}(\boldsymbol{x}^{\prime})\mbox{
and
}\boldsymbol{f}(\boldsymbol{x})\neq\boldsymbol{f}(\boldsymbol{x}^{\prime})$
(2)

in which the relation operators $\leq$ and $\neq$ are defined as:

$\displaystyle\boldsymbol{f}(\boldsymbol{a})$
$\displaystyle\leq\boldsymbol{f}(\boldsymbol{b})\Leftrightarrow
f_{i}(\boldsymbol{a})\leq
f_{i}(\boldsymbol{b}),~{}~{}~{}\forall~{}i\in\left\\{1,\ldots,m\right\\}$ (3)
$\displaystyle\boldsymbol{f}(\boldsymbol{a})$
$\displaystyle\neq\boldsymbol{f}(\boldsymbol{b})\Leftrightarrow\exists~{}i\in\left\\{1,\ldots,m\right\\}:f_{i}(\boldsymbol{a})\neq
f_{i}(\boldsymbol{b})$ (4)

in which $\boldsymbol{a}$ and $\boldsymbol{b}$ represent two different
decision vectors. In short, this dominance is usually expressed as
$\boldsymbol{f}(\boldsymbol{x})\prec\boldsymbol{f}(\boldsymbol{x}^{\prime})$.

All solutions that are not dominated by any other decision vector of a given
set are called nondominated regarding this set. In this way, the Pareto
optimal set $\mathcal{P}$ is defined as the set of nondominated solutions:

$\mathcal{P}=\left\\{\boldsymbol{x}^{*}\in\Omega\mid\nexists~{}\boldsymbol{x}\in\Omega:\boldsymbol{f}(\boldsymbol{x})\prec\boldsymbol{f}(\boldsymbol{x}^{*})\right\\}~{}.$
(5)

The image of this set in the objective space is called the Pareto front
$\mathcal{F}=\boldsymbol{f}(\mathcal{P})$, i.e.:

$\mathcal{F}=\left\\{\boldsymbol{y}=\boldsymbol{f}(\boldsymbol{x}):\boldsymbol{x}\in\mathcal{P}\right\\}~{}.$
(6)

##  3 Relaxed Dominance Criteria

A few years ago, Ikeda et al., (2001) proposed a relaxed dominance criterion
to deal with dominance resistant solutions, i.e., solutions that are extremely
inferior to others in at least one objective, but hardly-dominated. The idea
behind this approach, called $\alpha$-dominance, is to set lower and upper
bounds of trade-off rates between two objectives, such that solutions
characterized by a small improvement in some objectives (compared to the
amount of detriment in one objective), which would be nondominated according
to the standard Pareto dominance, would be rejected. Since there is no
explicit formula for the calculation of $\alpha$, its choice is usually left
to the designer, representing a difficult and problem-dependent task.
Furthermore, even though the selective pressure induced by the
$\alpha$-approach improves the convergence of a MOEA, there is no guarantee of
finding a representative well-spread estimation of the real Pareto front.

Following Ikeda’s work, Laumanns et al., (2002) proposed two relaxed dominance
methods: the additive and the multiplicative $\epsilon$-dominance schemes.
These mechanisms act as an archiving strategy to ensure both properties of
convergence towards the Pareto-optimal front and diversity among the solutions
found. These archiving approaches essentially promote the convergence to a
representative well distributed approximation of the global Pareto front,
always preserving the best solutions found at intermediate iterations of the
algorithm. These techniques guarantee that no two archived solutions can share
a given $\epsilon_{i}$ neighborhood on the $i^{th}$ objective (with
$i=1,\ldots,m$, $m$ being the number of objectives), with the $\epsilon$
values usually provided by the designer to control the size (resolution) of
the solution set. However, useful $\epsilon$ values are normally not known
before executing a MOEA, since the equations used to estimate the values of
$\epsilon$ are only valid for the case of linear Pareto fronts. This leads to
difficulties in computing appropriate values of $\epsilon$ to provide the
desired number of nondominated points, since the geometrical features of the
Pareto-optimal front are commonly unknown by the designer, and the
$\epsilon$-dominance strategy can lose a high number of viable solutions when
the $\epsilon$ values are badly estimated. Moreover, this approach tends to
neglect viable solutions since it does not allow two points with a difference
of $\epsilon_{i}$ in the $i^{th}$ objective to be mutually nondominated.
Because of this property, it is usually not possible to obtain solutions at
the corners of the estimated Pareto front, contributing negatively to the
spread of solutions along its extension (Hernández-Díaz et al.,, 2007). The
$\epsilon$-dominance concept has been used and adapted in a number of studies,
e.g. to improve the convergence characteristics of MOEAs (Zhao and Suganthan,,
2010; Li et al.,, 2011; Hernández-Díaz et al.,, 2011), to address specific
problems (Aguirre and Tanaka,, 2009), or to tackle many-objective optimization
problems (Pasia et al.,, 2011).

An alternative $\epsilon$-dominance strategy to overcome some of the
limitations of $\epsilon$-dominance was later proposed by Hernández-Díaz et
al., (2007). This technique, called Pareto adaptive $\epsilon$-dominance
(pa$\epsilon$-dominance for short), considers not only a different $\epsilon$
value for each objective, but a vector
$\boldsymbol{\epsilon}=\left(\boldsymbol{\epsilon}_{1},\ldots,\boldsymbol{\epsilon}_{m}\right)$
associated to the objectives $\boldsymbol{f}=\left(f_{1},\ldots,f_{m}\right)$
depending on the geometrical characteristics of the Pareto-optimal front. In
this way, different intensities of dominance are considered for each objective
according to the position of each solution along the Pareto front. The size of
the boxes, defined by the $\epsilon$ values, is adapted depending on the
corresponding position in the objective space, so that smaller boxes are
created where needed (e.g., at the extremes of the Pareto front), and larger
ones in other less problematic parts of the front. Unfortunately, this
approach also has some drawbacks. There are problems in which
pa$\epsilon$-dominance is not able to maintain a good distribution of
solutions at the extreme parts of the Pareto front. Moreover, the adaptive
vector $\boldsymbol{\epsilon}$ depends on the quality of the front used for
its estimation (Hernández-Díaz et al.,, 2007). In fact, this affects the
spread of solutions along the front and also has a negative effect on the
performance of the method.

Sato et al., (2007) proposed to control the dominance area of solutions in
order to induce appropriate selective pressure and ranking (in the objective
domain) in MOEAs. This strategy controls the degree of contraction or
expansion of the dominance area of solutions by modifying the fitness value of
each objective function, which is attained by changing a user-defined
parameter $S_{i}$, for all $i=1,\dots,m$. As shown in (Sato et al.,, 2007),
when $S_{i}<0.5$, the $i^{th}$ fitness value $f_{i}(\boldsymbol{x})$ is
increased to $f_{i}^{\prime}(\boldsymbol{x})>f_{i}(\boldsymbol{x})$. On the
other hand, when $S_{i}>0.5$, $f_{i}(\boldsymbol{x})$ is decreased to
$f_{i}^{\prime}(\boldsymbol{x})<f_{i}(\boldsymbol{x})$. At last, when
$S_{i}=0.5$, $f_{i}^{\prime}(\boldsymbol{x})=f_{i}(\boldsymbol{x})$, which is
equivalent to conventional dominance. However, since different rankings can be
produced, the optimum parameter $\boldsymbol{S}^{*}$ that yields maximum
search performance depends strongly on the problem. In addition, either
convergence or diversity can be emphasized by contracting or expanding the
dominance area, but not always both simultaneously (Sato et al.,, 2007).

To address some of the limitations discussed above, a relaxation of the strict
dominance concept based on an extension of the $\epsilon$-dominance criterion
and on the use of cones to control the dominance region of solutions was
suggested in (Batista et al., 2011b, ). The main idea of this relaxation was
to maintain the good convergence properties of $\epsilon$-dominance while
simultaneously improving the control over the diversity and resolution of the
estimated Pareto front, providing a dominance relation that is less sensitive
to the geometrical features of the front. As seen before, some relaxed forms
of dominance are sensitive to the loss of efficient solutions, mainly those
located in certain portions of the front, e.g., segments with degenerated
tradeoffs (small gains for one objective at the expense of large losses in
another) and the extremes of the front.

This loss of potentially interesting solutions is an undesired characteristic
in many cases, since it usually leads to degradation in the diversity of the
fronts obtained. In fact, the design of effective mechanisms to maintain
diversity remains as a key issue in a number of cases, for instance, when
extending particle swarm optimizers to solve multiobjective optimization
problems (Villalobos-Arias et al.,, 2005). Additionally, better diversity in
the sampling of the Pareto front might be useful for (i) some interactive
decision-making tools based on reconstructing the Pareto front (Chen et al.,,
2002; Chankong and Haimes,, 2008); (ii) the search for design rules and
principles, such as in the innovization concept advanced by Deb and
Srinivasan, 2006a ; or (iii) a posteriori analysis of the Pareto front to
identify the variables that are relevant for the tradeoffs (Deb and
Srinivasan, 2006b, ). As shown in (Batista et al., 2011b, ), since the
influence on the ordering of points performed by the cone of dominance can be
limited to a local neighborhood in the objective space, the
cone$\epsilon$-criterion enables the approximation of nondominated points in
some adjacent boxes that would otherwise be $\epsilon$-dominated, including
the extreme parts of the front, which contributes positively for the
performance of the method. In this work we expand and improve the concepts
initially introduced in that preliminary work, providing greater insight into
its characteristics and addressing some of the questions regarding the
diversity promoting properties of the cone$\epsilon$-dominance.

##  4 Diversity Preservation into MOEA’s Memory Archive

As indicated before, some techniques have been suggested over the years in
order to approach and maintain a good distribution of solutions in the space
of objectives. The ones that are used in this work are discussed in this
section.

###  4.1 Clustering

In optimization problems defined in continuous search domains, the Pareto-
optimal set can be extremely large or even infinite. Nevertheless, since
approximating a large set of nondominated solutions is useless from the
decision maker’s point of view, pruning the external archive while maintaining
some desirable characteristics is necessary or even mandatory. For instance,
the SPEA proposed by Zitzler and Thiele, (1999) uses a clustering technique
named _average linkage method_. This approach is employed to prune the
contents of the bounded memory archive of SPEA. In a general way, (i) each
external nondominated solution compose initially a distinct cluster; (ii) the
two closest ones amalgamate into a larger one; (iii) the second step is
repeated until the desired number is achieved; finally, (iv) the reduced
nondominated set is obtained by selecting a representative solution per
cluster.

Another example concerns to the clustering technique used in SPEA2 (Zitzler et
al.,, 2001). This mechanism presents similarities with the truncation method
of SPEA, however it does not loose boundary points. In short, when the
nondominated front exceeds the archive limit, the point which has the smallest
distance to the $k$ nearest ones in the archive is deleted, at each stage,
thus reducing the memory archive to its desired size. Also, if there are
several solutions with minimum distance, the tie is broken by considering the
second smallest distance and so forth.

Some different clustering techniques have also been adopted in MOEAs for
maintaining diversity in their memory archives and even in decision variable
space, e.g., (Pulido and Coello Coello,, 2004; Janson and Merkle,, 2005;
Padhye et al.,, 2009).

###  4.2 Crowding Distance

The crowding distance mechanism was suggested by Deb et al., (2002) to perform
diversity preservation in the NSGA-II. This algorithm uses the Pareto
dominance to rank the solutions into different fronts (rank based on dominance
depth). In each front, the crowding distance of a solution is estimated by
calculating the sum of the Euclidean distances among the two neighboring
solutions from either side of the solution along the objectives. As
illustrated in Fig. 1(a), the crowding distance of a solution $i$ is given by
the sum $f_{1}^{(i)}+f_{2}^{(i)}$.

Variations of the crowding comparison operator have been used by several
researchers to maintain diversity in a memory archive, e.g., (Raquel and
Naval,, 2005; Wang et al.,, 2010).

###  4.3 Pareto $\epsilon$-Dominance

Since the performance of the multiplicative $\epsilon$-approach proposed by
Laumanns et al., (2002) tends to be more sensitive to the geometrical
characteristics of a frontier (Hernández-Díaz et al.,, 2007), only the
additive scheme will be discussed hereinafter.

Formally, supposing that all objectives $f_{i}$,
$i\in\left\\{1,\ldots,m\right\\}$, are to be minimized, and also that $1\leq
f_{i}\leq K$, for all $i$, then, given a vector
$\boldsymbol{y}\in\mathbb{R}^{m}$ and $\epsilon>0$, $\boldsymbol{y}$ is said
to $\epsilon$-dominate $\boldsymbol{y}^{\prime}\in\mathbb{R}^{m}$, denoted as
$\boldsymbol{y}\stackrel{{\scriptstyle\epsilon}}{{\prec}}\boldsymbol{y}^{\prime}$,
if and only if,

$y_{i}-\epsilon\leq y_{i}^{\prime},\mbox{ for all
}i\in\left\\{1,\ldots,m\right\\}~{}.$ (7)

As introduced by Laumanns et al., (2002), a two-level selection mechanism is
implemented in the $\epsilon$-dominance approach. First, it creates a
hypergrid in the objective space where each box uniquely contains one solution
vector. Basically, a box-level dominance relation is used, so that the
algorithm always maintains a set of nondominated boxes, thus guaranteeing the
diversity property. Second, if two vectors share the same box, the usual
Pareto dominance relation is applied, so that the best one is selected and
convergence is guaranteed. However, if none of these two vectors dominates the
other, it is usual to keep the point closest to the origin of the box, i.e.,
to the corner where all objectives would have the lowest values within that
box (Deb et al.,, 2003; Deb et al., 2005a, ).

The $\epsilon$-dominance mechanism generates a hypergrid in the objective
space with $\left(\left(K-1\right)/\epsilon\right)^{m}$ boxes which
accommodate a maximum of $\left(\left(K-1\right)/\epsilon\right)^{m-1}$ non
$\epsilon$-dominated points. Then, supposing that the designer wants a maximum
of $T$ non $\epsilon$-dominated points in the archive, the $\epsilon$ value
can be easily calculated as:

$\epsilon=\frac{K-1}{T^{1/(m-1)}}~{}\cdot$ (8)

It should be noted, however, that the $\epsilon$-dominance strategy is only
able to obtain this number $T$ in cases where the Pareto-front is linear. For
other cases (see Fig. 1(b)) this value is merely an upper limit, with the
actual number of nondominated points found being much smaller (Hernández-Díaz
et al.,, 2007). Also, it is interesting to note that the given definitions can
be generalized by considering a different $\epsilon$ value for each objective.

(a) Crowding distance assignment.

(b) Diversity access by $\epsilon$-dominance.

Figure 1: (a) Crowding distance assignment used by the NSGA-II. For example,
the crowding measure for the solution $i$ is given by
$f_{1}^{(i)}+f_{2}^{(i)}$. (b) Distribution achieved by $\epsilon$-dominance.
It is normal to lose the extreme points of the Pareto front, as well as points
located in segments of the front that are almost parallel to the objective
axes. As shown, the points (⚫) are $\epsilon$-Pareto, and the points (⚪),
although Pareto-optimal, are $\epsilon$-dominated. This illustrates the fact
that a high number of nondominated solutions can be lost if the decison maker
does not take into account, or does not know beforehand, the geometric
features of the true Pareto front of the problem to be solved.

##  5 The Pareto Cone $\epsilon$-Dominance Strategy

Before we present the formal definition of the cone$\epsilon$-approach, a
conceptual interpretation is given (Batista et al., 2011b, ). Figure 2
contrasts both the $\epsilon$-dominance and the cone$\epsilon$-dominance
strategies, and emphasizes the grid generated in the objective space and the
different areas dominated by a solution $\boldsymbol{y}$. In order to get a
nondominated solution set, the cone$\epsilon$-dominance mechanism entails both
the area dominated by the cone and by the standard Pareto dominance, i.e., the
shaded region in Fig. 2(b). Indeed, the hypervolume dominated by
$\boldsymbol{y}$ using the cone$\epsilon$-dominance approach represents a
relaxation of that dominated by $\boldsymbol{y}$ when using the usual
dominance. As illustrated, the proposed relaxation enables the approximation
of nondominated points in some adjacent boxes that would be
$\epsilon$-dominated. Essentially, since the influence on the ordering of
points performed by the cone$\epsilon$-dominance is limited to a local
neighborhood in the objective space, its effects are less dramatic than those
of the $\epsilon$-dominance, which makes it possible to obtain a better
distribution of solutions.

Note that the cone$\epsilon$-approach differs in some aspects from all the
mechanisms reviewed in this paper: (i) the selection pressure is ensured by
using both the Pareto criterion and the cone of dominance; (ii) the
convergence and diversity preservation of the solutions found are
simultaneously guaranteed by means of the hypergrid adopted; and (iii) a well-
spread nondominated archive is expected due to the local effect of the cone in
the objective domain. This last feature is illustrated in the Fig. 3.

The cone$\epsilon$-dominance can also be seen as a hybrid between
$\epsilon$-dominance and the proper efficiency with respect to cones discussed
in (Miettinen,, 1998). The use of cones into the MOEA’s structure can also be
found in (Branke et al.,, 2001; Shukla et al.,, 2010; Hirsch et al.,, 2011).

(a) $\epsilon$-dominance.

(b) cone$\epsilon$-dominance.

Figure 2: Illustration of the region dominated by a solution $\boldsymbol{y}$
regarding the use of the $\epsilon$ and cone$\epsilon$ concepts for a two-
objective minimization problem.

(a) Diversity by cone$\epsilon$-dominance.

(b) Diversity by $\epsilon$-dominance.

Figure 3: (a) Distribution performed by the cone$\epsilon$-approach in a
connected Pareto front. (b) The $\epsilon$-relation is very sensitive to the
geometry of the frontier and, therefore, inefficient to estimate
representative solutions of the Pareto set.

###  5.1 Basic Definitions

In this section we present the basic definitions needed for the formal
presentation of the cone$\epsilon$-approach (Batista et al., 2011b, ).

Definition 1. (Cone) A set $\mathcal{C}$ is a cone if
$\lambda\boldsymbol{y}\in C$ for any $\boldsymbol{y}\in C$ and
$\forall~{}\lambda\geq 0$. $\Box$

Definition 2. (Generated cone) For two vectors $\boldsymbol{w}_{1}$ and
$\boldsymbol{w}_{2}$, the cone generated by these vectors is the set
$\mathcal{C}=\left\\{\boldsymbol{z}:\boldsymbol{z}=\lambda_{1}\boldsymbol{w}_{1}+\lambda_{2}\boldsymbol{w}_{2},~{}\forall~{}\lambda_{1},\lambda_{2}\geq
0\right\\}$. $\Box$

Extending the generated cone concept to $m$ dimensions, we have that the
hypercone generated by the vectors $\boldsymbol{w}_{i}$, for all
$i\in\left\\{1,\ldots,m\right\\}$, is the set
$\mathcal{C}=\left\\{\boldsymbol{z}:\boldsymbol{z}=\lambda_{1}\boldsymbol{w}_{1}+\ldots+\lambda_{i}\boldsymbol{w}_{i}+\ldots+\lambda_{m}\boldsymbol{w}_{m},~{}\forall~{}\lambda_{i}\geq
0\right\\}$.

The definitions above form the basis for a mechanism to control the
hypervolume dominated by a specific cone $\mathcal{C}$. For the 2D case (see
Fig. 4), it is easy to see that, with respect to the origin of the box,
$\boldsymbol{w}_{1}=\left[\epsilon_{1}~{}~{}\kappa\epsilon_{2}\right]^{T}$ and
$\boldsymbol{w}_{2}=\left[\kappa\epsilon_{1}~{}~{}\epsilon_{2}\right]^{T}$.
The cone $\mathcal{C}$ can, therefore, be rewritten as:

$\mathcal{C}=\Bigg{\\{}\boldsymbol{z}:\overbrace{\left[\begin{array}[]{c}z_{1}\\\
z_{2}\end{array}\right]}^{\boldsymbol{z}}=\overbrace{\left[\begin{array}[]{cc}\epsilon_{1}&\kappa\epsilon_{1}\\\
\kappa\epsilon_{2}&\epsilon_{2}\end{array}\right]}^{\Psi}\overbrace{\left[\begin{array}[]{c}\lambda_{1}\\\
\lambda_{2}\end{array}\right]}^{\boldsymbol{\lambda}},~{}\forall~{}\lambda_{1},\lambda_{2}\geq
0\Bigg{\\}}$ (9)

in which the parameter $\kappa\in[0,1)$ controls the opening of the cone
$\mathcal{C}$, and $\Psi$ is the cone-dominance matrix, which in fact controls
the hypervolume dominated by $\mathcal{C}$. Notice that the
cone$\epsilon$-dominance strategy tends toward the $\epsilon$-dominance
strategy when $\kappa\rightarrow 0$, and degenerates into the usual Pareto
dominance for $\kappa=1$.

Figure 4: Mechanism used in the cone$\epsilon$-dominance approach to control
the hypervolume dominated by a specific cone.

For the general case with $m$-objectives we have:

$\Psi(\epsilon_{i},\kappa)=\left[\begin{array}[]{ccccc}\epsilon_{1}&\kappa\epsilon_{1}&\ldots&\kappa\epsilon_{1}\\\
\kappa\epsilon_{2}&\epsilon_{2}&\ldots&\kappa\epsilon_{2}\\\
\vdots&\vdots&\ddots&\vdots\\\
\kappa\epsilon_{m}&\kappa\epsilon_{m}&\ldots&\epsilon_{m}\end{array}\right]~{}\cdot$
(10)

which can be used to define the cone$\epsilon$-dominance strategy.

Definition 3. (Cone $\epsilon$-dominance) Given two vectors $\boldsymbol{y}$,
$\boldsymbol{y}^{\prime}\in\mathbb{R}^{m}$, $\boldsymbol{y}$ is said to
cone$\epsilon$-dominate $\boldsymbol{y}^{\prime}$ if and only if,
$\boldsymbol{y}$ Pareto-dominates $\boldsymbol{y}^{\prime}$ or the solution of
the linear system $\Psi\boldsymbol{\lambda}=\boldsymbol{z}$, with
$\boldsymbol{z}=\boldsymbol{y}^{\prime}-[\boldsymbol{y}-\boldsymbol{\epsilon}]$,
and $\epsilon_{i}>0$, gives $\lambda_{i}\geq
0~{}\forall~{}i\in\left\\{1,\ldots,m\right\\}$. Equivalently, we say
$\boldsymbol{y}\stackrel{{\scriptstyle
cone~{}\epsilon}}{{\prec}}\boldsymbol{y}^{\prime}$ if and only if,

$\left(\boldsymbol{y}\prec\boldsymbol{y}^{\prime}\right)~{}\vee~{}\left(\Psi\boldsymbol{\lambda}=\boldsymbol{z}\mid\lambda_{i}\geq
0\mbox{ for all }i=\left\\{1,\ldots,m\right\\}\right)~{}.$ (11)

$\Box$

Observe that this in fact represents a relaxed form of Pareto dominance, since
the Pareto-dominated region of any given point is a subset of the
cone$\epsilon$-dominated region. Moreover, it is reasonably easy to verify
that both the properties of convergence and diversity are satisfied by this
criterion. The convergence property is ensured by the Pareto-dominance
characteristics of the cone$\epsilon$ criterion and by the storing of
nondominated solutions in the archive population of a MOEA. As for the
diversity, since each box accommodates only a single vector, this property is
also guaranteed.

###  5.2 Maintaining a Cone $\epsilon$-Pareto Front

Similarly to the $\epsilon$-dominance, the archive update function for the
cone$\epsilon$-dominance strategy also employs a two level concept. At the
first, a discretization of the objective space into boxes is used, with a
single solution within each box. Any algorithm based on the
cone$\epsilon$-dominance relation always maintains a set of nondominated
solutions according to this criterion, which ensures the diversity property.
To this end, every solution in the archive is assigned a box index
$\boldsymbol{b}\in\mathbb{R}^{m}$:

$\boldsymbol{b}_{i}(\boldsymbol{y})=\left\\{\begin{array}[]{ll}\epsilon_{i}\left\lfloor
y_{i}/\epsilon_{i}\right\rfloor&,~{}~{}i\in\left\\{1,\ldots,m\right\\},~{}~{}\mbox{for
minimizing }f_{i}\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
\epsilon_{i}\left\lceil
y_{i}/\epsilon_{i}\right\rceil&,~{}~{}i\in\left\\{1,\ldots,m\right\\},~{}~{}\mbox{for
maximizing }f_{i}\end{array}\right.$ (12)

where $\left\lfloor\cdot\right\rfloor$ and $\left\lceil\cdot\right\rceil$
return the closest lower and upper integer to their argument, respectively.

On the second level, if two vectors share the same box, the former solution is
only replaced if it is Pareto dominated within the box, or else by a point
closest to the origin of the box, thus guaranteeing convergence. Algorithm 1
summarizes these concepts.

Input: $\mathcal{H}$, $\Psi$, $\boldsymbol{y}$

1 begin

2 if  _$\boldsymbol{y}$ is cone$\epsilon$-dominated by any
$\boldsymbol{y}^{\prime}\in\mathcal{H}$_ then

3 Reject $\boldsymbol{y}$;

4

5 else if  _$\boldsymbol{y}$ shares the same box with an archive member
$\boldsymbol{y}^{\prime}$_ then

6 if  _$\boldsymbol{y}$ dominates $\boldsymbol{y}^{\prime}$ or
$\boldsymbol{y}$ is closer to the origin of the box than
$\boldsymbol{y}^{\prime}$_ then

7 Delete all of the cone$\epsilon$-dominated archive members;

8 Replace $\boldsymbol{y}^{\prime}$ by $\boldsymbol{y}$;

9

10 else

11 Reject $\boldsymbol{y}$;

12

13

14 else if  _$\boldsymbol{y}$ cone$\epsilon$-dominates any
$\boldsymbol{y}^{\prime}\in\mathcal{H}$_ then

15 Delete all of the cone$\epsilon$-dominated archive members;

16 Insert $\boldsymbol{y}$ into the archive;

17

18 else

19 Insert $\boldsymbol{y}$ into the archive;

20

21

22

Output: $\mathcal{H}^{\prime}$

Algorithm 1 Archive update function performed by cone$\epsilon$-dominance.

###  5.3 Evaluating the Archive Size

Assume that all objectives $f_{i}$, $i\in\left\\{1,\ldots,m\right\\}$, are to
be minimized, and also that $1\leq f_{i}\leq K$, for all $i$. As discussed
previously, the cone$\epsilon$-dominance approach divides the objective space
into $\left(\left(K-1\right)/\epsilon\right)^{m}$ boxes, each allowed to
contain a single solution in the archive population $\mathcal{H}$. Since the
usual dominance relation ensures a monotonic front between the extreme boxes
of the hypergrid, the maximum number of boxes that can be “touched” by any
front is limited. However, the estimation of feasible solutions inside these
touched boxes depends on the connectivity of the Pareto front and on the value
of $\kappa$.

(a) Connected front.

(b) Concave surface.

(c) Disconnected front.

Figure 5: Illustration of the relation between $\epsilon$ and the size of the
archive $\mathcal{H}$ (Batista et al., 2011b, ).

In general, the number of “touched” boxes is maximum if a connected monotonic
front exists between the extreme boxes of the hypergrid. Figures 5(a) and 5(b)
illustrate two possible situations in which the number of estimated
cone$\epsilon$-Pareto solutions is maximum, i.e., equal to the number of boxes
that are touched by the front which is equal to five and seven, respectively.
These values are calculated by a simple relation governed by the number of
objectives, the range of the objective values, and the values of
$\epsilon_{i}$:

$\left|\mathcal{H}\right|\leq
m\left[\left(\prod\limits_{i=1}^{m-1}\frac{K-1}{\epsilon_{i}}\right)-\left(\prod\limits_{i=1}^{m-2}\frac{K-1}{\epsilon_{i}}\right)\right]+1~{}.$
(13)

It must be emphasized, however, that the maximum size of $\mathcal{H}$ cannot
be reached for a problem where the Pareto front is disconnected. Nevertheless,
the cone$\epsilon$-dominance relation is still likely to estimate one solution
for each box touched by the front.

In the particular case where $\epsilon_{i}=\epsilon$, $\forall
i\in\left\\{1,\ldots,m\right\\}$, the previous expression can be rewritten as:

$\left|\mathcal{H}\right|\leq
m\left[\left(\frac{K-1}{\epsilon}\right)^{m-1}-\left(\frac{K-1}{\epsilon}\right)^{m-2}\right]+1~{}.$
(14)

Figure 5(c) illustrates a problem where the Pareto front is disconnected,
which means that the maximum size of $\mathcal{H}$ cannot be reached. However,
the cone$\epsilon$-dominance relation is still likely to estimate one solution
for each box touched by the front.

The $\epsilon$ value that provides a maximum of $T$ nondominated points in the
archive, according to the cone$\epsilon$-criterion, can be easily calculated
from (14) by solving (15):

$\left(\frac{T-1}{m}\right)\epsilon^{(m-1)}+\left(K-1\right)^{m-2}\epsilon-\left(K-1\right)^{m-1}=0~{}.$
(15)

Finally, observe that specific bounds on the objective values are not used in
the Alg. 1 and are not required to ensure the convergence. They are only
employed to demonstrate the relation between $\epsilon$ and the size of the
archive $\mathcal{H}$.

##  6 Experimental Design and Validation of the Proposed Approach

To validate the proposed cone$\epsilon$-dominance approach, six algorithms are
considered for the experimental study: the first two are the well-known NSGA-
II (Deb et al.,, 2002) and SPEA2 (Zitzler et al.,, 2001) methods, both based
on the standard Pareto dominance relation; the third is a modified NSGA-II, in
which the crowding distance is re-calculated after each point removal in order
to improve diversity; the fourth is a clustering version of the NSGA-II, in
which the crowding distance procedure is replaced with a clustering approach;
the fifth is a steady-state MOEA based on the $\epsilon$-dominance strategy;
and the sixth is implemented by modifying the last one, replacing the
$\epsilon$-dominance mechanism with the cone$\epsilon$-dominance approach. The
comparison of the two last methods will enable us to show the performance of
the same algorithm with and without cone$\epsilon$-dominance. Furthermore,
this study will contrast the ability of the cone$\epsilon$-criterion in
maintaining diversity against operators which handle diversity efficiently,
e.g., the clustering operator (Zitzler and Thiele,, 1999) used in the C-NSGA-
II and the truncation method of SPEA2 (Zitzler et al.,, 2001).

The multiobjective evolutionary algorithms are briefly described next:

  1. 1.

NSGA-II: This algorithm was proposed by Deb et al., (2002). In general terms,
the parent and offspring populations are combined and evaluated using the fast
nondominated sorting approach, an elitist approach, and an efficient crowding
mechanism to implement density estimation. When the memory archive set is
greater than its capability, only those solutions that are maximally apart
from their neighbors, according to the crowding measure, are maintained. In
the original NSGA-II the crowding measure is calculated only once, and this
information is used to remove a specific number of solutions with the smallest
crowding value. Considering a population with $N$ points, the computational
complexity is $\mathcal{O}\left(mN^{2}\right)$ for the nondominated sort,
$\mathcal{O}\left(mN\log N\right)$ for the crowding distance assignment, which
makes the algorithm computationally fast, and $\mathcal{O}\left(N\log
N\right)$ for the crowded comparison operator. In this way, the overall
complexity of the NSGA-II is dominated by the nondominated sorting procedure,
i.e., $\mathcal{O}\left(mN^{2}\right)$. Though the NSGA-II has been shown to
perform quite well for two or three-objective problems, its crowding operator
is not adequate in maintaining a good distribution of solutions in a problem
having many objective functions (Deb et al., 2005a, ).

  2. 2.

SPEA2: The SPEA2 (Zitzler et al.,, 2001) presents some differences with
respect to its predecessor, SPEA (Zitzler and Thiele,, 1999): (i) it employs a
fine-grained fitness assignment strategy; (ii) it uses a nearest neighbor
density estimation technique which guides the search more efficiently; (iii)
the archive size is fixed; (iv) it has an enhanced archive truncation method
that guarantees the preservation of boundary solutions; and (v) only members
of the archive participate in the mating selection process. The computational
complexity of the fitness assignment procedure is controlled by the density
estimator, $\mathcal{O}\left(mN^{2}\log N\right)$, but the worst run-time
complexity of the SPEA2 is dominated by the truncation operator, that is
$\mathcal{O}\left(mN^{3}\right)$. However, on average the complexity will be
$\mathcal{O}\left(mN^{2}\log N\right)$ since the sorting of the distances
governs the overall complexity (Zitzler et al.,, 2001).

  3. 3.

NSGA-II*: Instead of using the standard NSGA-II’s crowding distance mechanism,
the NSGA-II* performs an improved diversity preservation approach by re-
calculating the crowding distance after each point removal. In this situation,
the worst-case complexity of the crowding distance assignment is
$\mathcal{O}\left(m(2N)\log(2N)\right)$. This means that the overall
complexity of the NSGA-II* is $\mathcal{O}\left(mN^{2}\right)$, which is due
to the nondominated sorting procedure, as observed for the NSGA-II. However, a
better distributed set of solutions is supposed to be achieved with the NSGA-
II*.

  4. 4.

C-NSGA-II: This method was discussed in (Deb et al.,, 2003), and is a
straightforward replacement of the NSGA-II’s crowding mechanism with the
clustering approach used in SPEA (Zitzler and Thiele,, 1999), which has a
computational complexity of $\mathcal{O}\left(mN^{3}\right)$, involving
Euclidean distance calculations. In spite of the larger computational time
required by this method, the clustered NSGA-II is expected to find a better
distributed set of nondominated solutions than the original NSGA-II.

  5. 5.

$\epsilon$-MOEA: This approach was proposed by Deb et al., (2003); Deb et al.,
2005a , and consists of a steady-state MOEA based on the $\epsilon$-dominance
concept introduced in Laumanns et al., (2002). In this method, two
populations, evolutionary population and archive population, are evolved
simultaneously, and two offspring solutions are created by using one solution
from each population. Each offspring is then used to update both parent and
archive populations in a steady-state manner, thereby providing better chances
of creating good offspring solutions. Note, however, that the archive
population is updated based on the $\epsilon$-dominance concept, whereas an
usual domination concept is used to update the parent population. Since a
fitness assignment method is not implemented in the $\epsilon$-MOEA, its
complexity includes only the operations of maintaining and truncating the
memory population, being these techniques handled together in the archive
update process. The computational complexity of the $\epsilon$-MOEA, for the
entire population, requires comparisons of $\mathcal{O}\left(mN^{2}\right)$ in
the worst case. However, it is important to note that a computational time of
$\mathcal{O}\left(mN^{2}\right)$ is required in addition due to the
simultaneous evolution of two independent populations (Li et al.,, 2011). A
careful observation reveals that the $\epsilon$-MOEA procedure emphasizes
nondominated solutions, maintains the diversity in the archive by allowing
only one solution inside each hyperbox on the Pareto-optimal front, and
represents an elitist approach (Deb et al.,, 2003).

  6. 6.

cone$\epsilon$-MOEA: This algorithm is a modification of the $\epsilon$-MOEA,
in which we include cone$\epsilon$-dominance instead of the regular
$\epsilon$-dominance concept. Due to the small differences regarding the
archive update strategies of these methods, a similar discussion can be made
for the time complexity of the cone$\epsilon$-MOEA. In this way, considering
one updating operation of a single offspring (see Algorithm 1), the complexity
order of total comparisons required is $\mathcal{O}\left(mN\right)$ to check
the Pareto dominance and $\mathcal{O}\left(m^{2}N\right)$ to check the cone of
dominance, in which the $m^{2}$ time cost refers to a linear system solution
by factorization. However, since no further comparisons are required if the
candidate solution is dominated at the beginning of the process, the average
time complexity gets lower, and the best case run-time required is
$\mathcal{O}\left(m\right)$. In a straightforward manner, the computational
complexity of the cone$\epsilon$-MOEA, for the entire population, requires a
run-time of $\mathcal{O}\left((mN)^{2}\right)$ in the worst case.

Based on the previous discussion, and without considering the effect of
special data structures, it is possible to sort the algorithms with respect to
computational complexity. In this way, the sequence of run-time efficiency
from the lowest to highest is: C-NSGA-II, SPEA2, cone$\epsilon$-MOEA, NSGA-
II*, NSGA-II, and $\epsilon$-MOEA. These methods have been implemented in
Matlab, and some of them are based on the source codes available in (KanGAL,,
2011)111The implementation of the algorithms, the samples of the true Pareto
fronts used, and the raw and processed results of the experiments can be
retrieved from (Batista et al., 2011c, ).. Further information about test
problems, performance metrics, parameter settings, and statistical design are
presented in the following sections.

###  6.1 Benchmark Test Problems

In a previous work (Batista et al., 2011b, ), we performed a limited
simulation study in which the choice of problems was directed by the
geometrical characteristics of the Pareto fronts rather than the difficulty of
solving each test problem. Now, in order to carry out a more appropriate
analysis, we consider more complex test problems, each providing a different
kind of challenge for multiobjective evolutionary algorithms. First, we have
chosen two continuous and simple tests, Deb52 (Deb,, 1999) and Poloni’s
problem (Poloni,, 1995), both employed in the previous work (Batista et al.,
2011b, ). The five two-objective ZDT problems were also selected (Zitzler et
al.,, 2000). Since this paper deals only with continuous-parameter problems,
the ZDT5 problem was not considered as it is defined for binary strings. At
last, we consider nine three-objective DTLZ test problems (Deb et al., 2005b,
). In this study, the constraints were handled using a simple penalty method
($f_{i}+1000\sum_{j=1}^{n_{g}}\max\left(0,g_{j}\right)$) for all the
algorithms. Further details about these benchmark problems are presented in
Table 1 and throughout the paper.

Table 1: Analytical test problems adopted in the experimental study.

$\begin{array}[]{l|c|c|l|l}\hline\cr\mbox{Problem}&$n$&\mbox{Bounds}&\mbox{Objective
functions}&\mbox{Pareto front}\\\ \hline\cr\mbox{P1:
Deb52}&2&[0,1]&f_{1}(\boldsymbol{x})=1-\exp(-4x_{1})\sin^{4}(10\pi
x_{1})&\mbox{Concave}\\\ &&&f_{2}(\boldsymbol{x})=g(x_{2})h(x_{1})\mbox{ ;
}g(x_{2})=1+x_{2}^{2}&\\\
&&&h(x_{1})=\left\\{\begin{array}[]{ll}1-\left(\frac{f_{1}(\boldsymbol{x})}{g(x_{2})}\right)^{10}&\mbox{if
}f_{1}(\boldsymbol{x})\leq g(x_{2})\\\ 0&\mbox{otherwise.}\\\
\end{array}\right.&\\\ \hline\cr\mbox{P2:
Pol}&2&[-\pi,\pi]&f_{1}(\boldsymbol{x})=1+(A_{1}-B_{1})^{2}+(A_{2}-B_{2})^{2}&\mbox{Nonconvex
and}\\\
&&&f_{2}(\boldsymbol{x})=(x_{1}+3)^{2}+(x_{2}+1)^{2}&\mbox{disconnected}\\\
&&&A_{1}=0.5\sin 1-2\cos 1+\sin 2-1.5\cos 2&\\\ &&&A_{2}=1.5\sin 1-\cos
1+2\sin 2-0.5\cos 2&\\\ &&&B_{1}=0.5\sin x_{1}-2\cos x_{1}+\sin x_{2}-1.5\cos
x_{2}&\\\ &&&B_{2}=1.5\sin x_{1}-\cos x_{1}+2\sin x_{2}-0.5\cos x_{2}&\\\
\hline\cr\mbox{P3: ZDT1}&30&[0,1]&f_{1}(\boldsymbol{x})=x_{1}\mbox{ ;
}f_{2}(\boldsymbol{x})=1-\sqrt{x_{1}/g(\boldsymbol{x})}&\mbox{Convex}\\\
&&&g(\boldsymbol{x})=1+9\left(\sum_{i=2}^{n}x_{i}\right)/(n-1)&\\\
\hline\cr\mbox{P4: ZDT2}&30&[0,1]&f_{1}(\boldsymbol{x})=x_{1}\mbox{ ;
}f_{2}(\boldsymbol{x})=1-\left(x_{1}/g(\boldsymbol{x})\right)^{2}&\mbox{Nonconvex}\\\
&&&g(\boldsymbol{x})=1+9\left(\sum_{i=2}^{n}x_{i}\right)/(n-1)&\\\
\hline\cr\mbox{P5: ZDT3}&30&[0,1]&f_{1}(\boldsymbol{x})=x_{1}\mbox{ ;
}f_{2}(\boldsymbol{x})=1-\sqrt{h}-h\sin\left(10\pi
x_{1}\right)&\mbox{Disconnected}\\\
&&&g(\boldsymbol{x})=1+9\left(\sum_{i=2}^{n}x_{i}\right)/(n-1)\mbox{ and
}h=x_{1}/g(\boldsymbol{x})&\mbox{convex parts}\\\ \hline\cr\mbox{P6:
ZDT4}&10&x_{1}\in[0,1]&f_{1}(\boldsymbol{x})=x_{1}&\mbox{Nonconvex}\\\
&&x_{i}\in[-5,5],&f_{2}(\boldsymbol{x})=1-\sqrt{x_{1}/g(\boldsymbol{x})}&\mbox{(multimodal}\\\
&&i=2,...,n&g(\boldsymbol{x})=1+10\left(n-1\right)+\sum_{i=2}^{n}\left(x_{i}^{2}-10\cos\left(4\pi
x_{i}\right)\right)&\mbox{problem)}\\\ \hline\cr\mbox{P7:
ZDT6}&10&[0,1]&f_{1}(\boldsymbol{x})=1-\exp(-4x_{1})\sin^{6}(6\pi
x_{1})&\mbox{Nonconvex}\\\
&&&f_{2}(\boldsymbol{x})=1-(f_{1}/g(\boldsymbol{x}))^{2}&\mbox{(nonuniform}\\\
&&&g(\boldsymbol{x})=1+9\left[\left(\sum_{i=2}^{n}x_{i}\right)/(n-1)\right]^{0.25}&\mbox{search
space)}\\\ \hline\cr\mbox{P8:
DTLZ1}&7&[0,1]&f_{1}(\boldsymbol{x})=0.5x_{1}x_{2}(1+g(\boldsymbol{x}_{m}))&\mbox{Linear}\\\
&&&f_{2}(\boldsymbol{x})=0.5x_{1}(1-x_{2})(1+g(\boldsymbol{x}_{m}))&\mbox{hyperplane}\\\
&&&f_{3}(\boldsymbol{x})=0.5(1-x_{1})(1+g(\boldsymbol{x}_{m}))&\\\
&&&g(\boldsymbol{x}_{m})=100\left(\left|\boldsymbol{x}_{m}\right|+\sum_{x_{i}\in\boldsymbol{x}_{m}}h_{i}\right)&\\\
&&&h_{i}=\left(x_{i}-0.5\right)^{2}-\cos(20\pi(x_{i}-0.5))&\\\
\hline\cr\mbox{P9:
DTLZ2}&n&[0,1]&f_{1}(\boldsymbol{x})=(1+g(\boldsymbol{x}_{m}))\prod_{i=1}^{m-1}\cos(x_{i}^{\alpha}\pi/2)&\mbox{Concave
surface}\\\
&&&f_{j}(\boldsymbol{x})=(1+g(\boldsymbol{x}_{m}))\left(\prod_{i=1}^{m-j}\cos(x_{i}^{\alpha}\pi/2)\right)\sin(x_{M}^{\alpha}\pi/2)&\\\
&&&f_{m}(\boldsymbol{x})=(1+g(\boldsymbol{x}_{m}))\sin(x_{1}^{\alpha}\pi/2)\mbox{
; }M=m-j+1&\\\
&&&g(\boldsymbol{x}_{m})=\sum_{x_{i}\in\boldsymbol{x}_{m}}\left(x_{i}-0.5\right)^{2}\mbox{
; }\alpha=1&\\\ &&&n=m+k-1\mbox{ , }k=10\mbox{ and }j=2,\ldots,m-1&\\\
\hline\cr\mbox{P10: DTLZ3}&12&[0,1]&f_{1}\mbox{, }f_{2}\mbox{ and }f_{3}\mbox{
as in DTLZ2}&\mbox{Concave surface}\\\
&&&g(\boldsymbol{x}_{m})=100\left(\left|\boldsymbol{x}_{m}\right|+\sum_{x_{i}\in\boldsymbol{x}_{m}}h_{i}\right)&\\\
&&&h_{i}=\left(x_{i}-0.5\right)^{2}-\cos(20\pi(x_{i}-0.5))&\\\
\hline\cr\mbox{P11: DTLZ4}&12&[0,1]&f_{1}\mbox{, }f_{2}\mbox{ and }f_{3}\mbox{
as in DTLZ2 with }\alpha=100&\mbox{Concave surface}\\\ \hline\cr\mbox{P12:
DTLZ5}&12&[0,1]&f_{1}(\boldsymbol{x})=(1+g(\boldsymbol{x}_{m}))\cos(\theta_{1})\cos(\theta_{2})&\mbox{Concave
curve}\\\
&&&f_{2}(\boldsymbol{x})=(1+g(\boldsymbol{x}_{m}))\cos(\theta_{1})\sin(\theta_{2})&\\\
&&&f_{3}(\boldsymbol{x})=(1+g(\boldsymbol{x}_{m}))\sin(\theta_{1})&\\\
&&&\theta_{1}=x_{1}\pi/2\mbox{ and
}\theta_{2}=\frac{\pi}{4(1+g(r))}(1+2g(r)x_{2})&\\\
&&&g(\boldsymbol{x}_{m})=\sum_{x_{i}\in\boldsymbol{x}_{m}}\left(x_{i}-0.5\right)^{2}&\\\
\hline\cr\mbox{P13: DTLZ6}&12&[0,1]&f_{1}\mbox{, }f_{2}\mbox{ and }f_{3}\mbox{
as in DTLZ5}&\mbox{Concave curve}\\\
&&&g(\boldsymbol{x}_{m})=\sum_{x_{i}\in\boldsymbol{x}_{m}}x_{i}^{0.1}&\\\
\hline\cr\mbox{P14: DTLZ7}&22&[0,1]&f_{1}(\boldsymbol{x})=x_{1}\mbox{ ;
}f_{2}(\boldsymbol{x})=x_{2}\mbox{ ;
}f_{3}(\boldsymbol{x})=(1+g(\boldsymbol{x}_{m}))h&\mbox{Disconnected}\\\
&&&h=3-\sum_{i=1}^{2}\left[\frac{f_{i}}{1+g}\left(1+\sin(3\pi
f_{i})\right)\right]&\mbox{(noncontiguous}\\\
&&&g(\boldsymbol{x}_{m})=1+\frac{9}{20}\sum_{x_{i}\in\boldsymbol{x}_{m}}x_{i}&\mbox{convex
regions)}\\\ \hline\cr\mbox{P15:
DTLZ8}&30&[0,1]&f_{i}(\boldsymbol{x})=0.1\sum_{j=10(i-1)+1}^{10j}x_{i}\mbox{ ;
}i=1,2,3&\mbox{Combination of}\\\
&&&g_{i}(\boldsymbol{x})=f_{3}(\boldsymbol{x})+4f_{i}(\boldsymbol{x})-1\geq
0\mbox{ ; }i=1,2&\mbox{a straight line}\\\
&&&g_{3}(\boldsymbol{x})=2f_{3}(\boldsymbol{x})+f_{1}(\boldsymbol{x})+f_{2}(\boldsymbol{x})-1\geq
0&\mbox{and a hyperplane}\\\ \hline\cr\mbox{P16:
DTLZ9}&30&[0,1]&f_{i}(\boldsymbol{x})=\sum_{j=10(i-1)+1}^{10j}x_{i}^{0.1}\mbox{
; }i=1,2,3&\mbox{Concave curve}\\\
&&&g_{i}(\boldsymbol{x})=f_{3}^{2}(\boldsymbol{x})+f_{i}^{2}(\boldsymbol{x})-1\geq
0\mbox{ ; }i=1,2&\\\ \hline\cr\end{array}$

###  6.2 Performance Metrics

Evolutionary multiobjective optimization techniques are required to consider
two different goals, i.e., beyond acquiring convergence to the Pareto-optimal
front, an equally important task is to find and maintain a diverse set of
solutions. In addition, achieving these goals in a small computational time is
also an important issue in MOEAs. To consider this multi-criterion nature in
the evaluation of multi-objective algorithms, regarding the convergence and
diversity of the solutions found, we have used four different metrics.

In order to assess how near the solutions found are from the Pareto-optimal
front, the convergence metric $(\gamma)$, proposed by Deb et al., (2002), was
considered. This quality indicator, see Fig. 6(a), measures the distance
between the obtained nondominated front $\mathcal{H}$ and a detailed sampling
of the true Pareto-optimal front $\mathcal{H}^{*}$:

$\gamma=\frac{\sum_{i=1}^{\left|\mathcal{H}\right|}d_{i}}{\left|\mathcal{H}\right|}$
(16)

where $d_{i}$ is the Euclidean distance, in the objective space, between the
solution $i\in\mathcal{H}$ and the nearest member of $\mathcal{H}^{*}$, and
the operator $\left|~{}\cdot~{}\right|$ returns the cardinality of the set in
its argument. So, the lower the $\gamma$ value the better the convergence of
the solutions in $\mathcal{H}$. As illustrated in the Fig. 6(a), a result with
$\gamma=0$ means $\mathcal{H}\subseteq\mathcal{H}^{*}$, in which all the
estimated solutions are Pareto-optimal.

Because MOEAs are required to achieve an interesting trade-off between the
convergence and distribution of the solutions approximated, the diversity
metric $(\Delta)$ (Deb et al.,, 2002) was also used. This quality indicator
measures the extent of spread achieved among the obtained nondominated
solutions in $\mathcal{H}$. Considering that it is desirable to obtain a set
of solutions that spans the entire Pareto-optimal region, the $\Delta$ value
is defined as:

$\Delta=\frac{\sum_{i=1}^{m}d_{i}^{e}+\sum_{i=1}^{\left|\mathcal{H}\right|}\left|d_{i}-\bar{d}\right|}{\sum_{i=1}^{m}d_{i}^{e}+\left|\mathcal{H}\right|\bar{d}}$
(17)

where $d_{i}^{e}$ denotes the Euclidean distance between the extreme points in
$\mathcal{H}$ and $\mathcal{H}^{*}$ along the $i^{th}$ coordinate, and $d_{i}$
measures the Euclidean distance of each point in $\mathcal{H}$ to its closest
neighbor. Therefore, a lower $\Delta$ value indicates a better distribution of
solutions, as can be seen from Fig. 6(b). Notice that a result with $\Delta=0$
means the extreme points of $\mathcal{H}^{*}$ have been found and $d_{i}$
equals to $\bar{d}$ for all $i$.

The third metric, known as S-metric or Hypervolume $(HV)$ metric (Zitzler and
Thiele,, 1999), calculates the hypervolume enclosed by the estimated front
$\mathcal{H}$ and a reference point $\boldsymbol{y}_{ref}$ dominated by all
solutions in this front. Regarding a minimization MOP, the larger the
dominated hypervolume, the better the front is (see the Fig. 6(c)). For all
test problems, the reference point was defined as $10\%$ greater than the
upper boundaries of the real Pareto-optimal front. Formally, this metric is
described as the Lebesgue measure $\Lambda$ of the union of hypercubes $h_{i}$
defined by a nondominated point $\boldsymbol{y}_{i}$ and
$\boldsymbol{y}_{ref}$:

$HV(\mathcal{H})=\Lambda\left(\left\\{\bigcup\limits_{i}h_{i}\mid\boldsymbol{y}_{i}\in\mathcal{H}\right\\}\right)=\Lambda\left(\bigcup\limits_{\boldsymbol{y}\in\mathcal{H}}\left\\{\boldsymbol{y}^{\prime}\mid\boldsymbol{y}\prec\boldsymbol{y}^{\prime}\prec\boldsymbol{y}_{ref}\right\\}\right)$
(18)

Even though this metric estimates both convergence and diversity of the
solutions in $\mathcal{H}$, it is more sensitive to the convergence of the
points towards the Pareto-optimal front, and it prefers convex regions to non-
convex ones (Zitzler et al.,, 2007). For instance, as shown in Fig. 6(c), the
contribution of the solution $i^{\prime}$ to the $HV$ measure is greater than
that of the solution $i^{\prime\prime}$.

In order to evaluate the convergence of each algorithm in contrast to all
others, a generalization of the Coverage of Two Sets metric (Zitzler and
Thiele,, 1999) is proposed here. This metric, called here the Coverage of Many
Sets $(CS)$, quantifies the domination of the final population of one
algorithm over the union of the remaining ones. The proposed $CS$ function is
stated as:

$CS(X_{i},U_{i})=\frac{\left|a^{\prime\prime}\in
U_{i};~{}\exists~{}a^{\prime}\in X_{i}:~{}a^{\prime}\preceq
a^{\prime\prime}\right|}{\left|U_{i}\right|}$ (19)

where $X_{i}$ and $U_{i}$ are sets of objective vectors, and
$a^{\prime}\preceq a^{\prime\prime}$ means that $a^{\prime}$ covers
$a^{\prime\prime}$, that is, either $a^{\prime}\prec a^{\prime\prime}$ or
$a^{\prime}=a^{\prime\prime}$. Function $CS$ maps the pair
$\left(X_{i},U_{i}\right)$ to the interval $[0,1]$, in which $X_{i}$, for all
$i=1,\ldots,k$, denotes the final Pareto front resulting from algorithm $i$,
and $U_{i}$, defined as:

$U_{i}=\bigcup\limits_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{k}X_{j}$ (20)

represents the union of the final Pareto fronts of all the $k$ algorithms,
except $i$. The value $CS\left(X_{i},U_{i}\right)=1$ implies that all points
in $U_{i}$ are dominated by or equal to points in $X_{i}$. The opposite,
$CS\left(X_{i},U_{i}\right)=0$, represents the situation when none of the
points in $U_{i}$ are covered by the set $X_{i}$. Note that for $k=2$, the
$CS$ metric tends to the Coverage of Two Sets proposed in (Zitzler and
Thiele,, 1999). Although the measures $CS\left(X_{i},U_{i}\right)$ and
$CS\left(U_{i},X_{i}\right)$ are not complementary, by simplicity this study
considers only the comparisons $CS\left(X_{i},U_{i}\right)$, for all
$i=1,\ldots,k$. Figure 6(d) shows the $CS$ measures for a particular $k=3$
case. As illustrated, $CS\left(X_{1},U_{1}\right)=4/12$,
$CS\left(X_{2},U_{2}\right)=8/12$, and $CS\left(X_{3},U_{3}\right)=0$,
indicating the superior quality of the front $X_{2}$, followed by $X_{1}$ and
$X_{3}$, respectively.

For the first two metrics, a detailed sampling of the true Pareto-optimal
front of each problem must be known. Since we are dealing with test problems,
the true Pareto-optimal front is not difficult to obtain. In this work, we
have used uniformly spaced Pareto-optimal solutions as the approximation of
the true Pareto-optimal front. The reference fronts used here can be retrieved
online (Batista et al., 2011c, ).

(a) Convergence metric ($\gamma$).

(b) Diversity metric ($\Delta$).

(c) S-metric (HV).

(d) Coverage of many sets ($CS$).

Figure 6: Quality indicators adopted in the experimental study. In these
illustrations, the sets $\mathcal{H}^{*}$ and $\mathcal{H}$ represent the
global Pareto front and the estimated one, respectively. According to letter
(c), the contribution of the solution $i^{\prime}$ for the $HV$ measure is
greater than that of the solution $i^{\prime\prime}$. In (d), the coverage of
many sets measures are $CS\left(X_{1},U_{1}\right)=4/12$,
$CS\left(X_{2},U_{2}\right)=8/12$, and $CS\left(X_{3},U_{3}\right)=0$,
indicating the superior quality of the front $X_{2}$, followed by $X_{1}$ and
$X_{3}$.

###  6.3 Parameter Settings

With the objective of providing a comparison baseline for the performance of
the proposed cone$\epsilon$-dominance strategy, some parameter settings were
adopted for all the algorithms: population size $N=100$, and probabilities of
crossover and mutation $p_{xover}=1$ and $p_{mut}=1/n$, respectively. Since
all test problems involve real-valued decision variables, we have used the
simulated binary crossover (SBX) operator (Deb and Agrawal,, 1995) and the
polynomial mutation operator (Deb and Goyal,, 1996) to create the offspring
solutions. The crossover ($\eta_{xover}$) and mutation ($\eta_{mut}$)
distribution indices adopted for each problem are shown in Table 2.
Furthermore, for a fair comparison, all algorithms are iterated for a fixed
number of solution evaluations: 100,000 evaluations for ZDT4 and DTLZ8; 50,000
for DTLZ9; 30,000 for DTLZ3 and DTLZ4; and 20,000 for the remaining ones.

Table 2: Crossover and mutation distribution indices adopted for the algorithms. | Deb52, Pol, ZDT1, ZDT3 | ZDT4, DTLZ1, | ZDT2,
---|---|---|---
| ZDT6, DTLZ2, DTLZ4 | DTLZ3 | DTLZ5–DTLZ9
$\eta_{xover}$ | 15 | 2 | 2
$\eta_{mut}$ | 20 | 20 | 5
Table 3: Calculated (first line) and estimated (second line)
$\boldsymbol{\epsilon}$ values for the algorithms
$\underline{\boldsymbol{\epsilon}}$-MOEA and
$\underline{\boldsymbol{c}}$one$\epsilon$-MOEA. These values have been
obtained in order to get roughly 100 solutions in the archive at the end of
the run. The sign “–” means that it was not necessary to estimate values for
$\boldsymbol{\epsilon}$. In these cases, the specific calculation method has
provided adequate values.

Two-Objective Test Problems
---
Alg. | Deb52 | Pol | ZDT1 | ZDT2 | ZDT3 | ZDT4 | ZDT6
$\boldsymbol{\epsilon}$ | [0.0083,0.010] | [0.16,0.25] |  $\epsilon_{i}$ = 0.010 |  $\epsilon_{i}$ = 0.010 | [0.0085,0.018] |  $\epsilon_{i}$ = 0.010 | [0.0072,0.0093]
| [0.003,0.003] | [0.038,0.038] |  $\epsilon_{i}$ = 0.0075 |  $\epsilon_{i}$ = 0.0076 |  $\epsilon_{i}$ = 0.0026 |  $\epsilon_{i}$ = 0.0065 |  $\epsilon_{i}$ = 0.0067
$\boldsymbol{c}$ | [0.0164,0.0198] | [0.3168,0.4950] |  $\epsilon_{i}$ = 0.0198 |  $\epsilon_{i}$ = 0.0198 | [0.0168,0.0356] |  $\epsilon_{i}$ = 0.0198 | [0.0143,0.0184]
| – | [0.20,0.33] | – | – | [0.012,0.025] | – | –

Three-Objective Test Problems
---
Alg. | DTLZ1 | DTLZ2 | DTLZ3 | DTLZ4 | DTLZ5
$\boldsymbol{\epsilon}$ |  $\epsilon_{i}$ = 0.05 |  $\epsilon_{i}$ = 0.10 |  $\epsilon_{i}$ = 0.10 |  $\epsilon_{i}$ = 0.10 | [0.007,0.007,0.01]
| [0.02,0.02,0.05] | [0.06,0.06,0.066] | [0.06,0.06,0.066] |  $\epsilon_{i}$ = 0.062 |  $\epsilon_{i}$ = 0.005
$\boldsymbol{c}$ |  $\epsilon_{i}$ = 0.0798 |  $\epsilon_{i}$ = 0.1595 |  $\epsilon_{i}$ = 0.1595 |  $\epsilon_{i}$ = 0.1595 | [0.014,0.014,0.02]
| [0.05,0.05,0.0833] | – | – | – |  $\epsilon_{i}$ = 0.025

Alg. | DTLZ6 | DTLZ7 | DTLZ8 | DTLZ9
---|---|---|---|---
$\boldsymbol{\epsilon}$ | [0.007,0.007,0.01] | [0.086,0.086,0.3386] | [0.075,0.075,0.10] |  $\epsilon_{i}$ = 0.010
|  $\epsilon_{i}$ = 0.005 |  $\epsilon_{i}$ = 0.05 | [0.02,0.02,0.04] |  $\epsilon_{i}$ = 0.025
$\boldsymbol{c}$ | [0.014,0.014,0.02] | [0.1372,0.1372,0.5404] | [0.12,0.12,0.16] |  $\epsilon_{i}$ = 0.0198
|  $\epsilon_{i}$ = 0.017 | [0.12,0.12,0.30] |  $\epsilon_{i}$ = 0.03 | [0.27,0.27,0.25]

In order to obtain a final archive population size of $100$ solutions, we have
used the $\boldsymbol{\epsilon}$ values showed in Table 3. The calculated
$\boldsymbol{\epsilon}$ values have been obtained following the guidelines
provided by Laumanns et al., (2002) for the $\epsilon$-MOEA (8), and by
Batista et al., 2011b  for the cone$\epsilon$-MOEA (15). Since both
$\epsilon$-dominance and cone$\epsilon$-dominance may lose nondominated
points, we have employed estimated $\boldsymbol{\epsilon}$ values to get
roughly $100$ solutions in the final archive. For the $\epsilon$-approach, the
estimated values were obtained from (Deb et al.,, 2003; Deb et al., 2005a, ),
whereas for the cone$\epsilon$-approach the estimation was performed by
testing different $\boldsymbol{\epsilon}$ values to get roughly $100$
solutions in the archive at the end of the run.

Note that, for the cone$\epsilon$-MOEA, the estimated $\boldsymbol{\epsilon}$
values on Table 3 differ from the calculated ones only in eight of the sixteen
problems, among which three have a disconnected Pareto-optimal front, and a
non adequate calculated $\boldsymbol{\epsilon}$ value was expected. On the
other hand, in the $\epsilon$-MOEA all $\boldsymbol{\epsilon}$ values have
changed. This phenomenon suggests that the cone$\epsilon$-dominance approach
may be less susceptible to the loss of nondominated solutions than the
$\epsilon$-dominance. This situation is clear from Fig. 7, which shows the
distribution of solutions obtained on problems Deb52 and Pol in order to get
$T=30$ and $T=50$ solutions, respectively. The number of solutions
approximated by both approaches is indicated in the figure by considering the
cardinality of the final Pareto front $\left|\mathcal{H}\right|$. In this
simulation, the $\boldsymbol{\epsilon}$ values were calculated using (8) for
the $\epsilon$-MOEA and (15) for the cone$\epsilon$-MOEA. The boxes within
which a solution lies are also shown in the figure. It is interesting to
notice that all solutions are nondominated according to the
$\epsilon$-dominance or the cone$\epsilon$-dominance depending on the strategy
adopted, and each box is occupied by only one solution. Furthermore, the
cone$\epsilon$-approach has achieved a better spread of solutions in both
tests, even in the Pol problem, which has a disconnected Pareto-optimal front.
Because of the loss of efficient solutions by both strategies, only the
“estimated” $\boldsymbol{\epsilon}$ values are considered in the remainder of
the experimental study.

As the cone$\epsilon$-dominance criterion is influenced by the $\kappa$
parameter, we have performed some preliminary testing to observe the effect of
different $\kappa$ values on the performance of the cone$\epsilon$-MOEA. Table
4 shows the effect of this parameter on the values of the unary metrics
$\gamma$, $\Delta$ and $HV$ for the benchmark problems Deb52, ZDT1, and DTLZ2
with three and four objective functions. In this simple test, intermediate
values for $\kappa$ seem to yield reasonably good performance values for all
metrics, from which a value of $\kappa=0.5$ was chosen for all experiments
conducted in this work. However, as the number of objectives increases,
smaller values of $\kappa$ (i.e., $\kappa\rightarrow 0$) seem to be more
adequate, at the cost of the loss of some nondominated solutions since the
effect of the cone of dominance approximates that one of the
$\epsilon$-approach.

Table 4: Influence of different $\kappa$ values on the performance of the
cone$\epsilon$-MOEA on test problems Deb52, ZDT1, and DTLZ2 with three and
four objective functions. Median (M) and standard deviation (SD) over 30
independent runs are shown. Intermediate values for $\kappa$ seem to yield
reasonably good performance values for all metrics. A more appropriate study
is required in order to formally characterize the effect of this parameter.

$\begin{array}[]{clccccccccccc}\\\
\hline\cr\lx@intercol\mbox{Metric}\hfil\lx@intercol&\lx@intercol\hfil\mbox{$\kappa$
; Deb52}\hfil\lx@intercol\\\
\cline{3-13}\cr&&\mbox{0.0}&\mbox{0.1}&\mbox{0.2}&\mbox{0.3}&\mbox{0.4}&\mbox{0.5}&\mbox{0.6}&\mbox{0.7}&\mbox{0.8}&\mbox{0.9}&\mbox{0.99}\\\
\hline\cr\mbox{$\gamma$}&\mbox{M}&0.0006&0.0006&0.0005&0.0006&0.0005&0.0006&0.0006&0.0006&0.0006&0.0006&0.0006\\\
&\mbox{SD}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&0.0001&0.0001\\\
\mbox{$\Delta$}&\mbox{M}&0.6766&0.6813&0.5244&0.2991&0.2552&0.2432&0.2648&0.2892&0.3147&0.3194&0.3199\\\
&\mbox{SD}&0.0004&0.0021&0.0025&0.0027&0.0034&0.0039&0.0017&0.0019&0.0016&0.0042&0.0066\\\
\mbox{HV}&\mbox{M}&0.2735&0.2779&0.2794&0.2802&0.2806&0.2806&0.2806&0.2806&0.2806&0.2806&0.2806\\\
&\mbox{SD}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}&<10^{-4}\\\
\mbox{$\left|\mathcal{H}\right|$}&\mbox{M}&19.00&51.00&74.00&93.00&101.00&101.00&101.00&101.00&101.00&101.00&101.00\\\
&\mbox{SD}&<10^{-4}&<10^{-4}&0.2537&0.3457&0.4842&<10^{-4}&<10^{-4}&<10^{-4}&0.1826&0.1826&0.1826\\\
\hline\cr\\\
\hline\cr\lx@intercol\mbox{Metric}\hfil\lx@intercol&\lx@intercol\hfil\mbox{$\kappa$
; ZDT1}\hfil\lx@intercol\\\
\cline{3-13}\cr&&\mbox{0.0}&\mbox{0.1}&\mbox{0.2}&\mbox{0.3}&\mbox{0.4}&\mbox{0.5}&\mbox{0.6}&\mbox{0.7}&\mbox{0.8}&\mbox{0.9}&\mbox{0.99}\\\
\hline\cr\mbox{$\gamma$}&\mbox{M}&0.0103&0.0069&0.0055&0.0059&0.0074&0.0040&0.0042&0.0051&0.0053&0.0050&0.0038\\\
&\mbox{SD}&0.0072&0.0038&0.0057&0.0047&0.0049&0.0042&0.0060&0.0058&0.0040&0.0050&0.0034\\\
\mbox{$\Delta$}&\mbox{M}&0.3046&0.5543&0.3678&0.2084&0.1818&0.1812&0.1898&0.1937&0.1934&0.1956&0.1891\\\
&\mbox{SD}&0.0122&0.0607&0.0480&0.0408&0.0235&0.0220&0.0234&0.0251&0.0240&0.0232&0.0155\\\
\mbox{HV}&\mbox{M}&0.8435&0.8561&0.8602&0.8607&0.8598&0.8652&0.8650&0.8636&0.8633&0.8638&0.8657\\\
&\mbox{SD}&0.0115&0.0066&0.0094&0.0079&0.0082&0.0069&0.0099&0.0096&0.0066&0.0083&0.0057\\\
\mbox{$\left|\mathcal{H}\right|$}&\mbox{M}&37.00&63.00&84.50&98.00&100.00&101.00&101.00&101.00&101.00&101.00&101.00\\\
&\mbox{SD}&0.6397&5.7211&2.8730&5.0901&3.8201&0.5467&0.8584&0.9371&0.9377&1.3515&0.7112\\\
\hline\cr\\\
\hline\cr\lx@intercol\mbox{Metric}\hfil\lx@intercol&\lx@intercol\hfil\mbox{$\kappa$
; DTLZ2 (m = 3)}\hfil\lx@intercol\\\
\cline{3-13}\cr&&\mbox{0.0}&\mbox{0.1}&\mbox{0.2}&\mbox{0.3}&\mbox{0.4}&\mbox{0.5}&\mbox{0.6}&\mbox{0.7}&\mbox{0.8}&\mbox{0.9}&\mbox{0.99}\\\
\hline\cr\mbox{$\gamma$}&\mbox{M}&0.0062&0.0069&0.0072&0.0070&0.0074&0.0079&0.0074&0.0076&0.0074&0.0078&0.0072\\\
&\mbox{SD}&0.0002&0.0013&0.0015&0.0013&0.0012&0.0014&0.0010&0.0019&0.0007&0.0014&0.0009\\\
\mbox{$\Delta$}&\mbox{M}&0.0503&0.6066&0.3029&0.2411&0.2386&0.2308&0.2274&0.2175&0.2079&0.2173&0.1982\\\
&\mbox{SD}&0.0041&0.0422&0.0357&0.0302&0.0264&0.0219&0.0316&0.0275&0.0306&0.0295&0.0239\\\
\mbox{HV}&\mbox{M}&0.6731&0.7149&0.7383&0.7435&0.7458&0.7469&0.7467&0.7469&0.7470&0.7470&0.7471\\\
&\mbox{SD}&0.0066&0.0042&0.0023&0.0012&0.0007&0.0006&0.0005&0.0005&0.0005&0.0005&0.0003\\\
\mbox{$\left|\mathcal{H}\right|$}&\mbox{M}&21.00&69.00&88.00&93.00&94.50&95.00&95.00&95.00&95.00&95.00&94.00\\\
&\mbox{SD}&1.3047&3.1639&2.8367&2.0424&1.7750&1.9464&2.2894&2.0197&1.5643&2.2614&1.7100\\\
\hline\cr\\\
\hline\cr\lx@intercol\mbox{Metric}\hfil\lx@intercol&\lx@intercol\hfil\mbox{$\kappa$
; DTLZ2 (m = 4)}\hfil\lx@intercol\\\
\cline{3-13}\cr&&\mbox{0.0}&\mbox{0.1}&\mbox{0.2}&\mbox{0.3}&\mbox{0.4}&\mbox{0.5}&\mbox{0.6}&\mbox{0.7}&\mbox{0.8}&\mbox{0.9}&\mbox{0.99}\\\
\hline\cr\mbox{$\gamma$}&\mbox{M}&0.0001&0.0311&0.0385&0.0312&0.0449&0.0404&0.0445&0.0488&0.0590&0.0489&0.0534\\\
&\mbox{SD}&0.0001&0.0198&0.0218&0.0239&0.0283&0.0240&0.0369&0.0281&0.0284&0.0241&0.0304\\\
\mbox{$\Delta$}&\mbox{M}&0.1390&0.4700&0.3602&0.3296&0.3299&0.3429&0.3377&0.3258&0.3253&0.3304&0.3319\\\
&\mbox{SD}&0.1173&0.0304&0.0307&0.0226&0.0255&0.0263&0.0187&0.0262&0.0210&0.0254&0.0259\\\
\mbox{$\left|\mathcal{H}\right|$}&\mbox{M}&14.00&79.50&90.00&92.00&95.00&96.00&95.50&97.00&98.00&95.50&97.00\\\
&\mbox{SD}&1.9815&4.9642&4.8476&4.2372&4.6307&4.2129&5.8530&5.0496&4.5945&4.6233&4.3423\\\
\hline\cr\end{array}$

(a) cone$\epsilon$-MOEA distribution on test Deb52:
$\boldsymbol{\epsilon}=(0.0535,0.0645)$, $\left|\mathcal{H}\right|=31$,
$T=30$. (b) $\epsilon$-MOEA distribution on test Deb52:
$\boldsymbol{\epsilon}=(0.0277,0.0333)$, $\left|\mathcal{H}\right|=12$,
$T=30$. (c) cone$\epsilon$-MOEA distribution on test Pol:
$\boldsymbol{\epsilon}=(0.6180,0.9804)$, $\left|\mathcal{H}\right|=38$,
$T=50$. (d) $\epsilon$-MOEA distribution on test Pol:
$\boldsymbol{\epsilon}=(0.3151,0.5000)$, $\left|\mathcal{H}\right|=09$,
$T=50$.

Figure 7: Efficient solutions generated by cone$\epsilon$-MOEA and
$\epsilon$-MOEA considering calculated $\boldsymbol{\epsilon}$ values to
provide a maximum of $T$ solutions in the final archive. The number of
solutions approximated by each approach is indicated by the cardinality of the
final Pareto front $\left|\mathcal{H}\right|$. The presented fronts are the
outcome of a typical run.

###  6.4 Statistical Design

To evaluate the possible differences in the performance of the methods tested,
we have employed tests designed to detect statistically significant
differences and to estimate their sizes. The tests were performed
independently for each of the quality metrics described previously, and are
outlined in the following paragraphs.

The raw dataset used for this analysis was composed of the final Pareto-front
obtained on $50$ independent runs of each algorithm on each problem, from
which the quality criteria described earlier were calculated. The statistical
tests were then performed using these calculated values.

In order to detect whether any of the six algorithms presented a significant
overall difference (i.e., when integrated over the whole test set), we
employed a randomized complete block design (Montgomery,, 2008), with the
algorithms as levels of the experimental factor Algorithm, and the problems as
blocking levels. The statistical model considered by this experimental design
can be defined, e.g., for the $\Delta$ metric, as:

$y^{\Delta}_{ijk}=\mu^{\Delta}+\tau_{i}^{\Delta}+\beta_{j}^{\Delta}+\epsilon_{ijk}^{\Delta}$
(21)

in which $y^{\Delta}_{ijk}$ represents the value of this metric on the
$k^{th}$ run of the $i^{th}$ algorithm on the $j^{th}$ problem; $\mu^{\Delta}$
is the mean of $\Delta$ over all observations; $\tau_{i}^{\Delta}$ is the
contribution due to the $i^{th}$ algorithm; $\beta_{j}^{\Delta}$ is the
component due to the $j^{th}$ problem (block); and $\epsilon_{ijk}^{\Delta}$
is the residual.

For each metric considered, the ANOVA was applied222In all tests reported, the
assumptions of the parametric ANOVA were verified by means of residual
analysis. While small deviations from the normality assumption were detected,
these were deemed unimportant, since the test is robust to moderate departures
from normality in balanced designs with large sample sizes (Montgomery,,
2008). to test the null hypothesis of absence of difference among the six
algorithms evaluated over all problems against it two-sided alternative, i.e.
(again, using the $\Delta$ metric as an example):

$\begin{split}\begin{cases}H_{0}^{\Delta}:\tau_{i}^{\Delta}=0,\ \forall
i\vskip 6.0pt plus 2.0pt minus 2.0pt\\\ H_{1}^{\Delta}:\tau_{i}^{\Delta}\neq
0\ \mbox{for any }i\end{cases}\end{split}$ (22)

After the significance tests, estimators of the block effects were obtained
(Montgomery,, 2008) and subtracted from the samples, thus allowing for
problem-independent pairwise comparisons of the algorithms. Bootstrap
estimations of the distribution of the means (Davison and Hinkley,, 1997) for
each algorithm on each metric were then obtained, which enabled the use of
simple two-sample t-tests and parametric estimators of the differences of
means (Montgomery and Runger,, 2006) for the post-hoc pairwise comparisons.
These comparisons allow not only the determination of statistical
significance, but also the estimation of the magnitude of the differences,
which is important in assessing the practical significance of results (Ellis,,
2010).

Besides the overall analysis performed using the design above, statistical
tests were also performed for each problem independently, in order to assess
the relative performance of the algorithms on each of the test problems. For
each problem and performance metric, the null hypotheses used were those of
equality of mean values against two-sided alternative hypotheses. An one-way
ANOVA design was used in these tests, and effect size estimators were also
calculated using bootstrap, allowing post-hoc pairwise comparisons between
each pair of algorithms. In all cases, the significance level used was $99\%$.

##  7 Results and Discussion

The results obtained for the four performance metrics considered, using the
estimated values of $\boldsymbol{\epsilon}$, are summarized in this section.
The analysis of the results is divided in four parts. In the first we present
an overall analysis of the results, based on the observed general performance
of the algorithms tested on the four performance criteria defined in the
previous section. The second part consists of a detailed analysis of the
behavior of each algorithm for the individual problems, in terms of their
known properties and geometrical characteristics of the Pareto fronts. The
third part consists of a cardinality analysis regarding the number of
solutions estimated by both $\epsilon$ and cone$\epsilon$-approaches, and
their effects on the performance of the methods. Finally, possible extensions
are dicussed in the fourth part.

###  7.1 Overall Analysis

The results obtained in the overall analysis are summarized in Tables 5 and 6,
and on Fig. 8. From these results, a few considerations can be made over the
relative performance of the six algorithms on the test set used.

Table 5: Overall pairwise differences. The numbers indicate the point
estimator for the difference between the row algorithm and the column
algorithm (after removing problem effects). Differences marked with an
asterisk indicate significance at the $99\%$ familywise confidence level,
adjusted using Tukey’ HSD procedure (Montgomery and Runger,, 2006).

$\mathbf{\Delta}$ metric
---
Algorithm | NSGA-II |  $\epsilon$-MOEA | cone$\epsilon$-MOEA | C-NSGA-II | SPEA2
$\epsilon$-MOEA | -0.03* |  |  |  |
cone$\epsilon$-MOEA | -0.20* | -0.17* |  |  |
C-NSGA-II | -0.23* | -0.19* | -0.02* |  |
SPEA2 | -0.36* | -0.32* | -0.15* | -0.13* |
NSGA-II* | -0.08* | -0.04* | 0.12* | 0.15* | 0.28*

Ordering∗∗: | SPEA2 C-NSGA-II cone$\epsilon$-MOEA NSGA-II* $\epsilon$-MOEA NSGA-II
---|---

$\mathbf{\gamma}$ metric
---
Algorithm | NSGA-II |  $\epsilon$-MOEA | cone$\epsilon$-MOEA | C-NSGA-II | SPEA2
$\epsilon$-MOEA | -0.06* |  |  |  |
cone$\epsilon$-MOEA | -4e-3 | 0.06* |  |  |
C-NSGA-II | 0.06* | 0.12* | 0.07* |  |
SPEA2 | 0.04* | 0.10* | 0.04* | -0.02 |
NSGA-II* | 6e-3 | 0.07* | 0.01 | -0.06* | -0.03*

Ordering∗∗: |  $\epsilon$-MOEA $\overline{\mbox{cone$\epsilon$-MOEA\ \ NSGA-II\ \ NSGA-II*}}$ $\overline{\mbox{SPEA2\ \ C-NSGA-II}}$
---|---

CS metric
---
Algorithm | NSGA-II |  $\epsilon$-MOEA | cone$\epsilon$-MOEA | C-NSGA-II | SPEA2
$\epsilon$-MOEA | 1e-3 |  |  |  |
cone$\epsilon$-MOEA | 0.01 | 9e-3 |  |  |
C-NSGA-II | -0.04* | -0.04* | -0.05* |  |
SPEA2 | -0.05* | -0.06* | -0.06* | -0.01 |
NSGA-II* | -0.04* | -0.04* | -0.05* | -3e-3 | 9e-3

Ordering∗∗: |  $\overline{\mbox{cone$\epsilon$-MOEA\ \ $\epsilon$-MOEA\ \ NSGA-II}}$ $\overline{\mbox{C-NSGA-II\ \ NSGA-II*\ \ SPEA2}}$
---|---

HV metric
---
Algorithm | NSGA-II |  $\epsilon$-MOEA | cone$\epsilon$-MOEA | C-NSGA-II | SPEA2
$\epsilon$-MOEA | -4e-3 |  |  |  |
cone$\epsilon$-MOEA | 2e-3 | 6e-3 |  |  |
C-NSGA-II | -9e-3 | -5e-3 | -0.01* |  |
SPEA2 | -8e-3 | -4e-3 | -0.01* | 9e-4 |
NSGA-II* | -0.01* | -6e-3 | -0.01* | -1e-3 | -2e-3

Ordering∗∗: | $\overline{\mbox{cone$\epsilon$-MOEA\ }\underline{\mbox{\ NSGA-II\ }\overline{\mbox{\ $\epsilon$-MOEA\ }}}}\overline{\underline{\mbox{\ SPEA2\ \ C-NSGA-II\ }}\mbox{\ NSGA-II*}}$
---|---

** Ordering from best to worst. Horizontal bars indicate the absence of
statistically significant differences.

(a) $\Delta$ metric. (b) $\gamma$ metric. (c) $HV$ metric. (d) $CS$ metric.

Figure 8: Individual $99\%$ confidence intervals for the estimated average
performance of the six algorithms on each metric. The dotted line marks the
overall mean result for the metric.

First of all, the diversity metric $\Delta$ presented an interesting ordering
of performance, with the SPEA2 presenting the best performance, the C-NSGA-II
in second ($+0.13$)333The numbers in parentheses represent the difference in
performance between a given algorithm and the best average result., the
cone$\epsilon$-MOEA in third ($+0.15$), the NSGA-II* in fourth ($+0.28$), the
$\epsilon$-MOEA in fifth ($+0.32$) and the NSGA-II ranking last ($+0.36$).
Table 5 shows that the magnitude of the differences for this metric are
considerable from a practical point of view, with an absolute difference of
$36\%$ between the best (SPEA2) and worst (NSGA-II) algorithms. The difference
magnitudes reported in Table 5 also suggest that even though the
cone$\epsilon$-MOEA has occupied the third position in this rank, its
difference from the C-NSGA-II is of little effect from a practical point of
view. It is also relevant to have in mind that the computational complexity of
the cone$\epsilon$ approach is significantly inferior to those of the two
first ones (SPEA2 and C-NSGA-II). Figure 8 and Table 6 provide the estimated
absolute values and confidence intervals for each algorithm, which provide the
necessary information for the derivation of the significance-based ranks given
in Table 6.

The performance ordering of the convergence metric $\gamma$ shows that the
$\epsilon$-MOEA was able to significantly outperform all other approaches in
this metric, with the cone$\epsilon$-MOEA ($+0.06$), the NSGA-II ($+0.06$) and
the NSGA-II* ($+0.07$) tied in second, and the SPEA2 ($+0.10$) and C-NSGA-II
($+0.12$) tied in fifth. The effect sizes in this metric are also significant,
with an absolute difference between the best ($\epsilon$-MOEA) and the worst
(C-NSGA-II and SPEA2) of around $12\%$. Figure 8 and Tables 5 and 6 again
provide detailed information on the magnitude of the differences and their
statistical significance.

The generalized coverage metric $CS$ presented an ordering with the
cone$\epsilon$-MOEA, $\epsilon$-MOEA ($-9e-3$), and NSGA-II ($-0.01$) tied in
first, and the C-NSGA-II ($-0.05$), NSGA-II* ($-0.05$), and SPEA2 ($-0.06$)
tied in fourth. Here the magnitudes of the effects were relatively modest,
with an absolute difference of about $6\%$ between the best and worst
performing algorithms. Given the relatively low absolute average values for
this metric (see Table 6), it is reasonable to suppose that in most cases all
algorithms were finding at least some solutions very close to the true Pareto
front, which would contribute to lowering the value of this metric for all
others. However, the lack of other works to use as a comparison baseline for
this metric makes it difficult to make any conclusive statement about the
practical significance of these effect sizes.

Finally, the hypervolume metric $HV$ had a somewhat more complex ordering
structure, with no statistically significant difference among the
cone$\epsilon$-MOEA, NSGA-II ($-2e-3$), and $\epsilon$-MOEA ($-6e-3$); among
the NSGA-II, $\epsilon$-MOEA, SPEA2 ($-0.01$) and C-NSGA-II ($-0.01$)
algorithms; and also among the $\epsilon$-MOEA, SPEA2, C-NSGA-II, and NSGA-II*
($-0.01$) methods. The horizontal bars at the bottom of the HV section of
Table 5 illustrate this ordering of algorithms. In any case, the magnitudes of
the differences were of very little practical significance, with an average
absolute difference of only about $1\%$ between the best and worst-ranking
algorithms. It is also important to notice here that, as explained in the next
section, this quality metric is not an adequate indicator to be used with
problems such as DTLZ8 and DTLZ9, and was not calculated for these two
problems.

Table 6: Aggregation of overall results. The values indicate the estimated global mean and standard error for each algorithm (after removing problem effects), and the value in parentheses is the algorithm rank within each metric. Ties were treated using the mean rank method, so that the sum of ranks is always the same for each metric. The final ranking is then calculated as an ordering of the average ranks. Algorithm | $\Delta$ | $\gamma$ | HV | CS | Avg. Rank | Final Rank
---|---|---|---|---|---|---
NSGA-II | 0.64$\pm$4e-3 | 0.40$\pm$6e-3 | 0.95$\pm$2e-3 | 0.14$\pm$5e-3 | 3.5 | 3
(6) | (3) | (3) | (2)
$\epsilon$-MOEA | 0.61$\pm$8e-3 | 0.34$\pm$7e-3 | 0.94$\pm$2e-3 | 0.14$\pm$5e-3 | 2.88 | 2
(5) | (1) | (3.5) | (2)
cone$\epsilon$-MOEA | 0.44$\pm$5e-3 | 0.40$\pm$6e-3 | 0.95$\pm$2e-3 | 0.15$\pm$4e-3 | 2.5 | 1
(3) | (3) | (2) | (2)
C-NSGA-II | 0.42$\pm$5e-3 | 0.47$\pm$6e-3 | 0.94$\pm$2e-3 | 0.10$\pm$4e-3 | 4.13 | 5.5
(2) | (5.5) | (4) | (5)
SPEA2 | 0.29$\pm$4e-3 | 0.44$\pm$6e-3 | 0.94$\pm$2e-3 | 0.09$\pm$4e-3 | 3.88 | 4
(1) | (5.5) | (4) | (5)
NSGA-II* | 0.57$\pm$5e-3 | 0.41$\pm$5e-3 | 0.94$\pm$2e-3 | 0.10$\pm$4e-3 | 4.13 | 5.5
(4) | (3) | (4.5) | (5)

Figure 9: Compromise between convergence ($\gamma$) and diversity ($\Delta$)
for each algorithm. The bars represent $99\%$ confidence intervals.

Table 6 shows a rank analysis of the algorithms, in which the results of the
four metrics are ranked and aggregated using a simple averaging criterion, and
an overall rank is calculated for the algorithms. No statistical analysis of
significance was possible for the rank analysis, due to the non-independence
of the four performance metrics, which could bias the results of most
inferential procedures. However, the final ranks suggest that the
cone$\epsilon$-MOEA was capable of presenting an efficient and balanced
performance over all the performance metrics considered. The $\epsilon$-MOEA
comes in a close second, mainly due to its poor performance on the $\Delta$
metric. The NSGA-II ranked third, also mainly as a consequence of its very
poor diversity performance. The SPEA2 was ranked fourth, and the C-NSGA-II and
NSGA-II* were tied in the last place. From the observation of their
performance rankings on all metrics other than the $\Delta$, it is possible to
observe that both the SPEA2 and the C-NSGA-II seem to have attained good
diversity performance at the cost of their convergence properties. As for the
NSGA-II*, the modification of the crowding distance procedure seems to have
indeed improved its diversity when compared against the original NSGA-II, but
not enough to be competitive with the cone$\epsilon$-MOEA, C-NSGA-II or SPEA2
algorithms.

To conclude this overall analysis, Fig. 9 illustrates the compromise between
the convergence and diversity obtained for each algorithm. From this figure,
it is possible to notice that the best tradeoffs with respect to these two
performance indicators are obtained by the SPEA2, cone$\epsilon$-MOEA and
$\epsilon$-MOEA. As expected from the data shown in Table 6, the
cone$\epsilon$-MOEA is shown to provide an efficient balance between
convergence and diversity, which makes it an interesting and competitive
approach for the solutions of MOPs in which a good sampling of the Pareto
optimal front is desired.

###  7.2 Problemwise Analysis

The discussion presented in this section consists of an analysis of the
behavior of each algorithm for the individual problems. The results obtained
for each problem are summarized in Tables 7 and 8.

Table 7: Problemwise comparison of the algorithms on the four performance
metrics used, for problems Deb52, Pol and the ZDT family. The values reported
represent the mean and standard error obtained for each combination of
algorithm, problem and performance metric.

Problem | Metric | NSGA-II |  $\epsilon$-MOEA | cone$\epsilon$-MOEA | C-NSGA-II | SPEA2 | NSGA-II*
---|---|---|---|---|---|---|---
Deb52 | $\gamma$ | 0.58$\pm$7e-3 | 0.56$\pm$3e-3 | 0.56$\pm$2e-3 | 0.58$\pm$1e-2 | 0.55$\pm$7e-3 | 0.55$\pm$6e-3
$\Delta$ | 0.53$\pm$8e-3 | 0.99$\pm$6e-4 | 0.32$\pm$3e-4 | 0.33$\pm$6e-3 | 0.20$\pm$3e-3 | 0.41$\pm$8e-3
HV | 0.99$\pm$7e-5 | 0.99$\pm$6e-5 | 0.99$\pm$0 | 0.99$\pm$1e-4 | 0.99$\pm$3e-5 | 0.99$\pm$7e-5
CS | 0.02$\pm$7e-4 | 0.03$\pm$10e-4 | 0.03$\pm$8e-4 | 0.02$\pm$8e-4 | 0.03$\pm$9e-4 | 0.02$\pm$8e-4
Pol | $\gamma$ | 0.20$\pm$2e-2 | 0.13$\pm$6e-4 | 0.19$\pm$2e-2 | 0.19$\pm$9e-3 | 0.15$\pm$2e-3 | 0.16$\pm$2e-3
$\Delta$ | 0.58$\pm$1e-2 | 0.98$\pm$9e-4 | 0.29$\pm$6e-3 | 0.38$\pm$7e-3 | 0.24$\pm$3e-3 | 0.36$\pm$8e-3
HV | 1.00$\pm$1e-5 | 1.00$\pm$5e-6 | 1.00$\pm$4e-6 | 1.00$\pm$3e-5 | 1.00$\pm$4e-6 | 1.00$\pm$8e-6
CS | 0.04$\pm$1e-3 | 0.04$\pm$1e-3 | 0.04$\pm$2e-3 | 0.04$\pm$10e-4 | 0.06$\pm$1e-3 | 0.05$\pm$1e-3
Zdt1 | $\gamma$ | 0.16$\pm$2e-2 | 0.30$\pm$3e-2 | 0.23$\pm$3e-2 | 0.18$\pm$2e-2 | 0.30$\pm$2e-2 | 0.19$\pm$2e-2
$\Delta$ | 0.79$\pm$1e-2 | 0.70$\pm$4e-3 | 0.37$\pm$6e-3 | 0.50$\pm$8e-3 | 0.29$\pm$7e-3 | 0.56$\pm$1e-2
HV | 0.99$\pm$10e-4 | 0.98$\pm$2e-3 | 0.98$\pm$2e-3 | 0.98$\pm$9e-4 | 0.98$\pm$1e-3 | 0.98$\pm$1e-3
CS | 0.33$\pm$2e-2 | 0.20$\pm$3e-2 | 0.27$\pm$3e-2 | 0.30$\pm$2e-2 | 0.15$\pm$2e-2 | 0.27$\pm$2e-2
Zdt2 | $\gamma$ | 0.43$\pm$9e-3 | 0.65$\pm$8e-3 | 0.30$\pm$4e-3 | 0.80$\pm$1e-2 | 0.41$\pm$8e-3 | 0.41$\pm$8e-3
$\Delta$ | 0.76$\pm$1e-2 | 0.56$\pm$3e-3 | 0.38$\pm$4e-3 | 0.50$\pm$7e-3 | 0.28$\pm$4e-3 | 0.58$\pm$1e-2
HV | 0.99$\pm$1e-4 | 0.99$\pm$7e-5 | 0.99$\pm$4e-5 | 0.98$\pm$1e-4 | 0.99$\pm$8e-5 | 0.99$\pm$9e-5
CS | 0.07$\pm$3e-3 | 0.04$\pm$2e-3 | 0.13$\pm$3e-3 | 0.01$\pm$6e-4 | 0.08$\pm$3e-3 | 0.07$\pm$3e-3
Zdt3 | $\gamma$ | 0.16$\pm$2e-2 | 0.15$\pm$1e-2 | 0.17$\pm$10e-3 | 0.19$\pm$1e-2 | 0.35$\pm$3e-2 | 0.24$\pm$2e-2
$\Delta$ | 0.67$\pm$1e-2 | 0.85$\pm$1e-2 | 0.57$\pm$2e-2 | 0.49$\pm$1e-2 | 0.33$\pm$7e-3 | 0.50$\pm$2e-2
HV | 0.98$\pm$2e-3 | 0.97$\pm$3e-3 | 0.97$\pm$2e-3 | 0.97$\pm$2e-3 | 0.95$\pm$4e-3 | 0.97$\pm$3e-3
CS | 0.41$\pm$3e-2 | 0.36$\pm$3e-2 | 0.33$\pm$2e-2 | 0.29$\pm$2e-2 | 0.12$\pm$2e-2 | 0.24$\pm$3e-2
Zdt4 | $\gamma$ | 0.27$\pm$2e-2 | 0.32$\pm$2e-2 | 0.35$\pm$2e-2 | 0.47$\pm$2e-2 | 0.69$\pm$2e-2 | 0.57$\pm$2e-2
$\Delta$ | 0.61$\pm$9e-3 | 0.59$\pm$1e-2 | 0.48$\pm$1e-2 | 0.58$\pm$1e-2 | 0.53$\pm$1e-2 | 0.68$\pm$1e-2
HV | 0.84$\pm$9e-3 | 0.80$\pm$1e-2 | 0.79$\pm$1e-2 | 0.73$\pm$1e-2 | 0.62$\pm$1e-2 | 0.66$\pm$2e-2
CS | 0.73$\pm$3e-2 | 0.58$\pm$4e-2 | 0.57$\pm$3e-2 | 0.41$\pm$4e-2 | 0.15$\pm$2e-2 | 0.21$\pm$4e-2
Zdt6 | $\gamma$ | 0.06$\pm$3e-2 | 0.04$\pm$1e-2 | 0.04$\pm$2e-2 | 0.04$\pm$1e-2 | 0.01$\pm$3e-3 | 0.03$\pm$2e-2
$\Delta$ | 0.52$\pm$1e-2 | 0.24$\pm$1e-2 | 0.27$\pm$1e-2 | 0.37$\pm$9e-3 | 0.18$\pm$8e-3 | 0.37$\pm$2e-2
HV | 0.96$\pm$1e-2 | 0.98$\pm$7e-3 | 0.98$\pm$10e-3 | 0.97$\pm$7e-3 | 0.99$\pm$1e-3 | 0.98$\pm$10e-3
CS | 0.18$\pm$2e-2 | 0.11$\pm$2e-2 | 0.18$\pm$2e-2 | 0.10$\pm$1e-2 | 0.21$\pm$2e-2 | 0.21$\pm$2e-2

Table 8: Problemwise comparison of the algorithms on the four performance
metrics used, for the DTLZ family. The values reported represent the mean and
standard error obtained for each combination of algorithm, problem and
performance metric.

Problem | Metric | NSGA-II |  $\epsilon$-MOEA | cone$\epsilon$-MOEA | C-NSGA-II | SPEA2 | NSGA-II*
---|---|---|---|---|---|---|---
Dtlz1 | $\gamma$ | 0.23$\pm$7e-3 | 0.13$\pm$1e-3 | 0.17$\pm$2e-2 | 0.39$\pm$1e-2 | 0.18$\pm$2e-3 | 0.17$\pm$2e-3
$\Delta$ | 0.34$\pm$3e-3 | 0.12$\pm$2e-3 | 0.05$\pm$1e-2 | 0.20$\pm$2e-2 | 0.08$\pm$1e-3 | 0.34$\pm$4e-3
HV | 0.95$\pm$6e-4 | 0.92$\pm$4e-4 | 0.95$\pm$2e-4 | 0.96$\pm$5e-4 | 0.97$\pm$1e-4 | 0.96$\pm$5e-4
CS | 0.02$\pm$8e-4 | 0.01$\pm$9e-4 | 0.03$\pm$1e-3 | 0.00$\pm$5e-4 | 0.02$\pm$1e-3 | 0.02$\pm$10e-4
Dtlz2 | $\gamma$ | 0.62$\pm$10e-3 | 0.70$\pm$7e-3 | 0.48$\pm$8e-3 | 0.75$\pm$1e-2 | 0.55$\pm$8e-3 | 0.48$\pm$6e-3
$\Delta$ | 0.81$\pm$10e-3 | 0.42$\pm$4e-3 | 0.42$\pm$6e-3 | 0.29$\pm$5e-3 | 0.16$\pm$3e-3 | 0.83$\pm$9e-3
HV | 0.89$\pm$9e-4 | 0.92$\pm$4e-4 | 0.94$\pm$1e-4 | 0.90$\pm$7e-4 | 0.93$\pm$4e-4 | 0.89$\pm$9e-4
CS | 0.03$\pm$1e-3 | 0.02$\pm$8e-4 | 0.06$\pm$2e-3 | 0.01$\pm$6e-4 | 0.03$\pm$1e-3 | 0.04$\pm$1e-3
Dtlz3 | $\gamma$ | 0.35$\pm$2e-2 | 0.25$\pm$1e-2 | 0.42$\pm$3e-2 | 0.50$\pm$3e-2 | 0.32$\pm$2e-2 | 0.26$\pm$1e-2
$\Delta$ | 0.33$\pm$9e-3 | 0.19$\pm$1e-2 | 0.29$\pm$2e-2 | 0.22$\pm$3e-2 | 0.15$\pm$2e-2 | 0.34$\pm$8e-3
HV | 0.90$\pm$10e-4 | 0.91$\pm$2e-2 | 0.91$\pm$1e-2 | 0.91$\pm$1e-3 | 0.93$\pm$3e-4 | 0.90$\pm$9e-4
CS | 0.02$\pm$1e-3 | 0.03$\pm$2e-3 | 0.04$\pm$2e-3 | 0.00$\pm$5e-4 | 0.02$\pm$2e-3 | 0.02$\pm$2e-3
Dtlz4 | $\gamma$ | 0.32$\pm$8e-3 | 0.41$\pm$2e-2 | 0.53$\pm$3e-2 | 0.34$\pm$1e-2 | 0.33$\pm$1e-2 | 0.30$\pm$3e-3
$\Delta$ | 0.67$\pm$2e-2 | 0.37$\pm$3e-2 | 0.43$\pm$2e-2 | 0.36$\pm$4e-2 | 0.22$\pm$3e-2 | 0.66$\pm$8e-3
HV | 0.88$\pm$1e-2 | 0.86$\pm$2e-2 | 0.86$\pm$2e-2 | 0.84$\pm$2e-2 | 0.87$\pm$2e-2 | 0.90$\pm$7e-4
CS | 0.04$\pm$2e-3 | 0.02$\pm$1e-3 | 0.03$\pm$2e-3 | 0.02$\pm$2e-3 | 0.03$\pm$2e-3 | 0.03$\pm$1e-3
Dtlz5 | $\gamma$ | 0.14$\pm$2e-3 | 0.26$\pm$3e-3 | 0.56$\pm$3e-2 | 0.22$\pm$5e-3 | 0.15$\pm$2e-3 | 0.13$\pm$1e-3
$\Delta$ | 0.74$\pm$2e-2 | 0.78$\pm$5e-3 | 0.83$\pm$9e-3 | 0.43$\pm$6e-3 | 0.26$\pm$4e-3 | 0.61$\pm$1e-2
HV | 0.99$\pm$1e-4 | 0.99$\pm$7e-5 | 0.98$\pm$3e-5 | 0.98$\pm$1e-4 | 0.99$\pm$4e-5 | 0.99$\pm$1e-4
CS | 0.06$\pm$2e-3 | 0.02$\pm$1e-3 | 0.04$\pm$1e-3 | 0.03$\pm$1e-3 | 0.05$\pm$2e-3 | 0.07$\pm$2e-3
Dtlz6 | $\gamma$ | 0.84$\pm$8e-3 | 0.94$\pm$3e-3 | 0.83$\pm$3e-3 | 0.83$\pm$5e-3 | 0.84$\pm$7e-3 | 0.84$\pm$8e-3
$\Delta$ | 0.78$\pm$1e-2 | 0.99$\pm$6e-4 | 0.45$\pm$6e-4 | 0.51$\pm$7e-3 | 0.32$\pm$5e-3 | 0.62$\pm$1e-2
HV | 0.99$\pm$1e-4 | 0.99$\pm$4e-5 | 0.99$\pm$2e-5 | 0.99$\pm$9e-5 | 0.99$\pm$3e-5 | 0.99$\pm$1e-4
CS | 0.02$\pm$6e-4 | 0.03$\pm$8e-4 | 0.04$\pm$8e-4 | 0.03$\pm$8e-4 | 0.03$\pm$8e-4 | 0.01$\pm$7e-4
Dtlz7 | $\gamma$ | 0.74$\pm$9e-3 | 0.47$\pm$2e-3 | 0.74$\pm$10e-3 | 0.86$\pm$9e-3 | 0.67$\pm$8e-3 | 0.74$\pm$10e-3
$\Delta$ | 0.78$\pm$1e-2 | 0.55$\pm$9e-3 | 0.71$\pm$9e-3 | 0.56$\pm$1e-2 | 0.52$\pm$8e-3 | 0.78$\pm$8e-3
HV | 0.92$\pm$1e-3 | 0.91$\pm$1e-3 | 0.93$\pm$7e-4 | 0.91$\pm$1e-3 | 0.94$\pm$6e-4 | 0.92$\pm$9e-4
CS | 0.06$\pm$2e-3 | 0.08$\pm$2e-3 | 0.06$\pm$2e-3 | 0.02$\pm$1e-3 | 0.09$\pm$3e-3 | 0.07$\pm$2e-3
Dtlz8 | $\gamma$ | 0.45$\pm$2e-2 | 0.03$\pm$10e-4 | 0.24$\pm$5e-3 | 0.52$\pm$1e-2 | 0.83$\pm$9e-3 | 0.58$\pm$1e-2
$\Delta$ | 0.82$\pm$9e-3 | 0.69$\pm$9e-3 | 0.70$\pm$9e-3 | 0.49$\pm$8e-3 | 0.36$\pm$8e-3 | 0.80$\pm$1e-2
CS | 0.11$\pm$4e-3 | 0.20$\pm$3e-3 | 0.24$\pm$4e-3 | 0.10$\pm$3e-3 | 0.10$\pm$3e-3 | 0.09$\pm$3e-3
Dtlz9 | $\gamma$ | 0.90$\pm$5e-3 | 0.16$\pm$7e-3 | 0.57$\pm$5e-3 | 0.64$\pm$7e-3 | 0.77$\pm$7e-3 | 0.90$\pm$5e-3
$\Delta$ | 0.62$\pm$5e-3 | 0.71$\pm$1e-2 | 0.52$\pm$6e-3 | 0.49$\pm$6e-3 | 0.45$\pm$4e-3 | 0.61$\pm$6e-3
CS | 0.11$\pm$8e-3 | 0.51$\pm$1e-2 | 0.34$\pm$7e-3 | 0.23$\pm$8e-3 | 0.24$\pm$8e-3 | 0.13$\pm$7e-3

The Deb52 problem tests the ability of the MOEA to find non-convex Pareto-
optimal solutions (Deb,, 1999). As this simple and continuous problem has only
two decision variables, all algorithms presented a similar performance in
terms of the $\gamma$, $HV$ and $CS$ measures. On the other hand, the SPEA2
achieved the best $\Delta$ value, followed by cone$\epsilon$-MOEA, C-NSGA-II,
NSGA-II* and NSGA-II, respectively. Since the optimal frontier of this problem
is represented by a “strong” concave curve, and the $\epsilon$-mechanism is
sensitive to the loss of feasible solutions, the $\epsilon$-MOEA ranked last
in this metric.

The Pol problem, although simple, provides difficulties by introducing
discontinuities in the non-convex Pareto front (Poloni,, 1995). As the
parameter space is also mapped with two variables, very small differences were
observed for the $HV$ and $CS$ indicators. The cone$\epsilon$-MOEA found the
second best $\Delta$ value, and a convergence measure similar to that of the
approaches based on the NSGA-II. As an illustration, Fig. 7 gives a good idea
of the convergence and distribution of the solutions accomplished by both
$\epsilon$-approaches in a particular case.

The ZDT1 30-variable problem has a convex Pareto-optimal front (Zitzler et
al.,, 2000). The NSGA-II approaches achieved the best $\gamma$ and $CS$
values. The cone$\epsilon$-MOEA also ranked well in these metrics. Because of
the cone$\epsilon$-concept, this method achieved the second best spread of
solutions. The SPEA2 obtained the overall best value in this metric. The ZDT2
30-variable problem represents the counterpart of ZDT1, and tests the ability
of the MOEA to approximate non-convex Pareto-optimal solutions (Zitzler et
al.,, 2000). The cone$\epsilon$-MOEA has significantly outperformed the other
methods in the $\gamma$ and $CS$ indicators. Regarding the diversity found,
the SPEA2 and the cone$\epsilon$-MOEA were the two best ones. No significant
differences were observed for the $HV$ indicator in both problems.

As we have observed for the Pol test, the ZDT3 problem provides difficulties
by introducing discontinuities in the objective space, with a Pareto front
consisting of several disconnected convex parts. However, the Pareto-optimal
set in the 30-variable parameter space has no discontinuity (Zitzler et al.,,
2000). The cone$\epsilon$-MOEA ranked third in the $\gamma$ and $CS$ measures.
The $\epsilon$-MOEA and the SPEA2 presented either a good convergence or
diversity: the $\epsilon$-MOEA found the best $\gamma$ value and the worst
$\Delta$ measure; on the other hand, the SPEA2 ranked first for the diversity
measure at the cost of the worst $\gamma$ value.

The ZDT4 10-variable problem contains a number of local Pareto-optimal fronts
and, therefore, tests the ability of the MOEA to deal with multimodality,
providing challenges for a MOEA to converge to the global Pareto front
(Zitzler et al.,, 2000). In this problem, the cone$\epsilon$-MOEA has achieved
the best spread of solutions. On the other hand, the NSGA-II performed better
in terms of convergence, $HV$ and $CS$, followed by the $\epsilon$-MOEA and
cone$\epsilon$-MOEA in these metrics.

In the ZDT6 10-variable problem, the Pareto-optimal solutions are nonuniformly
distributed along the global Pareto front (Zitzler et al.,, 2000).
Essentially, very small differences were observed for both convergence and
hypervolume metrics. The SPEA2 presented the best diversity, followed by the
$\epsilon$-MOEA and the cone$\epsilon$-MOEA. For the $CS$ indicator, the
$\epsilon$-MOEA and the C-NSGA-II ranked last.

The three-objective DTLZ1 test problem has seven decision variables, and the
Pareto-optimal solutions lie on a plane satisfying $\sum_{i=1}^{m}f_{i}=0.5$,
such that $f_{i}\in[0,0.5]$ for all $i=\left\\{1,2,3\right\\}$ (Deb et al.,
2005b, ). Small differences were observed for the algorithms regarding the
$HV$ and $CS$ metrics. The $\epsilon$-MOEA has achieved the best convergence,
followed by the cone$\epsilon$-MOEA and the NSGA-II*, and the
cone$\epsilon$-method has outperformed the other algorithms in finding a well-
distributed set of solutions. The $\epsilon$-MOEA was better than the NSGA-II
-based methods in terms of diversity, possibly due to the use of the
$\epsilon$-dominance, which tends to avoid the clumping of solutions. Figure
10 presents a visual comparison between the outcome of a typical run of each
MOEA considered in this study. It can be seen that the distribution of
solutions obtained by the NSGA-II and NSGA-II* are poorer than the other
MOEAs.

(a) NSGA-II. (b) NSGA-II*. (c) C-NSGA-II. (d) SPEA2. (e) $\epsilon$-MOEA. (f)
cone$\epsilon$-MOEA.

Figure 10: Efficient solutions generated for problem DTLZ1, considering the
estimated $\boldsymbol{\epsilon}$ values provided in Table 3. The fronts
presented are the outcome of a typical run.

The DTLZ2 12-variable test problem has a spherical Pareto-optimal front
satisfying $\sum_{i=1}^{m}f_{i}^{2}=1$, such that $f_{i}\in[0,1]$ for all
$i=\left\\{1,2,3\right\\}$ (Deb et al., 2005b, ). For the hypervolume and
coverage of many sets quality indicators, the use of cone$\epsilon$-dominance
was able to provide a small advantage over the other approaches. Additionally,
the cone$\epsilon$-MOEA and NSGA-II* have achieved the best convergence, and
the clustering techniques employed in SPEA2 and C-NSGA-II provided the overall
best $\Delta$ values. Although no significant difference was found between the
cone$\epsilon$-MOEA and the $\epsilon$-MOEA for the diversity metric, it can
be seen from Fig. 11 that the border regions of the Pareto front were better
mapped by the cone$\epsilon$-method. Notice that, even though these fronts are
the outcome of a single run, both the $\epsilon$ and the
cone$\epsilon$-concepts yielded very similar distributions of points in all
runs. In the case of the $\epsilon$-MOEA, there seems to be a considerable gap
between the boundary solutions and their nearest neighbors. This happens
because of the gentle slope near the boundary solutions on a spherical
surface, and the fact that the $\epsilon$-dominance does not allow any point
to be nondominated within an $\epsilon_{i}$ in the $i^{th}$ objective (Deb et
al.,, 2003; Deb et al., 2005a, ).

(a) NSGA-II. (b) NSGA-II*. (c) C-NSGA-II. (d) SPEA2. (e) $\epsilon$-MOEA. (f)
cone$\epsilon$-MOEA.

Figure 11: Efficient solutions generated for problem DTLZ2, considering the
estimated $\boldsymbol{\epsilon}$ values provided in Table 3. The fronts
presented are the outcome of a typical run.

The DTLZ3 12-variable test problem has a large number of local Pareto-optimal
fronts and investigates the ability of the MOEA to converge to the global
Pareto front. As in the previous problem, the global front satisfies
$\sum_{i=1}^{m}f_{i}^{2}=1$, such that $f_{i}\in[0,1]$ for all
$i=\left\\{1,2,3\right\\}$ (Deb et al., 2005b, ). In this problem, the
cone$\epsilon$-MOEA ranked fourth in the $\Delta$ and fifth in the $\gamma$
metric. The SPEA2 and $\epsilon$-MOEA presented the best overall values in
these indicators, respectively. Small differences were observed for the other
two quality measures.

The DTLZ4 12-variable test problem introduces a nonuniform density of
solutions on the three-objective Pareto-optimal set, and so it investigates
the ability of the MOEA to maintain a good distribution of solutions on the
corresponding front. The global optimal front is the surface represented by
$\sum_{i=1}^{m}f_{i}^{2}=1$, in which $f_{i}\in[0,1]$ for all
$i=\left\\{1,2,3\right\\}$ (Deb et al., 2005b, ). As this problem has a
greater density of solutions near the $f_{1}$–$f_{2}$ and $f_{1}$–$f_{3}$
planes, some runs of all the algorithms achieved solutions only in these
regions. In this problem, the NSGA-II* was the best in terms of convergence,
followed by NSGA-II, SPEA2, C-NSGA-II, and $\epsilon$-MOEA, which presented an
advantage over the cone$\epsilon$-method. Concerning the spread of solutions,
the best methods were the SPEA2 and C-NSGA-II, followed by both
$\epsilon$-approaches. Although there is a difference in $\Delta$ in favor of
the $\epsilon$-MOEA when compared with the cone$\epsilon$-MOEA, it can be seen
from Fig. 12 that the border regions of the Pareto front were better mapped by
the cone$\epsilon$-method. Note again that, even though these fronts are the
outcome of a single run, both the $\epsilon$ and the cone$\epsilon$-concepts
yielded similar distributions of solutions in all runs. Regarding the $HV$ and
$CS$ metrics, only small differences were observed between the algorithms.

(a) NSGA-II. (b) NSGA-II*. (c) C-NSGA-II. (d) SPEA2. (e) $\epsilon$-MOEA. (f)
cone$\epsilon$-MOEA.

Figure 12: Efficient solutions generated for problem DTLZ4, considering the
estimated $\boldsymbol{\epsilon}$ values provided in Table 3. The fronts
presented are the outcome of a typical run.

The three-objective 12-variable DTLZ5 problem has a Pareto-optimal front
satisfying $\sum_{i=1}^{m}f_{i}^{2}=1$, in which $f_{i}\in[0,1]$ for all
$i=\left\\{1,2,3\right\\}$. This problem tests the ability of the MOEA to find
a lower-dimensional Pareto curve while dealing with a higher-dimensional
objective space (Deb et al., 2005b, ). In this problem, the NSGA-II* achieved
the best convergence, closely followed by the NSGA-II and SPEA2. For the
diversity indicator, the clustering-based methods (SPEA2 and C-NSGA-II)
provided the best overall $\Delta$ values. The cone$\epsilon$-MOEA performed
worst in these metrics. The DTLZ6 problem is basically a harder version of the
DTLZ5 (Deb et al., 2005b, ). In this case, the cone$\epsilon$-MOEA presented
the best and the second best performances in finding well-converged and well-
distributed solutions, respectively. Once more, the $HV$ and $CS$ metrics
presented very small differences between the algorithms tested.

The 22-variable DTLZ7 problem has a disconnected set of Pareto-optimal areas
in the search space, and tests the ability of the algorithm to maintain
subpopulations in different Pareto front regions (Deb et al., 2005b, ). The
$\epsilon$-MOEA has outperformed the other methods in the convergence
indicator, and the SPEA2 achieved the best $\Delta$ value, closely followed by
$\epsilon$-MOEA and C-NSGA-II. The cone$\epsilon$-method ranked fourth in this
measure. Again, the $HV$ and $CS$ metrics presented very small differences
between the algorithms, with positive gains in favor of the SPEA2 in both
cases. Figure 13 shows the outcome of a typical run of each algorithm
considered in this study. Since the influence on the ordering of points
performed by the cone of dominance was limited to a local neighborhood in the
objective space when using $\kappa=0.5$, the cone$\epsilon$-method was not
able to estimate a well-spread approximation of the global Pareto front (see
Fig. 13(f)), which may also have adversely affected the convergence of the
approach. Due to this feature, both the cone$\epsilon$-MOEA and the NSGA-II
presented a very similar distribution of solutions on this test. Note,
however, that the performance of the cone$\epsilon$-approach should be
improved by considering a smaller $\kappa$ value, and a more desirable quality
for the convergence and diversity of the solutions could be found.

(a) NSGA-II. (b) NSGA-II*.

(c) C-NSGA-II. (d) SPEA2.

(e) $\epsilon$-MOEA.

(f) cone$\epsilon$-MOEA.

Figure 13: Efficient solutions generated for problem DTLZ7, considering the
estimated $\boldsymbol{\epsilon}$ values provided in Table 3. The fronts
presented are the outcome of a typical run.

In the three-objective 30-variable DTLZ8 problem, the overall Pareto-optimal
front is a combination of a straight line and a hyperplane. This problem
involves $n_{g}=3$ inequality constraints: the straight line is the
intersection of the first $n_{g}-1$ constraints, and the hyperplane is
represented by the last one. MOEAs may find difficulty in achieving solutions
in both the Pareto regions and also in maintaining a good distribution of
solutions on the hyperplane (Deb et al., 2005b, ). In this study, the
constraints were handled using a simple penalty method
($f_{i}+1000\sum_{j=1}^{n_{g}}\max\left(0,g_{j}\right)$) for all the
algorithms. As discussed elsewhere (Deb et al.,, 2003; Deb et al., 2005a, ;
Deb et al., 2005b, ), the domination-based MOEAs suffer, in this particular
problem, from what is known as the “redundancy problem”. In general terms,
even though some solutions are found on the global Pareto front, there exist
many other weakly nondominated solutions in the final archive. In Fig. 14(a),
these “redundant solutions” are those that are on the adjoining surfaces of
the Pareto-optimal front, and their presence in the final nondominated set is
very difficult to eradicate in real-parameter MOEAs (Ikeda et al.,, 2001). Due
to this feature, the obtained set of solutions may wrongly find a higher-
dimensional surface than the Pareto-optimal front, although the global Pareto
front may be of smaller dimension. As a direct consequence, the S-metric
($HV$) is not appropriate for this problem, as the redundant solutions may
affect the hypervolume indicator, and was not considered in this case. Figure
14(e) shows that many of these redundant solutions get $\epsilon$-dominated by
the Pareto-optimal solutions found. In the case of the
cone$\epsilon$-approach, as the cone of dominance presented a local effect in
the objective space when using $\kappa=0.5$, and only few redundant solutions
were removed from the archive (see Fig. 14(f)). Again, it is important to
observe that a global effect can be performed by choosing a smaller $\kappa$
value (e.g, a $\kappa\rightarrow 0$ would achieve the same results as the
$\epsilon$-dominance). From Table 8, it can be seen that the clustering-based
methods achieved the best $\Delta$ values. Although no significant difference
was found between the $\epsilon$-methods in finding well-distributed
solutions, they have outperformed the NSGA-II and NSGA-II* in this indicator.
The $\epsilon$-MOEA found the best result for the $\gamma$ metric, followed by
the cone$\epsilon$-algorithm. The overall best $CS$ measure was obtained by
the cone$\epsilon$-MOEA.

(a) NSGA-II.

(b) NSGA-II*.

(c) C-NSGA-II. (d) SPEA2.

(e) $\epsilon$-MOEA.

(f) cone$\epsilon$-MOEA.

Figure 14: Efficient solutions generated for problem DTLZ8, considering the
estimated $\boldsymbol{\epsilon}$ values provided in Table 3. The fronts
presented are the outcome of a typical run.

The 30-variable DTLZ9 problem is also created using the constraint surface
approach employed by the DTLZ8 test. The Pareto-optimal front is a curve in a
three-dimensional objective space, and the density of solutions gets thinner
towards the Pareto-optimal region. This curve lies on the intersection of all
the inequality constraints. As in the previous problem, the DTLZ9 test
investigates the ability of the MOEA to deal with redundant solutions (Deb et
al., 2005b, ). Again, the constraints were handled using a penalty method
($f_{i}+1000\sum_{j=1}^{n_{g}}\max\left(0,g_{j}\right)$) for all the
algorithms. A two-dimensional plot of the Pareto-optimal curve on the
$f_{1}$–$f_{3}$ plane allows an easier visualization of the solutions obtained
(see Fig. 15). From Table 8, it can be seen that the clustering-based methods
and the cone$\epsilon$-MOEA have presented the best measures in terms of
solution diversity. However, the $\epsilon$-MOEA strongly outperformed all
other approaches in terms of $\gamma$ and $CS$ values, followed by the
cone$\epsilon$-MOEA in these metrics. Once again, the S-metric ($HV$) was not
considered in this problem since redundant solutions may affect the ability of
the hypervolume indicator to adequately quantify the quality of the fronts.

(a) NSGA-II. (b) NSGA-II*.

(c) C-NSGA-II. (d) SPEA2.

(e) $\epsilon$-MOEA.

(f) cone$\epsilon$-MOEA.

Figure 15: Efficient solutions generated for problem DTLZ9, considering the
estimated $\boldsymbol{\epsilon}$ values provided in Table 3. The fronts
presented are the outcome of a typical run.

###  7.3 Cardinality analysis

It is interesting to verify that the estimated values of $\epsilon$ used for
the $\epsilon$-MOEA and cone$\epsilon$-MOEA were able to provide final archive
sizes of approximately 100 solutions. To be more precise, the mean and
standard error for the cardinality of the $\epsilon$-MOEA was $96.6\pm 0.51$,
and for the cone$\epsilon$-MOEA it was $98.2\pm 0.40$. While both approaches
had mean values slightly smaller than the nominal 100 (which was achieved in
all runs by the NSGA-II, NSGA-II*, C-NSGA-II and SPEA2), the analysis of the
results obtained does not indicate that this had any deleterious effect on
these two algorithms. Figure 16 illustrates the distribution of the
cardinalities of the final fronts obtained by the $\epsilon$-MOEA and
cone$\epsilon$-MOEA for the 16 test problems considered.

Figure 16: Distribution of cardinalities of the final fronts returned by the
$\epsilon$-MOEA and cone$\epsilon$-MOEA for the test problems considered.

While the results above suggest no evidence that the mean final population
sizes of the algorithms are different, some questions remain open to question.
In particular the issue of whether the $\epsilon$\- or
cone$\epsilon$-dominance criteria are capable of producing competitive or
superior results when exactly the same number of points is considered is a
point that deserves some attention.

To address the issue above, a complementary experiment was performed in order
to isolate the possible effects of different cardinalities from the effects of
the dominance criteria themselves. The experiment takes advantage from the
fact that the test problems employed in this work have known, very finely-
sampled Pareto-optimal fronts, and was designed as follows.

A number of desired points was chosen, and each dominance criteria was allowed
to choose a set of exactly the prescribed size. In the case of the
$\epsilon$-dominance, the values of $\epsilon$ were calculated using (8), and
afterwards adjusted using small perturbations until the desired number of
solutions was reached. For the cone$\epsilon$-dominance, eleven values of
$\kappa$ were chosen, covering the entire range of this parameter. For each
value, the $\epsilon$ values were calculated using (15), and then adjusted to
match the desired cardinality. For the one-at-a-time crowding distance and the
clustering criteria from SPEA1/C-NSGA-II and SPEA2, no parameters other than
the desired set size had to be set.

After each of the dominance criteria returned the selected sample from the
Pareto set, the quality metrics $\Delta$ and $HV$ were calculated for the
sets444In this case comparisons of the $\gamma$ and $CS$ metrics would be
pointless, since all criteria are already working on samples of the Pareto-
optimal front., and compared using the criteria as the experimental factor
levels, and both Problem and Cardinality as blocking factors. The parameters
used in the experiment were:

  * •

Cardinality: twenty values of cardinality, from 181 to 200 points, were used;

  * •

Problems: six problems from the test set were used: DEB52, POL, DTLZ2, DTLZ7,
ZDT1 and ZDT3;

  * •

Criteria: the criteria tested were:

    * –

Crowding distance (CD), as implemented in the NSGA-II*;

    * –

Clustering (CL1), as implemented in the C-NSGA-II;

    * –

Clustering (CL2), as implemented in the SPEA2;

    * –

$\epsilon$-dominance ($\epsilon$);

    * –

cone$\epsilon$-dominance (cone);

  * •

$\kappa$ values: the kappa values were chosen so as to cover the entire
variation range of this parameter. Nine equally-spaced points were included
($0.1,\ 0.2,\ldots,\ 0.9$) together with two extreme points ($0.01$ and
$0.99$);

The analysis of this experiment was performed in two parts: first, the effects
of the $\kappa$ values on the performance of the cone$\epsilon$-criterion on
the $\Delta$ and $HV$ metrics was investigated. The best-, median-, and worst-
case instances of the cone$\epsilon$-dominance were then included in an
overall comparison against the performance of the remaining criteria. In both
cases, the randomized complete block design (Montgomery,, 2008) was used, with
both Problems and Desired Cardinality as blocking factors.

Figure 17 shows the results concerning the effects of $\kappa$ on the
performance of the cone$\epsilon$-dominance approach in this matched-
cardinality experiment. These results suggest that, for the test problems
employed, the performance of this approach regarding both metrics is
reasonably robust to variations in the value of $\kappa$ within a broad range,
with no significant differences for $\kappa\geq 0.3$. As the value of $\kappa$
is reduced below $0.3$, the performance concerning the $\Delta$ metric
significantly worsens, with a marked decrease in the effectiveness of the
criterion for $\kappa=0.01$. For the $HV$ metric no statistically significant
differences were detected.

(a) $\kappa$ effects on the $\Delta$ metric (b) $\kappa$ effects on the $HV$
metric

Figure 17: Effect of different $\kappa$ values on the performance of the
cone$\epsilon$-dominance, after removing problem and cardinality effects. The
bars represent $99\%$ confidence intervals, and the dotted line shows the
grand mean for each metric.

From the results of the $\kappa$ analysis, three cases were selected for the
overall analysis: worst, median, and best case performance. Since the $HV$
metric presented no significant differences, the point estimator of the
$\Delta$ metric was used as the ranking variable. The worst, median, and best
cases were, respectively, $\kappa=0.01$, $\kappa=0.5$ and $\kappa=0.8$. These
three variants of the cone$\epsilon$ approach were compared with the remaining
criteria, yielding the results shown in Figure 18.

(a) Comparison of mean $\Delta$ performance (b) Results for $HV$ metric

Figure 18: Performance comparison of the dominance criteria used in the MOEAs
employed in this work, after removing problem and cardinality effects. The
bars represent $99\%$ confidence intervals, and the dotted line shows the
grand mean for each metric.

These results suggest an interesting view of the ability of the dominance
criteria examined to maintain a diverse set of solutions. As expected, the
$HV$ metric did not present relevant differences555While the differences were
statistically significant, the practical effect is essentially zero, as the
magnitude of the differences for this metric is vanishingly small., as the
criteria were already working on samples taken from the Pareto-optimal set. As
for the $\Delta$ metric, it can be seen that the statistical ordering of the
results, at the $99\%$ confidence level, was:

$\overline{\mbox{CL2\ \ }\underline{\mbox{\ CL1\ \ \ CD\ }}}\underline{\mbox{\
Cone0.8\ \ \ Cone0.5\ }}\ \ \ \ \mbox{Cone0.01}\ \ \ \ \epsilon$

These results suggest that the only criterion significantly better than the
best- and median-case cone$\epsilon$-approach was the clustering mechanism
used in the SPEA2. The $\epsilon$-dominance had the worst diversity
performance in this experiment, followed by the worst-case cone$\epsilon$
which, with $\kappa=0.01$, was close to the condition where it effectively
becomes equal to the $\epsilon$ criterion. In general, the
cone$\epsilon$-dominance criterion seems to be a competitive alternative to
the clustering and crowding distance approaches as long as it employs its
usual values for the $\kappa$ parameter, e.g., the median case $\kappa=0.5$
used in the experiments reported throughout this paper. Moreover, the
$\epsilon$-based methods were not designed for the purpose of imposing an
ordering on a non-dominated set, while the clustering and crowding distance
methods allow for this ordering, making them less sensitive, in terms of
diversity metrics, to the removal of points. These characteristics may explain
the statistical ordering observed in the experiment.

###  7.4 Possible Extensions

From the results discussed previously, one can notice that while generally
effective, the cone$\epsilon$-method appears to have difficulties in dealing
with the redundant solutions (see Fig. 15(f)), which may negatively affect its
behavior. In problems where redundant solutions are known to exist, the
performance of the cone$\epsilon$-approach could be improved by using smaller
$\kappa$ values. More importantly, some authors argue that the probability of
occurrence of this “redundancy problem” becomes higher as the number of
objectives is increased (Deb et al., 2005b, ). Indeed, an extension of the
cone$\epsilon$-concept should be used to handle these classes of problems. A
simple and efficient strategy can be obtained by extending the area cone-
dominated by a solution, e.g., as in the $\alpha$-dominance (Ikeda et al.,,
2001). Since in this case the cone of dominance ensures the Pareto optimality,
only a single linear system should be solved, thus decreasing the
computational complexity of the method. Also, the adoption of a hypergrid in
the objective domain remains very useful in order to guarantee the convergence
and the diversity preservation during the evolution of a MOEA. However, a
deeper investigation on the characteristics and performance of the
cone$\epsilon$-dominance on many-objective optimization problems (MaOPs) is
left for a future work. Preliminary results on the ability of the $\epsilon$\-
and cone$\epsilon$-dominance criteria to adequately order sets of points in
multi and many-objective optimization problems can be found in (Batista et
al., 2011a, ).

##  8 Conclusions

In this paper we presented a relaxed form of dominance, named the
cone$\epsilon$-dominance, for improving the convergence and diversity
preservation capabilities of multiobjective evolutionary algorithms. This
alternative dominance criterion can be seen as a hybridization of the concepts
of $\epsilon$-dominance and of proper efficiency with respect to cones. After
presenting the mathematical definition of the cone$\epsilon$-dominance, we
incorporated this concept into a steady-state MOEA, and showed that the
computational complexity of this MOEA is $\mathcal{O}\left((mN)^{2}\right)$ in
the worst case.

A comprehensive evaluation of the cone$\epsilon$-dominance strategy was also
presented. We employed six MOEAs to evaluate the relative effectiveness of the
proposed methodology: NSGA-II, SPEA2, C-NSGA-II, $\epsilon$-MOEA, NSGA-II*,
and cone$\epsilon$-MOEA. The performance of these algorithms was quantified by
means of four different metrics, selected to evaluate the convergence and
diversity of the final fronts obtained by each algorithm: the convergence
metric ($\gamma$), the diversity metric ($\Delta$), the hypervolume or
s-metric ($HV$), and a proposed generalized form of the coverage of two sets
metric, named the coverage of many sets ($CS$). For the experimental design,
sixteen two- and three-objective test problems were selected, including the
classic DTLZ and ZDT benchmark families.

The results obtained for the four metrics were analyzed both in terms of
statistical significance and magnitude of the differences. From the analysis
of the effect sizes, calculated for each algorithm on each metric after
removing the problem effects, we could conclude that the cone$\epsilon$-MOEA
was capable of providing an efficient balance between convergence and
diversity. These results strongly support the conclusion that the proposed
algorithm is an interesting and competitive approach for the solutions of
multiobjective optimization problems.

## Acknowledgment

The authors would like to acknowledge the insightful comments provided by the
anonymous reviewers, who greatly contributed for the quality of this paper.
Support for this work was provided by the National Council for Research and
Development (CNPq), grants 306910/2006-3, 141819/2009-0, 305506/2010-2 and
472446/2010-0; and by the Research Foundation of the State of Minas Gerais
(FAPEMIG, Brazil), grants Pronex: TEC 01075/09 and Pronem: CEX APQ-04611-10.

## References

  * Aguirre and Tanaka, (2009) Aguirre, H. and Tanaka, K. (2009).  Adaptive $\epsilon$-ranking on MNK-landscapes.  In IEEE Symposium on Computational Intelligence in Multi-Criteria Decision-Making, MCDM 2009, pages 104–111.
  * (2) Batista, L., Campelo, F., Guimarães, F., and Ramírez, J. (2011a).  A comparison of dominance criteria in many-objective optimization problems.  In Proceedings of the IEEE Congress on Evolutionary Computation, CEC 2011, volume 1, pages 2359–2366, New Orleans. IEEE Press.
  * (3) Batista, L., Campelo, F., Guimarães, F., and Ramírez, J. (2011b).  Pareto cone $\epsilon$-dominance: Improving convergence and diversity in multiobjective evolutionary algorithms.  In Takahashi, R. H., Deb, K., Wanner, E. F., and Greco, S., editors, Evolutionary Multi-Criterion Optimization, EMO 2011, LNCS 6576, pages 76–90. Springer.
  * (4) Batista, L., Campelo, F., Guimarães, F., and Ramírez, J. (2011c).  Source codes, samples of the true Pareto fronts, and raw and processed results of the experiments.  Online: http://ppgee.ufmg.br/~fcampelo/files/EC-2011-01/.
  * Branke et al., (2001) Branke, J., Kau$\beta$ler, T., and Schmeck, H. (2001).  Guidance in evolutionary multi-objective optimization.  Advances in Engineering Software, 32(6):499–507.
  * Chankong and Haimes, (2008) Chankong, V. and Haimes, Y. Y. (2008).  Multiobjective Decision Making: Theory and Methodology.  Dover, New York.
  * Chen et al., (2002) Chen, Y.-L., Liao, W.-B., Yang, Y.-R., Shen, K.-Y., Wang, S.-C., and Chang, Y.-C. (2002).  The interactive surrogate worth trade-off method for multi-objective decision-making in reactive power sources planning.  In Proceedings of the International Conference on Power System Technology, volume 2, pages 863–866.
  * Coello Coello et al., (2002) Coello Coello, C. A., Van Veldhuizen, D. A., and Lamont, G. B. (2002).  Evolutionary Algorithms for Solving Multi-Objective Problems, volume 5 of Genetic Algorithms and Evolutionary Computation.  Kluwer Academic Publishers.
  * Davison and Hinkley, (1997) Davison, A. and Hinkley, D. (1997).  Bootstrap Methods and their Application.  Cambridge University Press, 1st edition.
  * Deb, (1999) Deb, K. (1999).  Multi-objective genetic algorithms: Problem difficulties and construction of test problems.  Evolutionary Computation, 7(3):205–230.
  * Deb, (2001) Deb, K. (2001).  Multiobjective Optimization Using Evolutionary Algorithms.  John Wiley & Sons, Chichester.
  * Deb and Agrawal, (1995) Deb, K. and Agrawal, R. (1995).  Simulated binary crossover for continuous search space.  Complex Systems, 9(2):115–148.
  * Deb and Goyal, (1996) Deb, K. and Goyal, M. (1996).  A combined genetic adaptive search (GeneAS) for engineering design.  Computer Science and Informatics, 26(4):30–45.
  * Deb et al., (2003) Deb, K., Mohan, M., and Mishra, S. (2003).  Towards a quick computation of well-spread Pareto-optimal solutions.  In Fonseca, C. M., Fleming, P. J., Zitzler, E., Thiele, L., and Deb, K., editors, Evolutionary Multi-Criterion Optimization, EMO 2003, LNCS 2632, pages 222–236. Springer.
  * (15) Deb, K., Mohan, M., and Mishra, S. (2005a).  Evaluating the $\epsilon$-dominance based multi-objective evolutionary algorithm for a quick computation of Pareto-optimal solutions.  Evolutionary Computation, 13(2):501–525.
  * Deb et al., (2002) Deb, K., Pratap, A., Agarwal, S., and Meyarivan, T. (2002).  A fast and elitist multiobjective genetic algorithm: NSGA-II.  IEEE Transactions on Evolutionary Computation, 6(2):182–197.
  * (17) Deb, K. and Srinivasan, A. (2006a).  Innovization: Innovative design principles through optimization.  In Proceedings of the Genetic and Evolutionary Computation Conference, GECCO 2006, pages 1629–1636.
  * (18) Deb, K. and Srinivasan, A. (2006b).  Monotonicity analysis, evolutionary multi-objective optimization, and discovery of design principles.  Technical Report 2006004, KanGAL.
  * (19) Deb, K., Thiele, L., Laumanns, M., and Zitzler, E. (2005b).  Scalable test problems for evolutionary multiobjective optimization.  In Abraham, A., Jain, L., and Goldberg, R., editors, Evolutionary Multiobjective Optimization, Theoretical Advances and Applications, pages 105–145. Springer.
  * Ellis, (2010) Ellis, P. D. (2010).  The Essential Guide to Effect Sizes: An Introduction to Statistical Power, Meta-Analysis and the Interpretation of Research Results.  Cambridge University Press.
  * Farhang-Mehr and Azarm, (2002) Farhang-Mehr, A. and Azarm, S. (2002).  Diversity assessment of Pareto optimal solution sets: An entropy approach.  In Proceedings of the 2002 Congress on Evolutionary Computation, CEC 2002, volume 1, pages 723–728.
  * Hernández-Díaz et al., (2007) Hernández-Díaz, A., Santana-Quintero, L., Coello Coello, C., and Molina, J. (2007).  Pareto-adaptive $\epsilon$-dominance.  Evolutionary Computation, 15(4):493–517.
  * Hernández-Díaz et al., (2011) Hernández-Díaz, A. G., Santana-Quintero, L. V., Coello Coello, C. A., Molina, J., and Caballero, R. (2011).  Improving the efficiency of $\epsilon$-dominance based grids.  Information Sciences, 181(15):3101–3129.
  * Hirsch et al., (2011) Hirsch, C., Shukla, P. K., and Schmeck, H. (2011).  Variable preference modeling using multi-objective evolutionary algorithms.  In Proceedings of the 6th International Conference on Evolutionary Multi-Criterion Optimization, EMO 2011, LNCS 6576, pages 91–105.
  * Ikeda et al., (2001) Ikeda, K., Kita, H., and Kobayashi, S. (2001).  Failure of Pareto-based MOEAs: Does non-dominated really mean near to optimal?  In Proceedings of the 2001 Congress on Evolutionary Computation, CEC 2001, volume 2, pages 957–962, Seoul, South Korea.
  * Janson and Merkle, (2005) Janson, S. and Merkle, D. (2005).  A new multi-objective particle swarm optimization algorithm using clustering applied to automated docking.  In Hybrid Metaheuristics, HM 2005, LNCS 3636, pages 128–141.