RSG, a Method for Pareto Front Approximation and Reference Set Generation

1. Introduction

Multi-objective optimization has become an integral part of the decision-making for many real-world problems. In a multi-objective optimization problem (MOP), one is faced with the issue of concurrently optimizing k individual objectives. The set of optimal solutions is called the Pareto set. The image of the Pareto set is called the Pareto front. The latter set is, in many cases, most important for the decision maker (DM) since it provides him/her with an overview of the optimal performances for his/her project. What makes MOPs hard to deal with is that one can expect that both Pareto set and front form – at least locally and under certain assumptions on the model – objects of dimension $k-1$ ([1]). For the numerical treatment of MOPs, specialized evolutionary algorithms, called multi-objective evolutionary algorithms (MOEAs), have caught the interest of many researchers and practitioners during the last three decades ([2]). MOEAs are population-based and hence allow for the obtaining of a finite approximation of the entire solution set in one run of the algorithm. For the performance assessment of the outcome sets, several different indicators have been proposed so far (e.g., [3,4,5,6,7,8]). Some of these performance indicators are distance-based and require a “suitable” finite-size representation of the Pareto front. While until now, a vast variety of different MOEAs has been proposed and analyzed, it is fair to say that the generation of suitable reference sets has played a rather minor role in the evolutionary multi-objective optimization (EMO) community. It is evident that such reference sets should be complete. Further, as we show in this work, a biased representation can lead to misleading or even false information.

To fill this gap, we propose in this work the Reference Set Generator (RSG). The main steps of RSG are as follows: (i) a first approximation $A_y$ of the Pareto front is either taken or generated. This set can be, in principle, of arbitrary size, and the elements can be non-uniformly distributed along the Pareto front (i.e., biased). However, all of these elements have to be "close enough" to the set of interest. In order to obtain a bias-free approximation, (ii) component detection and (iii) a filling step is applied to $A_y$. Finally, (iv) a reduction step is applied. The RSG is applicable to MOPs with Pareto fronts of, in principle, any shape and dimension. Further, the resulting reference set is of any desired magnitude. We will show the strength of the novel method on several benchmark problems.

The remainder of this paper is organized as follows: in Section 2, we briefly recall the background for the understanding of this work. We further discuss the related work and the performance indicators that benefit from our approach. In Section 3, we first motivate the need for bias-free complete finite-size Pareto front approximations, and propose the Reference Set Generator (RSG). In Section 4, we present some numerical results on selected benchmark problems and compare the RSG to related algorithms. Finally, we draw our conclusions in Section 5 and give possible paths for future research.

2. Background and Related Work

We consider multi-objective optimization problems (MOPs) that can be mathematically expressed via

$$\underset{x\in Q}{min}F\left(x\right).$$

(MOP)

Hereby, the map F is defined as

$$F:Q\to {\mathbb{R}}^k,F\left(x\right)={(f_1\left(x\right),\dots ,f_k\left(x\right))}^T,$$

(1)

where we assume each of the individual objectives $f_i:Q\to \mathbb{R}$, $i=1\dots ,k$, to be continuous. We stress, however, that the method we propose in the sequel, RSG, can, in principle, also be applied to discrete problems. Q is the domain of the objective functions that is typically expressed by equality and inequality constraints.

In order to define the optimality of an MOP, one can use the concept of dominance [9].

Definition 1.

(a) Let $v,w\in {\mathbb{R}}^k$. Then the vector v is less than w ($v{<}_{p}w$), if $v_i<w_i$ for all $i\in \{1,\dots ,k\}$. The relation ${\le }_p$ is defined analogously.

(b) 

$y\in Q$ is dominated by a point $x\in Q$ ($x\prec y$) with respect to (MOP) if $F\left(x\right){\le }_{p}F\left(y\right)$ and $F\left(x\right)\ne F\left(y\right)$.

(c) 

$x\in Q$ is called a Pareto point or Pareto optimal if there is no $y\in Q$ that dominates x.

(d) 

The set $P_Q$ of Pareto optimal solutions

$$P_Q:=\{x\in Q:\nexists y\in Qs.t.y\prec x\}$$

(2)

is called the Pareto set.

(e) 

The image $F\left(P_Q\right)$ of the Pareto set is called the Pareto front.

One can expect that both the Pareto set and the Pareto front form under certain conditions, at least locally objects of dimension $k-1$. For details, we refer to [1]. Due to this "curse of dimensionality" It is, hence, not possible for an evolutionary multi-objective optimization algorithm (MOEA) to keep all promising candidate solutions (e.g., all non-dominated ones) during the algorithm run. It is, hence inevitable – at least for continuous problems – to select which of the promising solutions should be kept in order to obtain a "suitable" approximation of the solution set (in most cases, the Pareto front of the given MOP). Within MOEAs, this process is termed "selection". Another term, which can be used synonymously, is "archiving". The latter is typically used when the MOEA is equipped with an external archive.

