New Findings on Indicator-based Multi-Objective Evolutionary Algorithms: A Brief Summary

Quality indicators (QIs) are functions that assign a real value to a set that represents the Pareto front approximation of a multi-objective optimization problem. In the evolutionary multi-objective optimization community, QIs have been mainly employed in two ways: (1) for the performance assessment of multi-objective evolutionary algorithms (MOEAs), which produce Pareto front approximations, and (2) to be adopted as the backbone of selection mechanisms of MOEAs. Regardless of the continuing advances on QIs and their utilization in MOEAs, there are currently a vast number of open questions in this research area. In this doctoral thesis, we have focused on two main research directions: the design of new selection mechanisms based on the competition and cooperation of multiple QIs, aiming to compensate for the weaknesses (in terms of convergence and diversity properties) of individual QIs with the strengths of the others. The second research axis is the generation of new QIs that are compliant with the Pareto dominance relation extended to sets. Such QIs have a direct impact on the type of conclusions that can be drawn about the performance of MOEAs. Our experimental results have shown that the use of multiple QIs either to design new selection mechanisms or to construct new Pareto-compliant QIs is a promising research direction that can improve the capa-bilities of MOEAs and that allows for a performance assessment of MOEAs with a higher degree of conﬁdence.


Introduction
In several engineering, scientific, and industrial fields, there is the necessity to solve problems that involve the simultaneous optimization of several objective functions. When the objective functions are in mutual conflict, these problems are commonly known as multi-objective optimization problems (MOPs). There are a number of mathematical programming techniques available to solve MOPs. However, these techniques are constrained by the mathematical properties of the objective functions, i.e., differentiability and continuity of the objective functions. If these properties are not present, mathematical programming methods will not succeed. In order to overcome this issue, bio-inspired metaheuristic techniques such as multi-objective evolutionary algorithms (MOEAs) [1] have become a popular choice, since they exhibit good results when tackling complex MOPs. MOEAs do not require the objective functions to fulfill specific mathematical properties and, in addition, they are population-based methods which allows the generation of multiple approximate solutions in a single run.
Currently, there are several paradigms to design MOEAs which have resulted in a plethora of MOEAs currently available [1,2]. Hence, an important question is how to know which MOEA has the best performance (in terms of convergence and diversity of solutions, among other aspects) on different MOPs. In the early days of evolutionary multi-objective optimization (EMO), performance assessment was based on qualitative observations of the Pareto front approximations produced by MOEAs. However, as the number of objective functions increases, a visual comparison is no longer possible. Consequently, quality indicators (QIs), which are set-based functions, started to be adopted to quantitatively assess Pareto front approxi- mations [3,4]. Due to the utilization of QIs, it is possible to assess MOEAs and draw conclusions which are mathematically sound. Since the late 1990s, several QIs have been proposed in the EMO field [4], each having different characteristics (see Figure 1). Throughout the years, QIs focused on assessing convergence and convergence-diversity have received particular attention because they indicate which MOEA produces Pareto front approximations closer to the true Pareto front. Among these convergence QIs, 1 a remarkable one is the hypervolume (HV) [5] indicator which measures the size of the objective space dominated by an approximation set bounded by an anti-optimal reference point. Moreover, HV and its variants are the only unary QIs that are known to be Pareto-compliant. In spite of the nice mathematical properties of the HV, its calculation is computationally expensive as the number of objective functions increases. Hence, some other less expensive but mathematically weaker QIs have been proposed in recent years. For instance, the R2 indicator [6], the inverted generational distance plus (IGD + ) [7], the additive indicator ( + ) [8], and the average Hausdorff distance indicator (∆ p ) [9].
The use of QIs is not only restricted to the assessment of MOEAs. Since a QI indicates how good an MOEA is, it is straightforward to use them to guide the evolutionary process of MOEAs. The underlying idea of the so-called indicator-based MOEAs (IB-MOEAs) is to employ a QI (or maybe more than one) as the backbone of a mechanism of an MOEA. Figure 2 shows the taxonomy of indicator-based mechanisms (IB-Mechanisms) proposed by Falcón-Cardona and Coello Coello [11]. The employment of IB-MOEAs has had a remarkable impact in the EMO field because these IB-Mechanisms increase the selection pressure of an MOEA, allowing it to tackle MOPs having a high number of objective functions (i.e., the so-called manyobjective optimization problems) and, furthermore, they introduce new preferences that bias the search behavior of an MOEA to discover new solutions. However, IB-MOEAs have been mainly designed using a single QI, inheriting its strengths but also its weaknesses.
Recently, Ishibuchi et al. [12] pointed out that some MOEAs are overspecialized on solving MOPs with regular Pareto front shapes. The geometries of these regular Pareto fronts are highly correlated with the form of a simplex. Mainly, MOEAs using a set of convex weight vectors 2 (which form a simplex) suffer from this overspecialization since if the Pareto front shape is not regular, these MOEAs do not behave properly. In this light, the authors emphasized that it is mandatory to design new MOEAs having a Pareto front shape invariance property.
The doctoral thesis [10] (in which this work is based on) is focused on two main problems in the EMO field. First, the design of Pareto front shape-invariant MOEAs through the use of selection mechanisms based on multiple QIs. To fulfill this first goal, we followed two directions: competition and cooperation of different QIs to compensate for the weaknesses of a given baseline QI with the strengths of the others. The use of multiple QIs is the underlying idea of the Multi-Indicator-based MOEAs (MIB-MOEAs) which had been scarcely explored by the EMO community [13,11]. Furthermore, our proposed MIB-MOEAs have demonstrated to be Pareto front shape-invariant, producing Pareto front approximation with good convergence and diversity properties. Our second research direction is the construction of new Paretocompliant QIs that exhibit preferences different from those of the HV. In order to accomplish this goal, we proposed a mathematical framework for the combination of multiple QIs (having specific mathematical properties) such that the resulting combined QI is Pareto-compliant. This mathematical framework has increased the number of currently available Pareto-compliant QIs.
The rest of this document is organized as follows. Section 2 provides some mathematical definitions required to make this paper self-contained. Section 3 is devoted to show our MIB-MOEAs based on a com-  petitive scheme. Section 4 describes our proposed approaches following a cooperative framework. Section 5 outlines how to design new Pareto-compliant QIs based on the combination of multiple indicators. Finally, our main conclusions and some possible paths for future work are shown in Section 6.

