# Obtaining Smoothly Navigable Approximation Sets in Bi-Objective Multi-Modal
Optimization

Renzo J. Scholman  Centrum Wiskunde & InformaticaAmsterdamThe Netherlands
Delft University of Technology, Delft, The NetherlandsDelftThe Netherlands
<EMAIL_ADDRESS>,  Anton Bouter  Centrum Wiskunde &
InformaticaAmsterdamThe Netherlands<EMAIL_ADDRESS>,  Leah R.M. Dickhoff
Leiden University Medical CenterLeidenThe Netherlands<EMAIL_ADDRESS>,
Tanja Alderliesten  Leiden University Medical CenterLeidenThe Netherlands
<EMAIL_ADDRESS>and  Peter A.N. Bosman  Centrum Wiskunde &
InformaticaAmsterdamThe Netherlands<EMAIL_ADDRESS>
(2022)

###### Abstract.

Even if a Multi-modal Multi-Objective Evolutionary Algorithm (MMOEA) is
designed to find all locally optimal approximation sets of a Multi-modal
Multi-objective Optimization Problem (MMOP), there is a risk that the found
approximation sets are not smoothly navigable because the solutions belong to
various niches, reducing the insight for decision makers. Moreover, when the
multi-modality of MMOPs increases, this risk grows and the trackability of
finding all locally optimal approximation sets decreases. To tackle these
issues, two new MMOEAs are proposed: Multi-Modal Bézier Evolutionary Algorithm
(MM-BezEA) and Set Bézier Evolutionary Algorithm (Set-BezEA). Both MMOEAs
produce approximation sets that cover individual niches and exhibit inherent
decision-space smoothness as they are parameterized by Bézier curves. MM-BezEA
combines the concepts behind the recently introduced BezEA and MO-HillVallEA
to find all locally optimal approximation sets. Set-BezEA employs a novel
multi-objective fitness function formulation to find limited numbers of
diverse, locally optimal, approximation sets for MMOPs of high multi-modality.
Both algorithms, but especially MM-BezEA, are found to outperform the MMOEAs
MO_Ring_PSO_SCD and MO-HillVallEA on MMOPs of moderate multi-modality with
linear Pareto sets. Moreover, for MMOPs of high multi-modality, Set-BezEA is
found to indeed be able to produce high-quality approximation sets, each
pertaining to a single niche.

Evolutionary algorithms, multi-modal multi-objective optimization, niching,
Bézier curve estimation

††doi: 10.1145/nnnnnnn.nnnnnnn††isbn: 978-x-xxxx-xxxx-x/YY/MM††conference: The
Genetic and Evolutionary Computation Conference 2022; July 9–13, 2022; Boston,
USA††journalyear: 2022††copyright: acmcopyright††booktitle: GECCO ’22: Genetic
and Evolutionary Computation Conference, July 09–13, 2022, Boston, USA††price:
††ccs: Mathematics of computing Evolutionary algorithms

##  1\. Introduction

Many real-world optimization problems have multiple conflicting objectives,
whereby improvement in one objective often results in the deterioration of
another. Multi-Objective Evolutionary Algorithms (MOEAs), like the Multi-
Objective Real-Valued Gene-pool Optimal Mixing (MO-RV-GOMEA) (Bouter et al.,
2017b) algorithm, are widely accepted to be well-suited to solve such Multi-
objective Optimization Problems (MOPs) (Deb, 2001). The aim is to obtain a set
of solutions, called the approximation set, such that all solutions are non-
dominated and the set itself is close to the set of Pareto-optimal solutions.

A more complex type of MOPs is that of Multi-modal MOPs (MMOPs), where the
goal is not to find one, but multiple, if not all, (local) Pareto sets (PSs).
Each of these pertains to a niche, a subset of the search space, where a
single mode resides. The PSs may, however, well map to the same Pareto Front
(PF) in objective space, similar to having multiple (locally) optimal
solutions of the same quality in a single-objective problem, e.g., the sine
function.

In order to have MOEAs solve MMOPs, they need additional tools that prevent
their convergence to a single niche of the landscape (Mahfoud, 1996). Niching
(Li et al., 2017) is one of many diversity-preserving tools used by Multi-
modal MOEAs (MMOEAs) to effectively and simultaneously search for solutions
near the (local) PS in each niche. Niching has been successfully applied to
established MMOEAs in the form of the multi-objective particle swarm
optimization using ring topology and special crowding distance
(MO_Ring_PSO_SCD) algorithm (Yue et al., 2017) and Omni-optimizer (Deb and
Tiwari, 2008) among others.

Figure 1. Approximation set and front with parallel coordinates plot as
produced by MO-RV-GOMEA on MinDist problem, where
$f_{0}(\boldsymbol{x})=min(||\boldsymbol{x}-[1,-1]||,||\boldsymbol{x}-[-1,1]||)$
and
$f_{1}(\boldsymbol{x})=min(||\boldsymbol{x}-[1,1]||,||\boldsymbol{x}-[-1,-1]||)$.
Shaded blue and red regions correspond to niches with global PSs.

Visualization of problem definition.

Most MMOEAs do not explicitly model multiple approximation sets, but include
diversity preserving techniques to ensure that solutions from multiple niches
are maintained. The result of these MMOEAs is usually given in the form of a
single approximation front, often derived from (a subset of) the elitist
archive. A decision maker can then investigate this front by traversing the
solutions for desired trade-offs. However, the underlying solutions are taken
from several distinct niches. This results in observing a counterintuitive
change in decision variable values when navigating the approximation front.
Figure 1 shows such an approximation front and set that contains solutions
from both modes on the MinDist problem (Maree et al., 2019) and demonstrates
the counterintuitive changes in decision variable values in the parallel
coordinates plot (i.e., if one were to traverse the front and inspect
solutions). It shows that in objective space a front is found that looks to
have approximated the PF to (near) optimality, but the solutions jump around
throughout decision space. Especially in highly multi-modal problems,
solutions can originate from many different modes, which may make it more
insightful for decision makers to be able to request approximation sets for a
specific number of modes.

The issue of counterintuitive navigation along the approximation front has
also been explored in recent works, which introduced a new indicator-based
MOEA for bi-objective optimization called BezEA (Maree et al., 2020). A new
problem formulation for population-based MOEAs was introduced whereby they
parameterized approximation sets as Bézier curves. This ensured the
navigational smoothness of an approximation set, whilst still being able to
find good approximation sets when using the HyperVolume indicator (HV)
(Zitzler and Thiele, 1999). By design, BezEA disallowed curves to dominate
parts of themselves to ensure that the approximation set constitutes a single
niche in the landscape.

Recent work that included the concept of niching showed promising results in
maintaining multiple approximation sets in a population-based MMOEA called the
MO-HillVallEA (Maree et al., 2019). The authors extended the concept of Hill-
Valley Clustering (HVC) (Maree et al., 2018) for MOPs to Multi-Objective HVC
(MO-HVC) for MMOPs. It is capable of finding and preserving approximation
sets, one for each niche, in parallel over time by considering Pareto
domination per niche. However, MO-HillVallEA produces approximation sets that
are not inherently smooth due to slight oscillation around the PS.
Furthermore, no limit on the number of approximation sets can be set, which
degrades the quality of the approximation sets in highly multi-modal problems
as the population is divided over all niches.

In this work, the notions of niching through HVC and Bézier curve
parameterizations will be combined. The use of niching allows to effectively
search the multi-modal landscape. Whereas the use of Bézier curve
parameterizations not only enforces the smooth and intuitive navigability that
is desired by decision makers, but also enforces each approximation set to be
within a single niche. Furthermore, the use of the HV indicator allows the new
algorithms to converge closer to the PS as compared to the Pareto dominance-
based algorithms (Berghammer et al., 2012). Two new algorithms are proposed
that combine the aforementioned techniques, called Multi-Modal Bézier
Evolutionary Algorithm (MM-BezEA) and Set Bézier Evolutionary Algorithm (Set-
BezEA). These respectively try to find all approximation sets for a given MMOP
or a specific number of approximation sets for MMOPs that are highly multi-
modal, where each approximation set consists of solutions from a single mode.

In order to combine the techniques of Bézier curve parameterizations and HVC
into the two proposed new algorithms, several contributions are made. First,
the problem of how to niche approximation sets in the form of Bézier curve
parameterizations is resolved. Second, initialization of approximation sets
within a single niche is enabled, as otherwise, clustering becomes ambiguous
if they span multiple niches. Last, a new fitness function formulation is
introduced to find a preset number of approximation sets for highly multi-
modal problems.