Most of the existing MOEAs can be divided into three main classes: (a) MOEAs that are based on the concept of dominance (e.g., [10,11,12,13]), (b) MOEAs that are based on decompositions (e.g., [14,15,16,17,18]), and (c) MOEAs that make use of an indicator function (e.g., [19,20,21,22,23]). The selection strategies of the first generation of MOEAs of class (a) are based on a combination of non-dominated sorting and niching (e.g., [24,25,26]). Later, elite preservation was included, leading to an increased overall performance (and, as a consequence, better Pareto front approximations). This holds, e.g., for SPEA ([27]), PAES ([13]), SPEA-II ([11]), and NSGA-II ([10]). Theoretical studies on selection mechanisms have been done by the groups of Rudolph ([28,29,30,31]) and Hanne ([32,33,34,35]). All of these selection mechanisms deal with the abilities of the populations to reach the Pareto set/front. On the other hand, the distributions of the individuals along the Pareto sets/fronts have not been considered. The selection mechanisms of the MOEAs within classes (b) and (c) follow directly from the construction of the algorithms: selection in a MOEA of class (b) is done by considering the values of the chosen scalarization functions. Analogously, the selection in a MOEA of class (c) is done by considering the indicator contributions.

Existing (external) archiving strategies can also be divided into three classes: (a) unbounded archivers, (b) implicitly bounded archivers, and (c) bounded archivers. Unbounded archivers store all promising solutions during the algorithm run. The magnitudes of such archives can exceed any given threshold if the algorithm is run long enough. Unbounded archivers have, e.g., been used and analyzed in [36,37,38,39,40,41,42]. $ϵ$-dominance ([43]) can be viewed as a weaker concept of dominance. This relation allows single solutions to “cove” entire parts of the Pareto front of a given MOP, which is the basis for most implicitly bounded archivers. Such strategies have first been considered in the context of evolutionary multi-objective optimization (EMO) by Laumanns et al. ([44]). Later works proposed and analyzed different approximations (such as gap-free approximations of the Pareto front) and deal with different sets of interest (e.g., the consideration of all nearly optimal or all $ϵ$-locally optimal solutions), ([40,41,45,46,47,48,49,50]). Finally, bounded archivers have, e.g., been proposed in [51,52] where adaptive grid selections have been utilized. Bounded archivers tailored to the use of the Hypervolume indicator have been suggested in [53,54]). Laumanns and Zenklusen have proposed two bounded archivers that aim for $ϵ$-approximations of the Pareto front ([55]).

All of the selection/archiving strategies mentioned above have in common that they aim for a "best approximation" of the set of interest out of the given (finite) data. A related but slightly different problem is to generate a "suitable" (in particular complete and bias-free) finite size approximation of the set of interest S for the sake of comparisons, even if S is known approximately or even analytically but not “trivia” (e.g., a line segment). Tian et al. [56] describe the procedure of PlatEMO [57], which uses uniform sampling on the simplex and then maps these points to the particular Pareto front of each problem. However, an analytical expression or characterization of the Pareto front is required, and in some cases, the obtained set is not completely uniform. A method related to RSG can be found in [58] which has the aim of guiding the iterates of a particular Newton method toward the Pareto front. In this work, we extend this idea for the purpose of generating complete and bias-free Pareto front approximations of relatively large magnitudes (in particular compared to population sizes used in EMO).

Finally, reference sets as the ones generated by RSG are helpful for the evaluation of the performance qualities of candidate sets (populations) in EMO. More precisely, such sets are required for all distance-based indicators. The earliest such indicators are the Generational Distance (GD, [3]) and the Inverted Generational Distance. (IGD, [5]). Later, the indicator ${\Delta }_p$ ([6,59]) has been proposed that is a combination of slight variants of GD and IGD, and that can be viewed as an averaged version of the Hausdorff distance $d_H$. So far, there exist several extensions of these performance indicators. For instance, the consideration of continuous sets – either only the Pareto front or also the candidate solution set – has been done in [8,60] leading to modifications of $IGD$ and ${\Delta }_p$. The indicators IGD+ ([7]) and DOA ([61]) are modifications of IGD that are Pareto compliant.

3. Reference Set Generator (RSG)

In this section, we first motivate the need for complete and bias-free finite size Pareto front approximations, and then propose the Reference Set Generator (RSG) that targets for such sets.

3.1. Motivation

Distance-based indicators require a “suitable” finite-size approximation of the Pareto front in order to give a "correct" value for the approximation quality of the considered candidate solution set. This holds in particular for the above-mentioned indicators GD, IGD, and ${\Delta }_p$, together with their variants. Such representations are ideally uniformly spread along the entire Pareto front ([60]). This, however, is a non-trivial task unless the Pareto front is given analytically and has a relatively simple form (e.g., linear or perfectly spherical). Regrettably, this is the case for only a few test problems (e.g., DTLZ1 and DTLZ2). On the other hand, there exist quite a few benchmark MOPs where the shape of the Pareto set is relatively simple. For such problems, it is tempting to choose uniform samples from the Pareto set (i.e., $X=\{s_1,\dots ,s_m\}$, where all $s_i\in P_Q$), and to use the respective image $Y=F\left(X\right)$ to represent of the Pareto front. The following discussion shows, however, that this approach has to be handled with care since it can induce unwanted biases in the approximations that, in turn, may result in misleading indicator values.

As the first example, consider the one-dimensional bi-objective problem

$$F\left(x\right)=\left(\begin{array}{c}1-{\displaystyle \frac{1}{x}}\\ {\displaystyle \frac{1}{x}}\end{array}\right).$$