Background
Throughout this paper, we focus, without loss of generality, on unconstrained multi-objective optimization problems for minimization [1], which are defined as follows: T is the n-dimensional vector of decision variables and X ⊆ R n ; f j : X → R, j = 1, . . . , m are the objective functions.
Definition 2.1 (Unary Quality Indicator [8]). A unary quality indicator I is a function I : Ψ → R, which assigns a real value to a Pareto front approximation. Ψ is the set of all approximation sets. Definition 2.2 (Hypervolume indicator [5]). Given an anti-optimal reference point r ∈ R m , the hypervolume indicator is defined as follows: where L(·) denotes the Lebesgue measure in R m . Definition 2.3 (Unary R2 indicator [6]). The unary R2 indicator is defined as follows: where W is a set of weight vectors and u w : R m → R is a scalarizing function defined by a weight vector w ∈ W that assigns a real value to each m-dimensional vector.
Definition 2.4 (IGD + indicator [7]). The IGD + for minimization, is defined as follows: Definition 2.5 (Unary + indicator [8]). The unary + indicator gives the minimum distance by which a Pareto front approximation needs to or can be translated in each dimension in objective space such that a reference set is weakly dominated. Mathematically, it is defined as follows: To define the averaged Hausdorff distance (∆ p ), it is first necessary to introduce a variant of the indicators Generational Distance (GD) [3] and Inverted Generational Distance (IGD) [14], denoted as GD p and IGD p , respectively. Definition 2.6 (GD p indicator [9]).
Definition 2.8 (Averaged Hausdorff Distance indicator (∆ p ) [9]). For a given p > 0, ∆ p is defined as follows: As with IGD, the ∆ p indicator requires an aspiration set. ∆ p was proposed to eliminate some shortcomings of IGD such as its sensitivity to the cardinality of sets [9]. Definition 2.9 (Indicator contribution). Let I be any indicator in the set {HV, R2, IGD + , + , ∆ p }. The individual contribution C of a solution a ∈ A to the indicator value is given as follows: Definition 2.10 (Riesz s-energy [15]). For a given s > 0, the Riesz s-energy indicator is defined as follows: where · represents the Euclidean distance. As s → ∞, E s prefers more uniform solutions. This indicator measures the even distribution of a set of points in d-dimensional manifolds.
Definition 2.11 (Riesz s-energy individual contribution). The individual contribution C of a solution a ∈ A to the Riesz s-energy indicator is as follows:

Competitive Paradigm
In the specialized literature, there are several IB-MOEAs. A common feature of these IB-MOEAs is that they are mainly based on a single QI, which implies that an IB-MOEA benefits from its strengths and is affected by its weaknesses. For instance, HV-based MOEAs generate a high number of solutions around the Pareto front's knee and on the boundaries for concave Pareto fronts, while they tend to generate evenly distributed solutions for convex and linear Pareto fronts and they also have good performance on degenerate ones [16,17]. In contrast, R2-based MOEAs generate well-distributed solutions on concave Pareto fronts but they do not do the same on convex, degenerate or disconnected Pareto fronts since the convex weight vectors intersect such geometries [18,12]. In this section and in the remaining ones, we will overview how to design MIB-MOEAs that take advantage of multiple IB-Mechanisms or IB-MOEAs to improve the convergence and diversity properties, aiming to generate robust MOEAs, i.e., MOEAs which are invariant to the Pareto front shape.
A competition usually involves a scramble for resources where these resources can be food, territory, water, minerals, or even money. From the point of view of the design of MIB-MOEAs, we defined a competition as the capacity of an IB-Mechanism to be selected and executed as long as it remains possible. Under this paradigm, we proposed two approaches. On the one hand, the Multi-Indicator Hyper-heuristic (MIHPS) establishes a competition between four indicator-based density estimators (IB-DEs) [19]. The core idea is that the IB-DEs compete to be selected as many times as possible, biasing the decision of a hyper-heuristc that uses a Markov chain, based on their capability to improve the R2-based convergence graph, related  to the evolutionary process, during a given time window of analysis. On the other hand, we proposed the ensemble indicator-based MOEA (EIB-MOEA) [20] where five IB-DEs compete to obtain a greater weight value which is employed by the AdaBoost algorithm (an ensemble method belonging to the large class of machine learning algorithms). The larger the weight value, the greater the influence of the IB-DE in the selection process. The setting of the weight values is based on an online learning process that for each IB-DE considers its effect on the overall quality of the population. A greater emphasis is given to those IB-DEs that are capable of improving the quality of the population in terms of convergence and diversity. Due to space limitations, in the following sections, MIHPS and EIB-MOEA are only briefly sketched along with a discussion of their performance results. Figure 3 presents in a graphical way how MIHPS works. MIHPS is a steady-state MOEA 3 that uses the Pareto dominance relation 4 as the backbone of its environmental selection and, in addition, it has four IB-DEs based on the indicators: R2, IGD + , + , and ∆ p . Instead of using all the IB-DEs simultaneously, MIHPS selects a specific IB-DE to be employed during the next T w iterations. The selection of an IB-DE depends on an analysis of the R2 convergence graph of the algorithm. The waveform of the R2 convergence graph reflects the influence of all the IB-DEs that have been selected. If during the time-window T w assigned to a given IB-DE, the R2 convergence graph presents an ascending behavior based on a linear regression of the R2 samples, then the IB-DE is rewarded by increasing its probability to be selected. These selection probabilities are modelled by a Markov chain which is shown in Figure 4. In case that a descending behavior is exhibited, the IB-DE is punished by decreasing its selection probability. Hence, the IB-DEs compete to obtain a greater selection mechanism that biases the decision of the heuristic selection mechanism shown in Figure 3. The inspiration of MIHPS comes from the idea that each IB-DE would have a good performance on MOPs where the underlying QI has strong preferences. Hence, depending on the search difficulties and the Pareto front shape of the MOP, the strongest IB-DE will be continuously selected, i.e., it would have greater selection probabilities which weights the transitions in the Markov chain. MIHPS was compared with three state-of-the-art MOEAs that employ a set of convex weight vectors in  Overall, we found that some IB-DEs are consistently employed at specific stages of the evolutionary process (taken from [19]). some of their mechanisms: the Non-dominated Sorting Genetic Algorithm III (NSGA-III) [21], the Multi-Objective Evolutionary Algorithm based on Decomposition (MOEA/D) [22], and the R2-Evolutionary Multi-Objective Algorithm (R2-EMOA) [18]. These four algorithms were tested under similar conditions to allow a fair comparison. Thus, the same common parameters settings were employed. Additionally, to analyze the convergence and diversity properties of the adopted MOEAs, we decided to use a test suite of MOPs denoted as the Walking-Fish-Group (WFG) problems [23]. The WFG test suite comprises 9 nine MOPs with distinct search difficulties and Pareto front shapes that, in addition, are scalable in the number of objective functions. In total, we defined 45 test instances since there are nine WFG test problems and we employed them with 2, 3, 5, 6, and 10 objective functions. Figure 5 presents a summary of the obtained results where HV and the Riesz s-energy indicators were employed for the assessment of the algorithms. By observing the figure, it is clear that MIHPS is ranked first, considering both QIs, in most of the test instances while it is ranked fourth in a minimum number of times. These results imply that MIHPS is performing well on a wide variety of problems. Finally, a remarkable result in presented in Figure 6. BY analyzing the full results, we found that some IB-DEs are consistently employed at very specific stages of the evolutionary process. This result supports the idea that it is possible to design an MIB-MOEA that can exploit the strengths of its baseline QIs (IB-Mechanisms) depending on the MOP being tackled.