##  2\. Bézier parameterizations

One of the key features of the newly proposed algorithms is that Bézier
parameterizations are used as approximation sets for bi-objective optimization
(Maree et al., 2020). This allows the algorithm to model the approximation set
as a smooth curve through decision space.

###  2.1. Definition of solution set

An $\ell$-dimensional Bézier curve $\mathbf{B}(t;C_{q})$ can be defined using
$q\geq 2$ control points $c_{i}$ in an ordered set
$C_{q}=\left\\{\boldsymbol{c}_{1},...,\boldsymbol{c}_{q}\right\\}$, where
$\ell$ is the problem dimensionality and
$\boldsymbol{c}_{i}\in\mathbb{R}^{\ell}$. The full notation is:

(1)
$\small\mathbf{B}(t;C_{q})=\sum_{i=1}^{q}\binom{q-1}{i-1}(1-t)^{q-i-1}t^{i-1}\boldsymbol{c}_{i}\text{
for }0\leq t\leq 1$

The endpoints of the Bézier curve are always defined by the first and last
control points, whilst the other control points are normally not located on
the curve. A solution set of given size $p$,
$S_{p,q}(C_{q})=\left\\{\boldsymbol{x}_{1},...,\boldsymbol{x}_{p}\right\\}$
with $\boldsymbol{x}_{i}\in\mathbb{R}^{\ell}$, can now be parameterized by a
Bézier curve by selecting an evenly spread set of $p$ points
$\boldsymbol{x}_{i}$. Figure 2 visualizes two solution sets $S_{p,q}(C_{q})$
parameterized by Bézier curves. The solution set $S_{p,q}(C_{q})$ is formally
defined as:

(2) $\small
S_{p,q}\left(C_{q}\right)=\left\\{\mathbf{B}\left(\dfrac{0}{p-1};C_{q}\right),\mathbf{B}\left(\dfrac{1}{p-1};C_{q}\right),...,\mathbf{B}\left(\dfrac{p-1}{p-1};C_{q}\right)\right\\}$

$S_{p,q}(C_{q})$ is parameterized for (M)MOEAs by taking the concatenation of
the decision variables in the set $C_{q}$ as a solution (Maree et al., 2020;
Beume et al., 2007). This results in a solution being of the form
$[\boldsymbol{c}_{1},...,\boldsymbol{c}_{q}]\in\mathbb{R}^{q*\ell}$.

Figure 2. Bézier curves with $q\in{2,3}$ control points in black. In red is
the interpolated curve with in blue the $p=11$ points evenly spread in the
domain of $t$ along the curve. (Maree et al., 2020)

Visualization of Bézier curves.

###  2.2. Evaluation

To evaluate a solution set $S_{p,q}(C_{q})$, a number of new functions were
previously introduced (Maree et al., 2020). These functions are briefly
explained in the following paragraphs. Figure 3 illustrates these functions to
give the reader a more graphical indication.

A new function $A^{nb}\left(S_{p,q}\right)$ has been introduced that
calculates a navigational Bézier (nb) order $o_{nb}$. This order is defined as
starting from the best solution for objective $f_{0}$ to the best solution in
$f_{1}$. All solutions that are dominated by others, which lie on the curve,
are omitted from the subset that defines the navigational order. An
approximation set $\mathcal{A}_{p,q,o_{nb}}$ is then the resulting subset of
$S_{p,q}(C_{q})$. The quality of the approximation set
$\mathcal{A}_{p,q,o_{nb}}$ can now be evaluated, e.g., using the uncrowded HV
indicator (Maree et al., 2021).

A new constraint function $C\left(\mathcal{S}_{p,q},o_{nb}\right)$ was also
introduced. It is employed in order to not only push all dominated solutions
on the curve towards the undominated region of the search space, but also to
prevent the curve from wrapping around itself in objective space. This may for
instance happen if a curve stretches across two local PSs, which is not
preferential. The constraint function uses the uncrowded distance metric
$ud_{f}(\boldsymbol{x}_{i},\mathcal{A})$ (Touré et al., 2019), which measures
the Euclidean distance from a dominated point $\boldsymbol{x}_{i}$ to the
approximation boundary $\delta f(\mathcal{A}_{p,q,o_{nb}})$ in objective
space. Furthermore, to further increase pressure towards the unfolding of
Bézier curves in objective space, all dominated solutions and those not in
$\mathcal{A}_{p,q,o_{nb}}$ are pulled towards their neighbouring solutions on
the Bézier curve by taking the Euclidean distance in objective space between
these solution and their neighbours as an additional constraint value. All
dominated solutions from $S_{p,q}(C_{q})$ now have their uncrowded distance
values and the Euclidean distances in objective space to neighbours of those
not in $\mathcal{A}_{p,q,o_{nb}}$ summed up as a constraint for the total
solution set. In combination with constraint domination (Deb, 2000), this
constraint pushes all solutions along the Bézier curve towards the undominated
region.

Figure 3. Evaluation of Bézier parameterizations. (Maree et al., 2020)

##  3\. Niching methods

To allow the two algorithms proposed in this paper, i.e., MM-BezEA and Set-
BezEA, to effectively search the multi-modal landscape, several previously
introduced niching methods are used and combined. These are employed in order
to extend on the uni-modal search that is originally performed by BezEA. As
the number of modes is usually unidentified beforehand, the algorithms need to
be able to adapt to the number of modes present in an MMOP, or the number a
user is after.

###  3.1. HVC and MO HVC

HVC is a so-called two-stage niching approach that clusters and evolves the
population for multi-modal single-objective optimization problems. In each
generation, the first stage is used to locate each of the distinct niches, for
each of which a core search algorithm is initialized in the second stage.

At the heart of the HVC approach is the Hill-Valley Test (HVT) (Ursem, 1999),
which can be utilized to determine whether two solutions reside in the same
niche. This test first determines an edge between two solutions
$\boldsymbol{x}_{i}$ and $\boldsymbol{x}_{j}$ in the search space. Along this
edge, $N_{t}$ evenly spread points are evaluated, determined by the distance
between the two solutions divided by the expected edge length. If any of these
$N_{t}$ test points have a fitness that is worse than that of
$\boldsymbol{x}_{i}$ and $\boldsymbol{x}_{j}$, the test detects that there is
a hill in between them. Consequently, the two solutions are to be put in
separate clusters. On the other hand, if all $N_{t}$ points have equal or
better fitness values than both $\boldsymbol{x}_{i}$ and $\boldsymbol{x}_{j}$,
these two solutions belong to the same valley and are to be clustered
together. A pseudocode example of this test can be found in Algorithm 1.

Input: Solutions $\boldsymbol{x}_{i}$, $\boldsymbol{x}_{j}$,
$N_{t}\in\mathbb{Z}$, objective function $f$

for  _$k=1,...,N_{t}$_ do

$\boldsymbol{x}_{test}=\boldsymbol{x}_{i}+\dfrac{k}{N_{t}+1}(\boldsymbol{x}_{j}-\boldsymbol{x}_{i})$

if  _$\text{max}\left\\{f(\boldsymbol{x}_{i}),f(\boldsymbol{x}_{j})\right\\}
<f(\boldsymbol{x}_{test})$_ then return false

return true

Algorithm 1 [Boolean] = HillValleyTest($\boldsymbol{x}_{i}$,
$\boldsymbol{x}_{j}$, $N_{t}$, $f$) (Ursem, 1999)

In order to determine in which order the solutions are to be clustered (i.e.,
undergo the HVT), the concept of the nearest better tree (Preuss, 2010) is
employed. It sorts the solutions in descending order of fitness, whereby the
first solution automatically forms the first cluster. Then, the population is
traversed in this descending order, starting at the second best solution, and
each solution $\boldsymbol{x}_{i}$ is tested with the HVT to decide if it
belongs to the same cluster as its nearest better neighbour
$\boldsymbol{x}_{j}$. If the test returns false, it is repeated at most
$\ell+1$ times with the next-nearest better neighbours $\boldsymbol{x}_{j}$,
before concluding that the solution $\boldsymbol{x}_{i}$ belongs to a new
cluster.

The MO-HillVallEA algorithm (Maree et al., 2019) expands on the previous HVC
approach in the form of MO-HVC. It uses the same concept of the HVT, but now
performs clustering for each of the $m$ objective functions separately, which
results in $m$ cluster sets of the population. To obtain a single cluster set,
the intersection of all clustering sets is taken, similar to the colored
regions of Figure 1.