(3)

Let the domain be given by $Q=[0.1,3]$, then the Pareto set is identical to Q, and the Pareto front is the line segment that connects the points $a={(-\mathrm{9,10})}^{\top }$ and $b={(2/\mathrm{3,1}/3)}^{\top }$. Figure 1 shows the result when using N equally spaced points along the Pareto set. As it can be seen for $N=50$ and $N=500$, there is a clear bias of the images toward the right lower end of the Pareto front. For $N=\mathrm{10,000}$, the Pareto front approximation is "complete" (at least from the practical point of view) and appears to be perfect. However, it possesses the same bias. To see the impact of the reference set on the performance indicators consider the two hypothetical outcomes (e.g., possible results from different MOP solvers):

$$\begin{array}{cc}\hfill A& :=\left\{a+{\displaystyle \frac{2i-1}{10}}(b-a)|i\in \{1,\dots ,5\}\right\},\mathrm{and}\hfill \\ \hfill B& :=\left\{a+{\displaystyle \frac{i}{10}}(b-a)|i\in \{6,\dots ,10\}\right\}.\hfill \end{array}$$

(4)

Figure 2 shows the two sets together with the Pareto front. Note that A is the perfect 5-element approximation of the Pareto front: the elements are equally distributed along the Pareto front, and the extreme points are shifted “half-way in” ([60]). The set B is certainly not perfect as it, e.g., misses to "cover" more than half of the front. Table 1 shows values $I(O,R)$ for different distance-based indicators, the outcomes $O\in \{A,B\}$, and different representations R of the Pareto front. For $R=R_{\mathrm{10,000}}^x$ (the one shown in Fig. 1 (c)), all indicators – except $d_H={\Delta }_{\infty }$ – yield lower values for B than for A, indicating (erroneously) that B is better than A. This is not the case for $d_H$ since the Hausdorff distance is determined by the maximum of the considered distances and not by an average of those (however, $d_H$ has other disadvantages in the context of EMO, most prominently that it punishes single outliers [6]). The situation changes when selecting $R=R_{\mathrm{10,000}}^y$ as a representation of the front. This representation also contains $N=\mathrm{10,000}$ elements, but these are chosen uniformly along the Pareto front. Now, there is a tie for the two GD variants, and for all other indicators, A leads to better values than B. These values are indeed very close to the `correct’ values: all exact GD values are equal to zero since A and B are contained in the Pareto front. In order to compute the exact IGD values, a particular integral has to be solved [60]. When using $R_{10000}^y$ as representation, the computation of the IGD values can be interpreted as a Riemann sum with $N=\mathrm{10,000}$ equally sized sub-intervals leading to perfect solutions, at least from the practical point of view.

We repeat the process, but now using only $N=100$ elements for the representation (see Table 1). We can see the same trend, i.e., that B appears to be better for $R_{100}^x$ while A appears to be better when using $R_{100}^y$. Further, we see that the indicator values for $R_{100}^y$ are already quite close to the exact values (i.e., when using $R_{10000}^y$). While the proper choice of N may not be an issue for bi-objective problems ($k=2$), this may get important for a larger number of objectives due to the “curse of dimensionality”: at least for continuous problems, one can expect that the Pareto front forms under certain (mild) assumptions a manifold of dimension $k-1$ ([1]).

We next consider further test problems. Figure 4, Figure 5 and Figure 6 show the Pareto front approximations for the commonly used test problems DTLZ1, ZDT1, and ZDT3, respectively, where the representation has been obtained via uniform sampling in decision variable space along the Pareto set. In all cases, the Pareto front representations appear to be perfect for large enough values of N. However, certain biases can be observed for lower values of N. For all test problems, we selected two points out of the representation and showed their nearest neighbors. These values differ by one order of magnitude which confirms that the solutions are not uniformly distributed along the fronts. Using such representations, the same issues can arise as discussed above.

To conclude, the suitable representation of the Pareto front of a given MOP is crucial when considering distance-based performance indicators that use an average of the distances considered. Such representations are ideally equally distributed over the front. If the representation contains a bias, this may result in misleading indicator values, leading, in turn, to a wrong evaluation of the obtained results. In particular, the approach to performing the sampling along the Pareto set is, though tempting, not appropriate for such indicator-based indicators. In the sequel, we will propose a method that aims to achieve a uniform Pareto front representations.

Figure 3. Distribution resulting from uniformly sampling the Pareto set of DTLZ1 with N points. In (a) $d(y_1,N{N}_{y_1})=0.068181240846835$, $d(y_2,N{N}_{y_2})=0.008729713347982$.

Figure 3. Distribution resulting from uniformly sampling the Pareto set of DTLZ1 with N points. In (a) $d(y_1,N{N}_{y_1})=0.068181240846835$, $d(y_2,N{N}_{y_2})=0.008729713347982$.

Figure 4. Distribution resulting from uniformly sampling the Pareto set of DTLZ2 with N points. In (a) $d(y_1,N{N}_{y_1})=0.174311485495316$, $d(y_2,N{N}_{y_2})=0.030268871802677$.

Figure 4. Distribution resulting from uniformly sampling the Pareto set of DTLZ2 with N points. In (a) $d(y_1,N{N}_{y_1})=0.174311485495316$, $d(y_2,N{N}_{y_2})=0.030268871802677$.

Figure 5. Distribution resulting from uniformly sampling the Pareto set of ZDT1 with N points. In (a) $d(y_1,N{N}_{y_1})=0.101010101010101$, $d(y_2,N{N}_{y_2})=0.011440746020984$.

Figure 5. Distribution resulting from uniformly sampling the Pareto set of ZDT1 with N points. In (a) $d(y_1,N{N}_{y_1})=0.101010101010101$, $d(y_2,N{N}_{y_2})=0.011440746020984$.

Figure 6. Distribution resulting from uniformly sampling the Pareto set of ZDT3 with N points. In (a) $d(y_1,N{N}_{y_1})=0.066834581606972$, $d(y_2,N{N}_{y_2})=0.004486990395919$.

Figure 6. Distribution resulting from uniformly sampling the Pareto set of ZDT3 with N points. In (a) $d(y_1,N{N}_{y_1})=0.066834581606972$, $d(y_2,N{N}_{y_2})=0.004486990395919$.

3.2. RSG

In the following, we assume we are interested in a Pareto front approximation of size N for a given MOP. Further, we assume that we are given a set $A_y$ of (in principle) arbitrary size ℓ of non-dominated, possibly non-uniformly distributed points that are "close enough" to the PF. Below, we will discuss different strategies to obtain $A_y$. Given this data, the Reference Set Generation (RSG) process consists of three main steps: component detection, filling, and reduction. The idea is to fill the gaps between the points within each connected point of $A_y$. This leads to a more complete set F with a higher cardinality than $A_y$, which can then be reduced to obtain a uniform reference set of size N. In general, PFs can be disconnected, and if we simply fill the gaps in $A_y$, we may introduce points that do not belong to the PF (Figure 7(b)). Therefore, component detection must be performed before applying the filling process to each detected component (Figure 7(a)).

The general procedure of RSG is presented in Algorithm 1, and each of the main steps will be explained in the following subsections. Figures 8, 12, 9, and 13 illustrate all the steps for the problems ZDT3, CONV3-4, DTLZ7 and CONV4-2F. It is important to note that the cases for $k=2$ and $k\ge 3$ objectives exhibit slight variations in the component detection step and follow completely different procedures for the filling step.

Algorithm 1 Reference Set Generation (RSG)

Input:  Starting set $A_y$, filling size $N_f$, output size N. Output:  PF reference Z. 1: $C=\{C^1,\dots ,C^{nc}\}\leftarrow$ ComponentDetection($A_y$) 2: for $i=1,\dots ,nc$ do               ▹ for each component 3:      $F^i\leftarrow$ Filling($C^i,\lceil {\displaystyle \frac{N_f·\left|C^i\right|}{|A_y|}}\rceil$) 4: end for 5: $F\leftarrow {\bigcup }_i{F}^i$ 6: $Z\leftarrow$ reduction($F,N$) 7: return Z

3.3. Component Detection

Since the PF might be disconnected, we first apply a component detection on $A_y$. We use DBSCAN [62] in the objective space for this purpose for three main reasons: (i) the number of components does not need to be known a priori, (ii) the method detects outliers, and (iii) we have observed that a density-based approach works better than a distance-based one (e.g., k-means) for component detection. DBSCAN has two parameters: $minpts$ and r. To make the component detection process “parameter-free”, we compute a small grid-search to find the optimal values of $minpts$ and r based on the weakest link function defined in [63]. The default parameter values are $minpts\in \left\{\mathrm{2,3}\right\}$ and $r\in \{0.10\overline{d},0.11\overline{d},\dots ,0.15\overline{d}\}$ for bi-objective problems, and $minpts\in \left\{\mathrm{3,4}\right\}$ and $r\in \{0.19\overline{d},0.20\overline{d},\dots ,0.23\overline{d}\}$ otherwise, where $\overline{d}$ is the average pairwise distance between all points. However, these parameters can be adjusted when calling RSG. A summary of the component detection process is presented in Algorithm 2. In the following, we describe the remaining steps for a single connected component. If multiple components exist, the procedures must be repeated analogously for each component $C^i$ identified by Algorithm 2.

Algorithm 2 Component Detection

Input:  starting set $A_y=\{p_1,\dots ,p_{\ell }\}$, number of objectives k Output:  Number of clusters $nc$, clusters $C=\{C^1,\dots ,C^{nc}\}$. 1: Set $\overline{d}\leftarrow {\displaystyle \frac{2}{\ell (\ell -1)}}{\sum }_{p_i\ne p_j}|p_i-p_j|$ 2: if $k=2$ then 3:     Set $r\text{\_}range\leftarrow \{0.1\overline{d},0.11\overline{d},\dots ,0.16\overline{d}$} 4:     Set $minpts\text{\_}range\leftarrow \{\mathrm{2,3}$} 5: else 6:     Set $r\text{\_}range\leftarrow \{0.19\overline{d},0.20\overline{d},\dots ,0.23\overline{d}$} 7:     Set $minpts\text{\_}range\leftarrow \{\mathrm{3,4}$} 8: end if 9: Set $wlmin\leftarrow \infty$ 10: for $minpts\mathbf{in}minpts\text{\_}range$ do 11:     for $r\mathbf{in}r\text{\_}range$ do 12:         $C_t,n{c}_t\leftarrow$ DBSCAN($A_y,r,minpts$) 13:         $wl\leftarrow$ WeakestLink($C_t$) 14:         if $wl\le wlmin$ then 15:            $C\leftarrow C_t$ 16:            $nc\leftarrow n{c}_t$ 17:            $wlmin\leftarrow wl$ 18:         end if 19:     end for 20: end for 21: return $C=\{C^1,\dots ,C^{nc}\},nc$

3.4. Filling

Even if we know the PS a priori, a uniform sampling of the PS will not result in a uniform sampling of the PF. We assume that we have a set of points $A_y$ that is not uniformly distributed. However, if we fill the gaps and select N points from the filled set, we can obtain a more uniform set, leading to better IGD approximations when selecting points from these filled sets. The idea behind the filling step is to create a set that is as uniform as possible, so that the reduction step (in particular, k-means) does not get stuck in non-uniform local optima, which would lead to non-uniform final sets. The next task is, therefore, to compute $N_f$ solutions that are ideally uniformly distributed along $A_y$.

This process is performed differently for $k=2$ and $k\ge 3$ objectives:

• For $k=2$, we sort the points of $A_y=\{p_1,\dots ,p_{\ell }\}$ in increasing order of $f_1$, i.e., the first objective. Then, we consider the piecewise linear curve formed by the segments between $p_1$ and $p_2$, $p_2$ and $p_3$, and so on. The total length of this curve is given by $\left|L\right|={\sum }_{i=1}^{\ell -1}\left|L_i\right|$, where $|L_i|=\parallel p_i-p_{i+1}{\parallel }_2$. To perform the filling, we arrange the $N_f$ desired points along the curve L such that the first point is $p_1$ and the subsequent points are distributed equidistantly along L. This is achieved by placing each point at a distance of ${\delta }_{\ell }=\left|L\right|/(N_f-1)$ from the previous one along L. See Algorithm 3 for details.
• The filling process for $k\ge 3$ consists of several intermediate steps that must be described first; see Algorithm 4 for a general outline of the procedure. The procedure is as follows: First, to better represent $A_y$ (particularly for the filling step), we triangulate this set in $k-1$ dimensional space. This is done because the PF for continuous MOPs forms a set whose dimension is at most $k-1$. To achieve this, we compute a “normal vector” $\eta$ to $A_y$ using equation 9, and then we project it onto the $k-1$ hyperplane normal to $\eta$, obtaining the projected set $P_{k-1}$. After this, we compute the Delaunay triangulation [64] of $P_{k-1}$, which provides a triangulation $T^{\prime }$ that can be used in the original k-dimensional space. For concave PFs, the triangulation may include triangles (or simplex for $k>3$) that extend beyond $A_y$ (Figure 12(d)), so a removal strategy is applied to eliminate these triangles and obtain the final triangulation T. Finally, each triangle $t_i\in T$ is uniformly filled at random with a number of points proportional to its area (or volume for $k>3$), resulting in the filled set F of size $N_f$.
  We will now describe each step in more detail in the following:

  −

  Computing “normal vector” $\left(\eta \right)$. Since the front is not known, we compute the normal direction $\eta$ orthogonal to the convex hull defined by the minimal elements of $A_y$. More precisely, we compute $\eta$ as follows: if $A_y=\{p_1,\dots ,p_{\ell }\}$, choose

  $$y_{\left(i\right)}^m\in arg\underset{j=1,\dots ,N^{\prime }}{min}{p}_{j,i},i=1,\dots ,k,$$

  (5)

  where $p_{j,i}$ denotes the i-th element of $p_j$, and set

  $$M:=(y_{\left(2\right)}^m-y_{\left(1\right)}^m,y_{\left(3\right)}^m-y_{\left(1\right)}^m,\dots ,y_{\left(k\right)}^m-y_{\left(1\right)}^m)\in {\mathbb{R}}^{k\times (k-1)}.$$

  (6)

  Next, compute a QR-factorization of M, i.e.,

  $$M=QR=(q_1,\dots ,q_k)R,$$

  (7)

  where $Q\in {\mathbb{R}}^{k\times k}$ is an orthogonal matrix with column vectors $q_i$, and $R\in {\mathbb{R}}^{k\times (k-1)}$ is a right upper triangular matrix. Then, the vector

  $$\eta =-\mathrm{sgn}\left(q_{k,1}\right){\displaystyle \frac{q_k}{\parallel q_k{\parallel }_2}}$$

  (8)

  is the desired shifting direction. Since Q is orthogonal, the vectors $v_1:=q_1,\dots ,v_{k-1}:=q_{k-1}$ form an orthonormal basis of the hyperplane that is orthogonal to $\eta$. That is, these vectors can be used for the construction of $P_{k-1}$.

  −

  $k-1$ Projection $\left(P_{k-1}\right)$. We use $\eta$ as the first axis of a new coordinate system $(\eta ,v_1,\dots ,v_{k-1})$, where the vectors $v_i$ are defined as above. In this coordinate system, the orthonormal vectors $v_1,\dots ,v_{k-1}$ form the basis of a hyperplane orthogonal to $\eta$. The projection of $A_y$ onto this hyperplane ($P_{k-1}$) is achieved by first expressing $A_y$ in this new coordinate system as $P=\beta \eta +{\beta }_1{v}_1+\dots +{\beta }_{k-1}{v}_{k-1}$, and then removing the first coordinate, yielding $P_{k-1}={\beta }_1{v}_1+\dots +{\beta }_{k-1}{v}_{k-1}$.

  −

  Delaunay Triangulation $\left(T^{\prime }\right)$. Compute the Delaunay triangulation of $P_{k-1}$. This returns $T^{\prime }$, a list of size ${\delta }^{\prime }$ containing the indices of $P_{k-1}$ that form the triangles (or simplices for $k>3$). The list $T^{\prime }$ serves as the triangulation for the k-dimensional set $A_y$, which is possible because $T^{\prime }$ consists of indices, making it independent of the dimension. We use ${\delta }^{\prime }$ to denote the number of triangles obtained, $\{T^{\prime }(i,1),\dots ,T^{\prime }(i,k)\}$ the indices of the vertices forming triangle i and $\{p_{T^{\prime }(i,1)},\dots ,p_{T^{\prime }(i,k)}\}$ to denote the corresponding vertices of triangle i.

  −

  Triangle Cleaning $\left(T\right)$. We identify three types of unwanted triangles: those with large sides, those with large areas, and those where the matrix containing the coordinates of the vertices has a large condition number. The type of cleaning applied depends on the problem; however, the procedure remains the same for any problematic triangle case and is outlined in Algorithm 5. First, the property ${\rho }_i$ (area, largest side, or condition number) is computed for all the triangles $i=1,\dots ,\delta$. Next, triangles i with ${\rho }_i>\tau$ are removed. The parameter $\tau$ is also problem-dependent, and the specific values used for each problem will be detailed in the results section.

  −

  Triangle Filling $\left(F\right)$. For each triangle $t_i\in T$ with area $a_i$, we generate $\lceil {\displaystyle \frac{a_i}{A}}{N}_f\rceil$ points uniformly at random inside triangle $t_i$, following the procedure described in [65]. That is, the number of points is proportional to the area (or volume) of each triangle (or simplex). Here, $A={\sum }_{i=1}^{\delta }{a}_i$ is the total area of the triangulation.

Algorithm 3 Filling ($k=2$ Objectives)

Input:  starting set $A_y=\{p_1,\dots ,p_{\ell }\}$, filling size $N_f$ Output:  Filled set $F=\{y_1,\dots ,y_{N_f}\}$ 1: $X=\{x_1,\dots ,x_{\ell }\}\leftarrow$ sort $A_y$ according to its first objective $f_1$ 2: $L_i\leftarrow \parallel x_i-x_{i+1}\parallel \forall i=1,\dots ,\ell -1$ 3: $L\leftarrow {\sum }_{i=1}^{\ell -1}{L}_i$ 4: ${\delta }_{\ell }\leftarrow L/(N_f-1)$ 5: $dist\text{\_}left\leftarrow (0,0,\dots ,0)\in {\mathbb{R}}^{\ell }$ 6: for $i=1:\ell -1$ do            ▹ compute number of points per segment 7:     $ratio\leftarrow \left(\right|L_i|+dist\text{\_}left\left(i\right))/{\delta }_{\ell }$ 8:     $points\text{\_}per\text{\_}segment\left(i\right)\leftarrow$$\lfloor ratio\rfloor$ 9:     $dist\text{\_}left(i+1)\leftarrow (ratio-\lfloor ratio\rfloor )·{\delta }_{\ell }$ 10: end for 11: $count=1$ 12: for $i=1:\ell -1$ do              ▹ for each line segment 13:     if $points\text{\_}per\text{\_}segment\left(i\right)>0$ then      ▹ check if a point lands in segment $L_i$ 14:         ${\nu }_i:=(x_{i+1}-x_i)/L_i$ 15:         $y_{count}=x_i+({\delta }_{\ell }-dist\text{\_}left\left(i\right))·{\nu }_i$ 16:         $count\leftarrow count+1$ 17:         for $j=2:points\text{\_}per\text{\_}segment\left(i\right)$ do        ▹ if $L_i$ has more than one point 18:            $y_{count}\leftarrow y_{count-1}+{\delta }_{\ell }·{\nu }_i$ 19:            $count\leftarrow count+1$ 20:         end for 21:     end if 22: end for 23: return $F=\{y_1,\dots ,y_{N_f}\}$

Algorithm 4 Filling ($k\ge 3$ Objectives)

Input:  starting set $A_y=\{p_1,\dots ,p_{\ell }\}$, filling size $N_f$ Output:  Filled set $F=\{y_1,\dots ,y_{N_f}\}$ 1: $\eta \leftarrow$ normal_vector$\left(A_y\right)$ 2: $P_{k-1}\leftarrow$ projection$(A_y,\eta )$ 3: $T^{\prime }\leftarrow$ DelaunayTriangulation$\left(P_{k-1}\right)$ 4: $T\leftarrow$ TriangleCleaning$\left(T^{\prime }\right)$ 5: $F\leftarrow$ TriangleFilling$(T,A_y)$ 6: return $F=\{y_1,\dots ,y_{N_f}\}$

Algorithm 5: Triangle Cleaning

Input:  starting set $A_y=\{p_1,\dots ,p_l\}$, triangulation $T^{\prime }\in {\mathbb{R}}^{k\times {\delta }^{\prime }}$, number of triangles ${\delta }^{\prime }$, parameter threshold $\tau$ Output:  cleaned triangulation T of size $\delta$ 1: if chosen property is area then 2:     ${\rho }_i\leftarrow {\displaystyle \frac{1}{k!}}det{\left[\left(\begin{array}{c}\begin{array}{c}{p}_{T^{\prime }(i,2)}^T-p_{T^{\prime }(i,1)}^T\\ p_{T^{\prime }(i,3)}^T-p_{T^{\prime }(i,1)}^T\\ ⋮\\ p_{T^{\prime }(i,k)}^T-p_{T^{\prime }(i,1)}^T\end{array}\end{array}\right){\left(\begin{array}{c}\begin{array}{c}{p}_{T^{\prime }(i,2)}^T-p_{T^{\prime }(i,1)}^T\\ p_{T^{\prime }(i,3)}^T-p_{T^{\prime }(i,1)}^T\\ ⋮\\ p_{T^{\prime }(i,k)}^T-p_{T^{\prime }(i,1)}^T\end{array}\end{array}\right)}^T\right]}^{1/2}\forall i=1,\dots ,{\delta }^{\prime }$ 3: else if chosen property is largest side then 4:     ${\rho }_i\leftarrow$ largest side of simplex with vertices $\{p_{T^{\prime }(i,1)},\dots ,p_{T^{\prime }(i,k)}\}$ 5: else if chosen property is condition number then 6:     ${\rho }_i\leftarrow \kappa \left(\begin{array}{c}{p}_{T^{\prime }(i,1)}\\ p_{T^{\prime }(i,2)}\\ ⋮\\ p_{T^{\prime }(i,k)}\end{array}\right)$        ▹$\kappa \left(A\right)$ is the condition number of matrix A 7: end if 8: $\rho \leftarrow {\sum }_{i=1}^{{\delta }^{\prime }}{\rho }_i$ 9: $T\leftarrow T^{\prime }$ 10: if  ${\rho }_i>\tau ·\rho$  then 11:     remove triangle i from T 12: end if 13: return T, number of triangles $\delta$

3.5. Reduction

Once we have computed the filled set F, we need to select N points that are ideally evenly distributed along F. To this end, we use k-means clustering with N clusters, as there is a strong relationship between k-means and the optimal $IGD$ subset selection ([66,67]). The resulting N cluster centroids form the PF reference set Z. Note that this reduction method can be modified if needed to obtain an ideal reference according to other types of indicators (i.e., those that are not distance-based).

3.6. Obtaining $A_y$

RSG requires an initial approximation $A_y$ of the Pareto front. Note that by construction of the algorithm, this set can have small imperfections (which can be removed by the filling step) and can also have biases in the approximation (reduction step). However, it is desired that $A_y$ "captures" the shape of the entire Pareto front. The computation of such an approximation is certainly problem-dependent. For our computations, we have used the following three main procedures to get $A_y$:

• sampling: For some benchmark problems, either the Pareto set or the Pareto front is given in analytic form. If a sampling can be performed in objective space (e.g., for linear fronts, the remaining steps of the RSG may not be needed to further improve the quality of the solution set. If the sampling is performed in decision variable space, the elements of the resulting image $A_y$ may not be uniformly distributed along the Pareto front as discussed above. However, in that case, the filling and reduction step will help to remove biases.
• archiving: The result of an MOEA or any other MOP solver can, of course, be taken. This could be either the final archive of the population, via merging several populations of the same or several runs ([58]), or via using external (unbounded) archives ([50]). Note that this includes taking a reference set from a given repository. We have used archiving, e.g., for the test problems WFG3-9, DTLZ1-4, DTLZ7, ZDT1-6, CONV3, CONV3-4, and CONV4-2F.
• continuation: An alternative to the above mentioned techniques is to make use of multi-objective continuation methods, probably in combination with the use of several different starting points. In particular, we have used the Pareto Tracer (PT, [68,69,70]), a state-of-the-art continuation method that is able to treat problems of in principle any dimensions (both n and k), can handle general constraints and that can even detect local degeneration of the solution set. We have used PT, e.g., for the test problem WFG1, WFG2, DTLZ5, and DTLZ6.

3.7. Complexity Analysis

The overall complexity is $\mathcal{O}(\gamma {\ell }^2+{\tau }_{2}Nk{N}_f)$ for $k=2$ (regardless of the number of components), and $\mathcal{O}(\gamma {\ell }^2+\delta k^2+{\tau }_{2}Nk{N}_f)$ for $k\ge 3$, where ℓ is the size of the initial approximation $A_y$, k is the number of objectives, $\delta$ is the number of triangles in the Delaunay triangulation, ${\tau }_2$ is the number of iterations of k-means (bounded to 500 in this work), N is the desired size of the reference set Z, and $N_f$ is the size of the filling. This assumes $\ell \gg k$ and that the triangle-cleaning method used is based on the longest side (which was the method applied to all the references presented in this work). Typically, obtaining a decent approximation requires a large value of $N_f$, making the clustering step the dominant one and thus reducing the overall complexity to $\mathcal{O}\left({\tau }_{2}Nk{N}_f\right)$ for any k. We now present the complexity analysis in detail for each step separately, considering a single component. We now present the complexity analysis in detail for each step separately, considering a single component.

• Component Detection. The time complexity is $\mathcal{O}({(\ell )}^2+\gamma (\ell log(\ell )+{(\ell )}^2))$ which accounts for the computation of the average distance, plus the size of the grid search ($\gamma$) multiplied by the sum of the complexities of DBSCAN and the WeakestLink computation. Here, ℓ is the size of $A_y$, and $\gamma$ represents the number of parameter combinations of the grid search, with $\gamma =14$ for $k=2$ and $\gamma =10$ for $k\ge 3$ using the default values. If it is previously known that the Pareto front is connected, then the parameters of DBSCAN can be correctly adjusted, and $\gamma$ can be set to 1.
• Filling. The time complexity depends on the number of objectives:

  −

  For $k=2$ the time complexity is $\mathcal{O}(\ell +k{N}_f)$, which accounts for sorting and placing the $N_f$ points along the line segments.

  −

  For $k\ge 3$ the time complexity is $\mathcal{O}([k\ell +k{(k-1)}^2]+\ell k^2+\ell log\ell +[(k-1)\delta +k{N}_f])$ due to the computations involved in determining the normal vector $\eta$, changing coordinates and projecting, performing the Delaunay triangulation, and filling the triangles. Here, $\delta$ represents the size of the cleaned Delaunay triangulation, i.e., the number of triangles. Additionally, triangle cleaning must be considered, though its complexity depends on the method used. It is given by $\mathcal{O}\left(\delta k^3\right)$ when cleaning is based on area or the condition number (due to determinant computation), or $\mathcal{O}\left(\delta k^2\right)$ when cleaning based on the longest side.

• Select Reference Set T. The time complexity is $\mathcal{O}\left({\tau }_{2}Nk{N}_f\right)$ due to the k-means clustering algorithm. Here, ${\tau }_2$ is the number of iterations of k-means.

4. Results

In this section, we show the strength of the novel approach to selected test problems. We further show – as far as possible – comparisons to related methods.

First, we show all the steps of the RSG on five representative problems: Figure 8 shows these steps for ZDT3 ($k=2$ and disconnected), Figure 9 for DTLZ7 ($k=3$ and disconnected), Figure 11 for CONV3 (defined in Appendix A, $k=3$), Figure 12 for CONV3-4 (defined in Appendix A, $k=3$) and Figure 13 for CONV4-2F (defined in Appendix A, $k=4$ and disconnected).

In Figure 8 and Figure 9, we have used starting sets $A_y$ that have a slight bias. The last subfigures show the results of RSG for different values of N. For Figure 11, Figure 12 and Figure 13, we demonstrate the universality of RSG, as it does not need any analytical information of the PF since it utilizes the result of a MOEA and an archiver to generate a reference of any desired size.

Additionally, a comparison between RSG, PlatEMO [57], Pymoo, and the filling step is presented in Figure 14 and Figure 15 for the WFG2, ZDT1, ZDT3, DTLZ2, Convex DTLZ2 (CDTLZ2), C2-DTLZ2, and DTLZ7 test problems. Although the reference sets provided by PlatEMO and Pymoo are of high quality, they still exhibit some bias in certain bi-objective problems (such as ZDT1 and ZDT3), and especially in three-objective problems such as WFG2, and CDTLZ2. Furthermore, for problems like WFG2, and DTLZ7, the number of points in the reference sets of PlatEMO and Pymoo is limited to a fixed set, in contrast to RSG, which can generate any desired number of points.

Upon acceptance of this paper, we will provide reference sets on our website http://neo.cinvestav.mx:3005/ for the problems considered in this work as well as other problems (in the future, we intend to provide data for all commonly used benchmark test functions). The website will allow users to request reference sets of any desired size for any of the aforementioned problems. Further, the codes of RSG will be found on our GitHub page https://github.com/aerfangel/RSG.

5. Conclusions and Future Work

In this paper, we have addressed the problem of obtaining bias-free and complete finite-size approximations of the solution sets (Pareto fronts) of multi-objective optimization problems (MOPs). Such approximations are, in particular, required for the fair usage of distance-based performance indicators, which are frequently used in evolutionary multi-objective optimization (EMO). If the Pareto front approximations are biased or incomplete, the use of these performance indicators can lead to misleading or false information. To address this issue, we have proposed the Reference Set Generator (RSG). This method starts with an initial (probably biased) approximation of the Pareto front. An unbiased approximation is then computed via component detection, filling, and a reduction to the desired size. The RSG can be applied to Pareto fronts of any shape and dimension. We have finally demonstrated the strength of the novel approach on several benchmark problems.

In the future, we intend to use the RSG on the Pareto fronts of all commonly used continuous test problems. Special attention has to be paid to problems with degenerated fronts (i.e., problems where the Pareto front does locally not form an object of dimension $k-1$). In the current approach, we handle degeneracy using the PT to obtain a filled set and, from there, use the reduction step. For future work, we will explore if the projection can be modified to fill such sets, which may lead to a more general approach to handling degeneration.