Ensemble Indicator-based Multi-Objective Evolutionary Algorithm
Currently, a hot research area is machine learning. In this research area, there are several algorithms focused on regression, classification and clustering, among other tasks. In order to promote a competition Ensemble−indicator−based density estimator among IB-Mechanisms, we took advantage of the ensemble methods which are employed in machine learning to construct meta-models. A remarkable aspect of ensemble methods is that they construct a strong learner based on the weaker learners. EIB-MOEA defines a competition among five IB-DEs (using HV, R2, IGD + , + , and ∆ p as baseline QIs) by using the AdaBoost ensemble method [24]. Under this approach, we consider the individual IB-DEs as the weaker learners that are to be combined to construct a stronger multi-indicator-based density estimator (MIB-DE). AdaBoost constructs a stronger learner on the basis of a linear combination of the weaker learners. Figure 7 shows the linear combination that EIB-MOEA defines. To this aim, given an approximation set A = { a 1 , a 2 , . . . , a N } where each a i ∈ R m , the stronger learner is defined as follows: In this equation, w j ≥ 0 is the weight value (relative importance) of the j th indicator I j , j = 1, 2, . . . , k. Furthermore, rank Ij ( a, A) denotes the rank assigned to the individual contribution of a to the indicator I j where rank 1 is related to the worst solution and rank N is given to the best one. Hence, the underlying idea is to make a consensus between all the IB-DEs to delete the worst-contributing solution to the overall quality of the population. The calculation of the weight values, which are the core of the linear combination, is based on an online learning approach that biases the reward or punishment assigned by the AdaBoost algorithm, using an exponential-loss approach. The online learning procedure performs a simulation of the effect of each IB-DE on the overall quality of the main population. If the decision drawn by the IB-DE helps to improve the quality of the population, this is accounted into a record structure. Similarly, if the effect of the IB-DE is negative, this is also registered in the record structure. AdaBoost employs the record structures of all the IB-DEs to update the weight values where a grater value is assigned to the best-performing IB-DEs. Due to the competition scheme of EIB-MOEA, it is expected that the best IB-DEs biased the decisions made throughout the whole evolutionary process.
To test the performance of EIB-MOEA, we compared it with five steady-state IB-MOEAs that use each baseline QI (HV, R2, IGD + , + , and ∆ p ) in an individual density estimator, giving rise to the algorithms: S-Metric-Selection Evolutionary Multi-Objective Algorithm (SMS-EMOA) [17], R2-EMOA [18], IGD + -based Many-Objective Evolutionary Algorithm (IGD + -MaOEA) [25], + -MaOEA, and ∆ p -MaOEA (the last two are defined similarly to IGD + -MaOEA), respectively. This experiment aims to determine if the competition scheme is responsible of better convergence and diversity properties than using individually the IB-Mechanisms. We adopted the test suites Deb-Thiele-Laumanns-Zitzler (DTLZ) [26], WFG, and their inverted versions defined by Ishibuchi et al. [12] which are denoted as DTLZ −1 and WFG −1 . We employed two-and three-objective versions of the MOPs of the adopted test suites. Regarding the performance assessment of the algorithms, we decided to utilize seven QIs to get a better idea of the improvements obtained by EIB-MOEA from which five of them (HV, R2, IGD + , + , and ∆ p ) are focused on measuring convergence, and the remaining ones (Riesz s-energy and Solow-Polasky Diversity) assess diversity. Figure 8 summarizes the experimental results by means of boxplots that represent the statistical ranks associated to the MOEAs.
In this case, the closer the mean to zero and the tighter the boxplot (implying a lower variance), the better the algorithm. Overall, EIB-MOEA is the best-performing algorithm in most of the seven QIs, especially for the convergence QIs. Regarding the diversity QIs, EIB-MOEA presents a competitive behavior. Some examples of the Pareto front approximations generated by EIB-MOEA and the adopted steady-state IB-MOEAs are presented in Figure 9. From this figure, it is possible to observe that EIB-MOEA has some drawbacks to produce uniformly distributed solutions which supports its competitive performance in terms of Riesz s-energy and the Solow-Polasky Diversity indicator (SPD) [27]. In spite of this lack of uniformity, EIB-MOEA succeeds to exploit the strengths of the individual IB-DEs to construct a stronger MIB-DE, regarding the convergence of solutions.