###  3.2. Restart scheme with elitist archive

Various algorithms implemented a form of a restart scheme whereby the
population size is increased over time. Examples of such schemes are the
interleaved multistart scheme (Bouter et al., 2017b) and the restart-
Covariance Matrix Adaptation Evolution Strategy with Increased Population
(IPOP-CMA-ES) algorithm (Auger and Hansen, 2005). In HillVallEA (Maree et al.,
2018), an elitist archive is combined with a restart scheme, where the
population size is doubled after each restart as in IPOP-CMA-ES (Auger and
Hansen, 2005). By employing the HVT to check if a solution resides in another
niche, the elitist archive of HillVallEA is capable of holding the elites for
each of the modes, despite it being developed for single objective problems.

To prevent HillVallEA from revisiting already searched modes, it makes use of
the elitist archive, which is inspired by the repelling subpopulations (RS-
CMSA) algorithm (Ahrari et al., 2017) that defines taboo regions close to
elites. The steps taken to discard the regions of the search space, for which
an elite was already found in one of the earlier populations, start with
adding the elites to the population of the current restart. Then, all
solutions are clustered using HVC, followed by discarding all clusters that
have one of the elites as their best solution. As a result of discarding these
regions of the search space, more attention is given to undiscovered parts of
the search space after each restart.

###  3.3. Fitness sharing

Fitness sharing (Holland, 1975) is another niching technique aimed at keeping
the population diverse. It does so by degrading certain individuals fitness
values if these solutions are close to one another. This is to be done using a
niche radius parameter $\sigma_{share}$ defining which solutions are close,
together with a scaling factor $\alpha$ that degrades the fitness (Goldberg et
al., 1992). However, it has been shown non-trivial to set correct values for
these parameters (Darwen and Yao, 1995).

##  4\. Multi Modal-Bézier Evolutionary Algorithm

In this section the Multi-Modal Bézier Evolutionary Algorithm, abbreviated as
MM-BezEA, is described. It has been created through the combination and
modification of techniques described in the previous sections. The most
notable of the modifications are the adjustments implemented in HVC in order
to apply it to Bézier curve parameterizations, as well as the initialization
of approximation sets within niches.

###  4.1. Clustering approximation sets

The Bézier curves are evaluated using the uncrowded HV measure (Maree et al.,
2021). Since this is a scalar, the HVC approach is essentially applicable in a
rather straightforward manner because it intuitively allows the clustering of
single-objective problems. However, the approximation set
$\mathcal{A}_{p,q,o_{nb}}$ that is used in the HV calculation only considers
the undominated indices of the Bézier solution set $S_{p,q}(C_{q})$ as defined
in $o_{nb}$. Hence, the objective value given to the solution set
$S_{p,q}(C_{q})$ seems highly dependent on how many dominated solutions there
are on the Bézier curve.

Input: Solutions sets $\mathbf{S_{i}}$, $\mathbf{S_{j}}$, int $N_{t}$,
objective functions $f_{0},...,f_{m-1}$

Output : Whether $\mathbf{S_{i}}$ and $\mathbf{S_{j}}$ belong to the same
niche

for  _$l=1,...,q$_ do

$\boldsymbol{c}_{i,l}=$ control point $l$ of $\mathbf{S_{i}}$

$\boldsymbol{c}_{j,l}=$ control point $l$ of $\mathbf{S_{j}}$

// check $\boldsymbol{c}_{i,l}$ and $\boldsymbol{c}_{j,l}$ are in same niche
for all m objectives

for  _$k=0,...,m-1$_ do

if  _HillValleyTest( $\boldsymbol{c}_{i,l}$, $\boldsymbol{c}_{j,l}$, $N_{t}$,
$f_{k}$)_ then return false

return true

Algorithm 2 [$B$] = Bezier-HillValleyTest($\mathbf{S_{i}}$, $\mathbf{S_{j}}$,
$N_{t}$, $f$)

To enable the clustering of Bézier solution sets $S_{p,q}(C_{q})$, the idea
behind MO-HVC can be used on the set of control points $C_{q}$ because each of
these is a single solution as normally defined in MO optimization. Also, since
a solution set is defined to be deteriorating in $f_{0}$ and improving in
$f_{1}$ according to $o_{nb}$, the order of the control points is inverted if
$f_{0}(\boldsymbol{c}_{1})<f_{0}(\boldsymbol{c}_{q})$ does not hold (Maree et
al., 2020). Accordingly, the $i$-th Bézier solution can be designated to be in
the same niche as the $j$-th Bézier solution if their control points
$\boldsymbol{c}_{l}^{i}\in C_{q}^{i}$ and $\boldsymbol{c}_{l}^{j}\in
C_{q}^{j}$ for $l=\left\\{1,...,q\right\\}$ are in the same niche. In a
general sense, the same HVC approach as used in HillVallEA is used, but
inspiration has been taken from the MO-HVC approach to produce a new test for
Bézier solution sets, which is shown in Algorithm 2.

###  4.2. Initialization within niches

The original BezEA algorithm initializes all solution sets by sampling from a
uniform distribution over the search space. As it is an MOEA that was not
designed for multi-modal optimization, the uniform initialization allows
solution sets to be initialized within or in between any niche(s). Clustering
these solutions with the newly introduced Bézier HVT will result in finding a
large number of separate niches as each control point has to be in the same
mode.

To prevent this, a new initialization method for Bezier solution sets is
proposed to enforce their initialization within a niche. First, an iteration
of MO-HVC is run on a set of $q*N$ solutions, $N$ being the population size,
that is sampled from a uniform distribution over the search space, where the
resulting clusters include all of the $\boldsymbol{x}_{test}$ solutions
resulting from applying Algorithm 1. Second, selection is performed on the
clusters proportional to their size in order to reduce their combined size,
with test solutions, down to $q*N$. Lastly, Bézier solution sets
$S_{p,q}(C_{q})$ are initialized by randomly choosing $q$ solutions as the
control points from one single cluster $\mathcal{C}$ as produced by MO-HVC. A
visualization of the result can be seen in Figure 4.

Figure 4. Initialization of Bézier solutions ($q=2$) for MinDist.

###  4.3. Algorithm overview

MM-BezEA has a similar structure as the restart scheme in HillVallEA (Maree et
al., 2018) that is described in section 3.2. Every iteration, the combination
of taboo regions and subsequent initialization of Bézier curves takes place.
For each of the resulting niches, a core search algorithm is run for one
generation, which in the case of MM-BezEA is the RV-GOMEA algorithm (Bouter et
al., 2017a) that is also used in BezEA. At the end of each generation, the
Bézier HVT of Algorithm 2 is used in the HVC step. This step takes all
solutions originating from all clusters and clusters them again for the next
generation. In between generations, the notion of cluster registration
(Bosman, 2010) is used on the cluster mean closest in decision space to
transfer the parameters for RV-GOMEA between the clusters of each generation.
An overview of the algorithm in the form of pseudocode is given in Algorithm
3.

Input: MO function $f$, population size $N$, number of test points $p$, number
of control points $q$, evaluation or runtime budget

Output : Elitist archive $\mathbb{E}=[\mathcal{E}_{0},\mathcal{E}_{1},...]$

$\mathcal{P}_{mo}=$ UniformSampling($q\cdot N$, $f$)

$\mathbb{C}_{mo}=$ MO-HillValleyClustering($\mathcal{P}_{mo}$, $f$)

$\mathbb{C}=$ InitializeBezierSolutions($\mathbb{C}_{mo}$, $q$, $p$, $f$)

$\mathbb{E}=$ InitializeElitistArchive($\mathbb{C}$)

while  _budget remaining_ do

$\mathcal{P}=\mathbb{E}$

for  _$\mathcal{C}_{i}\in\mathbb{C}$_ do

$\mathcal{O}_{i}=$ core_search_algorithm$\left(\mathcal{C}_{i}\right)$

$\mathcal{P}=\mathcal{P}\cup\mathcal{O}_{i}$

$\mathbb{C}_{prev}=\mathbb{C}$

$\mathbb{C}=$ BezierHillValleyClustering$\left(\mathcal{P},f\right)$

$\mathbb{E}=$ ConstructElitistArchive$\left(\mathbb{C},\mathbb{E}\right)$