Cooperative Paradigm
In the previous section, we investigated competitive schemes to design MIB-MOEAs. Both MIHPS and EIB-MOEA share a common premise: to select or to give more importance to the IB-Mechanism that provides the best convergence improvement. Hence, depending on the search difficulties and the Pareto front shape, a specific IB-DE is expected to prevail over others although these ones help during some stages of the search to achieve good convergence or diversity results. In contrast, in this section, we explore a cooperative approach. In other words, we aim that all the IB-Mechanisms or IB-MOEAs work together to obtain Pareto front approximations with higher convergence and diversity properties. In this doctoral thesis, we proposed two cooperative approaches. First, we designed an MOEA in which two IB-Mechanisms (IB-DEs especifically) cooperate. This approach, denoted as CRI-EMOA [28], employs an IGD + -DE to promote convergence, and a density estimator based on the Riesz s-energy to improve the diversity of the population, aiming to obtain a Pareto front shape-invariant performance. Unlike MIHPS, where multiple IB-DEs strictly focused on convergence compete, the IB-DEs of CRI-EMOA are complementary since one is oriented to convergence and the other one only promotes diversity. On the other hand, we proposed an MIB-MOEA, called Cooperative Multi-Indicator-based MOEA (cMIB-MOEA) [29], where multiple IB-MOEAs cooperate through the island model. The island model is a well-known model to design parallel evolutionary algorithms. The underlying idea is to let multiple steady-state MOEAs based on a single QI to evolve independently during a given number of epochs. Then, they cooperate by using a migration process where some selected solutions of a given island (IB-MOEA) are sent to other islands. In the following sections, we summarize the properties of CRI-EMOA and cMIB-MOEA.

EMOA based on the Cooperation between Riesz s-energy and IGD +
The Evolutionary Multi-Objective Algorithm based on the Cooperation of Riesz s-energy and IGD + (CRI-EMOA) is a steady-state MIB-MOEA that uses a Pareto-based environmental selection and promotes the cooperation between two IB-DEs. Since IGD + is a convergence QI and the Riesz s-energy is a diversity indicator, instead of competing, they actually cooperate, aiming to improve simultaneously these two quality aspects of a Pareto front approximation. Instead of analyzing an R2-based convergence graph as MIHPS does, CRI-EMOA controls the switching between the two IB-DEs by analyzing a convergence graph composed of quality values of the population, using an approximation to HV (to reduce the computational cost associated to calculate the exact HV value). For a given time window (T w ), CRI-EMOA creates a linear model of the samples and looks at the angle θ of the linear model and the coefficient of variation β. If θ ∈ [−θ,θ] ∧ β ≤β, whereθ andβ are user-supplied thresholds (sufficiently small), and the integrated non-dominated sorting algorithm [30] has created a single layer, then this indicates that the overall convergence behavior of CRI-EMOA is stagnated. In consequence, it is necessary to execute the Riesz s-energy-based density estimator to increase the diversity of the population, exploring new regions of the search space. Otherwise, the IGD + -DE is executed to promote convergence towards the Pareto front. The switching between the IB-DEs, where each one is focused on improving a characteristic of the Pareto front approximation, is the core of this cooperative scheme of IB-Mechanisms.
Similarly to EIB-MOEA, CRI-EMOA was tested using the DTLZ, WFG, DTLZ −1 , and WFG −1 test suites, varying the number of objective functions. The performance of CRI-EMOA was compared with MOEAs using and not using a set of convex weight vectors as part of their mechanisms. NSGA-III, MOEA/D, and MOMBI2 [31] use convex weight vectors while GDE-MOEA and ∆ p -MOEA do not utilize them. Tables 1 and 2 summarize the performance comparisons regarding the HV and SPD indicators, respectively. These tables show the number of times CRI-EMOA and the adopted MOEAs obtained the best statistical values, considering all the test instances. Moreover, each table breaks down the results by the number of objective functions (3, 5, and 10) and if the problems belong to the original DTLZ and WFG test suites or to their inverted versions. According to these experimental results, CRI-EMOA is the best-performing algorithm for both HV and SPD. Moreover, it is worth noting that CRI-EMOA is consistently the best algorithm when tackling the DTLZ −1 and WFG −1 problems. This implies that CRI-EMOA is not affected by irregular Pareto front shapes. To support this conclusion, we can observe Figure 10 that shows the Pareto front approximations, produced by CRI-EMOA, of the three-objective DTLZ1, DTLZ2, DTLZ7, WFG1, and WFG2, and their corresponding inverted versions. For both the original and the inverted problems, CRI-EMOA generates Pareto front approximations with good convergence, diversity, and coverage. Hence, this indicates that CRI-EMOA effectively exploits the strengths of its cooperative IB-DEs and, what is more, the experimental results present some insights about the Pareto front shape invariance of CRI-EMOA. Consequently, in constrast to the competitive paradigm, the inner cooperation of IB-Mechanisms of CRI-EMOA is a stronger step towards the generation of a more general MOEA whose performance is not affected by the Pareto front shapes.