$\mathbb{C}=$ RemoveElitesFrom$\left(\mathbb{C}\right)$

$\mathbb{C}=$ ClusterRegistration$\left(\mathbb{C},\mathbb{C}_{prev}\right)$

Algorithm 3 [$\mathbb{E}$] = MM-BezEA$\left(f,N,p,q\right)$

##  5\. Set Bézier Evolutionary Algorithm

As MM-BezEA is a niching method that divides its attention over all of the
niches found by HVC, it is intractable to use MM-BezEA in highly multi-modal
landscapes. In order to overcome this limitation, a second algorithm is
introduced that uses MO-RV-GOMEA (Bouter et al., 2017b) with an adapted
fitness function definition and encoding. Rather than encoding each solution
to be one Bézier curve, and clustering the population to find different
niches, this algorithm, named Set-BezEA, encodes and searches for a set of $b$
Bézier curves. It does so without further niching the population, which
prevents this approach from falling apart in highly multi-modal landscapes.

In the original BezEA, solutions represent Bézier curves consisting of the $q$
control points $C_{q}$. Given a problem dimensionality $\ell$, a Bézier
solution thus consists of $q\cdot\ell$ parameters. In Set-BezEA, a solution
represents $b$ Bézier curves consisting of $b\cdot q$ control points $C_{q}$
concatenated into a set of control points denoted as $C_{q,b}$. The control
points for the $i$-th Bézier solution are denoted as $C_{q}^{i}$. Hence,
solutions in Set-BezEA take the form
$[C_{q}^{1},...,C_{q}^{b}]\in\mathbb{R}^{b\cdot q\cdot\ell}$. The notation
$S_{p,q}^{i}$, equal to the definition of $S_{p,q}$ in Equation 2, is used to
denote the $i$-th Bézier solution of the set $S_{p,q,b}$ with $b$ Bézier
solutions.

###  5.1. Fitness function formulation

In order to evaluate the newly formulated solutions, a new fitness function is
proposed that can be seen in mathematical form in Equation 3. The original
fitness function formulation from BezEA is used in order to determine the
constraint $L(S_{p,q,b}(\phi))$ and $\text{ HV}_{f,S_{p,q,b}(\phi)}$ values
for all of the $b$ Bézier solutions combined. Fitness sharing (Holland, 1975)
is incorporated with a niching radius based on the expected edge length (Maree
et al., 2018), which is checked between each of the $b$ Bézier solutions
$S_{p,q}^{i}$.

However, fitness sharing only describes a minimal separation criterion.
Typically, in practice, decision makers would like to be able to find Bézier
solutions that are diverse, yet all of high quality, while limited in number.
For this reason, the problem is defined as an MO formulation. Herein, the
second objective is a distance objective $\Delta_{S_{p,q,b}}$ that is to be
maximized. It calculates the sum of cosine distances between all pairs of the
$b$ Bézier solutions to move them apart in decision space and promote finding
diverse approximation sets.

(3) maximize $f_{0}$: $\displaystyle\text{HV}_{f,S_{p,q,b}(\phi)}$
$\displaystyle=\sum_{i=0}^{b-1}\dfrac{\text{HV}_{f}\left(A^{nb}\left(S_{p,q}^{i}(\phi)\right)\right)}{n_{i}}$
maximize $f_{1}$: $\displaystyle\Delta_{S_{p,q,b}(\phi)}$
$\displaystyle=\sum_{i=0}^{b-1}\sum_{j=i+1}^{b-1}\dfrac{S_{p,q}^{i}(\phi)\cdot
S_{p,q}^{j}(\phi)}{||S_{p,q}^{i}(\phi)||\hskip 2.0pt||S_{p,q}^{j}(\phi)||}$
with constraint: $\displaystyle L\left(S_{p,q,b}(\phi)\right)$
$\displaystyle=\sum_{i=0}^{b-1}L_{b}\left(S_{p,q}^{i}(\phi),o_{nb}^{i}(\phi)\right)$

$\small\begin{split}\text{where
}&f:\mathbb{R}^{\ell}\rightarrow\mathbb{R}^{m},\phi\in\mathbb{R}^{b*q*\ell},\\\
&A^{nb}:\text{solution set }S_{p,q}^{i}(\phi)\rightarrow\text{approximation
set }\mathcal{A}_{p,q,o_{nb}}^{i}\\\ &n_{i}\text{ number of Bézier curves in
same niche as }S_{p,q}^{i}\\\ &o_{nb}^{i}(\phi)\text{ the navigational order
of the }i\text{-th Bézier solution}\\\
&L_{b}(S_{p,q}^{i}(\phi),o_{nb}^{i}(\phi))\text{ the Bézier constraint
function}\\\ \end{split}$

###  5.2. Adaptive steering

The MO formulation of the previous section creates a trade-off between finding
far apart solutions at the edges of the search space and those closer to the
PSs. However, especially assuming the global optimum cannot be obtained, the
primary goal is to obtain distinct solutions that are of high quality. Most
(M)MOEAs try to find improvements spread over the entire PF, whereas in some
situation only a part of the PF is interesting. In order to steer Set-BezEA
towards the interesting parts of the PF, the notion of adaptive steering
(Alderliesten et al., 2015) can be used on the first objective value. Starting
at allowing Bézier set solutions to be within $50\%$ of the best achieved
cumulative HV it scales linearly in 50 generations to $99\%$. This will force
all solutions towards higher HV values when combined with constraint
domination (Deb, 2000).

###  5.3. Linkage structure and partial evaluations

To optimize for the new fitness formulation, Set-BezEA uses MO-RV-GOMEA, which
uses a linkage model (Thierens, 2010) to determine which variables are to be
changed at the same time. It is often taken to be a linkage tree that is
constructed using the unweighted pair group method with arithmetic mean
(UPGMA) (Gronau and Moran, 2007) to iteratively cluster decision variables. A
new linkage structure is designed for Set-BezEA to guide MO-RV-GOMEA and allow
it to exploit the underlying structure of a problem. For this new linkage
structure, the ideas behind the linkage structures defined for BezEA and UHV-
GOMEA were combined (Maree et al., 2020).

When the $j$-th decision variable of a control point in Bézier curve changes,
the $j$-th decision variable of all test solutions of the approximation set
$S_{p,q}^{i}(C_{q})$ changes, which makes it beneficial to change the same
decision variable for all $q$ control points of a Bézier curve at once. This
is the first layer of the new linkage structure, where, for each of the $b$
Bézier solutions, their $j$-th decision variables of all $q$ control points
are grouped. Subsequently, the $\ell$ clusters of $q$ decision variables for
one Bézier solution are hierarchically grouped using UPGMA. The second layer
combines the $b$ resulting groups of the decision variables of all $b$ Bézier
solutions by concatenating them into one large linkage tree.

The newly defined linkage structure enables MO-RV-GOMEA to perform partial
evaluations in all situations. Even if the underlying problem setting is
black-box in nature, partial evaluations are enabled by the second layer, as
each of these Bézier curves can be evaluated separately for its HV
contribution to the first objective. The second distance objective always has
to be recalculated, but it does not require any evaluations of the original
fitness function.

###  5.4. MO-RV-GOMEA adaptations

The approximation sets $S_{p,q}^{i}(C_{q})$, defined by the Bézier curves, are
evaluated with the HV indicator that makes use of a reference point. But, the
choice of reference point highly influences whether or not the extremes of the
PF can be found (Auger et al., 2009).

The control points $\boldsymbol{c}_{1}^{i}$ and $\boldsymbol{c}_{q}^{i}$ are
likely to come as close as possible to the constraint region, defined by being
dominated by the reference point, to maximize the HV. However, if one of these
two control points falls just slightly in the constraint region, the
constraint function $C(S_{p,q,b}(\phi))$ combined with constraint domination
will make it one of the worst performing solutions. Furthermore, MO-RV-GOMEA
makes use of the Gene-pool Optimal Mixing (GOM) variation operator, thereby
disallowing deterioration of fitness. These two effects combined, quickly
reduce the variance in the population, since very slight adjustments can
significantly deteriorate the individual. In order to overcome this downside,
the GOM part of MO-RV-GOMEA has been adapted to allow deterioration of the
objective values, i.e., accepting a modification that would otherwise have
been rejected, with a probability of $5\%$, as was previously done in RV-GOMEA
(Bouter et al., 2017a). This adaptation enables MO-RV-GOMEA to effectively
search the outer regions of the search space close to the constraint space by
preserving more variance.

##  6\. Experiments

MM-BezEA and Set-BezEA are empirically benchmarked on several test problems.
The results are compared to MO-HillVallEA (Maree et al., 2019), MO-RV-GOMEA
(Bouter et al., 2017b), and MO_Ring_PSO_SCD (Yue et al., 2017).
MO_Ring_PSO_SCD is implemented through the PlatEMO framework (Tian et al.,
2017), where all of the used metrics and problems are manually implemented as
well. For the other algorithms, original C++ implementations are used.

###  6.1. Test problems

Several test problems are employed. First of these is the MinDist problem
(Maree et al., 2019) that was described in the introduction, where linear PSs
are to be found. The other employed functions are frequently used in
literature, namely OmniTest (Deb and Tiwari, 2008), Two on One (Preuss et al.,
2006), and Sympart {1,2,3} (Rudolph et al., 2006). Lastly, several problems
are taken from the Multi-modal Multi-objective test Function (MMF) benchmark
suite (Liang et al., 2018) in the form of MMF {1,2,12,14,15}.

A mix of PS and PF shapes have been chosen to determine the capabilities of
the two new algorithms on different problem types. Table 1 shows some of the
important characteristics for each of the problems. For all problems with a
configurable number of PSs, it is set to 2, likewise the problem dimension
$\ell$ is fixed to 2. In order to determine the values of the performance
indicators, the reference PSs will be made using 5000 points that adhere to
the analytical formulas describing the PSs. In the case of Two on One, a very
close approximation is used, as no analytical formula is available (Preuss et
al., 2006).

Table 1. Bi-objective problem instances and characteristics. Problem | $\ell$ | #PS | PS Shape | PF Shape
---|---|---|---|---
MinDist | $[2...\infty)\in\mathbb{Z}$ | n | Linear | Convex
Omni Test | $[2...\infty)\in\mathbb{Z}$ | $3^{\ell}$ | Linear | Convex
Two on One | 2 | 2 | Linear | Convex
Sympart 1 & 2 | 2 | 9 | Linear | Convex
Sympart 3 | 2 | 9 | Non-linear | Convex
MMF 1 & 2 | 2 | 2 | Non-linear | Convex
MMF 12 | $[2...\infty)\in\mathbb{Z}$ | n | Linear | Disconnected
MMF 14 & 15 | $[2...\infty)\in\mathbb{Z}$ | n | Linear | Concave

###  6.2. Benchmark setup

In order to get a fair comparison, each of the algorithms will be given an
equally sized budget of 200k function evaluations for each of the problems.
This removes the influence of the used programming languages, as the
computation time is not limited. The parameters of MO-RV-GOMEA,
MO_Ring_PSO_SCD, and MO-HillVallEA are set to the values reported in relevant
literature. Furthermore, for each problem and metric, the average over 31 runs
will be taken.

The elitist archives sizes $N_{\mathbb{E}}$ are set to be 1250 for MO-RV-GOMEA
and MO-HillVallEA. The population size $N$ is set to 96 for MO-RV-GOMEA and
250 for MO-HillVallEA. MO-RV-GOMEA uses a linkage tree as its linkage model,
with a total of 5 clusters. For MO_Ring_PSO_SCD the population size is 800
(Yue et al., 2017). For the MM-BezEA and Set-BezEA algorithms, the number of
control points $q$ for each approximation set is set to 2. Just like for the
original BezEA algorithm, these two algorithms are both given population sizes
of 76 (Maree et al., 2020). The number of test points $p$ is set to 7. Lastly,
Set-BezEA is always configured to produce $b=2$ Bézier solutions as
approximation sets.

###  6.3. Performance indicators

The HV indicator (Zitzler and Thiele, 1999) is used to see how well the
algorithms perform in getting close to the PF. As a result of the use of test
points in the MM-BezEA and Set-BezEA algorithms, these algorithms have a
limited amount of points in the approximation set that can be used to
calculate the HV values. Therefore, a subset of the approximation set will be
taken for the other algorithms to allow a fair comparison based on the HV
indicator. Specifically, the same number of test points is selected for a fair
comparison by means of greedy Hypervolume Subset Selection (gHSS) (Guerreiro
et al., 2016).

We further use a relatively new performance indicator for multi-modal multi
objective optimization, named Pareto Set Proximity (PSP) (Yue et al., 2017).
It is an indicator that determines how well all PSs are approximated by taking
the Cover Rate (CR), that shows how well the extremes of all PSs are captured,
divided by the Inverted Generational Distance in decision space (IGDX) (Zhang
and Li, 2007), which can be used to determine how close the approximation sets
are to the PSs. For the IGDX measure, the approximation sets as produced by
MM-BezEA and Set-BezEA are interpolated by taking 1000 intermediate points
before determining the IGDX value. This can be performed relatively easily as
interpolating these parameterizations does not require any extra fitness
evaluations.

Finally, we use a performance indicator regarding smoothness, for which we
follow the definition as introduced in the work on BezEA (Maree et al., 2020).
It captures how smooth an approximation set can be navigated in terms of
decision variables by measuring the detour length in decision space when
traversing the approximation set from one solution to the next via an
intermediate solution, as compared to going to the next solution directly. The
smoothness approaches its maximum value of 1 if all solutions would be
colinear in decision space, where the lowest possible value is 0. In cases
where the approximation set is clustered, like in MO-HillVallEA, MM-BezEA and
Set-BezEA, the average smoothness over all clusters will be taken. In the
other cases the smoothness over the entire approximation set is taken.

###  6.4. Results

The results for all problems and algorithms per indicator can be seen in Table
2.

The HV results clearly show that all algorithms are capable of performing
nearly equally in obtaining a good approximation front. However, The MM-BezEA
and Set-BezEA algorithms do deteriorate in performance on the MMF1 and 2
problems that have non-linear PSs. The deterioration is inherently caused by
the chosen parameterizations that create approximation sets which are linear
in shape. Another problem instance where a smaller HV for the new algorithms
is obtained, is that of MMF12. Here, despite MM-BezEA obtaining the best PSP
values, the approximation sets did not fully approximate the actual PSs and
did not cover the endpoints.

The PSP indicator shows similar results, except that MO-RV-GOMEA performs
worse as it is not an MMOEA and therefore does not explicitly search for
multiple niches. Again promising results for MM-BezEA and Set-BezEA are shown
in cases where linear PSs can be found. In the problems with non-linear PSs,
MO-HillVallEA and MO_Ring_PSO_SCD manage to find better approximations.

The smoothness results show, as intended, that the chosen parameterizations
inherently cause smooth approximation sets. Therefore, MM-BezEA and Set-BezEA
always obtain a perfect smoothness of 1.0. Other algorithms do not obtain
this, except for MO-RV-GOMEA on 2 of the 11 problems. A visualization of the
results of MM-BezEA on the MinDist problem is given in Figure 5. This figure
depicts the smooth course of the decision variables in the parallel
coordinates plot for the rightmost approximation set in decision space. It
contrasts sharply to the parallel coordinates plot of Figure 1.

Figure 5. Approximation sets and front with parallel coordinates plot for one
of the approximation sets as produced by MM-BezEA on MinDist problem (Maree et
al., 2019).

Visualization of MM-BezEA results.