Cooperative Multi-Indicator-based Multi-Objective Evolutionary Algorithm
Unlike the previous section where we investigated the cooperation of multiple IB-DEs, in this section we propose to perform the cooperation of multiple steady-state IB-MOEAs, mainly, SMS-EMOA, R2-EMOA, IGD + -MaOEA, + -MaOEA, and ∆ p -MaOEA. The cooperation is carried out through the island model in which a number of subpopulations isolatedly evolve, each one in an island. After some epochs, a master   island (in our case) collects all subpopulations to keep the best individuals, according to some criterion, and, then, performs a migration process to increase the diversity of the islands. Our proposed approach, denoted as cMIB-MOEA, is composed, as shown in Figure 11, of five islands, each one having a steady-state IB-MOEA, and a master island that manages a Riesz s-energy-based archive to maintain uniformly distributed solutions coming from the IB-MOEAs. It is worth noting that each IB-MOEA evolves a micro-population of size µ/k, where µ is the maximum archive size and k is the number of islands. The reason to use micro-populations is to keep low the computational cost of each island, especially for the SMS-EMOA island. Hernández-Gómez et al. [32] empirically demonstrated that the overall computational cost of SMS-EMOA remains constant as long as the number of objective functions increases when using micro-populations, i.e., using population sizes ranging from 5 to 30 individuals. The micro-populations evolve independently during f mig epochs (i.e., generations). After this period, each island sends its whole micro-population to the master island where all of them are merged with the current contents of the master island to form a temporal population. This temporal population is filtered using the Pareto dominance Table 2: Number of times CRI-EMOA and the adopted MOEAs are ranked as the best-performing algorithm, regarding SPD. Master island Riesz s-energy-based archive Figure 11: Design of cMIB-MOEA using the island model parallel approach (taken from [29]). relation and, if the cardinality of the resulting set is greater than µ, a Riesz s-energy-based selection is performed to iteratively reduce the number of solutions until reaching µ individuals. Finally, the migration process is executed. Each island will receive n mig randomly chosen solutions (whose origin is different to this island) from the archive. The immigrant solutions will replace the worst n mig solutions of the target island, according to the contributions to the related QI. This migration process aims to ensure that each island maintains a well-diversified population to explore more regions of the search space. Finally, cMIB-MOEA returns the contents of the archive as its Pareto front approximation.
To investigate what is the effect of the cooperation between the baseline IB-MOEAs, we decided to compare cMIB-MOEA with panmictic versions of its baseline IB-MOEAs, i.e., each one using a population size of µ individuals. We selected MOPs from the DTLZ, WFG, DTLZ −1 , WFG −1 test suites, as well as the Viennet problems [33], and the Lamé superspheres [34], using two and three objective functions. We considered all these test suites because they represent a wide range of search difficulties and Pareto front shapes. We assessed the performance of cMIB-MOEA and the adopted IB-MOEAs using HV, R2, IGD + , + , ∆ p , and SPD to identify if cMIB-MOEA maintains a robust performance over several QIs preferences. Figure 12 shows a summary of the performance results, using a heat map. As the color gets darker, it means that the related algorithm is ranked first or second a greater number of times. Hence, cMIB-MOEA is clearly the best algorithm for HV, R2, ∆ p , and SPD if we look at the heat map related to the first ranks. HV, R2, and ∆ p indicate that cMIB-MOEA exhibits a good convergence behavior while SPD shows that cMIB-MOEA produces Pareto front approximations with a high degree of diversity. This last argument is supported by Figure 13 that compares some Pareto front approximations generated by cMIB-MOEA and the panmictic IB-MOEAs. For the five three-objective problems shown in the figure, cMIB-MOEA is the one that always presents good diversity and coverage. By summarizing the results, we have strong insights that cMIB-MOEA performs better than the panmictic versions of its baseline IB-MOEAs. In other words, the cooperative scheme where the weaknesses of some IB-MOEAs are compensated for the strengths of the others is better than using IB-MOEAs that use a single QI. Furthermore, similarly to CRI-EMOA, cMIB-MOEA exhibits a Pareto front shape-invariant behavior. Hence, the cooperative paradigm tends to be superior than the competitive one.

Combination of Quality Indicators
The Pareto dominance relation is a binary order relation applied on m-dimensional vectors. This binary relation can be extended to compare pairs of sets (i.e., Pareto front approximations). Let A and B be two sets, containing m-dimensional vectors, we say that A is better than B (denoted as A B) if ∀ b ∈ B, ∃ a ∈ A such that a b and A = B. The binary relation imposes a partial order in the set of all approximation sets (Ψ), thus, we can use it to compare Pareto front approximations. On the other hand, QIs also impose an order in Ψ but sometimes this order is not equal to that of which is the more general optimality criterion among Pareto front approximations. These QIs are known as not Pareto-compliant indicators. A drawback of these QIs is that they can draw misleading conclusions when comparing MOEAs. In contrast, there are QIs with stronger properties known as Pareto compliance and weak Pareto compliance.  Pareto-compliant QIs to increase the possibilities to compare MOEAs, using QIs with preferences different to those of HV and its variants.
To construct new Pareto-compliant QIs based on multiple indicators, we require to define a vector of indicators I CIV = (I 1 , I 2 , . . . , I k ) that we denote as compliant indicator vector (CIV) if ∀j = 1, 2, . . . , k, I j is weakly Pareto-compliant and there is at least an index t ∈ {1, . . . , k} such that I t is Pareto-compliant. In addition, we say that I : R k → R is a combined indicator. Consequently, we claim that I( I CIV ) is Paretocompliant if I is an order-preserving function. The following Theorem formalizes the previous discussion. The previous Theorem establishes a mathematical framework for the construction of new Pareto-compliant QIs. However, it is clear the need of a baseline Pareto-compliant QI (i.e., HV) for the construction of I. Nevertheless, our experimental results where we constructed Pareto-compliant versions of R2, IGD + , and + show that they have preferences different to those of HV. Hence, our framework increases the richness of Pareto-compliant indicators for the assessment of MOEA without drawing misleanding conclusions.

Conclusions and Future Work
In this doctoral thesis, we have tackled the design of MOEAs whose performance is not degraded by the Pareto front shapes of the underlying MOPs and the design of new Pareto-compliant QIs. Regarding the first problem, we proposed two paradigms: competition and cooperation of multiple IB-Mechanisms (and, in some cases, IB-MOEAs) to construct a more robust MOEA that presents a Pareto front shapeinvariant property. When analyzing both paradigms, we discovered that the cooperative scheme presents some advantages over the competitive approach. Our proposed approaches CRI-EMOA and cMIB-MOEA (that follow the cooperative paradigm) generate Pareto front approximations with good convergence and diversity properies, and, more remarkably, they are invariant to the Pareto front shapes. This represents a remarkable advancement in the design of MOEAs since the use of multiple QIs has shown to be a promising research direction to create robust MOEAs. Our second research axis was the design of new Pareto-compliant QIs. Through the use of an order-preserving function that assigns a real value to a compliant indicator vector, we demonstrated that the Pareto compliance property is present in the resulting combined QI. The impact of this result is a wider variety of Pareto-compliant QIs with preferences different to those of HV. As part of our future work, we aim to mathematically study the properties of the combined indicators and, also, investigate what is the practical effect of Pareto compliance in the search ability of IB-MOEAs. Additionally, we desire to directly compare the competitive and cooperative schemes to design better MIB-MOEAs.