Table 2. Results (avg. ($\pm$ st.dev.)) per problem and algorithm over 31 runs, bold identifies best result with statistical significance (Wilcoxon rank-sum test with $\alpha=0.05$ and Holm-Bonferroni correction). Indicator | Problem | MM-BezEA | Set-BezEA | MO-HillVallEA | MO_Ring_PSO_SCD | MO-RV-GOMEA
---|---|---|---|---|---|---
HV | MinDist | 1.17e+02 (±9.43e-05) | 1.17e+02 (±3.02e-04) | 1.17e+02 (±1.62e-02) | 1.17e+02 (±5.95e-03) | 1.17e+02 (±1.36e-03)
OmniTest | 8.48e+00 (±2.18e-06) | 8.43e+00 (±2.34e-01) | 8.47e+00 (±1.96e-03) | 8.47e+00 (±7.34e-04) | 8.47e+00 (±4.82e-04)
Sympart 1 | 1.17e+02 (±2.16e-05) | 1.16e+02 (±1.67e+00) | 1.17e+02 (±1.42e-02) | 1.17e+02 (±7.64e-03) | 1.17e+02 (±1.30e-03)
Sympart 2 | 1.17e+02 (±2.91e-05) | 1.15e+02 (±1.97e+00) | 1.17e+02 (±7.77e-03) | 1.17e+02 (±8.75e-03) | 1.17e+02 (±4.47e-03)
Sympart 3 | 1.17e+02 (±9.61e-05) | 1.16e+02 (±1.76e+00) | 1.17e+02 (±1.61e-02) | 1.17e+02 (±9.21e-03) | 1.17e+02 (±4.91e-03)
TwoOnOne | 1.13e+02 (±1.33e-04) | 1.13e+02 (±9.27e-04) | 1.13e+02 (±2.39e-04) | 1.13e+02 (±1.82e-04) | 1.13e+02 (±1.10e-04)
MMF 1 | 6.04e-01 (±3.86e-02) | 6.29e-01 (±5.63e-02) | 8.05e-01 (±2.37e-04) | 8.05e-01 (±8.69e-05) | 8.05e-01 (±6.60e-05)
MMF 2 | 6.34e-01 (±2.01e-04) | 5.88e-01 (±9.42e-02) | 8.04e-01 (±6.90e-04) | 8.04e-01 (±9.59e-04) | 8.05e-01 (±1.75e-04)
MMF 12 | 1.78e+00 (±2.02e-06) | 1.38e+00 (±2.02e-01) | 2.06e+00 (±2.57e-03) | 2.06e+00 (±2.05e-03) | 2.06e+00 (±1.49e-04)
MMF 14 | 5.63e+00 (±1.33e-05) | 5.63e+00 (±6.40e-03) | 5.63e+00 (±7.26e-04) | 5.63e+00 (±1.92e-03) | 5.63e+00 (±2.23e-04)
MMF 15 | 5.56e+00 (±2.03e-02) | 5.56e+00 (±4.47e-03) | 5.57e+00 (±6.52e-04) | 5.56e+00 (±1.54e-03) | 5.57e+00 (±1.79e-04)
PSP | MinDist | 3.26e+02 (±7.31e+01) | 3.77e+02 (±9.41e+01) | 5.02e+01 (±2.74e+00) | 6.97e+01 (±6.73e+00) | 1.21e+00 (±1.65e+00)
OmniTest | 2.36e+02 (±8.81e+01) | 1.71e-02 (±3.34e-02) | 7.13e+01 (±2.76e+00) | 6.90e+01 (±9.56e+00) | 1.36e-01 (±2.43e-01)
Sympart 1 | 2.67e+02 (±1.21e+02) | 1.22e-03 (±2.89e-03) | 3.56e+01 (±1.57e+00) | 2.79e+01 (±3.69e+00) | 1.16e-02 (±2.70e-02)
Sympart 2 | 3.09e+02 (±8.46e+01) | 2.76e-04 (±3.54e-04) | 3.60e+01 (±7.42e-01) | 2.38e+01 (±2.46e+00) | 1.02e-02 (±1.71e-02)
Sympart 3 | 6.41e+01 (±7.77e+01) | 1.27e-03 (±2.17e-03) | 4.33e+01 (±2.10e+00) | 2.64e+01 (±5.62e+00) | 8.09e-03 (±1.33e-02)
TwoOnOne | 3.04e+02 (±2.18e+02) | 2.69e+02 (±1.15e+02) | 4.50e+01 (±7.21e-01) | 2.45e+01 (±1.03e+01) | 2.68e+00 (±1.11e+00)
MMF 1 | 7.22e+00 (±2.41e+00) | 4.56e-01 (±3.55e-01) | 3.17e+01 (±6.84e-01) | 3.80e+01 (±6.79e+00) | 1.02e+00 (±2.89e-01)
MMF 2 | 4.00e+00 (±2.49e+00) | 4.29e-02 (±5.23e-02) | 1.17e+02 (±1.21e+01) | 5.04e+01 (±1.51e+01) | 2.18e+00 (±1.41e+00)
MMF 12 | 2.42e+01 (±7.72e+00) | 4.83e-01 (±6.81e-01) | 1.94e+01 (±6.46e+00) | 1.53e+01 (±1.46e-01) | 8.67e+00 (±1.36e+00)
MMF 14 | 2.68e+03 (±4.82e+02) | 1.72e+03 (±1.60e+03) | 3.70e+02 (±1.03e+01) | 2.31e+02 (±2.16e+01) | 1.08e+00 (±2.84e+00)
MMF 15 | 2.73e+02 (±4.17e+01) | 2.96e+02 (±1.41e+02) | 2.65e+02 (±4.58e+00) | 2.44e+02 (±7.25e+00) | 2.24e+01 (±1.40e-02)
Smoothness | MinDist | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 8.09e-01 (±3.70e-02) | 7.63e-01 (±5.91e-03) | 8.94e-01 (±1.81e-01)
OmniTest | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 9.23e-01 (±5.65e-03) | 7.28e-01 (±9.76e-03) | 7.16e-01 (±1.96e-01)
Sympart 1 | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 8.76e-01 (±2.52e-02) | 6.84e-01 (±9.76e-03) | 7.01e-01 (±2.09e-01)
Sympart 2 | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 8.78e-01 (±2.14e-02) | 5.73e-01 (±9.30e-03) | 7.87e-01 (±1.68e-01)
Sympart 3 | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 8.70e-01 (±3.10e-02) | 5.29e-01 (±1.18e-02) | 8.30e-01 (±1.92e-01)
TwoOnOne | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 7.77e-01 (±1.21e-02) | 7.47e-01 (±9.36e-03) | 7.51e-01 (±1.59e-01)
MMF 1 | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 9.01e-01 (±3.98e-02) | 7.82e-01 (±8.49e-03) | 5.86e-01 (±9.13e-02)
MMF 2 | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 9.38e-01 (±1.26e-02) | 5.01e-01 (±8.92e-03) | 8.46e-01 (±1.60e-01)
MMF 12 | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 8.35e-01 (±2.73e-02) | 6.35e-01 (±9.86e-03) | 1.00e+00 (±0.00e+00)
MMF 14 | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 9.31e-01 (±1.17e-02) | 8.71e-01 (±1.11e-02) | 9.62e-01 (±1.05e-01)
MMF 15 | 1.00e+00 (±0.00e+00) | 1.00e+00 (±0.00e+00) | 9.19e-01 (±1.43e-02) | 8.63e-01 (±8.54e-03) | 1.00e+00 (±0.00e+00)

###  6.5. Scalability

Set-BezEA has been created to overcome scalability issues with MM-BezEA in
highly multi-modal problems. To evaluate this, the OmniTest problem is used,
that is scalable in its problem dimensions whilst simultaneously scaling the
number of PSs. MM-BezEA and MO-RV-GOMEA will be compared to Set-BezEA with
$b\in{2,3,4}$, to show the impact of requesting more approximation sets.
Problem dimensions $\ell\in\left\\{2,3,5,10,20,30,50,100\right\\}$ have been
considered to see how the algorithms behave. The evaluations budget is set to
$100000*\ell$, with the rest of the parameters equal to before.

The PSP metric is omitted, as OmniTest scales very fast in the number of PSs
if the problem dimensionality is increased. In case $\ell\geq 8$, there are
over 5000 PSs, when only 5000 points of the PSs are analytically determined
for the indicators. Thus, only the smoothness and HV indicators are used.
Furthermore, as MM-BezEA and Set-BezEA produce approximation sets that are
each within a single niche, the approximation set of MO-RV-GOMEA will be post-
processed with MO-HVC to compare the average HV per niche.

Figure 6 shows the results for the HV and smoothness indicators. The results
show that on average MO-RV-GOMEA seems to obtain reasonable approximation sets
within each niche. However, many of these niches have only one or a few
solutions, which is also depicted by the lower smoothness values of the entire
approximation set on the right. Furthermore, the worst performance in terms of
HV per niche can be nearly two orders of magnitude lower than Set-BezEA when
$\ell=100$. MM-BezEA is, like MO-RV-GOMEA, capable of finding a good
approximation set for at least one niche as seen by the fact that the best
observed performance of all algorithms coincide, but MM-BezEA also has the
worst observed performance. For this problem, it is caused by the fact that
there are so many niches ($3^{\ell}$), that the algorithm cannot keep track of
all of them.

In the case of Set-BezEA, the results indicate that it is capable of finding
decent and distinct approximation sets. When more approximation sets are
optimized for, the average HV can decline slightly. However, even HV per niche
remains high in comparison to the other algorithms, whilst the smoothness
stays at its maximum of 1, showing that Set-BezEA exhibits the intended
behavior.

Figure 6. Results of dimensionality comparison. The left diagram shows the HV
results, with the approximation set smoothness on the right. The averages are
depicted as a line and the min. and max. values as shaded areas over 31 runs.
MM-BezEA and Set-BezEA smoothness results overlap.

##  7\. Discussion

MM-BezEA and Set-BezEA did not cover the endpoints of the PSs in the case of
MMF12, which resulted in lower HV values. This can be caused by the fact that
the Bézier fitness function will constrain a solution when one of its control
points is dominated in objective space by one of the test points. As the
endpoints of each part of the discontinuous PF are close to being dominated,
i.e., close to the constraint space, it can lead to not entirely capturing the
discontinuous pieces of the PF and thus resulting in a lower HV.

The use of the HV indicator to score individuals in MM-BezEA and Set-BezEA is
useful in the sense that it is a Pareto compliant indicator (Fleischer, 2003;
Zitzler et al., 2003), but it does suffer from a downside. In some situations
the endpoints of the approximation sets cannot reach the endpoints of the
Pareto set (Auger et al., 2009). Even in cases where the number of test points
will be set to infinity, the reference point can never be set so that the
extremes are captured.

Even though the smoothness indicator tries to determine whether an
approximation set is smooth by measuring the detour length, it comes down to
determining the angle between solutions. This implies that it only considers
linear curves to be perfectly smooth. If an approximation set oscillates
slightly around a PS, the angles between solutions will therefore result in a
lower smoothness. In cases where the niches can be separated in a good manner,
another definition of smoothness that considers the oscillation around the PS
might be more useful.

Future work could investigate the further use of the Bézier parameterizations
with more control points to approximate non-linear Pareto sets. Furthermore,
Bézier simplexes (Kobayashi et al., 2019) might be usable for problems with
more than two objectives.

##  8\. Conclusion

We proposed two algorithms, MM-BezEA and Set-BezEA, to search for
parameterizations of approximation sets that define smooth curves in the
decision space for bi-objective optimization problems. In highly multi-modal
problems, good navigability of approximation sets is desired by decision
makers, which can be further improved if approximation sets for a specific
number of modes can be requested. The results show that the smoothness is
indeed enforced by the Bézier parameterizations. Furthermore, the encoding of
$b$ curves in Set-BezEA produces an exact number of approximation sets that
are of high quality for highly multi-modal problems. Both algorithms were
outperformed in problems with non-linear Pareto sets, but only low-order
Bézier curves were researched. If the non-linearity of a mode surpasses the
Bézier order, it can no longer be captured by the Bézier curve. All together,
the algorithms are still potentially very useful in practice, which will be
targeted in future work.

###### Acknowledgements.

Leah R.M. Dickhoff was supported by the Dutch Cancer Society (Sponsor KWF
kankerbestrijding https://www.kwf.nl/, Project N.Grant #12183) and Sponsor
Elekta https://www.elekta.com/.

## References

  * (1)
  * Ahrari et al. (2017) Ali Ahrari, Kalyanmoy Deb, and Mike Preuss. 2017.  Multimodal Optimization by Covariance Matrix Self-Adaptation Evolution Strategy with Repelling Subpopulations.  _Evolutionary Computation_ 25, 3 (Sep 2017), 439–471.  https://doi.org/10.1162/evco_a_00182
  * Alderliesten et al. (2015) Tanja Alderliesten, Peter A. N. Bosman, and Arjan Bel. 2015\.  Getting the most out of additional guidance information in deformable image registration by leveraging multi-objective optimization. In _Medical Imaging 2015: Image Processing_ , Vol. 9413. SPIE, Bellingham, WA, USA, 469 – 475.  https://doi.org/10.1117/12.2081438
  * Auger et al. (2009) Anne Auger, Johannes Bader, Dimo Brockhoff, and Eckart Zitzler. 2009.  Theory of the Hypervolume Indicator: Optimal $\mu$-Distributions and the Choice of the Reference Point. In _Proceedings of the tenth ACM SIGEVO workshop on Foundations of genertic algorithms_ (Orlando, Florida, USA) _(FOGA ’09)_. Association for Computing Machinery, New York, NY, USA, 87–102.  https://doi.org/10.1145/1527125.1527138
  * Auger and Hansen (2005) Anne Auger and Nikolaus Hansen. 2005.  A restart CMA evolution strategy with increasing population size. In _2005 IEEE Congress on Evolutionary Computation_ , Vol. 2. IEEE, New York, NY, USA, 1769–1776 Vol. 2.  https://doi.org/10.1109/CEC.2005.1554902
  * Berghammer et al. (2012) Rudolf Berghammer, Tobias Friedrich, and Frank Neumann. 2012\.  Convergence of set-based multi-objective optimization, indicators and deteriorative cycles.  _Theoretical Computer Science_ 456 (Oct 2012), 2–17.  https://doi.org/10.1016/J.TCS.2012.05.036
  * Beume et al. (2007) Nicola Beume, Boris Naujoks, and Michael Emmerich. 2007\.  SMS-EMOA: Multiobjective selection based on dominated hypervolume.  _European Journal of Operational Research_ 181, 3 (Sep 2007), 1653–1669.  https://doi.org/10.1016/j.ejor.2006.08.008
  * Bosman (2010) Peter A.N. Bosman. 2010\.  The Anticipated Mean Shift and Cluster Registration in Mixture-Based EDAs for Multi-Objective Optimization. In _Proceedings of the 12th Annual Conference on Genetic and Evolutionary Computation_ (Portland, Oregon, USA) _(GECCO ’10)_. Association for Computing Machinery, New York, NY, USA, 351–358.  https://doi.org/10.1145/1830483.1830549
  * Bouter et al. (2017a) Anton Bouter, Tanja Alderliesten, Cees Witteveen, and Peter A. N. Bosman. 2017a.  Exploiting Linkage Information in Real-Valued Optimization with the Real-Valued Gene-Pool Optimal Mixing Evolutionary Algorithm. In _Proceedings of the Genetic and Evolutionary Computation Conference_ (Berlin, Germany) _(GECCO ’17)_. Association for Computing Machinery, New York, NY, USA, 705–712.  https://doi.org/10.1145/3071178.3071272
  * Bouter et al. (2017b) Anton Bouter, Ngoc Hoang Luong, Cees Witteveen, Tanja Alderliesten, and Peter A. N. Bosman. 2017b.  The Multi-Objective Real-Valued Gene-Pool Optimal Mixing Evolutionary Algorithm. In _Proceedings of the Genetic and Evolutionary Computation Conference_ (Berlin, Germany) _(GECCO ’17)_. Association for Computing Machinery, New York, NY, USA, 537–544.  https://doi.org/10.1145/3071178.3071274
  * Darwen and Yao (1995) Paul Darwen and Xin Yao. 1995.  A dilemma for fitness sharing with a scaling function. In _Proceedings of 1995 IEEE International Conference on Evolutionary Computation_ , Vol. 1. IEEE, New York, NY, USA, 166–171.  https://doi.org/10.1109/ICEC.1995.489138
  * Deb (2000) Kalyanmoy Deb. 2000\.  An efficient constraint handling method for genetic algorithms.  _Computer Methods in Applied Mechanics and Engineering_ 186, 2 (Jun 2000), 311–338.  https://doi.org/10.1016/S0045-7825(99)00389-8
  * Deb (2001) Kalyanmoy Deb. 2001\.  _Multi-Objective Optimization Using Evolutionary Algorithms_.  John Wiley & Sons, Inc., USA.
  * Deb and Tiwari (2008) Kalyanmoy Deb and Santosh Tiwari. 2008.  Omni-optimizer: A generic evolutionary algorithm for single and multi-objective optimization.  _European Journal of Operational Research_ 185, 3 (Mar 2008), 1062–1087.  https://doi.org/10.1016/j.ejor.2006.06.042
  * Fleischer (2003) Mark Fleischer. 2003\.  The Measure of Pareto Optima Applications to Multi-Objective Metaheuristics. In _Proceedings of the 2nd International Conference on Evolutionary Multi-Criterion Optimization_ (Faro, Portugal) _(EMO’03)_. Springer-Verlag, Berlin, Heidelberg, 519–533.
  * Goldberg et al. (1992) David E. Goldberg, Kalyanmoy Deb, and Jeffrey Horn. 1992\.  Massive Multimodality, Deception, and Genetic Algorithms. In _Parallel Problem Solving from Nature 2_ (Brussels, Belgium) _(PPSN-2)_. Elsevier Science Publishers B. V., Amsterdam, The Netherlands, The Netherlands, 37–48.
  * Gronau and Moran (2007) Ilan Gronau and Shlomo Moran. 2007.  Optimal implementations of UPGMA and other common clustering algorithms.  _Inform. Process. Lett._ 104, 6 (Dec 2007), 205–210.  https://doi.org/10.1016/j.ipl.2007.07.002
  * Guerreiro et al. (2016) Andreia P. Guerreiro, Carlos M. Fonseca, and Luís Paquete. 2016\.  Greedy Hypervolume Subset Selection in Low Dimensions.  _Evolutionary Computation_ 24 (Sep 2016), 521–544.  Issue 3. https://doi.org/10.1162/EVCO_a_00188
  * Holland (1975) John H. Holland. 1975\.  _Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence_.  University of Michigan Press, Ann Arbor, MI, USA.
  * Kobayashi et al. (2019) Ken Kobayashi, Naoki Hamada, Akiyoshi Sannai, Akinori Tanaka, Kenichi Bannai, and Masashi Sugiyama. 2019\.  Bézier simplex fitting: Describing pareto fronts of simplicial problems with small samples in multi-objective optimization. In _Proceedings of the 33rd AAAI Conference on Artificial Intelligence, AAAI 2019, the 31st Innovative Applications of Artificial Intelligence Conference, IAAI 2019 and the 9th AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2019_ (Honolulu, Hawai, USA). AAAI press, Palo Alto, CA, USA, 2304–2313.
  * Li et al. (2017) Xiaodong Li, Michael G. Epitropakis, Kalyanmoy Deb, and Andries Engelbrecht. 2017.  Seeking Multiple Solutions: An Updated Survey on Niching Methods and Their Applications.  _IEEE Transactions on Evolutionary Computation_ 21, 4 (Aug 2017), 518–538.  https://doi.org/10.1109/TEVC.2016.2638437
  * Liang et al. (2018) Jing Liang, Boyang Qu, Dunwei Gong, and Cai Yue. 2018\.  _Problem Definitions and Evaluation Criteria for the CEC 2019 Special Session on Multimodal Multiobjective Optimization_.  Technical Report.  https://doi.org/10.13140/RG.2.2.33423.64164
  * Mahfoud (1996) Samir W. Mahfoud. 1996\.  _Niching Methods for Genetic Algorithms_.  Ph. D. Dissertation. University of Illinois at Urbana-Champaign, USA.  UMI Order No. GAX95-43663.
  * Maree et al. (2019) Stefanus C. Maree, Tanja Alderliesten, and Peter A. N. Bosman. 2019\.  Real-Valued Evolutionary Multi-Modal Multi-Objective Optimization by Hill-Valley Clustering. In _Proceedings of the Genetic and Evolutionary Computation Conference_ (Prague, Czech Republic) _(GECCO ’19)_. Association for Computing Machinery, New York, NY, USA, 568–576.  https://doi.org/10.1145/3321707.3321759
  * Maree et al. (2020) Stefanus C. Maree, Tanja Alderliesten, and Peter A. N. Bosman. 2020\.  Ensuring Smoothly Navigable Approximation Sets by Bézier Curve Parameterizations in Evolutionary Bi-objective Optimization. In _Parallel Problem Solving from Nature – PPSN XVI_. Springer International Publishing, Cham, 215–228.
  * Maree et al. (2021) Stefanus C. Maree, Tanja Alderliesten, and Peter A. N. Bosman. 2021\.  Uncrowded Hypervolume-based Multi-objective Optimization with Gene-pool Optimal Mixing.  _Evolutionary Computation_ (Dec 2021), 1–24.  https://doi.org/10.1162/evco_a_00303
  * Maree et al. (2018) Stefanus C. Maree, Tanja Alderliesten, Dirk Thierens, and Peter A. N. Bosman. 2018.  Real-Valued Evolutionary Multi-Modal Optimization Driven by Hill-Valley Clustering. In _Proceedings of the Genetic and Evolutionary Computation Conference_ (Kyoto, Japan) _(GECCO ’18)_. Association for Computing Machinery, New York, NY, USA, 857–864.  https://doi.org/10.1145/3205455.3205477
  * Preuss (2010) Mike Preuss. 2010\.  Niching the CMA-ES via Nearest-Better Clustering. In _Proceedings of the 12th Annual Conference Companion on Genetic and Evolutionary Computation_ (Portland, Oregon, USA) _(GECCO ’10)_. Association for Computing Machinery, New York, NY, USA, 1711–1718.  https://doi.org/10.1145/1830761.1830793
  * Preuss et al. (2006) Mike Preuss, Boris Naujoks, and Günter Rudolph. 2006\.  Pareto Set and EMOA Behavior for Simple Multimodal Multiobjective Functions. In _Parallel Problem Solving from Nature – PPSN IX_ (Reykjavik, Iceland) _(PPSN’06)_. Springer-Verlag, Berlin, Heidelberg, 513–522.  https://doi.org/10.1007/11844297_52
  * Rudolph et al. (2006) Günter Rudolph, Boris Naujoks, and Mike Preuss. 2006\.  Capabilities of EMOA to Detect and Preserve Equivalent Pareto Subsets. In _Evolutionary Multi-Criterion Optimization, 4th International Conference, EMO 2007, Matsushima, Japan, March 5-8, 2007, Proceedings_ _(Lecture Notes in Computer Science, Vol. 4403)_. Springer, Berlin, Heidelberg, 36–50.  https://doi.org/10.1007/978-3-540-70928-2_7
  * Thierens (2010) Dirk Thierens. 2010\.  The Linkage Tree Genetic Algorithm. In _Proceedings of Parallel Problem Solving from Nature, PPSN XI_ (Berlin, Heidelberg) _(Lecture Notes in Computer Science)_. Springer, Berlin, Heidelberg, 264–273.
  * Tian et al. (2017) Ye Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin. 2017\.  PlatEMO: A MATLAB platform for evolutionary multi-objective optimization.  _IEEE Computational Intelligence Magazine_ 12, 4 (Nov 2017), 73–87.
  * Touré et al. (2019) Cheikh Touré, Nikolaus Hansen, Anne Auger, and Dimo Brockhoff. 2019.  Uncrowded Hypervolume Improvement: COMO-CMA-ES and the Sofomore Framework. In _Proceedings of the Genetic and Evolutionary Computation Conference_ (Prague, Czech Republic) _(GECCO ’19)_. Association for Computing Machinery, New York, NY, USA, 638–646.  https://doi.org/10.1145/3321707.3321852
  * Ursem (1999) R.K. Ursem. 1999\.  Multinational evolutionary algorithms. In _Proceedings of the 1999 Congress on Evolutionary Computation-CEC99_ , Vol. 3. IEEE, New York, NY, USA, 1633–1640 Vol. 3.  https://doi.org/10.1109/CEC.1999.785470
  * Yue et al. (2017) Cai Yue, Boyang Qu, and Jing Liang. 2017.  A Multi-objective Particle Swarm Optimizer Using Ring Topology for Solving Multimodal Multi-objective Problems.  _IEEE Transactions on Evolutionary Computation_ 22 (Sep 2017), 805–817.  https://doi.org/10.1109/TEVC.2017.2754271
  * Zhang and Li (2007) Qingfu Zhang and Hui Li. 2007.  MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition.  _IEEE Transactions on Evolutionary Computation_ 11 (12 2007), 712–731.  Issue 6. https://doi.org/10.1109/TEVC.2007.892759
  * Zitzler and Thiele (1999) Eckart Zitzler and Lothar Thiele. 1999.  Multiobjective evolutionary algorithms: a comparative case study and the strength Pareto approach.  _IEEE Transactions on Evolutionary Computation_ 3, 4 (Nov 1999), 257–271.  https://doi.org/10.1109/4235.797969
  * Zitzler et al. (2003) Eckart Zitzler, Lothar Thiele, Marco Laumanns, Carlos M. Fonseca, and Viviane G. da Fonseca. 2003\.  Performance assessment of multiobjective optimizers: an analysis and review.  _IEEE Transactions on Evolutionary Computation_ 7, 2 (May 2003), 117–132.  https://doi.org/10.1109/TEVC.2003.810758